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