tools/clitkVFResampleGenericFilter.cxx

   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 CLITKVFRESAMPLEGENERICFILTER_CXX
  19 #define CLITKVFRESAMPLEGENERICFILTER_CXX
  20
  21 #include "clitkVFResampleGenericFilter.h"
  22
  23 //--------------------------------------------------------------------
  24 clitk::VFResampleGenericFilter::VFResampleGenericFilter():
  25   clitk::ImageToImageGenericFilter<Self>("VFResample")
  26 {
  27   InitializeImageType<2>();
  28   InitializeImageType<3>();
  29   //  InitializeImageType<4>();
  30   mApplyGaussianFilterBefore = false;
  31   mDefaultPixelValue = 0.0;
  32   mInterpolatorName = "NN";
  33   mBSplineOrder=3;
  34 }
  35 //--------------------------------------------------------------------
  36
  37
  38 //--------------------------------------------------------------------
  39 template<unsigned int Dim>
  40 void clitk::VFResampleGenericFilter::InitializeImageType()
  41 {
  42   //typedef itk::Vector<float,Dim> v3f;
  43   //ADD_IMAGE_TYPE(Dim, v3f);
  44   ADD_DEFAULT_VEC_IMAGE_TYPES
  45 }
  46 //--------------------------------------------------------------------
  47
  48
  49 //--------------------------------------------------------------------
  50 template<class ImageType>
  51 void clitk::VFResampleGenericFilter::UpdateWithInputImageType()
  52 {
  53
  54   if (m_NbOfComponents == 1) {
  55     std::cerr << "Error, only one components ? Use clitkImageResample instead." << std::endl;
  56     exit(0);
  57   }
  58   typedef typename ImageType::PixelType PixelType;
  59   if (m_NbOfComponents == 2) Update_WithDimAndPixelTypeAndComponent<ImageType::ImageDimension,PixelType,2>();
  60   if (m_NbOfComponents == 3) Update_WithDimAndPixelTypeAndComponent<ImageType::ImageDimension,PixelType,3>();
  61   if (m_NbOfComponents == 4) Update_WithDimAndPixelTypeAndComponent<ImageType::ImageDimension,PixelType,4>();
  62 }
  63 //--------------------------------------------------------------------
  64
  65
  66 //--------------------------------------------------------------------
  67 template<unsigned int Dim, class PixelType, unsigned int DimCompo>
  68 void clitk::VFResampleGenericFilter::Update_WithDimAndPixelTypeAndComponent()
  69 {
  70   // Reading input
  71   //  typedef itk::Vector<PixelType, DimCompo> DisplacementType;
  72   typedef PixelType DisplacementType;
  73   typedef itk::Image< DisplacementType, Dim > ImageType;
  74
  75   typename ImageType::Pointer input = clitk::readImage<ImageType>(m_InputFilenames, m_IOVerbose);
  76
  77   // Main filter
  78   typename ImageType::Pointer outputImage = ComputeImage<ImageType>(input);
  79
  80   // Write results
  81   SetNextOutput<ImageType>(outputImage);
  82 }
  83 //--------------------------------------------------------------------
  84
  85 //--------------------------------------------------------------------
  86 template<class ImageType>
  87 typename ImageType::Pointer
  88 clitk::VFResampleGenericFilter::ComputeImage(typename ImageType::Pointer inputImage)
  89 {
  90
  91   // Check options
  92   static unsigned int dim = ImageType::ImageDimension;
  93   if (mOutputSize.size() != dim) {
  94     std::cerr << "Please set size with " << dim << " dimensions." << std::endl;
  95     return ITK_NULLPTR;
  96   }
  97   if (mOutputSpacing.size() != dim) {
  98     std::cerr << "Please set spacing with " << dim << " dimensions." << std::endl;
  99     return ITK_NULLPTR;
 100   }
 101   mOutputOrigin.resize(dim);
 102
 103   if (mApplyGaussianFilterBefore && mSigma.size() != dim) {
 104     std::cerr << "Please set sigma with " << dim << " dimensions." << std::endl;
 105     return ITK_NULLPTR;
 106   }
 107
 108   // Some typedefs
 109   typedef typename ImageType::SizeType SizeType;
 110   typedef typename ImageType::SpacingType SpacingType;
 111   typedef typename ImageType::PointType PointType;
 112
 113   // Create Image Filter
 114 #if ( ITK_VERSION_MAJOR < 5 )
 115   typedef itk::VectorResampleImageFilter<ImageType,ImageType> FilterType;
 116 #else
 117   typedef itk::ResampleImageFilter<ImageType,ImageType> FilterType;
 118 #endif
 119   typename FilterType::Pointer filter = FilterType::New();
 120
 121   // Instance of the transform object to be passed to the resample
 122   // filter. By default, identity transform is applied
 123   typedef itk::AffineTransform<double, ImageType::ImageDimension> TransformType;
 124   typename TransformType::Pointer transform =  TransformType::New();
 125   filter->SetTransform(transform);
 126
 127   // Set filter's parameters
 128   SizeType outputSize;
 129   SpacingType outputSpacing;
 130   PointType outputOrigin;
 131   for(unsigned int i=0; i<ImageType::ImageDimension; i++) {
 132     outputSize[i] = mOutputSize[i];
 133     outputSpacing[i] = mOutputSpacing[i];
 134     outputOrigin[i] = inputImage->GetOrigin()[i];
 135   }
 136
 137   filter->SetSize(outputSize);
 138   filter->SetOutputSpacing(outputSpacing);
 139   filter->SetOutputOrigin(outputOrigin);
 140   filter->SetDefaultPixelValue(static_cast<typename ImageType::PixelType>(mDefaultPixelValue));
 141
 142   // Select interpolator
 143   if (mInterpolatorName == "nn") {
 144 #if ( ITK_VERSION_MAJOR < 5 )
 145     typedef itk::VectorNearestNeighborInterpolateImageFunction<ImageType, double> InterpolatorType;
 146 #else
 147     typedef itk::NearestNeighborInterpolateImageFunction<ImageType, double> InterpolatorType;
 148 #endif
 149     typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
 150     filter->SetInterpolator(interpolator);
 151   } else {
 152     if (mInterpolatorName == "linear") {
 153 #if ( ITK_VERSION_MAJOR < 5 )
 154       typedef itk::VectorLinearInterpolateImageFunction<ImageType, double> InterpolatorType;
 155 #else
 156       typedef itk::LinearInterpolateImageFunction<ImageType, double> InterpolatorType;
 157 #endif
 158       typename InterpolatorType::Pointer interpolator =  InterpolatorType::New();
 159       filter->SetInterpolator(interpolator);
 160     } else {
 161       std::cerr << "Sorry, I do not know the interpolator (for vector field) '" << mInterpolatorName
 162                 << "'. Known interpolators are :  nn, linear" << std::endl;
 163       exit(0);
 164     }
 165   }
 166
 167   // Build initial Gaussian bluring (if needed)
 168   typedef itk::RecursiveGaussianImageFilter<ImageType, ImageType> GaussianFilterType;
 169   std::vector<typename GaussianFilterType::Pointer> gaussianFilters;
 170   if (mApplyGaussianFilterBefore) {
 171     for(unsigned int i=0; i<ImageType::ImageDimension; i++) {
 172       // Create filter
 173       gaussianFilters.push_back(GaussianFilterType::New());
 174       // Set options
 175       gaussianFilters[i]->SetDirection(i);
 176       gaussianFilters[i]->SetOrder(GaussianFilterType::ZeroOrder);
 177       gaussianFilters[i]->SetNormalizeAcrossScale(false);
 178       gaussianFilters[i]->SetSigma(mSigma[i]); // in millimeter !
 179       // Set input
 180       if (i==0) gaussianFilters[i]->SetInput(inputImage);
 181       else gaussianFilters[i]->SetInput(gaussianFilters[i-1]->GetOutput());
 182     }
 183     filter->SetInput(gaussianFilters[ImageType::ImageDimension-1]->GetOutput());
 184   } else {
 185     filter->SetInput(inputImage);
 186   }
 187
 188   // Go !
 189   try {
 190     filter->Update();
 191   } catch( itk::ExceptionObject & err ) {
 192     std::cerr << "Error while filtering " << m_InputFilenames[0].c_str()
 193               << " " << err << std::endl;
 194     exit(0);
 195   }
 196
 197   // Return result
 198   return filter->GetOutput();
 199
 200 }
 201 //--------------------------------------------------------------------
 202
 203
 204 //--------------------------------------------------------------------
 205 void clitk::VFResampleGenericFilter::SetOutputSize(const std::vector<int> & size)
 206 {
 207   mOutputSize.resize(size.size());
 208   std::copy(size.begin(), size.end(), mOutputSize.begin());
 209 }
 210 //--------------------------------------------------------------------
 211
 212 //--------------------------------------------------------------------
 213 void clitk::VFResampleGenericFilter::SetOutputSpacing(const std::vector<double> & spacing)
 214 {
 215   mOutputSpacing.resize(spacing.size());
 216   std::copy(spacing.begin(), spacing.end(), mOutputSpacing.begin());
 217 }
 218 //--------------------------------------------------------------------
 219
 220 //--------------------------------------------------------------------
 221 void clitk::VFResampleGenericFilter::SetInterpolationName(const std::string & inter)
 222 {
 223   mInterpolatorName = inter;
 224 }
 225 //--------------------------------------------------------------------
 226
 227 //--------------------------------------------------------------------
 228 void clitk::VFResampleGenericFilter::SetGaussianSigma(const std::vector<double> & sigma)
 229 {
 230   mApplyGaussianFilterBefore = true;
 231   mSigma.resize(sigma.size());
 232   std::copy(sigma.begin(), sigma.end(), mSigma.begin());
 233 }
 234 //--------------------------------------------------------------------
 235
 236 #endif
 237