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