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