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 "clitkConvertBLUTCoeffsToVFFilter.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     typedef ConvertBLUTCoeffsToVFFilter<DeformationFieldType> FilterType;
  94     typename FilterType::Pointer filter = FilterType::New();
  95     filter->SetInputFileName(m_ArgsInfo.coeff_arg);
  96     filter->SetLikeFileName(m_InputFileName);
  97     filter->SetVerbose(m_Verbose);
  98     filter->Update();
  99     deformationField = filter->GetOutput();
 100   }
 101   else {
 102     //Read the deformation field
 103     typedef itk::ImageFileReader<DeformationFieldType> DeformationFieldReaderType;
 104     typename  DeformationFieldReaderType::Pointer deformationFieldReader= DeformationFieldReaderType::New();
 105     deformationFieldReader->SetFileName(m_ArgsInfo.vf_arg);
 106     deformationFieldReader->Update();
 107     deformationField =deformationFieldReader->GetOutput();
 108   }
 109
 110   // Intensity interpolator
 111   typedef clitk::GenericVectorInterpolator<args_info_clitkWarpImage, DeformationFieldType, double> GenericInterpolatorType;
 112   typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
 113   genericInterpolator->SetArgsInfo(m_ArgsInfo);
 114
 115
 116   // -------------------------------------------
 117   // Spacing like DVF
 118   // -------------------------------------------
 119   if (m_ArgsInfo.spacing_arg == 0) {
 120     // Calculate the region
 121     typename DeformationFieldType::SizeType newSize;
 122     for (unsigned int i=0 ; i <Dimension; i++)
 123       newSize[i]=(unsigned int) (input->GetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i]/deformationField->GetSpacing()[i]);
 124
 125     // Get the interpolator
 126     typedef clitk::GenericVectorInterpolator<args_info_clitkWarpImage, DeformationFieldType, double> GenericInterpolatorType;
 127     typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
 128     genericInterpolator->SetArgsInfo(m_ArgsInfo);
 129
 130     // Resample to match the extent of the input
 131     typename itk::VectorResampleImageFilter<DeformationFieldType, DeformationFieldType >::Pointer
 132     resampler =itk::VectorResampleImageFilter<DeformationFieldType, DeformationFieldType >::New();
 133     resampler->SetInput(deformationField);
 134     resampler->SetOutputSpacing(deformationField->GetSpacing());
 135     resampler->SetOutputDirection(deformationField->GetDirection());
 136     resampler->SetSize(newSize);
 137     resampler->SetOutputOrigin(input->GetOrigin());
 138     resampler->SetInterpolator(genericInterpolator->GetInterpolatorPointer());
 139
 140     // Update
 141     if (m_Verbose) std::cout<< "Resampling the VF..." <<std::endl;
 142     try {
 143       resampler->Update();
 144     } catch( itk::ExceptionObject & excp ) {
 145       std::cerr << "Problem resampling the input vector field file" << std::endl;
 146       std::cerr << excp << std::endl;
 147       return;
 148     }
 149     deformationField= resampler->GetOutput();
 150
 151   }
 152
 153   // -------------------------------------------
 154   // Spacing like input
 155   // -------------------------------------------
 156   else if (!m_ArgsInfo.coeff_given) {
 157     // Get size
 158     typename DeformationFieldType::SizeType newSize;
 159     for (unsigned int i=0 ; i <Dimension; i++)
 160       newSize[i]=input->GetLargestPossibleRegion().GetSize()[i];
 161
 162     // Resample to match the extent of the input
 163     typename itk::VectorResampleImageFilter<DeformationFieldType, DeformationFieldType >::Pointer
 164     resampler =itk::VectorResampleImageFilter<DeformationFieldType, DeformationFieldType >::New();
 165     resampler->SetInput(deformationField);
 166     resampler->SetOutputSpacing(input->GetSpacing());
 167     resampler->SetOutputDirection(deformationField->GetDirection());
 168     resampler->SetSize(newSize);
 169     resampler->SetOutputOrigin(input->GetOrigin());
 170     resampler->SetInterpolator(genericInterpolator->GetInterpolatorPointer());
 171
 172     // Update
 173     if (m_Verbose) std::cout<< "Resampling the VF..." <<std::endl;
 174     try {
 175       resampler->Update();
 176     } catch( itk::ExceptionObject & excp ) {
 177       std::cerr << "Problem resampling the input vector field file" << std::endl;
 178       std::cerr << excp << std::endl;
 179       return;
 180     }
 181     deformationField= resampler->GetOutput();
 182   }
 183
 184
 185   // -------------------------------------------
 186   // Forward Warp
 187   // -------------------------------------------
 188   typename    itk::ImageToImageFilter<InputImageType, InputImageType>::Pointer warpFilter;
 189   if (m_ArgsInfo.forward_flag) {
 190     //Forward warping: always linear
 191     typedef clitk::ForwardWarpImageFilter<InputImageType, InputImageType, DeformationFieldType> ForwardWarpFilterType;
 192     typename ForwardWarpFilterType::Pointer forwardWarpFilter= ForwardWarpFilterType::New();
 193     forwardWarpFilter->SetDeformationField( deformationField );
 194     forwardWarpFilter->SetEdgePaddingValue( static_cast<PixelType>(m_ArgsInfo.pad_arg) );
 195     warpFilter=forwardWarpFilter;
 196   }
 197
 198   // -------------------------------------------
 199   // Backward Warp
 200   // -------------------------------------------
 201   else {
 202     // Get the interpolator
 203     typedef clitk::GenericInterpolator<args_info_clitkWarpImage, InputImageType, double> GenericInterpolatorType;
 204     typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
 205     genericInterpolator->SetArgsInfo(m_ArgsInfo);
 206
 207     //Backward mapping
 208     typedef itk::WarpImageFilter<InputImageType, InputImageType, DeformationFieldType> BackwardWarpFilterType;
 209     typename BackwardWarpFilterType::Pointer backwardWarpFilter= BackwardWarpFilterType::New();
 210     backwardWarpFilter->SetDisplacementField( deformationField );
 211     backwardWarpFilter->SetEdgePaddingValue( static_cast<PixelType>(m_ArgsInfo.pad_arg) );
 212     backwardWarpFilter->SetOutputSpacing( deformationField->GetSpacing() );
 213     backwardWarpFilter->SetOutputOrigin( input->GetOrigin() );
 214     backwardWarpFilter->SetOutputSize( deformationField->GetLargestPossibleRegion().GetSize() );
 215     backwardWarpFilter->SetOutputDirection( input->GetDirection() );
 216     typename itk::VectorResampleImageFilter<DeformationFieldType, DeformationFieldType >::Pointer
 217     resampler =itk::VectorResampleImageFilter<DeformationFieldType, DeformationFieldType >::New();
 218     backwardWarpFilter->SetInterpolator(genericInterpolator->GetInterpolatorPointer());
 219     warpFilter=backwardWarpFilter;
 220   }
 221
 222
 223   // -------------------------------------------
 224   // Update
 225   // -------------------------------------------
 226   warpFilter->SetInput(input);
 227   if (m_Verbose) std::cout<< "Warping the input..." <<std::endl;
 228   try {
 229     warpFilter->Update();
 230   } catch( itk::ExceptionObject & excp ) {
 231     std::cerr << "Problem warping the input image" << std::endl;
 232     std::cerr << excp << std::endl;
 233     return;
 234   }
 235
 236
 237   // Output
 238   typedef itk::ImageFileWriter<OutputImageType> WriterType;
 239   typename WriterType::Pointer writer = WriterType::New();
 240   writer->SetFileName(m_ArgsInfo.output_arg);
 241   writer->SetInput(warpFilter->GetOutput());
 242   writer->Update();
 243
 244 }
 245
 246
 247 }//end clitk
 248
 249 #endif //#define clitkWarpImageGenericFilter_txx