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 ::Initialize( void ) throw (itk::ExceptionObject)
195 //Change the origin of the 2D input
196 typename InputImageType::ConstPointer inputPtr=this->GetInput();
197 m_ModifiedInput=InputImageType::New();
198 m_ModifiedInput->CopyInformation(inputPtr);
199 InputSizeType size=inputPtr->GetLargestPossibleRegion().GetSize();
200 InputSpacingType spacing=inputPtr->GetSpacing();
203 InputPointType newOrigin;
204 newOrigin[0] =-static_cast<double>(size[0]-1)*spacing[0]/2.0;//-static_cast<double>(size[0]-1)*spacing[0]/2.0;
205 newOrigin[1] =-static_cast<double>(size[1]-1)*spacing[1]/2.0;//-static_cast<double>(size[1]-1)*spacing[1]/2.0;
206 m_ModifiedInput->SetOrigin(newOrigin);
208 // Calculate the projection matrix
209 this->CalculateProjectionMatrix();
211 // Calculate the coordinate increments
212 this->CalculateCoordinateIncrements();
214 m_IsInitialized=true;
218 //-----------------------------------------------------------------------
219 // Calculate Projection Matrix
220 //-----------------------------------------------------------------------
221 template <class InputImageType, class OutputImageType>
223 BackProjectImageFilter<InputImageType, OutputImageType>
224 ::CalculateProjectionMatrix( void )
226 //InputSpacingType inputSpacing=this->GetInput()->GetSpacing();
228 // Projection on YZ plane+pixelTrans
229 itk::Matrix<double,3,4> temp;
230 temp.Fill(itk::NumericTraits<double>::Zero);
231 // temp(0,0)=-0.5*inputSpacing[0];
232 // temp(1,0)=-0.5*inputSpacing[1];
234 temp(0,1)=m_SourceToScreen;//Invert Y/X? for correspondence to real projections
235 temp(1,2)=m_SourceToScreen;
236 // temp(0,3)=m_SourceToAxis*0.5*inputSpacing[0];
237 // temp(1,3)=m_SourceToAxis*0.5*inputSpacing[1];
238 temp(2,3)=-m_SourceToAxis;//-m_IsoCenter[0]-m_SourceToAxis;
240 // Get the rotation parameter array
241 itk::Array<double> rotationParameters(3);
242 const double dtr = ( atan(1.0) * 4.0 ) / 180.0; //constant for converting degrees into radians
243 rotationParameters[0]= 0.;
244 rotationParameters[1]= 0.;
245 rotationParameters[2]= -dtr*(double) m_ProjectionAngle;
247 // We first perform rigid transform (of source and panel), then a centered rotation around the transformed center
248 itk::Matrix<double,3,3> rigidRotation = GetRotationalPartMatrix3D(m_RigidTransformMatrix);
249 itk::Vector<double,3> rigidTranslation = GetTranslationPartMatrix3D(m_RigidTransformMatrix);
250 OutputPointType transformedIsoCenter = rigidRotation * m_IsoCenter + rigidTranslation;
252 // Calculate the centered rotation matrix
253 itk::Matrix<double,4,4> centeredRotationMatrix=GetCenteredRotationMatrix3D(rotationParameters,transformedIsoCenter);
255 // Define the voxel pixel transforms (shift half a pixel/voxel)
256 itk::Matrix<double,4,4> voxelTrans; voxelTrans.SetIdentity();
257 for (unsigned int i=0;i<3;i++)voxelTrans(i,3)=0.5*m_OutputSpacing[i];
259 // Compose the rotation with the rigid transform matrix and the voxel transform
260 itk::Matrix<double,4,4> finalTransform = centeredRotationMatrix * m_RigidTransformMatrix;
261 itk::Matrix<double,4,4> invFinalTransform ( finalTransform.GetInverse());
262 invFinalTransform=invFinalTransform*voxelTrans;
264 // Apply rigid transform to the projection matrix
265 for (unsigned int i=0; i<3;i++)
266 for (unsigned int j=0; j<4;j++)
268 m_ProjectionMatrix(i,j)=itk::NumericTraits<double>::Zero;
269 for (unsigned int k=0; k<4;k++)
270 m_ProjectionMatrix(i,j)+=temp(i,k)*invFinalTransform(k,j);
275 //-----------------------------------------------------------------------
276 // Calculate the coordinate increments
277 //-----------------------------------------------------------------------
278 template <class InputImageType, class OutputImageType>
280 BackProjectImageFilter<InputImageType, OutputImageType>
281 ::CalculateCoordinateIncrements( void )
283 //Compute Line increment
284 m_LineInc[0]=m_ProjectionMatrix[0][0]*m_OutputSpacing[0];
285 m_LineInc[1]=m_ProjectionMatrix[1][0]*m_OutputSpacing[0];
286 m_LineInc[2]=m_ProjectionMatrix[2][0]*m_OutputSpacing[0];
288 //Compute column increment (and restart at the beginning of the line)
289 m_ColInc[0]=m_ProjectionMatrix[0][1]*m_OutputSpacing[1]-m_OutputSize[0]*m_LineInc[0];
290 m_ColInc[1]=m_ProjectionMatrix[1][1]*m_OutputSpacing[1]-m_OutputSize[0]*m_LineInc[1];
291 m_ColInc[2]=m_ProjectionMatrix[2][1]*m_OutputSpacing[1]-m_OutputSize[0]*m_LineInc[2];
293 //Compute plane increment (and restart at the beginning of the columns)
294 m_PlaneInc[0]=m_ProjectionMatrix[0][2]*m_OutputSpacing[2]-m_ProjectionMatrix[0][1]*m_OutputSpacing[1]*m_OutputSize[1];
295 m_PlaneInc[1]=m_ProjectionMatrix[1][2]*m_OutputSpacing[2]-m_ProjectionMatrix[1][1]*m_OutputSpacing[1]*m_OutputSize[1];
296 m_PlaneInc[2]=m_ProjectionMatrix[2][2]*m_OutputSpacing[2]-m_ProjectionMatrix[2][1]*m_OutputSpacing[1]*m_OutputSize[1];
300 //-----------------------------------------------------------------------
301 // Threaded generate data
302 //-----------------------------------------------------------------------
303 template <class InputImageType, class OutputImageType>
305 BackProjectImageFilter<InputImageType, OutputImageType>
306 ::ThreadedGenerateData( const OutputImageRegionType & outputRegionForThread, int threadId )
309 InputImageConstPointer inputPtr=this->GetInput();
312 OutputImagePointer outputPTr= this->GetOutput();
314 //Iterator for the thread region
315 typedef itk::ImageRegionIterator<OutputImageType> OutputIteratorType;
316 OutputIteratorType iterator(outputPTr, outputRegionForThread);
317 OutputSizeType outputSizeForThread=outputRegionForThread.GetSize();
319 //Some temp variables
320 OutputPixelType value;
321 HomogeneOutputPointType homOutputPoint;
322 HomogeneInputPointType homInputPoint;
323 OutputPointType oPoint;
324 InputPointType iPoint;
325 iPoint.Fill(itk::NumericTraits<double>::Zero);
326 OutputIndexType oIndex;
327 ContinuousInputIndexType iIndex;
329 //Get the first output coordinate
330 oIndex=iterator.GetIndex();//costly but only once a thread
331 outputPTr->TransformIndexToPhysicalPoint(oIndex,oPoint);
333 for (unsigned int i=0;i<OutputImageDimension;i++) homOutputPoint[i]=oPoint[i];
336 //Compute the first input coordinate (invert Y/X)
337 homInputPoint= (m_ProjectionMatrix * homOutputPoint);
338 iPoint[0]=-homInputPoint[0]/homInputPoint[2] + m_PanelShift[0];
339 iPoint[1]=homInputPoint[1]/homInputPoint[2] + m_PanelShift[1];
341 typedef itk::LinearInterpolateImageFunction< InputImageType, double > InterpolatorType;
342 typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
343 interpolator->SetInputImage(this->GetInput());
345 //Run over all output voxels
346 for (unsigned int i=0; i<outputSizeForThread[2]; i++)
348 for (unsigned int j=0; j<outputSizeForThread[1]; j++)
350 for (unsigned int k=0; k<outputSizeForThread[0]; k++)
352 iPoint[0]=-homInputPoint[0]/homInputPoint[2] + m_PanelShift[0];
353 iPoint[1]=homInputPoint[1]/homInputPoint[2] + m_PanelShift[1];
355 //Check wether inside, convert to index (use modified with correct origin)
356 if (m_ModifiedInput->TransformPhysicalPointToContinuousIndex(iPoint, iIndex) && interpolator->IsInsideBuffer(iIndex))
357 value = interpolator->EvaluateAtContinuousIndex(iIndex);
358 //Outside: padding value
360 value=m_EdgePaddingValue;
364 //Advance one step in the line
365 homInputPoint+=m_LineInc;
370 homInputPoint+=m_ColInc;
374 homInputPoint+=m_PlaneInc;