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
  24 namespace clitk
  25 {
  26
  27   //-----------------------------------------------------------------------
  28   //     Constructor
  29   //-----------------------------------------------------------------------
  30
  31   template<class InputImageType, class OutputImageType>
  32   BackProjectImageFilter< InputImageType, OutputImageType >
  33   ::BackProjectImageFilter()
  34   {
  35
  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();
  42
  43     //Parameters for output
  44     this->m_OutputSpacing.Fill(1);
  45     this->m_OutputOrigin.Fill(0.0);
  46     this->m_OutputSize.Fill( 512 );
  47
  48     this->m_IsInitialized=false;
  49   }
  50
  51
  52   //-----------------------------------------------------------------------
  53   //     PrintSelf
  54   //-----------------------------------------------------------------------
  55   template<class InputImageType, class OutputImageType>
  56   void
  57   BackProjectImageFilter< InputImageType, OutputImageType >
  58   ::PrintSelf(std::ostream& os, itk::Indent indent) const
  59   {
  60     this->Superclass::PrintSelf(os,indent);
  61
  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;
  67   }
  68
  69
  70   //-----------------------------------------------------------------------
  71   //     Set output info from image
  72   //-----------------------------------------------------------------------
  73   template <class InputImageType, class OutputImageType>
  74   void
  75   BackProjectImageFilter<InputImageType, OutputImageType>
  76   ::SetOutputParametersFromImage ( const OutputImagePointer image )
  77   {
  78     this->SetOutputOrigin( image->GetOrigin() );
  79     this->SetOutputSpacing( image->GetSpacing() );
  80     this->SetOutputSize( image->GetLargestPossibleRegion().GetSize() );
  81   }
  82
  83
  84   //-----------------------------------------------------------------------
  85   //     Set output info from image
  86   //-----------------------------------------------------------------------
  87   template <class InputImageType, class OutputImageType>
  88   void
  89   BackProjectImageFilter<InputImageType, OutputImageType>
  90   ::SetOutputParametersFromImage (const OutputImageConstPointer image )
  91   {
  92     this->SetOutputOrigin( image->GetOrigin() );
  93     this->SetOutputSpacing( image->GetSpacing() );
  94     this->SetOutputSize( image->GetLargestPossibleRegion().GetSize() );
  95   }
  96
  97
  98   //-----------------------------------------------------------------------
  99   //     Generate output information
 100   //-----------------------------------------------------------------------
 101   template<class InputImageType, class OutputImageType>
 102   void
 103   BackProjectImageFilter< InputImageType, OutputImageType >
 104   ::GenerateOutputInformation( void )
 105   {
 106     // Don not call the superclass' implementation of this method (!=Dimensions)
 107     // Superclass::GenerateOutputInformation();
 108
 109     // get pointer to the output
 110     OutputImagePointer outputPtr = this->GetOutput();
 111     InputImageConstPointer  inputPtr  = this->GetInput();
 112
 113     if ( !outputPtr || ! inputPtr)
 114       {
 115         return;
 116       }
 117
 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());
 123
 124
 125     //OutputImageRegionType region;
 126     outputLargestPossibleRegion.SetSize( m_OutputSize );
 127     outputPtr->SetLargestPossibleRegion( outputLargestPossibleRegion );
 128     outputPtr->SetSpacing( m_OutputSpacing );
 129     outputPtr->SetOrigin( m_OutputOrigin );
 130
 131   }
 132
 133
 134   //-----------------------------------------------------------------------
 135   //     Generate input Requested region
 136   //-----------------------------------------------------------------------
 137   template <class InputImageType, class OutputImageType>
 138   void
 139   BackProjectImageFilter<InputImageType,OutputImageType>
 140   ::GenerateInputRequestedRegion()
 141   {
 142     //Call superclass
 143     Superclass::GenerateInputRequestedRegion();
 144
 145     if (!this->GetOutput())
 146       {
 147         return;
 148       }
 149     OutputImageRegionType outputRegion = this->GetOutput()->GetRequestedRegion();
 150     InputImagePointer inputPtr =  const_cast<InputImageType *>(this->GetInput());
 151
 152     if ( !inputPtr )
 153       {
 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());
 161         throw e;
 162       }
 163
 164     InputImageRegionType inputRegion;
 165     inputRegion=inputPtr->GetLargestPossibleRegion();
 166     inputPtr->SetRequestedRegion(inputRegion);
 167
 168   }
 169
 170
 171
 172   //-----------------------------------------------------------------------
 173   //     Before threaded data
 174   //-----------------------------------------------------------------------
 175   template <class InputImageType, class OutputImageType>
 176   void
 177   BackProjectImageFilter<InputImageType, OutputImageType>
 178   ::BeforeThreadedGenerateData( void )
 179   {
 180     if (!m_IsInitialized) this->Initialize();
 181   }
 182
 183
 184  //-----------------------------------------------------------------------
 185   //     Initialize
 186   //-----------------------------------------------------------------------
 187   template <class InputImageType, class OutputImageType>
 188   void
 189   BackProjectImageFilter<InputImageType, OutputImageType>
 190   ::Initialize( void ) throw (itk::ExceptionObject)
 191   {
 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();
 198
 199     //Change the origin
 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);
 204
 205     // Calculate the projection matrix
 206     this->CalculateProjectionMatrix();
 207
 208     // Calculate the coordinate increments
 209     this->CalculateCoordinateIncrements();
 210
 211     m_IsInitialized=true;
 212   }
 213
 214
 215   //-----------------------------------------------------------------------
 216   //     Calculate Projection Matrix
 217   //-----------------------------------------------------------------------
 218   template <class InputImageType, class OutputImageType>
 219   void
 220   BackProjectImageFilter<InputImageType, OutputImageType>
 221   ::CalculateProjectionMatrix( void )
 222   {
 223     InputSpacingType inputSpacing=this->GetInput()->GetSpacing();
 224
 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];
 230     temp(2,0)=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;
 236
 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;
 243
 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;
 248
 249     // Calculate the centered rotation matrix
 250     itk::Matrix<double,4,4> centeredRotationMatrix=GetCenteredRotationMatrix3D(rotationParameters,transformedIsoCenter);
 251
 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];
 255
 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;
 260
 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++)
 264         {
 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);
 268         }
 269   }
 270
 271
 272   //-----------------------------------------------------------------------
 273   //     Calculate the coordinate increments
 274   //-----------------------------------------------------------------------
 275   template <class InputImageType, class OutputImageType>
 276   void
 277   BackProjectImageFilter<InputImageType, OutputImageType>
 278   ::CalculateCoordinateIncrements( void )
 279   {
 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];
 284
 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];
 289
 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];
 294   }
 295
 296
 297   //-----------------------------------------------------------------------
 298   //     Threaded generate data
 299   //-----------------------------------------------------------------------
 300   template <class InputImageType, class OutputImageType>
 301   void
 302   BackProjectImageFilter<InputImageType, OutputImageType>
 303   ::ThreadedGenerateData( const OutputImageRegionType & outputRegionForThread,  int threadId )
 304   {
 305     //Projection pointer
 306     InputImageConstPointer inputPtr=this->GetInput();
 307     InputPixelType * beginPtr=const_cast<InputPixelType *>(this->GetInput()->GetBufferPointer());
 308     InputPixelType * pp;
 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     InputSizeType inputSize=inputPtr->GetLargestPossibleRegion().GetSize();
 328     double dx,dy,dxm,dym;
 329     int lx, ly;
 330
 331     //Get the first output coordinate
 332     oIndex=iterator.GetIndex();//costly but only once a thread
 333     outputPTr->TransformIndexToPhysicalPoint(oIndex,oPoint);
 334
 335     for (unsigned int i=0;i<OutputImageDimension;i++) homOutputPoint[i]=oPoint[i];
 336     homOutputPoint[3]=1;
 337
 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];
 342
 343     //Run over all output voxels
 344     for (unsigned int i=0; i<outputSizeForThread[2]; i++)
 345       {
 346         for (unsigned int j=0; j<outputSizeForThread[1]; j++)
 347           {
 348             for (unsigned int k=0; k<outputSizeForThread[0]; k++)
 349               {
 350                 iPoint[0]=homInputPoint[0]/homInputPoint[2];
 351                 iPoint[1]=homInputPoint[1]/homInputPoint[2];
 352
 353                 //Check wether inside, convert to index (use modified with correct origin)
 354                 if( m_ModifiedInput->TransformPhysicalPointToContinuousIndex(iPoint, iIndex) )
 355                   {
 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))) );
 364
 365                   }
 366                 //Outside: padding value
 367                 else value=m_EdgePaddingValue;
 368
 369                 //Set it
 370                 iterator.Set(value);
 371
 372                 //Advance one step in the line
 373                 homInputPoint+=m_LineInc;
 374                 ++iterator;
 375               }
 376
 377             //Advance one line
 378             homInputPoint+=m_ColInc;
 379           }
 380
 381         //Advance one plane
 382         homInputPoint+=m_PlaneInc;
 383       }
 384   }
 385
 386
 387 } // namespace clitk
 388
 389
 390 #endif