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
  44     //Parameters for output
  45     this->m_OutputSpacing.Fill(1);
  46     this->m_OutputOrigin.Fill(0.0);
  47     this->m_OutputSize.Fill( 512 );
  48
  49     this->m_IsInitialized=false;
  50   }
  51
  52
  53   //-----------------------------------------------------------------------
  54   //     PrintSelf
  55   //-----------------------------------------------------------------------
  56   template<class InputImageType, class OutputImageType>
  57   void
  58   BackProjectImageFilter< InputImageType, OutputImageType >
  59   ::PrintSelf(std::ostream& os, itk::Indent indent) const
  60   {
  61     this->Superclass::PrintSelf(os,indent);
  62
  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;
  68   }
  69
  70
  71   //-----------------------------------------------------------------------
  72   //     Set output info from image
  73   //-----------------------------------------------------------------------
  74   template <class InputImageType, class OutputImageType>
  75   void
  76   BackProjectImageFilter<InputImageType, OutputImageType>
  77   ::SetOutputParametersFromImage ( const OutputImagePointer image )
  78   {
  79     this->SetOutputOrigin( image->GetOrigin() );
  80     this->SetOutputSpacing( image->GetSpacing() );
  81     this->SetOutputSize( image->GetLargestPossibleRegion().GetSize() );
  82   }
  83
  84
  85   //-----------------------------------------------------------------------
  86   //     Set output info from image
  87   //-----------------------------------------------------------------------
  88   template <class InputImageType, class OutputImageType>
  89   void
  90   BackProjectImageFilter<InputImageType, OutputImageType>
  91   ::SetOutputParametersFromImage (const OutputImageConstPointer image )
  92   {
  93     this->SetOutputOrigin( image->GetOrigin() );
  94     this->SetOutputSpacing( image->GetSpacing() );
  95     this->SetOutputSize( image->GetLargestPossibleRegion().GetSize() );
  96   }
  97
  98
  99   //-----------------------------------------------------------------------
 100   //     Generate output information
 101   //-----------------------------------------------------------------------
 102   template<class InputImageType, class OutputImageType>
 103   void
 104   BackProjectImageFilter< InputImageType, OutputImageType >
 105   ::GenerateOutputInformation( void )
 106   {
 107     // Don not call the superclass' implementation of this method (!=Dimensions)
 108     // Superclass::GenerateOutputInformation();
 109
 110     // get pointer to the output
 111     OutputImagePointer outputPtr = this->GetOutput();
 112     InputImageConstPointer  inputPtr  = this->GetInput();
 113
 114     if ( !outputPtr || ! inputPtr)
 115       {
 116         return;
 117       }
 118
 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());
 124
 125
 126     //OutputImageRegionType region;
 127     outputLargestPossibleRegion.SetSize( m_OutputSize );
 128     outputPtr->SetLargestPossibleRegion( outputLargestPossibleRegion );
 129     outputPtr->SetSpacing( m_OutputSpacing );
 130     outputPtr->SetOrigin( m_OutputOrigin );
 131
 132   }
 133
 134
 135   //-----------------------------------------------------------------------
 136   //     Generate input Requested region
 137   //-----------------------------------------------------------------------
 138   template <class InputImageType, class OutputImageType>
 139   void
 140   BackProjectImageFilter<InputImageType,OutputImageType>
 141   ::GenerateInputRequestedRegion()
 142   {
 143     //Call superclass
 144     Superclass::GenerateInputRequestedRegion();
 145
 146     if (!this->GetOutput())
 147       {
 148         return;
 149       }
 150     OutputImageRegionType outputRegion = this->GetOutput()->GetRequestedRegion();
 151     InputImagePointer inputPtr =  const_cast<InputImageType *>(this->GetInput());
 152
 153     if ( !inputPtr )
 154       {
 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());
 162         throw e;
 163       }
 164
 165     InputImageRegionType inputRegion;
 166     inputRegion=inputPtr->GetLargestPossibleRegion();
 167     inputPtr->SetRequestedRegion(inputRegion);
 168
 169   }
 170
 171
 172
 173   //-----------------------------------------------------------------------
 174   //     Before threaded data
 175   //-----------------------------------------------------------------------
 176   template <class InputImageType, class OutputImageType>
 177   void
 178   BackProjectImageFilter<InputImageType, OutputImageType>
 179   ::BeforeThreadedGenerateData( void )
 180   {
 181     if (!m_IsInitialized) this->Initialize();
 182   }
 183
 184
 185  //-----------------------------------------------------------------------
 186   //     Initialize
 187   //-----------------------------------------------------------------------
 188   template <class InputImageType, class OutputImageType>
 189   void
 190   BackProjectImageFilter<InputImageType, OutputImageType>
 191   ::Initialize( void ) throw (itk::ExceptionObject)
 192   {
 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();
 199
 200     //Change the origin
 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);
 205
 206     // Calculate the projection matrix
 207     this->CalculateProjectionMatrix();
 208
 209     // Calculate the coordinate increments
 210     this->CalculateCoordinateIncrements();
 211
 212     m_IsInitialized=true;
 213   }
 214
 215
 216   //-----------------------------------------------------------------------
 217   //     Calculate Projection Matrix
 218   //-----------------------------------------------------------------------
 219   template <class InputImageType, class OutputImageType>
 220   void
 221   BackProjectImageFilter<InputImageType, OutputImageType>
 222   ::CalculateProjectionMatrix( void )
 223   {
 224     InputSpacingType inputSpacing=this->GetInput()->GetSpacing();
 225
 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];
 231     temp(2,0)=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;
 237
 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;
 244
 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;
 249
 250     // Calculate the centered rotation matrix
 251     itk::Matrix<double,4,4> centeredRotationMatrix=GetCenteredRotationMatrix3D(rotationParameters,transformedIsoCenter);
 252
 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];
 256
 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;
 261
 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++)
 265         {
 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);
 269         }
 270   }
 271
 272
 273   //-----------------------------------------------------------------------
 274   //     Calculate the coordinate increments
 275   //-----------------------------------------------------------------------
 276   template <class InputImageType, class OutputImageType>
 277   void
 278   BackProjectImageFilter<InputImageType, OutputImageType>
 279   ::CalculateCoordinateIncrements( void )
 280   {
 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];
 285
 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];
 290
 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];
 295   }
 296
 297
 298   //-----------------------------------------------------------------------
 299   //     Threaded generate data
 300   //-----------------------------------------------------------------------
 301   template <class InputImageType, class OutputImageType>
 302   void
 303   BackProjectImageFilter<InputImageType, OutputImageType>
 304   ::ThreadedGenerateData( const OutputImageRegionType & outputRegionForThread,  int threadId )
 305   {
 306     //Projection pointer
 307     InputImageConstPointer inputPtr=this->GetInput();
 308
 309     //Volume pointer
 310     OutputImagePointer outputPTr= this->GetOutput();
 311
 312     //Iterator for the thread region
 313     typedef itk::ImageRegionIterator<OutputImageType> OutputIteratorType;
 314     OutputIteratorType iterator(outputPTr, outputRegionForThread);
 315     OutputSizeType outputSizeForThread=outputRegionForThread.GetSize();
 316
 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();
 327
 328     //Get the first output coordinate
 329     oIndex=iterator.GetIndex();//costly but only once a thread
 330     outputPTr->TransformIndexToPhysicalPoint(oIndex,oPoint);
 331
 332     for (unsigned int i=0;i<OutputImageDimension;i++) homOutputPoint[i]=oPoint[i];
 333     homOutputPoint[3]=1;
 334
 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];
 339
 340     typedef itk::LinearInterpolateImageFunction< InputImageType, double > InterpolatorType;
 341     typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
 342     interpolator->SetInputImage(this->GetInput());
 343
 344     //Run over all output voxels
 345     for (unsigned int i=0; i<outputSizeForThread[2]; i++)
 346       {
 347         for (unsigned int j=0; j<outputSizeForThread[1]; j++)
 348           {
 349             for (unsigned int k=0; k<outputSizeForThread[0]; k++)
 350               {
 351                 iPoint[0]=-homInputPoint[0]/homInputPoint[2] + m_PanelShift[0];
 352                 iPoint[1]=homInputPoint[1]/homInputPoint[2] + m_PanelShift[1];
 353
 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
 358                 else
 359                   value=m_EdgePaddingValue;
 360                 //Set it
 361                 iterator.Set(value);
 362
 363                 //Advance one step in the line
 364                 homInputPoint+=m_LineInc;
 365                 ++iterator;
 366               }
 367
 368             //Advance one line
 369             homInputPoint+=m_ColInc;
 370           }
 371
 372         //Advance one plane
 373         homInputPoint+=m_PlaneInc;
 374       }
 375   }
 376
 377
 378 } // namespace clitk
 379
 380
 381 #endif