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->SetSize(newSize);
 136     resampler->SetOutputOrigin(input->GetOrigin());
 137     resampler->SetInterpolator(genericInterpolator->GetInterpolatorPointer());
 138
 139     // Update
 140     if (m_Verbose) std::cout<< "Resampling the VF..." <<std::endl;
 141     try {
 142       resampler->Update();
 143     } catch( itk::ExceptionObject & excp ) {
 144       std::cerr << "Problem resampling the input vector field file" << std::endl;
 145       std::cerr << excp << std::endl;
 146       return;
 147     }
 148     deformationField= resampler->GetOutput();
 149
 150   }
 151
 152   // -------------------------------------------
 153   // Spacing like input
 154   // -------------------------------------------
 155   else if (!m_ArgsInfo.coeff_given) {
 156     // Get size
 157     typename DeformationFieldType::SizeType newSize;
 158     for (unsigned int i=0 ; i <Dimension; i++)
 159       newSize[i]=input->GetLargestPossibleRegion().GetSize()[i];
 160
 161     // Resample to match the extent of the input
 162     typename itk::VectorResampleImageFilter<DeformationFieldType, DeformationFieldType >::Pointer
 163     resampler =itk::VectorResampleImageFilter<DeformationFieldType, DeformationFieldType >::New();
 164     resampler->SetInput(deformationField);
 165     resampler->SetOutputSpacing(input->GetSpacing());
 166     resampler->SetSize(newSize);
 167     resampler->SetOutputOrigin(input->GetOrigin());
 168     resampler->SetInterpolator(genericInterpolator->GetInterpolatorPointer());
 169
 170     // Update
 171     if (m_Verbose) std::cout<< "Resampling the VF..." <<std::endl;
 172     try {
 173       resampler->Update();
 174     } catch( itk::ExceptionObject & excp ) {
 175       std::cerr << "Problem resampling the input vector field file" << std::endl;
 176       std::cerr << excp << std::endl;
 177       return;
 178     }
 179     deformationField= resampler->GetOutput();
 180   }
 181
 182
 183   // -------------------------------------------
 184   // Forward Warp
 185   // -------------------------------------------
 186   typename    itk::ImageToImageFilter<InputImageType, InputImageType>::Pointer warpFilter;
 187   if (m_ArgsInfo.forward_flag) {
 188     //Forward warping: always linear
 189     typedef clitk::ForwardWarpImageFilter<InputImageType, InputImageType, DeformationFieldType> ForwardWarpFilterType;
 190     typename ForwardWarpFilterType::Pointer forwardWarpFilter= ForwardWarpFilterType::New();
 191     forwardWarpFilter->SetDeformationField( deformationField );
 192     forwardWarpFilter->SetEdgePaddingValue( static_cast<PixelType>(m_ArgsInfo.pad_arg) );
 193     warpFilter=forwardWarpFilter;
 194   }
 195
 196   // -------------------------------------------
 197   // Backward Warp
 198   // -------------------------------------------
 199   else {
 200     // Get the interpolator
 201     typedef clitk::GenericInterpolator<args_info_clitkWarpImage, InputImageType, double> GenericInterpolatorType;
 202     typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
 203     genericInterpolator->SetArgsInfo(m_ArgsInfo);
 204
 205     //Backward mapping
 206     typedef itk::WarpImageFilter<InputImageType, InputImageType, DeformationFieldType> BackwardWarpFilterType;
 207     typename BackwardWarpFilterType::Pointer backwardWarpFilter= BackwardWarpFilterType::New();
 208     backwardWarpFilter->SetDisplacementField( deformationField );
 209     backwardWarpFilter->SetEdgePaddingValue( static_cast<PixelType>(m_ArgsInfo.pad_arg) );
 210     backwardWarpFilter->SetOutputSpacing( deformationField->GetSpacing() );
 211     backwardWarpFilter->SetOutputOrigin( input->GetOrigin() );
 212     backwardWarpFilter->SetOutputSize( deformationField->GetLargestPossibleRegion().GetSize() );
 213     typename itk::VectorResampleImageFilter<DeformationFieldType, DeformationFieldType >::Pointer
 214     resampler =itk::VectorResampleImageFilter<DeformationFieldType, DeformationFieldType >::New();
 215     backwardWarpFilter->SetInterpolator(genericInterpolator->GetInterpolatorPointer());
 216     warpFilter=backwardWarpFilter;
 217   }
 218
 219
 220   // -------------------------------------------
 221   // Update
 222   // -------------------------------------------
 223   warpFilter->SetInput(input);
 224   if (m_Verbose) std::cout<< "Warping the input..." <<std::endl;
 225   try {
 226     warpFilter->Update();
 227   } catch( itk::ExceptionObject & excp ) {
 228     std::cerr << "Problem warping the input image" << std::endl;
 229     std::cerr << excp << std::endl;
 230     return;
 231   }
 232
 233
 234   // Output
 235   typedef itk::ImageFileWriter<OutputImageType> WriterType;
 236   typename WriterType::Pointer writer = WriterType::New();
 237   writer->SetFileName(m_ArgsInfo.output_arg);
 238   writer->SetInput(warpFilter->GetOutput());
 239   writer->Update();
 240
 241 }
 242
 243
 244 }//end clitk
 245
 246 #endif //#define clitkWarpImageGenericFilter_txx