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