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   ::Initialize( void ) throw (itk::ExceptionObject)
 194   {
 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();
 201
 202     //Change the origin
 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);
 207
 208     // Calculate the projection matrix
 209     this->CalculateProjectionMatrix();
 210
 211     // Calculate the coordinate increments
 212     this->CalculateCoordinateIncrements();
 213
 214     m_IsInitialized=true;
 215   }
 216
 217
 218   //-----------------------------------------------------------------------
 219   //     Calculate Projection Matrix
 220   //-----------------------------------------------------------------------
 221   template <class InputImageType, class OutputImageType>
 222   void
 223   BackProjectImageFilter<InputImageType, OutputImageType>
 224   ::CalculateProjectionMatrix( void )
 225   {
 226     //InputSpacingType inputSpacing=this->GetInput()->GetSpacing();
 227
 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];
 233     temp(2,0)=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;
 239
 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;
 246
 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;
 251
 252     // Calculate the centered rotation matrix
 253     itk::Matrix<double,4,4> centeredRotationMatrix=GetCenteredRotationMatrix3D(rotationParameters,transformedIsoCenter);
 254
 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];
 258
 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;
 263
 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++)
 267         {
 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);
 271         }
 272   }
 273
 274
 275   //-----------------------------------------------------------------------
 276   //     Calculate the coordinate increments
 277   //-----------------------------------------------------------------------
 278   template <class InputImageType, class OutputImageType>
 279   void
 280   BackProjectImageFilter<InputImageType, OutputImageType>
 281   ::CalculateCoordinateIncrements( void )
 282   {
 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];
 287
 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];
 292
 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];
 297   }
 298
 299
 300   //-----------------------------------------------------------------------
 301   //     Threaded generate data
 302   //-----------------------------------------------------------------------
 303   template <class InputImageType, class OutputImageType>
 304   void
 305   BackProjectImageFilter<InputImageType, OutputImageType>::ThreadedGenerateData( const OutputImageRegionType & outputRegionForThread,  itk::ThreadIdType threadId )
 306   {
 307     //Projection pointer
 308     InputImageConstPointer inputPtr=this->GetInput();
 309
 310     //Volume pointer
 311     OutputImagePointer outputPTr= this->GetOutput();
 312
 313     //Iterator for the thread region
 314     typedef itk::ImageRegionIterator<OutputImageType> OutputIteratorType;
 315     OutputIteratorType iterator(outputPTr, outputRegionForThread);
 316     OutputSizeType outputSizeForThread=outputRegionForThread.GetSize();
 317
 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
 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] + m_PanelShift[0];
 338     iPoint[1]=homInputPoint[1]/homInputPoint[2] + m_PanelShift[1];
 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