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://oncora1.lyon.fnclcc.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"
27 //-----------------------------------------------------------------------
29 //-----------------------------------------------------------------------
31 template<class InputImageType, class OutputImageType>
32 BackProjectImageFilter< InputImageType, OutputImageType >
33 ::BackProjectImageFilter()
36 //Parameters for backprojection
37 this->m_IsoCenter.Fill(0.0);
38 this->m_SourceToScreen = 1536.0;
39 this->m_SourceToAxis = 1000.0;
40 this->m_EdgePaddingValue = itk::NumericTraits<OutputPixelType>::Zero;//density images
41 this->m_RigidTransformMatrix.SetIdentity();
43 //Parameters for output
44 this->m_OutputSpacing.Fill(1);
45 this->m_OutputOrigin.Fill(0.0);
46 this->m_OutputSize.Fill( 512 );
48 this->m_IsInitialized=false;
52 //-----------------------------------------------------------------------
54 //-----------------------------------------------------------------------
55 template<class InputImageType, class OutputImageType>
57 BackProjectImageFilter< InputImageType, OutputImageType >
58 ::PrintSelf(std::ostream& os, itk::Indent indent) const
60 this->Superclass::PrintSelf(os,indent);
62 os << indent << "IsoCenter: " << m_IsoCenter << std::endl;
63 os << indent << "SourceToScreen: " << m_SourceToScreen << std::endl;
64 os << indent << "SourceToAxis: " << m_SourceToAxis << std::endl;
65 os << indent << "Edge Padding Value: " << m_EdgePaddingValue << std::endl;
66 os << indent << "Rigid Transform matrix: " << m_EdgePaddingValue << std::endl;
70 //-----------------------------------------------------------------------
71 // Set output info from image
72 //-----------------------------------------------------------------------
73 template <class InputImageType, class OutputImageType>
75 BackProjectImageFilter<InputImageType, OutputImageType>
76 ::SetOutputParametersFromImage ( const OutputImagePointer image )
78 this->SetOutputOrigin( image->GetOrigin() );
79 this->SetOutputSpacing( image->GetSpacing() );
80 this->SetOutputSize( image->GetLargestPossibleRegion().GetSize() );
84 //-----------------------------------------------------------------------
85 // Set output info from image
86 //-----------------------------------------------------------------------
87 template <class InputImageType, class OutputImageType>
89 BackProjectImageFilter<InputImageType, OutputImageType>
90 ::SetOutputParametersFromImage (const OutputImageConstPointer image )
92 this->SetOutputOrigin( image->GetOrigin() );
93 this->SetOutputSpacing( image->GetSpacing() );
94 this->SetOutputSize( image->GetLargestPossibleRegion().GetSize() );
98 //-----------------------------------------------------------------------
99 // Generate output information
100 //-----------------------------------------------------------------------
101 template<class InputImageType, class OutputImageType>
103 BackProjectImageFilter< InputImageType, OutputImageType >
104 ::GenerateOutputInformation( void )
106 // Don not call the superclass' implementation of this method (!=Dimensions)
107 // Superclass::GenerateOutputInformation();
109 // get pointer to the output
110 OutputImagePointer outputPtr = this->GetOutput();
111 InputImageConstPointer inputPtr = this->GetInput();
113 if ( !outputPtr || ! inputPtr)
118 // Set the output image largest possible region. Use a RegionCopier
119 // so that the input and output images can be different dimensions.
120 OutputImageRegionType outputLargestPossibleRegion;
121 this->CallCopyInputRegionToOutputRegion(outputLargestPossibleRegion,
122 inputPtr->GetLargestPossibleRegion());
125 //OutputImageRegionType region;
126 outputLargestPossibleRegion.SetSize( m_OutputSize );
127 outputPtr->SetLargestPossibleRegion( outputLargestPossibleRegion );
128 outputPtr->SetSpacing( m_OutputSpacing );
129 outputPtr->SetOrigin( m_OutputOrigin );
134 //-----------------------------------------------------------------------
135 // Generate input Requested region
136 //-----------------------------------------------------------------------
137 template <class InputImageType, class OutputImageType>
139 BackProjectImageFilter<InputImageType,OutputImageType>
140 ::GenerateInputRequestedRegion()
143 Superclass::GenerateInputRequestedRegion();
145 if (!this->GetOutput())
149 OutputImageRegionType outputRegion = this->GetOutput()->GetRequestedRegion();
150 InputImagePointer inputPtr = const_cast<InputImageType *>(this->GetInput());
154 // Because DataObject::PropagateRequestedRegion() allows only
155 // InvalidRequestedRegionError, it's impossible to write simply:
156 // itkExceptionMacro(<< "Missing input " << idx);
157 itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
158 e.SetLocation(ITK_LOCATION);
159 e.SetDescription("Missing input.");
160 e.SetDataObject(this->GetOutput());
164 InputImageRegionType inputRegion;
165 inputRegion=inputPtr->GetLargestPossibleRegion();
166 inputPtr->SetRequestedRegion(inputRegion);
172 //-----------------------------------------------------------------------
173 // Before threaded data
174 //-----------------------------------------------------------------------
175 template <class InputImageType, class OutputImageType>
177 BackProjectImageFilter<InputImageType, OutputImageType>
178 ::BeforeThreadedGenerateData( void )
180 if (!m_IsInitialized) this->Initialize();
184 //-----------------------------------------------------------------------
186 //-----------------------------------------------------------------------
187 template <class InputImageType, class OutputImageType>
189 BackProjectImageFilter<InputImageType, OutputImageType>
190 ::Initialize( void ) throw (itk::ExceptionObject)
192 //Change the origin of the 2D input
193 typename InputImageType::ConstPointer inputPtr=this->GetInput();
194 m_ModifiedInput=InputImageType::New();
195 m_ModifiedInput->CopyInformation(inputPtr);
196 InputSizeType size=inputPtr->GetLargestPossibleRegion().GetSize();
197 InputSpacingType spacing=inputPtr->GetSpacing();
200 InputPointType newOrigin;
201 newOrigin[0] =-static_cast<double>(size[0]-1)*spacing[0]/2.0;//-static_cast<double>(size[0]-1)*spacing[0]/2.0;
202 newOrigin[1] =-static_cast<double>(size[1]-1)*spacing[1]/2.0;//-static_cast<double>(size[1]-1)*spacing[1]/2.0;
203 m_ModifiedInput->SetOrigin(newOrigin);
205 // Calculate the projection matrix
206 this->CalculateProjectionMatrix();
208 // Calculate the coordinate increments
209 this->CalculateCoordinateIncrements();
211 m_IsInitialized=true;
215 //-----------------------------------------------------------------------
216 // Calculate Projection Matrix
217 //-----------------------------------------------------------------------
218 template <class InputImageType, class OutputImageType>
220 BackProjectImageFilter<InputImageType, OutputImageType>
221 ::CalculateProjectionMatrix( void )
223 InputSpacingType inputSpacing=this->GetInput()->GetSpacing();
225 // Projection on YZ plane+pixelTrans
226 itk::Matrix<double,3,4> temp;
227 temp.Fill(itk::NumericTraits<double>::Zero);
228 // temp(0,0)=-0.5*inputSpacing[0];
229 // temp(1,0)=-0.5*inputSpacing[1];
231 temp(0,1)=m_SourceToScreen;//Invert Y/X? for correspondence to real projections
232 temp(1,2)=m_SourceToScreen;
233 // temp(0,3)=m_SourceToAxis*0.5*inputSpacing[0];
234 // temp(1,3)=m_SourceToAxis*0.5*inputSpacing[1];
235 temp(2,3)=-m_SourceToAxis;//-m_IsoCenter[0]-m_SourceToAxis;
237 // Get the rotation parameter array
238 itk::Array<double> rotationParameters(3);
239 const double dtr = ( atan(1.0) * 4.0 ) / 180.0; //constant for converting degrees into radians
240 rotationParameters[0]= 0.;
241 rotationParameters[1]= 0.;
242 rotationParameters[2]= -dtr*(double) m_ProjectionAngle;
244 // We first perform rigid transform (of source and panel), then a centered rotation around the transformed center
245 itk::Matrix<double,3,3> rigidRotation = GetRotationalPartMatrix3D(m_RigidTransformMatrix);
246 itk::Vector<double,3> rigidTranslation = GetTranslationPartMatrix3D(m_RigidTransformMatrix);
247 OutputPointType transformedIsoCenter = rigidRotation * m_IsoCenter + rigidTranslation;
249 // Calculate the centered rotation matrix
250 itk::Matrix<double,4,4> centeredRotationMatrix=GetCenteredRotationMatrix3D(rotationParameters,transformedIsoCenter);
252 // Define the voxel pixel transforms (shift half a pixel/voxel)
253 itk::Matrix<double,4,4> voxelTrans; voxelTrans.SetIdentity();
254 for (unsigned int i=0;i<3;i++)voxelTrans(i,3)=0.5*m_OutputSpacing[i];
256 // Compose the rotation with the rigid transform matrix and the voxel transform
257 itk::Matrix<double,4,4> finalTransform = centeredRotationMatrix * m_RigidTransformMatrix;
258 itk::Matrix<double,4,4> invFinalTransform ( finalTransform.GetInverse());
259 invFinalTransform=invFinalTransform*voxelTrans;
261 // Apply rigid transform to the projection matrix
262 for (unsigned int i=0; i<3;i++)
263 for (unsigned int j=0; j<4;j++)
265 m_ProjectionMatrix(i,j)=itk::NumericTraits<double>::Zero;
266 for (unsigned int k=0; k<4;k++)
267 m_ProjectionMatrix(i,j)+=temp(i,k)*invFinalTransform(k,j);
272 //-----------------------------------------------------------------------
273 // Calculate the coordinate increments
274 //-----------------------------------------------------------------------
275 template <class InputImageType, class OutputImageType>
277 BackProjectImageFilter<InputImageType, OutputImageType>
278 ::CalculateCoordinateIncrements( void )
280 //Compute Line increment
281 m_LineInc[0]=m_ProjectionMatrix[0][0]*m_OutputSpacing[0];
282 m_LineInc[1]=m_ProjectionMatrix[1][0]*m_OutputSpacing[0];
283 m_LineInc[2]=m_ProjectionMatrix[2][0]*m_OutputSpacing[0];
285 //Compute column increment (and restart at the beginning of the line)
286 m_ColInc[0]=m_ProjectionMatrix[0][1]*m_OutputSpacing[1]-m_OutputSize[0]*m_LineInc[0];
287 m_ColInc[1]=m_ProjectionMatrix[1][1]*m_OutputSpacing[1]-m_OutputSize[0]*m_LineInc[1];
288 m_ColInc[2]=m_ProjectionMatrix[2][1]*m_OutputSpacing[1]-m_OutputSize[0]*m_LineInc[2];
290 //Compute plane increment (and restart at the beginning of the columns)
291 m_PlaneInc[0]=m_ProjectionMatrix[0][2]*m_OutputSpacing[2]-m_ProjectionMatrix[0][1]*m_OutputSpacing[1]*m_OutputSize[1];
292 m_PlaneInc[1]=m_ProjectionMatrix[1][2]*m_OutputSpacing[2]-m_ProjectionMatrix[1][1]*m_OutputSpacing[1]*m_OutputSize[1];
293 m_PlaneInc[2]=m_ProjectionMatrix[2][2]*m_OutputSpacing[2]-m_ProjectionMatrix[2][1]*m_OutputSpacing[1]*m_OutputSize[1];
297 //-----------------------------------------------------------------------
298 // Threaded generate data
299 //-----------------------------------------------------------------------
300 template <class InputImageType, class OutputImageType>
302 BackProjectImageFilter<InputImageType, OutputImageType>
303 ::ThreadedGenerateData( const OutputImageRegionType & outputRegionForThread, int threadId )
306 InputImageConstPointer inputPtr=this->GetInput();
307 InputPixelType * beginPtr=const_cast<InputPixelType *>(this->GetInput()->GetBufferPointer());
311 OutputImagePointer outputPTr= this->GetOutput();
313 //Iterator for the thread region
314 typedef itk::ImageRegionIterator<OutputImageType> OutputIteratorType;
315 OutputIteratorType iterator(outputPTr, outputRegionForThread);
316 OutputSizeType outputSizeForThread=outputRegionForThread.GetSize();
318 //Some temp variables
319 OutputPixelType value;
320 HomogeneOutputPointType homOutputPoint;
321 HomogeneInputPointType homInputPoint;
322 OutputPointType oPoint;
323 InputPointType iPoint;
324 iPoint.Fill(itk::NumericTraits<double>::Zero);
325 OutputIndexType oIndex;
326 ContinuousInputIndexType iIndex;
327 InputSizeType inputSize=inputPtr->GetLargestPossibleRegion().GetSize();
328 double dx,dy,dxm,dym;
331 //Get the first output coordinate
332 oIndex=iterator.GetIndex();//costly but only once a thread
333 outputPTr->TransformIndexToPhysicalPoint(oIndex,oPoint);
335 for (unsigned int i=0;i<OutputImageDimension;i++) homOutputPoint[i]=oPoint[i];
338 //Compute the first input coordinate (invert Y/X)
339 homInputPoint= (m_ProjectionMatrix * homOutputPoint);
340 iPoint[0]=-homInputPoint[0]/homInputPoint[2];
341 iPoint[1]=homInputPoint[1]/homInputPoint[2];
343 //Run over all output voxels
344 for (unsigned int i=0; i<outputSizeForThread[2]; i++)
346 for (unsigned int j=0; j<outputSizeForThread[1]; j++)
348 for (unsigned int k=0; k<outputSizeForThread[0]; k++)
350 iPoint[0]=homInputPoint[0]/homInputPoint[2];
351 iPoint[1]=homInputPoint[1]/homInputPoint[2];
353 //Check wether inside, convert to index (use modified with correct origin)
354 if( m_ModifiedInput->TransformPhysicalPointToContinuousIndex(iPoint, iIndex) )
356 //Own (fast bilinear) interpolation
357 lx = (int)floor(iIndex[0]); dx = iIndex[0]-lx; dxm = 1.-dx;
358 ly = (int)floor(iIndex[1]); dy = iIndex[1]-ly; dym = 1.-dy;
359 pp = beginPtr + ly*inputSize[0]+lx;
360 value =static_cast<OutputPixelType>( ( dxm * dym*(double)(*pp)
361 + dx * dym*(double)(*(pp+1))
362 + dxm* dy *(double)(*(pp + inputSize[0]))
363 + dx * dy *(double)(*(pp + inputSize[0]+1))) );
366 //Outside: padding value
367 else value=m_EdgePaddingValue;
372 //Advance one step in the line
373 homInputPoint+=m_LineInc;
378 homInputPoint+=m_ColInc;
382 homInputPoint+=m_PlaneInc;