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