filters/clitkVFResampleGenericFilter.cxx

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