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();
44 //Parameters for output
45 this->m_OutputSpacing.Fill(1);
46 this->m_OutputOrigin.Fill(0.0);
47 this->m_OutputSize.Fill( 512 );
49 this->m_IsInitialized=false;
53 //-----------------------------------------------------------------------
55 //-----------------------------------------------------------------------
56 template<class InputImageType, class OutputImageType>
58 BackProjectImageFilter< InputImageType, OutputImageType >
59 ::PrintSelf(std::ostream& os, itk::Indent indent) const
61 this->Superclass::PrintSelf(os,indent);
63 os << indent << "IsoCenter: " << m_IsoCenter << std::endl;
64 os << indent << "SourceToScreen: " << m_SourceToScreen << std::endl;
65 os << indent << "SourceToAxis: " << m_SourceToAxis << std::endl;
66 os << indent << "Edge Padding Value: " << m_EdgePaddingValue << std::endl;
67 os << indent << "Rigid Transform matrix: " << m_EdgePaddingValue << std::endl;
71 //-----------------------------------------------------------------------
72 // Set output info from image
73 //-----------------------------------------------------------------------
74 template <class InputImageType, class OutputImageType>
76 BackProjectImageFilter<InputImageType, OutputImageType>
77 ::SetOutputParametersFromImage ( const OutputImagePointer image )
79 this->SetOutputOrigin( image->GetOrigin() );
80 this->SetOutputSpacing( image->GetSpacing() );
81 this->SetOutputSize( image->GetLargestPossibleRegion().GetSize() );
85 //-----------------------------------------------------------------------
86 // Set output info from image
87 //-----------------------------------------------------------------------
88 template <class InputImageType, class OutputImageType>
90 BackProjectImageFilter<InputImageType, OutputImageType>
91 ::SetOutputParametersFromImage (const OutputImageConstPointer image )
93 this->SetOutputOrigin( image->GetOrigin() );
94 this->SetOutputSpacing( image->GetSpacing() );
95 this->SetOutputSize( image->GetLargestPossibleRegion().GetSize() );
99 //-----------------------------------------------------------------------
100 // Generate output information
101 //-----------------------------------------------------------------------
102 template<class InputImageType, class OutputImageType>
104 BackProjectImageFilter< InputImageType, OutputImageType >
105 ::GenerateOutputInformation( void )
107 // Don not call the superclass' implementation of this method (!=Dimensions)
108 // Superclass::GenerateOutputInformation();
110 // get pointer to the output
111 OutputImagePointer outputPtr = this->GetOutput();
112 InputImageConstPointer inputPtr = this->GetInput();
114 if ( !outputPtr || ! inputPtr)
119 // Set the output image largest possible region. Use a RegionCopier
120 // so that the input and output images can be different dimensions.
121 OutputImageRegionType outputLargestPossibleRegion;
122 this->CallCopyInputRegionToOutputRegion(outputLargestPossibleRegion,
123 inputPtr->GetLargestPossibleRegion());
126 //OutputImageRegionType region;
127 outputLargestPossibleRegion.SetSize( m_OutputSize );
128 outputPtr->SetLargestPossibleRegion( outputLargestPossibleRegion );
129 outputPtr->SetSpacing( m_OutputSpacing );
130 outputPtr->SetOrigin( m_OutputOrigin );
135 //-----------------------------------------------------------------------
136 // Generate input Requested region
137 //-----------------------------------------------------------------------
138 template <class InputImageType, class OutputImageType>
140 BackProjectImageFilter<InputImageType,OutputImageType>
141 ::GenerateInputRequestedRegion()
144 Superclass::GenerateInputRequestedRegion();
146 if (!this->GetOutput())
150 OutputImageRegionType outputRegion = this->GetOutput()->GetRequestedRegion();
151 InputImagePointer inputPtr = const_cast<InputImageType *>(this->GetInput());
155 // Because DataObject::PropagateRequestedRegion() allows only
156 // InvalidRequestedRegionError, it's impossible to write simply:
157 // itkExceptionMacro(<< "Missing input " << idx);
158 itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
159 e.SetLocation(ITK_LOCATION);
160 e.SetDescription("Missing input.");
161 e.SetDataObject(this->GetOutput());
165 InputImageRegionType inputRegion;
166 inputRegion=inputPtr->GetLargestPossibleRegion();
167 inputPtr->SetRequestedRegion(inputRegion);
173 //-----------------------------------------------------------------------
174 // Before threaded data
175 //-----------------------------------------------------------------------
176 template <class InputImageType, class OutputImageType>
178 BackProjectImageFilter<InputImageType, OutputImageType>
179 ::BeforeThreadedGenerateData( void )
181 if (!m_IsInitialized) this->Initialize();
185 //-----------------------------------------------------------------------
187 //-----------------------------------------------------------------------
188 template <class InputImageType, class OutputImageType>
190 BackProjectImageFilter<InputImageType, OutputImageType>
191 ::Initialize( void ) throw (itk::ExceptionObject)
193 //Change the origin of the 2D input
194 typename InputImageType::ConstPointer inputPtr=this->GetInput();
195 m_ModifiedInput=InputImageType::New();
196 m_ModifiedInput->CopyInformation(inputPtr);
197 InputSizeType size=inputPtr->GetLargestPossibleRegion().GetSize();
198 InputSpacingType spacing=inputPtr->GetSpacing();
201 InputPointType newOrigin;
202 newOrigin[0] =-static_cast<double>(size[0]-1)*spacing[0]/2.0;//-static_cast<double>(size[0]-1)*spacing[0]/2.0;
203 newOrigin[1] =-static_cast<double>(size[1]-1)*spacing[1]/2.0;//-static_cast<double>(size[1]-1)*spacing[1]/2.0;
204 m_ModifiedInput->SetOrigin(newOrigin);
206 // Calculate the projection matrix
207 this->CalculateProjectionMatrix();
209 // Calculate the coordinate increments
210 this->CalculateCoordinateIncrements();
212 m_IsInitialized=true;
216 //-----------------------------------------------------------------------
217 // Calculate Projection Matrix
218 //-----------------------------------------------------------------------
219 template <class InputImageType, class OutputImageType>
221 BackProjectImageFilter<InputImageType, OutputImageType>
222 ::CalculateProjectionMatrix( void )
224 InputSpacingType inputSpacing=this->GetInput()->GetSpacing();
226 // Projection on YZ plane+pixelTrans
227 itk::Matrix<double,3,4> temp;
228 temp.Fill(itk::NumericTraits<double>::Zero);
229 // temp(0,0)=-0.5*inputSpacing[0];
230 // temp(1,0)=-0.5*inputSpacing[1];
232 temp(0,1)=m_SourceToScreen;//Invert Y/X? for correspondence to real projections
233 temp(1,2)=m_SourceToScreen;
234 // temp(0,3)=m_SourceToAxis*0.5*inputSpacing[0];
235 // temp(1,3)=m_SourceToAxis*0.5*inputSpacing[1];
236 temp(2,3)=-m_SourceToAxis;//-m_IsoCenter[0]-m_SourceToAxis;
238 // Get the rotation parameter array
239 itk::Array<double> rotationParameters(3);
240 const double dtr = ( atan(1.0) * 4.0 ) / 180.0; //constant for converting degrees into radians
241 rotationParameters[0]= 0.;
242 rotationParameters[1]= 0.;
243 rotationParameters[2]= -dtr*(double) m_ProjectionAngle;
245 // We first perform rigid transform (of source and panel), then a centered rotation around the transformed center
246 itk::Matrix<double,3,3> rigidRotation = GetRotationalPartMatrix3D(m_RigidTransformMatrix);
247 itk::Vector<double,3> rigidTranslation = GetTranslationPartMatrix3D(m_RigidTransformMatrix);
248 OutputPointType transformedIsoCenter = rigidRotation * m_IsoCenter + rigidTranslation;
250 // Calculate the centered rotation matrix
251 itk::Matrix<double,4,4> centeredRotationMatrix=GetCenteredRotationMatrix3D(rotationParameters,transformedIsoCenter);
253 // Define the voxel pixel transforms (shift half a pixel/voxel)
254 itk::Matrix<double,4,4> voxelTrans; voxelTrans.SetIdentity();
255 for (unsigned int i=0;i<3;i++)voxelTrans(i,3)=0.5*m_OutputSpacing[i];
257 // Compose the rotation with the rigid transform matrix and the voxel transform
258 itk::Matrix<double,4,4> finalTransform = centeredRotationMatrix * m_RigidTransformMatrix;
259 itk::Matrix<double,4,4> invFinalTransform ( finalTransform.GetInverse());
260 invFinalTransform=invFinalTransform*voxelTrans;
262 // Apply rigid transform to the projection matrix
263 for (unsigned int i=0; i<3;i++)
264 for (unsigned int j=0; j<4;j++)
266 m_ProjectionMatrix(i,j)=itk::NumericTraits<double>::Zero;
267 for (unsigned int k=0; k<4;k++)
268 m_ProjectionMatrix(i,j)+=temp(i,k)*invFinalTransform(k,j);
273 //-----------------------------------------------------------------------
274 // Calculate the coordinate increments
275 //-----------------------------------------------------------------------
276 template <class InputImageType, class OutputImageType>
278 BackProjectImageFilter<InputImageType, OutputImageType>
279 ::CalculateCoordinateIncrements( void )
281 //Compute Line increment
282 m_LineInc[0]=m_ProjectionMatrix[0][0]*m_OutputSpacing[0];
283 m_LineInc[1]=m_ProjectionMatrix[1][0]*m_OutputSpacing[0];
284 m_LineInc[2]=m_ProjectionMatrix[2][0]*m_OutputSpacing[0];
286 //Compute column increment (and restart at the beginning of the line)
287 m_ColInc[0]=m_ProjectionMatrix[0][1]*m_OutputSpacing[1]-m_OutputSize[0]*m_LineInc[0];
288 m_ColInc[1]=m_ProjectionMatrix[1][1]*m_OutputSpacing[1]-m_OutputSize[0]*m_LineInc[1];
289 m_ColInc[2]=m_ProjectionMatrix[2][1]*m_OutputSpacing[1]-m_OutputSize[0]*m_LineInc[2];
291 //Compute plane increment (and restart at the beginning of the columns)
292 m_PlaneInc[0]=m_ProjectionMatrix[0][2]*m_OutputSpacing[2]-m_ProjectionMatrix[0][1]*m_OutputSpacing[1]*m_OutputSize[1];
293 m_PlaneInc[1]=m_ProjectionMatrix[1][2]*m_OutputSpacing[2]-m_ProjectionMatrix[1][1]*m_OutputSpacing[1]*m_OutputSize[1];
294 m_PlaneInc[2]=m_ProjectionMatrix[2][2]*m_OutputSpacing[2]-m_ProjectionMatrix[2][1]*m_OutputSpacing[1]*m_OutputSize[1];
298 //-----------------------------------------------------------------------
299 // Threaded generate data
300 //-----------------------------------------------------------------------
301 template <class InputImageType, class OutputImageType>
303 BackProjectImageFilter<InputImageType, OutputImageType>
304 ::ThreadedGenerateData( const OutputImageRegionType & outputRegionForThread, int threadId )
307 InputImageConstPointer inputPtr=this->GetInput();
310 OutputImagePointer outputPTr= this->GetOutput();
312 //Iterator for the thread region
313 typedef itk::ImageRegionIterator<OutputImageType> OutputIteratorType;
314 OutputIteratorType iterator(outputPTr, outputRegionForThread);
315 OutputSizeType outputSizeForThread=outputRegionForThread.GetSize();
317 //Some temp variables
318 OutputPixelType value;
319 HomogeneOutputPointType homOutputPoint;
320 HomogeneInputPointType homInputPoint;
321 OutputPointType oPoint;
322 InputPointType iPoint;
323 iPoint.Fill(itk::NumericTraits<double>::Zero);
324 OutputIndexType oIndex;
325 ContinuousInputIndexType iIndex;
326 InputSizeType inputSize=inputPtr->GetLargestPossibleRegion().GetSize();
328 //Get the first output coordinate
329 oIndex=iterator.GetIndex();//costly but only once a thread
330 outputPTr->TransformIndexToPhysicalPoint(oIndex,oPoint);
332 for (unsigned int i=0;i<OutputImageDimension;i++) homOutputPoint[i]=oPoint[i];
335 //Compute the first input coordinate (invert Y/X)
336 homInputPoint= (m_ProjectionMatrix * homOutputPoint);
337 iPoint[0]=-homInputPoint[0]/homInputPoint[2];
338 iPoint[1]=homInputPoint[1]/homInputPoint[2];
340 typedef itk::LinearInterpolateImageFunction< InputImageType, double > InterpolatorType;
341 typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
342 interpolator->SetInputImage(this->GetInput());
344 //Run over all output voxels
345 for (unsigned int i=0; i<outputSizeForThread[2]; i++)
347 for (unsigned int j=0; j<outputSizeForThread[1]; j++)
349 for (unsigned int k=0; k<outputSizeForThread[0]; k++)
351 iPoint[0]=-homInputPoint[0]/homInputPoint[2] + m_PanelShift[0];
352 iPoint[1]=homInputPoint[1]/homInputPoint[2] + m_PanelShift[1];
354 //Check wether inside, convert to index (use modified with correct origin)
355 if (m_ModifiedInput->TransformPhysicalPointToContinuousIndex(iPoint, iIndex) && interpolator->IsInsideBuffer(iIndex))
356 value = interpolator->EvaluateAtContinuousIndex(iIndex);
357 //Outside: padding value
359 value=m_EdgePaddingValue;
363 //Advance one step in the line
364 homInputPoint+=m_LineInc;
369 homInputPoint+=m_ColInc;
373 homInputPoint+=m_PlaneInc;