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