filters/clitkVFResampleGenericFilter.cxx

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