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