filters/clitkImageResampleGenericFilter.cxx

   1 #ifndef CLITKIMAGERESAMPLEGENERICFILTER2_CXX
   2 #define CLITKIMAGERESAMPLEGENERICFILTER2_CXX
   3
   4 /**
   5  -------------------------------------------------------------------
   6  * @file   clitkImageResampleGenericFilter.cxx
   7  * @author David Sarrut <David.Sarrut@creatis.insa-lyon.fr>
   8  * @date   23 Feb 2008 08:37:53
   9
  10  * @brief
  11
  12  -------------------------------------------------------------------*/
  13
  14 #include "clitkImageResampleGenericFilter.h"
  15
  16 // itk include
  17 #include "itkImage.h"
  18 #include "itkImageFileReader.h"
  19 #include "itkImageSeriesReader.h"
  20 #include "itkImageFileWriter.h"
  21 #include "itkRecursiveGaussianImageFilter.h"
  22 #include "itkResampleImageFilter.h"
  23 #include "itkAffineTransform.h"
  24 #include "itkNearestNeighborInterpolateImageFunction.h"
  25 #include "itkLinearInterpolateImageFunction.h"
  26 #include "itkBSplineInterpolateImageFunction.h"
  27 #include "itkBSplineInterpolateImageFunctionWithLUT.h"
  28 #include "itkCommand.h"
  29
  30 //--------------------------------------------------------------------
  31 clitk::ImageResampleGenericFilter::ImageResampleGenericFilter():
  32   ImageToImageGenericFilter<Self>("ImageResample") {
  33   mApplyGaussianFilterBefore = false;
  34   mDefaultPixelValue = 0.0;
  35   mInterpolatorName = "NN";
  36   mBSplineOrder=3;
  37   InitializeImageTypeWithDim<2>();
  38   InitializeImageTypeWithDim<3>();
  39   InitializeImageTypeWithDim<4>();
  40 }
  41 //--------------------------------------------------------------------
  42
  43
  44 //--------------------------------------------------------------------
  45 template<unsigned int Dim>
  46 void clitk::ImageResampleGenericFilter::InitializeImageTypeWithDim() {
  47   ADD_IMAGE_TYPE(Dim, char);
  48   ADD_IMAGE_TYPE(Dim, uchar);
  49   ADD_IMAGE_TYPE(Dim, short);
  50   ADD_IMAGE_TYPE(Dim, int);
  51   ADD_IMAGE_TYPE(Dim, float);
  52 }
  53 //--------------------------------------------------------------------
  54
  55
  56 //--------------------------------------------------------------------
  57 template<class InputImageType>
  58 void clitk::ImageResampleGenericFilter::UpdateWithInputImageType() {
  59
  60   // Some typedefs
  61   typedef typename InputImageType::SizeType    SizeType;
  62   typedef typename InputImageType::SpacingType SpacingType;
  63   typedef typename InputImageType::PointType   PointType;
  64   typedef typename InputImageType::PixelType   PixelType;
  65   static unsigned int dim = InputImageType::ImageDimension;
  66
  67   // Reading input
  68   typename InputImageType::Pointer input = this->GetInput<InputImageType>(0);
  69
  70   // Warning
  71   if (!std::numeric_limits<PixelType>::is_signed) {
  72     if ((mInterpolatorName == "bspline") || (mInterpolatorName == "blut")) {
  73       std::cerr << "Warning : input pixel type is not signed, use bspline interpolation at your own risk ..." << std::endl;
  74     }
  75   }
  76
  77   // Check options
  78   if (mOutputSize.size() != dim) {
  79     std::cerr << "Please set size with " << dim << " dimensions." << std::endl;
  80     return;
  81   }
  82   if (mOutputSpacing.size() != dim) {
  83     std::cerr << "Please set spacing with " << dim << " dimensions." << std::endl;
  84     return;
  85   }
  86   mOutputOrigin.resize(dim);
  87
  88   if (mApplyGaussianFilterBefore && mSigma.size() != dim) {
  89     std::cerr << "Please set sigma with " << dim << " dimensions." << std::endl;
  90     return;
  91   }
  92
  93   // Create Image Filter
  94   typedef itk::ResampleImageFilter<InputImageType,InputImageType> FilterType;
  95   typename FilterType::Pointer filter = FilterType::New();
  96
  97   // Instance of the transform object to be passed to the resample
  98   // filter. By default, identity transform is applied
  99   typedef itk::AffineTransform<double, InputImageType::ImageDimension> TransformType;
 100   typename TransformType::Pointer transform =  TransformType::New();
 101   filter->SetTransform(transform);
 102
 103   // Set filter's parameters
 104   SizeType outputSize;
 105   SpacingType outputSpacing;
 106   PointType outputOrigin;
 107   for(unsigned int i=0; i<InputImageType::ImageDimension; i++) {
 108     outputSize[i] = mOutputSize[i];
 109     outputSpacing[i] = mOutputSpacing[i];
 110     outputOrigin[i] = input->GetOrigin()[i];
 111   }
 112
 113   filter->SetSize(outputSize);
 114   filter->SetOutputSpacing(outputSpacing);
 115   filter->SetOutputOrigin(outputOrigin);
 116   filter->SetDefaultPixelValue(static_cast<PixelType>(mDefaultPixelValue));//DS TODO//JV comme ça?
 117
 118   // Select interpolator
 119   if (mInterpolatorName == "nn") {
 120     typedef itk::NearestNeighborInterpolateImageFunction<InputImageType, double> InterpolatorType;
 121     typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
 122     filter->SetInterpolator(interpolator);
 123   }
 124   else {
 125     if (mInterpolatorName == "linear") {
 126       typedef itk::LinearInterpolateImageFunction<InputImageType, double> InterpolatorType;
 127       typename InterpolatorType::Pointer interpolator =  InterpolatorType::New();
 128       filter->SetInterpolator(interpolator);
 129     }
 130     else {
 131       if (mInterpolatorName == "bspline") {
 132         typedef itk::BSplineInterpolateImageFunction<InputImageType, double> InterpolatorType;
 133         typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
 134         interpolator->SetSplineOrder(mBSplineOrder);
 135         filter->SetInterpolator(interpolator);
 136       }
 137       else {
 138         if (mInterpolatorName == "blut") {
 139           typedef itk::BSplineInterpolateImageFunctionWithLUT<InputImageType, double> InterpolatorType;
 140           typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
 141           interpolator->SetSplineOrder(mBSplineOrder);
 142           interpolator->SetLUTSamplingFactor(mSamplingFactors[0]);
 143           filter->SetInterpolator(interpolator);
 144         }
 145         else {
 146           std::cerr << "Sorry, I do not know the interpolator '" << mInterpolatorName
 147                     << "'. Known interpolators are :  nn, linear, bspline, blut" << std::endl;
 148           exit(0);
 149         }
 150       }
 151     }
 152   }
 153
 154   // Build initial Gaussian bluring (if needed)
 155   typedef itk::RecursiveGaussianImageFilter<InputImageType, InputImageType> GaussianFilterType;
 156   std::vector<typename GaussianFilterType::Pointer> gaussianFilters;
 157   if (mApplyGaussianFilterBefore) {
 158     for(unsigned int i=0; i<InputImageType::ImageDimension; i++) {
 159       // Create filter
 160       gaussianFilters.push_back(GaussianFilterType::New());
 161       // Set options
 162       gaussianFilters[i]->SetDirection(i);
 163       gaussianFilters[i]->SetOrder(GaussianFilterType::ZeroOrder);
 164       gaussianFilters[i]->SetNormalizeAcrossScale(false);
 165       gaussianFilters[i]->SetSigma(mSigma[i]); // in millimeter !
 166       // Set input
 167       if (i==0) gaussianFilters[i]->SetInput(input);
 168       else gaussianFilters[i]->SetInput(gaussianFilters[i-1]->GetOutput());
 169     }
 170     filter->SetInput(gaussianFilters[InputImageType::ImageDimension-1]->GetOutput());
 171   }
 172   else {
 173     filter->SetInput(input);
 174   }
 175
 176   // Go !
 177   try {
 178     filter->Update();
 179   }
 180   catch( itk::ExceptionObject & err ) {
 181     std::cerr << "Error while filtering " << mInputFilenames[0].c_str()
 182               << " " << err << std::endl;
 183     exit(0);
 184   }
 185
 186   // Get result
 187    typename InputImageType::Pointer outputImage = filter->GetOutput();
 188
 189   // Write/save results
 190   this->SetNextOutput<InputImageType>(outputImage);
 191
 192 }
 193 //--------------------------------------------------------------------
 194
 195
 196 //--------------------------------------------------------------------
 197 void clitk::ImageResampleGenericFilter::SetOutputSize(const std::vector<int> & size) {
 198   mOutputSize.resize(size.size());
 199   std::copy(size.begin(), size.end(), mOutputSize.begin());
 200 }
 201 //--------------------------------------------------------------------
 202
 203 //--------------------------------------------------------------------
 204 void clitk::ImageResampleGenericFilter::SetOutputSpacing(const std::vector<double> & spacing) {
 205   mOutputSpacing.resize(spacing.size());
 206   std::copy(spacing.begin(), spacing.end(), mOutputSpacing.begin());
 207 }
 208 //--------------------------------------------------------------------
 209
 210 //--------------------------------------------------------------------
 211 void clitk::ImageResampleGenericFilter::SetInterpolationName(const std::string & inter) {
 212   mInterpolatorName = inter;
 213 }
 214 //--------------------------------------------------------------------
 215
 216 //--------------------------------------------------------------------
 217 void clitk::ImageResampleGenericFilter::SetGaussianSigma(const std::vector<double> & sigma) {
 218   mApplyGaussianFilterBefore = true;
 219   mSigma.resize(sigma.size());
 220   std::copy(sigma.begin(), sigma.end(), mSigma.begin());
 221 }
 222 //--------------------------------------------------------------------
 223
 224 #endif
 225