itk/clitkBackProjectImageFilter.txx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   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
   8
   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.
  12
  13   It is distributed under dual licence
  14
  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"
  24
  25 namespace clitk
  26 {
  27
  28   //-----------------------------------------------------------------------
  29   //     Constructor
  30   //-----------------------------------------------------------------------
  31
  32   template<class InputImageType, class OutputImageType>
  33   BackProjectImageFilter< InputImageType, OutputImageType >
  34   ::BackProjectImageFilter()
  35   {
  36
  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.;
  45
  46     //Parameters for output
  47     this->m_OutputSpacing.Fill(1);
  48     this->m_OutputOrigin.Fill(0.0);
  49     this->m_OutputSize.Fill( 512 );
  50
  51     this->m_IsInitialized=false;
  52   }
  53
  54
  55   //-----------------------------------------------------------------------
  56   //     PrintSelf
  57   //-----------------------------------------------------------------------
  58   template<class InputImageType, class OutputImageType>
  59   void
  60   BackProjectImageFilter< InputImageType, OutputImageType >
  61   ::PrintSelf(std::ostream& os, itk::Indent indent) const
  62   {
  63     this->Superclass::PrintSelf(os,indent);
  64
  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;
  70   }
  71
  72
  73   //-----------------------------------------------------------------------
  74   //     Set output info from image
  75   //-----------------------------------------------------------------------
  76   template <class InputImageType, class OutputImageType>
  77   void
  78   BackProjectImageFilter<InputImageType, OutputImageType>
  79   ::SetOutputParametersFromImage ( const OutputImagePointer image )
  80   {
  81     this->SetOutputOrigin( image->GetOrigin() );
  82     this->SetOutputSpacing( image->GetSpacing() );
  83     this->SetOutputSize( image->GetLargestPossibleRegion().GetSize() );
  84   }
  85
  86
  87   //-----------------------------------------------------------------------
  88   //     Set output info from image
  89   //-----------------------------------------------------------------------
  90   template <class InputImageType, class OutputImageType>
  91   void
  92   BackProjectImageFilter<InputImageType, OutputImageType>
  93   ::SetOutputParametersFromImage (const OutputImageConstPointer image )
  94   {
  95     this->SetOutputOrigin( image->GetOrigin() );
  96     this->SetOutputSpacing( image->GetSpacing() );
  97     this->SetOutputSize( image->GetLargestPossibleRegion().GetSize() );
  98   }
  99
 100
 101   //-----------------------------------------------------------------------
 102   //     Generate output information
 103   //-----------------------------------------------------------------------
 104   template<class InputImageType, class OutputImageType>
 105   void
 106   BackProjectImageFilter< InputImageType, OutputImageType >
 107   ::GenerateOutputInformation( void )
 108   {
 109     // Don not call the superclass' implementation of this method (!=Dimensions)
 110     // Superclass::GenerateOutputInformation();
 111
 112     // get pointer to the output
 113     OutputImagePointer outputPtr = this->GetOutput();
 114     InputImageConstPointer  inputPtr  = this->GetInput();
 115
 116     if ( !outputPtr || ! inputPtr)
 117       {
 118         return;
 119       }
 120
 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());
 126
 127
 128     //OutputImageRegionType region;
 129     outputLargestPossibleRegion.SetSize( m_OutputSize );
 130     outputPtr->SetLargestPossibleRegion( outputLargestPossibleRegion );
 131     outputPtr->SetSpacing( m_OutputSpacing );
 132     outputPtr->SetOrigin( m_OutputOrigin );
 133
 134   }
 135
 136
 137   //-----------------------------------------------------------------------
 138   //     Generate input Requested region
 139   //-----------------------------------------------------------------------
 140   template <class InputImageType, class OutputImageType>
 141   void
 142   BackProjectImageFilter<InputImageType,OutputImageType>
 143   ::GenerateInputRequestedRegion()
 144   {
 145     //Call superclass
 146     Superclass::GenerateInputRequestedRegion();
 147
 148     if (!this->GetOutput())
 149       {
 150         return;
 151       }
 152     OutputImageRegionType outputRegion = this->GetOutput()->GetRequestedRegion();
 153     InputImagePointer inputPtr =  const_cast<InputImageType *>(this->GetInput());
 154
 155     if ( !inputPtr )
 156       {
 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());
 164         throw e;
 165       }
 166
 167     InputImageRegionType inputRegion;
 168     inputRegion=inputPtr->GetLargestPossibleRegion();
 169     inputPtr->SetRequestedRegion(inputRegion);
 170
 171   }
 172
 173
 174
 175   //-----------------------------------------------------------------------
 176   //     Before threaded data
 177   //-----------------------------------------------------------------------
 178   template <class InputImageType, class OutputImageType>
 179   void
 180   BackProjectImageFilter<InputImageType, OutputImageType>
 181   ::BeforeThreadedGenerateData( void )
 182   {
 183     if (!m_IsInitialized) this->Initialize();
 184   }
 185
 186
 187  //-----------------------------------------------------------------------
 188   //     Initialize
 189   //-----------------------------------------------------------------------
 190   template <class InputImageType, class OutputImageType>
 191   void
 192   BackProjectImageFilter<InputImageType, OutputImageType>
 193 #if ( ( ITK_VERSION_MAJOR == 4 ) && ( ITK_VERSION_MINOR > 12 ) || ( ITK_VERSION_MAJOR > 4 ))
 194   ::Initialize( void )
 195 #else
 196   ::Initialize( void ) throw (itk::ExceptionObject)
 197 #endif
 198   {
 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();
 205
 206     //Change the origin
 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);
 211
 212     // Calculate the projection matrix
 213     this->CalculateProjectionMatrix();
 214
 215     // Calculate the coordinate increments
 216     this->CalculateCoordinateIncrements();
 217
 218     m_IsInitialized=true;
 219   }
 220
 221
 222   //-----------------------------------------------------------------------
 223   //     Calculate Projection Matrix
 224   //-----------------------------------------------------------------------
 225   template <class InputImageType, class OutputImageType>
 226   void
 227   BackProjectImageFilter<InputImageType, OutputImageType>
 228   ::CalculateProjectionMatrix( void )
 229   {
 230     //InputSpacingType inputSpacing=this->GetInput()->GetSpacing();
 231
 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];
 237     temp(2,0)=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;
 243
 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;
 250
 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;
 255
 256     // Calculate the centered rotation matrix
 257     itk::Matrix<double,4,4> centeredRotationMatrix=GetCenteredRotationMatrix3D(rotationParameters,transformedIsoCenter);
 258
 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];
 262
 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;
 267
 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++)
 271         {
 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);
 275         }
 276   }
 277
 278
 279   //-----------------------------------------------------------------------
 280   //     Calculate the coordinate increments
 281   //-----------------------------------------------------------------------
 282   template <class InputImageType, class OutputImageType>
 283   void
 284   BackProjectImageFilter<InputImageType, OutputImageType>
 285   ::CalculateCoordinateIncrements( void )
 286   {
 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];
 291
 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];
 296
 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];
 301   }
 302
 303
 304   //-----------------------------------------------------------------------
 305   //     Threaded generate data
 306   //-----------------------------------------------------------------------
 307   template <class InputImageType, class OutputImageType>
 308   void
 309   BackProjectImageFilter<InputImageType, OutputImageType>::ThreadedGenerateData( const OutputImageRegionType & outputRegionForThread,  itk::ThreadIdType threadId )
 310   {
 311     //Projection pointer
 312     InputImageConstPointer inputPtr=this->GetInput();
 313
 314     //Volume pointer
 315     OutputImagePointer outputPTr= this->GetOutput();
 316
 317     //Iterator for the thread region
 318     typedef itk::ImageRegionIterator<OutputImageType> OutputIteratorType;
 319     OutputIteratorType iterator(outputPTr, outputRegionForThread);
 320     OutputSizeType outputSizeForThread=outputRegionForThread.GetSize();
 321
 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;
 331
 332     //Get the first output coordinate
 333     oIndex=iterator.GetIndex();//costly but only once a thread
 334     outputPTr->TransformIndexToPhysicalPoint(oIndex,oPoint);
 335
 336     for (unsigned int i=0;i<OutputImageDimension;i++) homOutputPoint[i]=oPoint[i];
 337     homOutputPoint[3]=1;
 338
 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];
 343
 344     typedef itk::LinearInterpolateImageFunction< InputImageType, double > InterpolatorType;
 345     typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
 346     interpolator->SetInputImage(this->GetInput());
 347
 348     //Run over all output voxels
 349     for (unsigned int i=0; i<outputSizeForThread[2]; i++)
 350       {
 351         for (unsigned int j=0; j<outputSizeForThread[1]; j++)
 352           {
 353             for (unsigned int k=0; k<outputSizeForThread[0]; k++)
 354               {
 355                 iPoint[0]=-homInputPoint[0]/homInputPoint[2] + m_PanelShift[0];
 356                 iPoint[1]=homInputPoint[1]/homInputPoint[2] + m_PanelShift[1];
 357
 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
 362                 else
 363                   value=m_EdgePaddingValue;
 364                 //Set it
 365                 iterator.Set(value);
 366
 367                 //Advance one step in the line
 368                 homInputPoint+=m_LineInc;
 369                 ++iterator;
 370               }
 371
 372             //Advance one line
 373             homInputPoint+=m_ColInc;
 374           }
 375
 376         //Advance one plane
 377         homInputPoint+=m_PlaneInc;
 378       }
 379   }
 380
 381
 382 } // namespace clitk
 383
 384
 385 #endif