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   typedef itk::VectorResampleImageFilter<ImageType,ImageType> FilterType;
 115   typename FilterType::Pointer filter = FilterType::New();
 116
 117   // Instance of the transform object to be passed to the resample
 118   // filter. By default, identity transform is applied
 119   typedef itk::AffineTransform<double, ImageType::ImageDimension> TransformType;
 120   typename TransformType::Pointer transform =  TransformType::New();
 121   filter->SetTransform(transform);
 122
 123   // Set filter's parameters
 124   SizeType outputSize;
 125   SpacingType outputSpacing;
 126   PointType outputOrigin;
 127   for(unsigned int i=0; i<ImageType::ImageDimension; i++) {
 128     outputSize[i] = mOutputSize[i];
 129     outputSpacing[i] = mOutputSpacing[i];
 130     outputOrigin[i] = inputImage->GetOrigin()[i];
 131   }
 132
 133   filter->SetSize(outputSize);
 134   filter->SetOutputSpacing(outputSpacing);
 135   filter->SetOutputOrigin(outputOrigin);
 136   filter->SetDefaultPixelValue(static_cast<typename ImageType::PixelType>(mDefaultPixelValue));
 137
 138   // Select interpolator
 139   if (mInterpolatorName == "nn") {
 140     typedef itk::VectorNearestNeighborInterpolateImageFunction<ImageType, double> InterpolatorType;
 141     typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
 142     filter->SetInterpolator(interpolator);
 143   } else {
 144     if (mInterpolatorName == "linear") {
 145       typedef itk::VectorLinearInterpolateImageFunction<ImageType, double> InterpolatorType;
 146       typename InterpolatorType::Pointer interpolator =  InterpolatorType::New();
 147       filter->SetInterpolator(interpolator);
 148     } else {
 149       std::cerr << "Sorry, I do not know the interpolator (for vector field) '" << mInterpolatorName
 150                 << "'. Known interpolators are :  nn, linear" << std::endl;
 151       exit(0);
 152     }
 153   }
 154
 155   // Build initial Gaussian bluring (if needed)
 156   typedef itk::RecursiveGaussianImageFilter<ImageType, ImageType> GaussianFilterType;
 157   std::vector<typename GaussianFilterType::Pointer> gaussianFilters;
 158   if (mApplyGaussianFilterBefore) {
 159     for(unsigned int i=0; i<ImageType::ImageDimension; i++) {
 160       // Create filter
 161       gaussianFilters.push_back(GaussianFilterType::New());
 162       // Set options
 163       gaussianFilters[i]->SetDirection(i);
 164       gaussianFilters[i]->SetOrder(GaussianFilterType::ZeroOrder);
 165       gaussianFilters[i]->SetNormalizeAcrossScale(false);
 166       gaussianFilters[i]->SetSigma(mSigma[i]); // in millimeter !
 167       // Set input
 168       if (i==0) gaussianFilters[i]->SetInput(inputImage);
 169       else gaussianFilters[i]->SetInput(gaussianFilters[i-1]->GetOutput());
 170     }
 171     filter->SetInput(gaussianFilters[ImageType::ImageDimension-1]->GetOutput());
 172   } else {
 173     filter->SetInput(inputImage);
 174   }
 175
 176   // Go !
 177   try {
 178     filter->Update();
 179   } catch( itk::ExceptionObject & err ) {
 180     std::cerr << "Error while filtering " << m_InputFilenames[0].c_str()
 181               << " " << err << std::endl;
 182     exit(0);
 183   }
 184
 185   // Return result
 186   return filter->GetOutput();
 187
 188 }
 189 //--------------------------------------------------------------------
 190
 191
 192 //--------------------------------------------------------------------
 193 void clitk::VFResampleGenericFilter::SetOutputSize(const std::vector<int> & size)
 194 {
 195   mOutputSize.resize(size.size());
 196   std::copy(size.begin(), size.end(), mOutputSize.begin());
 197 }
 198 //--------------------------------------------------------------------
 199
 200 //--------------------------------------------------------------------
 201 void clitk::VFResampleGenericFilter::SetOutputSpacing(const std::vector<double> & spacing)
 202 {
 203   mOutputSpacing.resize(spacing.size());
 204   std::copy(spacing.begin(), spacing.end(), mOutputSpacing.begin());
 205 }
 206 //--------------------------------------------------------------------
 207
 208 //--------------------------------------------------------------------
 209 void clitk::VFResampleGenericFilter::SetInterpolationName(const std::string & inter)
 210 {
 211   mInterpolatorName = inter;
 212 }
 213 //--------------------------------------------------------------------
 214
 215 //--------------------------------------------------------------------
 216 void clitk::VFResampleGenericFilter::SetGaussianSigma(const std::vector<double> & sigma)
 217 {
 218   mApplyGaussianFilterBefore = true;
 219   mSigma.resize(sigma.size());
 220   std::copy(sigma.begin(), sigma.end(), mSigma.begin());
 221 }
 222 //--------------------------------------------------------------------
 223
 224 #endif
 225