tools/clitkConeBeamProjectImageFilter.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://oncora1.lyon.fnclcc.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 __clitkConeBeamProjectImageFilter_txx
  19 #define __clitkConeBeamProjectImageFilter_txx
  20 namespace clitk
  21 {
  22
  23
  24   //=========================================================================================================================
  25   //Constructor
  26   template <class InputImageType, class OutputImageType>
  27   ConeBeamProjectImageFilter<InputImageType, OutputImageType>::
  28   ConeBeamProjectImageFilter()
  29   {
  30     // Initialize with the default values (Elekta Synergie)
  31     m_Verbose=false;
  32     m_IsInitialized=false;
  33     m_NumberOfThreadsIsGiven=false;
  34
  35     // Geometry
  36     m_IsoCenter.Fill(0.0);
  37     m_SourceToScreen=1536.;
  38     m_SourceToAxis=1000.;
  39     m_ProjectionAngle=0.;
  40     m_RigidTransformMatrix.SetIdentity();
  41     m_EdgePaddingValue=itk::NumericTraits<OutputPixelType>::Zero;//density images
  42
  43     // Pipe
  44     m_Transform=TransformType::New();
  45     m_Resampler=ResampleImageFilterType::New();
  46     m_Interpolator= InterpolatorType::New();
  47     m_ExtractImageFilter=ExtractImageFilterType::New();
  48     m_FlipImageFilter=FlipImageFilterType::New();
  49
  50     // Parameters for output
  51     this->m_OutputSpacing.Fill(0.8);
  52     this->m_OutputOrigin.Fill(0.0);
  53     this->m_OutputSize.Fill( 512 );
  54
  55   }
  56
  57
  58   //-----------------------------------------------------------------------
  59   //     Set output info from image
  60   //-----------------------------------------------------------------------
  61   template <class InputImageType, class OutputImageType>
  62   void
  63   ConeBeamProjectImageFilter<InputImageType, OutputImageType>
  64   ::SetOutputParametersFromImage ( OutputImagePointer image )
  65   {
  66     this->SetOutputOrigin( image->GetOrigin() );
  67     this->SetOutputSpacing( image->GetSpacing() );
  68     this->SetOutputSize( image->GetLargestPossibleRegion().GetSize() );
  69   }
  70
  71
  72   //-----------------------------------------------------------------------
  73   //     Set output info from image
  74   //-----------------------------------------------------------------------
  75   template <class InputImageType, class OutputImageType>
  76   void
  77   ConeBeamProjectImageFilter<InputImageType, OutputImageType>
  78   ::SetOutputParametersFromImage ( OutputImageConstPointer image )
  79   {
  80     this->SetOutputOrigin( image->GetOrigin() );
  81     this->SetOutputSpacing( image->GetSpacing() );
  82     this->SetOutputSize( image->GetLargestPossibleRegion().GetSize() );
  83   }
  84
  85
  86   //=========================================================================================================================
  87   // Initialize
  88   template <class InputImageType, class OutputImageType>
  89   void
  90   ConeBeamProjectImageFilter<InputImageType, OutputImageType>
  91   ::Initialize(void) throw (itk::ExceptionObject)
  92   {
  93
  94     //====================================================================
  95     // Define the transform: composition of rigid transfo around origin  and a centered rotation
  96     // The center of rotation should be placed at the isocenter!
  97
  98     if (m_Verbose) std::cout<<"The isocenter is at "<< m_IsoCenter <<"..." <<std::endl;
  99
 100     // Get the rotation parameter array
 101     itk::Array<double> rotationParameters(3);
 102     const double dtr = ( atan(1.0) * 4.0 ) / 180.0; //constant for converting degrees into radians
 103     if (m_Verbose)std::cout<<"The projection angle is "<< m_ProjectionAngle <<"° (0° being lateral left)..."<< std::endl;
 104     rotationParameters[0]= 0.;
 105     rotationParameters[1]= 0.;
 106     rotationParameters[2]= -dtr*(double) m_ProjectionAngle;
 107
 108     // We first perform rigid transform (of source and panel), then a centered rotation around the transformed center
 109     itk::Matrix<double,3,3> rigidRotation = GetRotationalPartMatrix3D(m_RigidTransformMatrix);
 110     itk::Vector<double,3> rigidTranslation = GetTranslationPartMatrix3D(m_RigidTransformMatrix);
 111     typename InputImageType::PointType transformedCenter = rigidRotation * m_IsoCenter + rigidTranslation;
 112
 113     // Calculate the centered rotation matrix
 114     itk::Matrix<double,4,4> centeredRotationMatrix = GetCenteredRotationMatrix3D(rotationParameters,transformedCenter);
 115
 116     // Compose this rotation with the rigid transform matrix
 117     itk::Matrix<double,4,4> finalTransform = m_RigidTransformMatrix * centeredRotationMatrix ;
 118
 119     // Set the rotation
 120     itk::Matrix<double,3,3> finalRotation = GetRotationalPartMatrix3D(finalTransform);
 121     m_Transform->SetMatrix(finalRotation);
 122     if (m_Verbose)std::cout<<"The final rotation is "<<finalRotation<<"..."<<std::endl;
 123
 124     // Set the translation
 125     itk::Vector<double,3> finalTranslation = GetTranslationPartMatrix3D(finalTransform);
 126     m_Transform->SetTranslation(finalTranslation);
 127     if (m_Verbose)std::cout<<"The final translation is "<<finalTranslation<<"..."<<std::endl;
 128
 129
 130     //====================================================================
 131     //Define a resample filter for the projection
 132     if (m_NumberOfThreadsIsGiven) m_Resampler->SetNumberOfThreads(m_NumberOfThreads);
 133     m_Resampler->SetDefaultPixelValue(m_EdgePaddingValue);
 134     if (m_Verbose)std::cout<<"The edge padding value is "<<m_EdgePaddingValue<<"..."<<std::endl;
 135
 136     //JV Original raycast does not take into account origin, but does implicit shift to center of volume
 137     //JV patched in RayCastInterpolateImageFunctionWithOrigin
 138     m_Interpolator->SetTransform(m_Transform);
 139     m_Interpolator->SetThreshold(0.);
 140
 141     // Focalpoint: we presume that for an angle of 0° the kV is lateral left (for the patient on his back), or on the positive X-axis.
 142     typename  InterpolatorType::InputPointType focalpoint;
 143     focalpoint[0]= m_IsoCenter[0] + m_SourceToAxis;
 144     focalpoint[1]= m_IsoCenter[1];
 145     focalpoint[2]= m_IsoCenter[2];
 146     m_Interpolator->SetFocalPoint(focalpoint);
 147     if (m_Verbose)std::cout<<"The focalpoint is at "<< focalpoint <<"..."<< std::endl;
 148
 149     // Connect the interpolator and the transform with the m_Resampler
 150     m_Resampler->SetInterpolator( m_Interpolator );
 151     m_Resampler->SetTransform( m_Transform );
 152
 153     // Describe the output projection image
 154     typename  InputImageType::SizeType   sizeOuput;
 155     sizeOuput[0] = 1;    // number of pixels along X of the 2D DRR image
 156     sizeOuput[1] = m_OutputSize[0];  // number of pixels along Y of the 2D DRR image
 157     sizeOuput[2] = m_OutputSize[1];  // number of pixels along Z of the 2D DRR image
 158     m_Resampler->SetSize( sizeOuput );
 159     if (m_Verbose)std::cout<<"The output size is "<< m_OutputSize <<"..."<< std::endl;
 160
 161     typename  InputImageType::SpacingType spacingOutput;
 162     spacingOutput[0] = 1;    // pixel spacing along X of the 2D DRR image [mm]
 163     spacingOutput[1] = m_OutputSpacing[0];  // pixel spacing along Y of the 2D DRR image [mm]
 164     spacingOutput[2] = m_OutputSpacing[1];  // pixel spacing along Y of the 2D DRR image [mm]
 165     m_Resampler->SetOutputSpacing( spacingOutput );
 166     if (m_Verbose)std::cout<<"The output size is "<< m_OutputSpacing <<"..."<< std::endl;
 167
 168     // The position of the DRR is specified, we presume that for an angle of 0° the flatpanel is located at the negative x-axis
 169     // JV -1 seems to correspond better with shearwarp of Simon Rit
 170     typename  InterpolatorType::InputPointType originOutput;
 171     originOutput[0] = m_IsoCenter[0]- (m_SourceToScreen - m_SourceToAxis);
 172     originOutput[1] = m_IsoCenter[1]-static_cast<double>(sizeOuput[1]-1)*spacingOutput[1]/2.0 - m_PanelShift;
 173     originOutput[2] = m_IsoCenter[2]-static_cast<double>(sizeOuput[2]-1)*spacingOutput[2]/2.0;
 174     m_Resampler->SetOutputOrigin( originOutput );
 175     if (m_Verbose)std::cout<<"The origin of the flat panel is at "<< originOutput <<",..."<< std::endl;
 176
 177     // We define the region to be extracted. Projection was in the YZ plane, X should be set to zero
 178     typename InternalImageType::RegionType::SizeType sizeTemp=sizeOuput;
 179     sizeTemp[0]=0;
 180     typename InternalImageType::IndexType startTemp; //=temp->GetLargestPossibleRegion().GetIndex();
 181     startTemp.Fill(0);
 182
 183     typename InternalImageType::RegionType desiredRegion;
 184     desiredRegion.SetSize( sizeTemp );
 185     desiredRegion.SetIndex( startTemp );
 186     m_ExtractImageFilter->SetExtractionRegion( desiredRegion );
 187
 188     // Prepare the Flip Image filter
 189     typename  FlipImageFilterType::FlipAxesArrayType flipArray;
 190     flipArray[0] = 0;
 191     flipArray[1] = 1;
 192     m_FlipImageFilter->SetFlipAxes( flipArray );
 193
 194     // Initialization complete
 195     m_IsInitialized=true;
 196
 197   }
 198
 199   //=========================================================================================================================
 200   // Update
 201   template <class InputImageType, class OutputImageType>
 202   void
 203   ConeBeamProjectImageFilter<InputImageType, OutputImageType>
 204   ::Update(void)
 205   {
 206
 207     //==================================================================
 208     // If geometry changed reinitialize
 209     if (! m_IsInitialized) this->Initialize();
 210
 211     // Input
 212     m_Resampler->SetInput( m_Input );
 213
 214     //==================================================================
 215     // Execute the filter
 216     if (m_Verbose)std::cout<<"Starting the projection..."<<std::endl;
 217     try {
 218       m_Resampler->Update();
 219     }
 220     catch( itk::ExceptionObject & err )
 221       {
 222         std::cerr << "ExceptionObject caught! Projection failed!" << std::endl;
 223         std::cerr << err << std::endl;
 224       }
 225
 226     //==================================================================
 227     // Make a 2D image out of it
 228     typename InternalImageType::Pointer temp=m_Resampler->GetOutput();
 229     m_ExtractImageFilter->SetInput(temp);
 230
 231     // We should flip the Y axis
 232     m_FlipImageFilter->SetInput(m_ExtractImageFilter->GetOutput());
 233     m_FlipImageFilter->Update();
 234
 235     // Get output
 236     m_Output= m_FlipImageFilter->GetOutput();
 237     m_Output->Update();
 238     m_Output->SetOrigin(m_OutputOrigin);
 239   }
 240
 241   //=========================================================================================================================
 242   // GetOutput
 243   template <class InputImageType, class OutputImageType>
 244   typename OutputImageType::Pointer
 245   ConeBeamProjectImageFilter<InputImageType, OutputImageType>
 246   ::GetOutput(void)
 247   {
 248     //JV it is not safe to repeatedly call getoutput/ always call Update first
 249     return m_Output;
 250   }
 251
 252 }
 253
 254
 255 #endif