tools/clitkWarpImageGenericFilter.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 clitkWarpImageGenericFilter_txx
  19 #define clitkWarpImageGenericFilter_txx
  20
  21 /* =================================================
  22  * @file   clitkWarpImageGenericFilter.txx
  23  * @author
  24  * @date
  25  *
  26  * @brief
  27  *
  28  ===================================================*/
  29
  30 #include "itkVectorResampleImageFilter.h"
  31 #include "clitkCoeffsToDVF.h"
  32
  33 namespace clitk
  34 {
  35
  36 //-------------------------------------------------------------------
  37 // Update with the number of dimensions
  38 //-------------------------------------------------------------------
  39 template<unsigned int Dimension>
  40 void
  41 WarpImageGenericFilter::UpdateWithDim(std::string PixelType)
  42 {
  43   if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D and "<< PixelType<<"..."<<std::endl;
  44
  45   if(PixelType == "short") {
  46     if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed short..." << std::endl;
  47     UpdateWithDimAndPixelType<Dimension, signed short>();
  48   }
  49   //    else if(PixelType == "unsigned_short"){
  50   //       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl;
  51   //       UpdateWithDimAndPixelType<Dimension, unsigned short>();
  52   //     }
  53
  54   else if (PixelType == "unsigned_char") {
  55     if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
  56     UpdateWithDimAndPixelType<Dimension, unsigned char>();
  57   }
  58
  59   //     else if (PixelType == "char"){
  60   //       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and signed_char..." << std::endl;
  61   //       UpdateWithDimAndPixelType<Dimension, signed char>();
  62   //     }
  63   else {
  64     if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
  65     UpdateWithDimAndPixelType<Dimension, float>();
  66   }
  67 }
  68
  69
  70 //-------------------------------------------------------------------
  71 // Update with the number of dimensions and the pixeltype
  72 //-------------------------------------------------------------------
  73 template <unsigned int Dimension, class  PixelType>
  74 void
  75 WarpImageGenericFilter::UpdateWithDimAndPixelType()
  76 {
  77
  78   // ImageTypes
  79   typedef itk::Image<PixelType, Dimension> InputImageType;
  80   typedef itk::Image<PixelType, Dimension> OutputImageType;
  81   typedef itk::Vector<float, Dimension> DisplacementType;
  82   typedef itk::Image<DisplacementType, Dimension> DeformationFieldType;
  83
  84   // Read the input
  85   typedef itk::ImageFileReader<InputImageType> InputReaderType;
  86   typename InputReaderType::Pointer reader = InputReaderType::New();
  87   reader->SetFileName( m_InputFileName);
  88   reader->Update();
  89   typename InputImageType::Pointer input= reader->GetOutput();
  90
  91   typename DeformationFieldType::Pointer deformationField;
  92   if (m_ArgsInfo.coeff_given) {
  93     deformationField = BLUTCoeffsToDVF<DeformationFieldType>(m_ArgsInfo.coeff_arg, m_InputFileName, m_Verbose);
  94   }
  95   else {
  96     //Read the deformation field
  97     typedef itk::ImageFileReader<DeformationFieldType> DeformationFieldReaderType;
  98     typename  DeformationFieldReaderType::Pointer deformationFieldReader= DeformationFieldReaderType::New();
  99     deformationFieldReader->SetFileName(m_ArgsInfo.vf_arg);
 100     deformationFieldReader->Update();
 101     deformationField =deformationFieldReader->GetOutput();
 102   }
 103
 104   // Intensity interpolator
 105   typedef clitk::GenericVectorInterpolator<args_info_clitkWarpImage, DeformationFieldType, double> GenericInterpolatorType;
 106   typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
 107   genericInterpolator->SetArgsInfo(m_ArgsInfo);
 108
 109
 110   // -------------------------------------------
 111   // Spacing like DVF
 112   // -------------------------------------------
 113   if (m_ArgsInfo.spacing_arg == 0) {
 114     // Calculate the region
 115     typename DeformationFieldType::SizeType newSize;
 116     for (unsigned int i=0 ; i <Dimension; i++)
 117       newSize[i]=(unsigned int) (input->GetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i]/deformationField->GetSpacing()[i]);
 118
 119     // Get the interpolator
 120     typedef clitk::GenericVectorInterpolator<args_info_clitkWarpImage, DeformationFieldType, double> GenericInterpolatorType;
 121     typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
 122     genericInterpolator->SetArgsInfo(m_ArgsInfo);
 123
 124     // Resample to match the extent of the input
 125     typename itk::VectorResampleImageFilter<DeformationFieldType, DeformationFieldType >::Pointer
 126     resampler =itk::VectorResampleImageFilter<DeformationFieldType, DeformationFieldType >::New();
 127     resampler->SetInput(deformationField);
 128     resampler->SetOutputSpacing(deformationField->GetSpacing());
 129     resampler->SetSize(newSize);
 130     resampler->SetOutputOrigin(input->GetOrigin());
 131     resampler->SetInterpolator(genericInterpolator->GetInterpolatorPointer());
 132
 133     // Update
 134     if (m_Verbose) std::cout<< "Resampling the VF..." <<std::endl;
 135     try {
 136       resampler->Update();
 137     } catch( itk::ExceptionObject & excp ) {
 138       std::cerr << "Problem resampling the input vector field file" << std::endl;
 139       std::cerr << excp << std::endl;
 140       return;
 141     }
 142     deformationField= resampler->GetOutput();
 143
 144   }
 145
 146   // -------------------------------------------
 147   // Spacing like input
 148   // -------------------------------------------
 149   else if (!m_ArgsInfo.coeff_given) {
 150     // Get size
 151     typename DeformationFieldType::SizeType newSize;
 152     for (unsigned int i=0 ; i <Dimension; i++)
 153       newSize[i]=input->GetLargestPossibleRegion().GetSize()[i];
 154
 155     // Resample to match the extent of the input
 156     typename itk::VectorResampleImageFilter<DeformationFieldType, DeformationFieldType >::Pointer
 157     resampler =itk::VectorResampleImageFilter<DeformationFieldType, DeformationFieldType >::New();
 158     resampler->SetInput(deformationField);
 159     resampler->SetOutputSpacing(input->GetSpacing());
 160     resampler->SetSize(newSize);
 161     resampler->SetOutputOrigin(input->GetOrigin());
 162     resampler->SetInterpolator(genericInterpolator->GetInterpolatorPointer());
 163
 164     // Update
 165     if (m_Verbose) std::cout<< "Resampling the VF..." <<std::endl;
 166     try {
 167       resampler->Update();
 168     } catch( itk::ExceptionObject & excp ) {
 169       std::cerr << "Problem resampling the input vector field file" << std::endl;
 170       std::cerr << excp << std::endl;
 171       return;
 172     }
 173     deformationField= resampler->GetOutput();
 174   }
 175
 176
 177   // -------------------------------------------
 178   // Forward Warp
 179   // -------------------------------------------
 180   typename    itk::ImageToImageFilter<InputImageType, InputImageType>::Pointer warpFilter;
 181   if (m_ArgsInfo.forward_flag) {
 182     //Forward warping: always linear
 183     typedef clitk::ForwardWarpImageFilter<InputImageType, InputImageType, DeformationFieldType> ForwardWarpFilterType;
 184     typename ForwardWarpFilterType::Pointer forwardWarpFilter= ForwardWarpFilterType::New();
 185     forwardWarpFilter->SetDeformationField( deformationField );
 186     forwardWarpFilter->SetEdgePaddingValue( static_cast<PixelType>(m_ArgsInfo.pad_arg) );
 187     warpFilter=forwardWarpFilter;
 188   }
 189
 190   // -------------------------------------------
 191   // Backward Warp
 192   // -------------------------------------------
 193   else {
 194     // Get the interpolator
 195     typedef clitk::GenericInterpolator<args_info_clitkWarpImage, InputImageType, double> GenericInterpolatorType;
 196     typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
 197     genericInterpolator->SetArgsInfo(m_ArgsInfo);
 198
 199     //Backward mapping
 200     typedef itk::WarpImageFilter<InputImageType, InputImageType, DeformationFieldType> BackwardWarpFilterType;
 201     typename BackwardWarpFilterType::Pointer backwardWarpFilter= BackwardWarpFilterType::New();
 202 #if ITK_VERSION_MAJOR >= 4
 203     backwardWarpFilter->SetDisplacementField( deformationField );
 204 #else
 205     backwardWarpFilter->SetDeformationField( deformationField );
 206 #endif
 207     backwardWarpFilter->SetEdgePaddingValue( static_cast<PixelType>(m_ArgsInfo.pad_arg) );
 208     backwardWarpFilter->SetOutputSpacing( deformationField->GetSpacing() );
 209     backwardWarpFilter->SetOutputOrigin( input->GetOrigin() );
 210     backwardWarpFilter->SetOutputSize( deformationField->GetLargestPossibleRegion().GetSize() );
 211     typename itk::VectorResampleImageFilter<DeformationFieldType, DeformationFieldType >::Pointer
 212     resampler =itk::VectorResampleImageFilter<DeformationFieldType, DeformationFieldType >::New();
 213     backwardWarpFilter->SetInterpolator(genericInterpolator->GetInterpolatorPointer());
 214     warpFilter=backwardWarpFilter;
 215   }
 216
 217
 218   // -------------------------------------------
 219   // Update
 220   // -------------------------------------------
 221   warpFilter->SetInput(input);
 222   if (m_Verbose) std::cout<< "Warping the input..." <<std::endl;
 223   try {
 224     warpFilter->Update();
 225   } catch( itk::ExceptionObject & excp ) {
 226     std::cerr << "Problem warping the input image" << std::endl;
 227     std::cerr << excp << std::endl;
 228     return;
 229   }
 230
 231
 232   // Output
 233   typedef itk::ImageFileWriter<OutputImageType> WriterType;
 234   typename WriterType::Pointer writer = WriterType::New();
 235   writer->SetFileName(m_ArgsInfo.output_arg);
 236   writer->SetInput(warpFilter->GetOutput());
 237   writer->Update();
 238
 239 }
 240
 241
 242 }//end clitk
 243
 244 #endif //#define clitkWarpImageGenericFilter_txx