tools/clitkWarpImageGenericFilter.txx

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