1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
5 - University of LYON http://www.universite-lyon.fr/
6 - Léon Bérard cancer center http://www.centreleonberard.fr
7 - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
9 This software is distributed WITHOUT ANY WARRANTY; without even
10 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11 PURPOSE. See the copyright notices for more information.
13 It is distributed under dual licence
15 - BSD See included LICENSE.txt file
16 - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ===========================================================================**/
18 #ifndef __clitkBackProjectImageFilter_txx
19 #define __clitkBackProjectImageFilter_txx
20 #include "clitkBackProjectImageFilter.h"
21 #include "itkContinuousIndex.h"
22 #include "vnl/vnl_math.h"
23 #include "itkLinearInterpolateImageFunction.h"
28 //-----------------------------------------------------------------------
30 //-----------------------------------------------------------------------
32 template<class InputImageType, class OutputImageType>
33 BackProjectImageFilter< InputImageType, OutputImageType >
34 ::BackProjectImageFilter()
37 //Parameters for backprojection
38 this->m_IsoCenter.Fill(0.0);
39 this->m_SourceToScreen = 1536.0;
40 this->m_SourceToAxis = 1000.0;
41 this->m_EdgePaddingValue = itk::NumericTraits<OutputPixelType>::Zero;//density images
42 this->m_RigidTransformMatrix.SetIdentity();
43 this->m_PanelShift[0] = 0.;
44 this->m_PanelShift[1] = 0.;
46 //Parameters for output
47 this->m_OutputSpacing.Fill(1);
48 this->m_OutputOrigin.Fill(0.0);
49 this->m_OutputSize.Fill( 512 );
51 this->m_IsInitialized=false;
55 //-----------------------------------------------------------------------
57 //-----------------------------------------------------------------------
58 template<class InputImageType, class OutputImageType>
60 BackProjectImageFilter< InputImageType, OutputImageType >
61 ::PrintSelf(std::ostream& os, itk::Indent indent) const
63 this->Superclass::PrintSelf(os,indent);
65 os << indent << "IsoCenter: " << m_IsoCenter << std::endl;
66 os << indent << "SourceToScreen: " << m_SourceToScreen << std::endl;
67 os << indent << "SourceToAxis: " << m_SourceToAxis << std::endl;
68 os << indent << "Edge Padding Value: " << m_EdgePaddingValue << std::endl;
69 os << indent << "Rigid Transform matrix: " << m_EdgePaddingValue << std::endl;
73 //-----------------------------------------------------------------------
74 // Set output info from image
75 //-----------------------------------------------------------------------
76 template <class InputImageType, class OutputImageType>
78 BackProjectImageFilter<InputImageType, OutputImageType>
79 ::SetOutputParametersFromImage ( const OutputImagePointer image )
81 this->SetOutputOrigin( image->GetOrigin() );
82 this->SetOutputSpacing( image->GetSpacing() );
83 this->SetOutputSize( image->GetLargestPossibleRegion().GetSize() );
87 //-----------------------------------------------------------------------
88 // Set output info from image
89 //-----------------------------------------------------------------------
90 template <class InputImageType, class OutputImageType>
92 BackProjectImageFilter<InputImageType, OutputImageType>
93 ::SetOutputParametersFromImage (const OutputImageConstPointer image )
95 this->SetOutputOrigin( image->GetOrigin() );
96 this->SetOutputSpacing( image->GetSpacing() );
97 this->SetOutputSize( image->GetLargestPossibleRegion().GetSize() );
101 //-----------------------------------------------------------------------
102 // Generate output information
103 //-----------------------------------------------------------------------
104 template<class InputImageType, class OutputImageType>
106 BackProjectImageFilter< InputImageType, OutputImageType >
107 ::GenerateOutputInformation( void )
109 // Don not call the superclass' implementation of this method (!=Dimensions)
110 // Superclass::GenerateOutputInformation();
112 // get pointer to the output
113 OutputImagePointer outputPtr = this->GetOutput();
114 InputImageConstPointer inputPtr = this->GetInput();
116 if ( !outputPtr || ! inputPtr)
121 // Set the output image largest possible region. Use a RegionCopier
122 // so that the input and output images can be different dimensions.
123 OutputImageRegionType outputLargestPossibleRegion;
124 this->CallCopyInputRegionToOutputRegion(outputLargestPossibleRegion,
125 inputPtr->GetLargestPossibleRegion());
128 //OutputImageRegionType region;
129 outputLargestPossibleRegion.SetSize( m_OutputSize );
130 outputPtr->SetLargestPossibleRegion( outputLargestPossibleRegion );
131 outputPtr->SetSpacing( m_OutputSpacing );
132 outputPtr->SetOrigin( m_OutputOrigin );
137 //-----------------------------------------------------------------------
138 // Generate input Requested region
139 //-----------------------------------------------------------------------
140 template <class InputImageType, class OutputImageType>
142 BackProjectImageFilter<InputImageType,OutputImageType>
143 ::GenerateInputRequestedRegion()
146 Superclass::GenerateInputRequestedRegion();
148 if (!this->GetOutput())
152 OutputImageRegionType outputRegion = this->GetOutput()->GetRequestedRegion();
153 InputImagePointer inputPtr = const_cast<InputImageType *>(this->GetInput());
157 // Because DataObject::PropagateRequestedRegion() allows only
158 // InvalidRequestedRegionError, it's impossible to write simply:
159 // itkExceptionMacro(<< "Missing input " << idx);
160 itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
161 e.SetLocation(ITK_LOCATION);
162 e.SetDescription("Missing input.");
163 e.SetDataObject(this->GetOutput());
167 InputImageRegionType inputRegion;
168 inputRegion=inputPtr->GetLargestPossibleRegion();
169 inputPtr->SetRequestedRegion(inputRegion);
175 //-----------------------------------------------------------------------
176 // Before threaded data
177 //-----------------------------------------------------------------------
178 template <class InputImageType, class OutputImageType>
180 BackProjectImageFilter<InputImageType, OutputImageType>
181 ::BeforeThreadedGenerateData( void )
183 if (!m_IsInitialized) this->Initialize();
187 //-----------------------------------------------------------------------
189 //-----------------------------------------------------------------------
190 template <class InputImageType, class OutputImageType>
192 BackProjectImageFilter<InputImageType, OutputImageType>
193 #if ( ( ITK_VERSION_MAJOR == 4 ) && ( ITK_VERSION_MINOR > 12 ) || ( ITK_VERSION_MAJOR > 4 ))
196 ::Initialize( void ) throw (itk::ExceptionObject)
199 //Change the origin of the 2D input
200 typename InputImageType::ConstPointer inputPtr=this->GetInput();
201 m_ModifiedInput=InputImageType::New();
202 m_ModifiedInput->CopyInformation(inputPtr);
203 InputSizeType size=inputPtr->GetLargestPossibleRegion().GetSize();
204 InputSpacingType spacing=inputPtr->GetSpacing();
207 InputPointType newOrigin;
208 newOrigin[0] =-static_cast<double>(size[0]-1)*spacing[0]/2.0;//-static_cast<double>(size[0]-1)*spacing[0]/2.0;
209 newOrigin[1] =-static_cast<double>(size[1]-1)*spacing[1]/2.0;//-static_cast<double>(size[1]-1)*spacing[1]/2.0;
210 m_ModifiedInput->SetOrigin(newOrigin);
212 // Calculate the projection matrix
213 this->CalculateProjectionMatrix();
215 // Calculate the coordinate increments
216 this->CalculateCoordinateIncrements();
218 m_IsInitialized=true;
222 //-----------------------------------------------------------------------
223 // Calculate Projection Matrix
224 //-----------------------------------------------------------------------
225 template <class InputImageType, class OutputImageType>
227 BackProjectImageFilter<InputImageType, OutputImageType>
228 ::CalculateProjectionMatrix( void )
230 //InputSpacingType inputSpacing=this->GetInput()->GetSpacing();
232 // Projection on YZ plane+pixelTrans
233 itk::Matrix<double,3,4> temp;
234 temp.Fill(itk::NumericTraits<double>::Zero);
235 // temp(0,0)=-0.5*inputSpacing[0];
236 // temp(1,0)=-0.5*inputSpacing[1];
238 temp(0,1)=m_SourceToScreen;//Invert Y/X? for correspondence to real projections
239 temp(1,2)=m_SourceToScreen;
240 // temp(0,3)=m_SourceToAxis*0.5*inputSpacing[0];
241 // temp(1,3)=m_SourceToAxis*0.5*inputSpacing[1];
242 temp(2,3)=-m_SourceToAxis;//-m_IsoCenter[0]-m_SourceToAxis;
244 // Get the rotation parameter array
245 itk::Array<double> rotationParameters(3);
246 const double dtr = ( atan(1.0) * 4.0 ) / 180.0; //constant for converting degrees into radians
247 rotationParameters[0]= 0.;
248 rotationParameters[1]= 0.;
249 rotationParameters[2]= -dtr*(double) m_ProjectionAngle;
251 // We first perform rigid transform (of source and panel), then a centered rotation around the transformed center
252 itk::Matrix<double,3,3> rigidRotation = GetRotationalPartMatrix3D(m_RigidTransformMatrix);
253 itk::Vector<double,3> rigidTranslation = GetTranslationPartMatrix3D(m_RigidTransformMatrix);
254 OutputPointType transformedIsoCenter = rigidRotation * m_IsoCenter + rigidTranslation;
256 // Calculate the centered rotation matrix
257 itk::Matrix<double,4,4> centeredRotationMatrix=GetCenteredRotationMatrix3D(rotationParameters,transformedIsoCenter);
259 // Define the voxel pixel transforms (shift half a pixel/voxel)
260 itk::Matrix<double,4,4> voxelTrans; voxelTrans.SetIdentity();
261 for (unsigned int i=0;i<3;i++)voxelTrans(i,3)=0.5*m_OutputSpacing[i];
263 // Compose the rotation with the rigid transform matrix and the voxel transform
264 itk::Matrix<double,4,4> finalTransform = centeredRotationMatrix * m_RigidTransformMatrix;
265 itk::Matrix<double,4,4> invFinalTransform ( finalTransform.GetInverse());
266 invFinalTransform=invFinalTransform*voxelTrans;
268 // Apply rigid transform to the projection matrix
269 for (unsigned int i=0; i<3;i++)
270 for (unsigned int j=0; j<4;j++)
272 m_ProjectionMatrix(i,j)=itk::NumericTraits<double>::Zero;
273 for (unsigned int k=0; k<4;k++)
274 m_ProjectionMatrix(i,j)+=temp(i,k)*invFinalTransform(k,j);
279 //-----------------------------------------------------------------------
280 // Calculate the coordinate increments
281 //-----------------------------------------------------------------------
282 template <class InputImageType, class OutputImageType>
284 BackProjectImageFilter<InputImageType, OutputImageType>
285 ::CalculateCoordinateIncrements( void )
287 //Compute Line increment
288 m_LineInc[0]=m_ProjectionMatrix[0][0]*m_OutputSpacing[0];
289 m_LineInc[1]=m_ProjectionMatrix[1][0]*m_OutputSpacing[0];
290 m_LineInc[2]=m_ProjectionMatrix[2][0]*m_OutputSpacing[0];
292 //Compute column increment (and restart at the beginning of the line)
293 m_ColInc[0]=m_ProjectionMatrix[0][1]*m_OutputSpacing[1]-m_OutputSize[0]*m_LineInc[0];
294 m_ColInc[1]=m_ProjectionMatrix[1][1]*m_OutputSpacing[1]-m_OutputSize[0]*m_LineInc[1];
295 m_ColInc[2]=m_ProjectionMatrix[2][1]*m_OutputSpacing[1]-m_OutputSize[0]*m_LineInc[2];
297 //Compute plane increment (and restart at the beginning of the columns)
298 m_PlaneInc[0]=m_ProjectionMatrix[0][2]*m_OutputSpacing[2]-m_ProjectionMatrix[0][1]*m_OutputSpacing[1]*m_OutputSize[1];
299 m_PlaneInc[1]=m_ProjectionMatrix[1][2]*m_OutputSpacing[2]-m_ProjectionMatrix[1][1]*m_OutputSpacing[1]*m_OutputSize[1];
300 m_PlaneInc[2]=m_ProjectionMatrix[2][2]*m_OutputSpacing[2]-m_ProjectionMatrix[2][1]*m_OutputSpacing[1]*m_OutputSize[1];
304 //-----------------------------------------------------------------------
305 // Threaded generate data
306 //-----------------------------------------------------------------------
307 template <class InputImageType, class OutputImageType>
309 BackProjectImageFilter<InputImageType, OutputImageType>::ThreadedGenerateData( const OutputImageRegionType & outputRegionForThread, itk::ThreadIdType threadId )
312 InputImageConstPointer inputPtr=this->GetInput();
315 OutputImagePointer outputPTr= this->GetOutput();
317 //Iterator for the thread region
318 typedef itk::ImageRegionIterator<OutputImageType> OutputIteratorType;
319 OutputIteratorType iterator(outputPTr, outputRegionForThread);
320 OutputSizeType outputSizeForThread=outputRegionForThread.GetSize();
322 //Some temp variables
323 OutputPixelType value;
324 HomogeneOutputPointType homOutputPoint;
325 HomogeneInputPointType homInputPoint;
326 OutputPointType oPoint;
327 InputPointType iPoint;
328 iPoint.Fill(itk::NumericTraits<double>::Zero);
329 OutputIndexType oIndex;
330 ContinuousInputIndexType iIndex;
332 //Get the first output coordinate
333 oIndex=iterator.GetIndex();//costly but only once a thread
334 outputPTr->TransformIndexToPhysicalPoint(oIndex,oPoint);
336 for (unsigned int i=0;i<OutputImageDimension;i++) homOutputPoint[i]=oPoint[i];
339 //Compute the first input coordinate (invert Y/X)
340 homInputPoint= (m_ProjectionMatrix * homOutputPoint);
341 iPoint[0]=-homInputPoint[0]/homInputPoint[2] + m_PanelShift[0];
342 iPoint[1]=homInputPoint[1]/homInputPoint[2] + m_PanelShift[1];
344 typedef itk::LinearInterpolateImageFunction< InputImageType, double > InterpolatorType;
345 typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
346 interpolator->SetInputImage(this->GetInput());
348 //Run over all output voxels
349 for (unsigned int i=0; i<outputSizeForThread[2]; i++)
351 for (unsigned int j=0; j<outputSizeForThread[1]; j++)
353 for (unsigned int k=0; k<outputSizeForThread[0]; k++)
355 iPoint[0]=-homInputPoint[0]/homInputPoint[2] + m_PanelShift[0];
356 iPoint[1]=homInputPoint[1]/homInputPoint[2] + m_PanelShift[1];
358 //Check wether inside, convert to index (use modified with correct origin)
359 if (m_ModifiedInput->TransformPhysicalPointToContinuousIndex(iPoint, iIndex) && interpolator->IsInsideBuffer(iIndex))
360 value = interpolator->EvaluateAtContinuousIndex(iIndex);
361 //Outside: padding value
363 value=m_EdgePaddingValue;
367 //Advance one step in the line
368 homInputPoint+=m_LineInc;
373 homInputPoint+=m_ColInc;
377 homInputPoint+=m_PlaneInc;