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