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