filters/clitkImageResampleGenericFilter.txx

   1 #ifndef CLITKIMAGERESAMPLEGENERICFILTER_TXX
   2 #define CLITKIMAGERESAMPLEGENERICFILTER_TXX
   3
   4 /**
   5    ------------------------------------------------=
   6    * @file   clitkImageResampleGenericFilter.txx
   7    * @author David Sarrut <david.sarrut@creatis.insa-lyon.fr>
   8    * @date   23 Feb 2008 08:40:11
   9    *
  10    * @brief
  11    *
  12    *
  13    ------------------------------------------------=*/
  14
  15 //--------------------------------------------------------------------
  16 template<unsigned int Dim>
  17 void ImageResampleGenericFilter::Update_WithDim() {
  18
  19 #define TRY_TYPE(TYPE)                                                  \
  20   if (IsSameType<TYPE>(mPixelTypeName)) { Update_WithDimAndPixelType<Dim, TYPE>(); return; }
  21   TRY_TYPE(signed char);
  22   TRY_TYPE(uchar);
  23   TRY_TYPE(short);
  24   TRY_TYPE(ushort);
  25   TRY_TYPE(int);
  26   TRY_TYPE(unsigned int);
  27   TRY_TYPE(float);
  28   TRY_TYPE(double);
  29 #undef TRY_TYPE
  30
  31   std::string list = CreateListOfTypes<char, uchar, short, ushort, int, float, double>();
  32   std::cerr << "Error, I don't know the type '" << mPixelTypeName << "' for the input image '"
  33             << mInputFilenames[0] << "'." << std::endl << "Known types are " << list << std::endl;
  34   exit(0);
  35 }
  36 //--------------------------------------------------------------------
  37
  38 //--------------------------------------------------------------------
  39 template<unsigned int Dim, class PixelType>
  40 void ImageResampleGenericFilter::Update_WithDimAndPixelType() {
  41   // Reading input
  42   typedef itk::Image<PixelType,Dim> ImageType;
  43   //typename ImageType::Pointer input = clitk::readImage<ImageType>(mInputFilenames, mIOVerbose);
  44   typename ImageType::Pointer input = this->GetInput<ImageType>(0);
  45
  46   // Main filter
  47   typename ImageType::Pointer outputImage = ComputeImage<ImageType>(input);
  48
  49   // Write results
  50   this->SetNextOutput<ImageType>(outputImage);
  51   //clitk::writeImage<ImageType>(outputImage, mOutputFilename, mIOVerbose);
  52 }
  53 //--------------------------------------------------------------------
  54
  55 //--------------------------------------------------------------------
  56 template<class ImageType>
  57 typename ImageType::Pointer
  58 ImageResampleGenericFilter::ComputeImage(typename ImageType::Pointer inputImage) {
  59
  60   // Warning
  61   if (!std::numeric_limits<typename ImageType::PixelType>::is_signed) {
  62     if ((mInterpolatorName == "bspline") || (mInterpolatorName == "blut")) {
  63       std::cerr << "Warning : input pixel type is not signed, use bspline interpolation at your own risk ..." << std::endl;
  64     }
  65   }
  66
  67   // Check options
  68   static unsigned int dim = ImageType::ImageDimension;
  69   if (mOutputSize.size() != dim) {
  70     std::cerr << "Please set size with " << dim << " dimensions." << std::endl;
  71     return NULL;
  72   }
  73   if (mOutputSpacing.size() != dim) {
  74     std::cerr << "Please set spacing with " << dim << " dimensions." << std::endl;
  75     return NULL;
  76   }
  77   mOutputOrigin.resize(dim);
  78
  79   if (mApplyGaussianFilterBefore && mSigma.size() != dim) {
  80     std::cerr << "Please set sigma with " << dim << " dimensions." << std::endl;
  81     return NULL;
  82   }
  83
  84   // Some typedefs
  85   typedef typename ImageType::SizeType SizeType;
  86   typedef typename ImageType::SpacingType SpacingType;
  87   typedef typename ImageType::PointType PointType;
  88
  89
  90   // Create Image Filter
  91   typedef itk::ResampleImageFilter<ImageType,ImageType> FilterType;
  92   typename FilterType::Pointer filter = FilterType::New();
  93
  94   // Instance of the transform object to be passed to the resample
  95   // filter. By default, identity transform is applied
  96   typedef itk::AffineTransform<double, ImageType::ImageDimension> TransformType;
  97   typename TransformType::Pointer transform =  TransformType::New();
  98   filter->SetTransform(transform);
  99
 100   // Set filter's parameters
 101   SizeType outputSize;
 102   SpacingType outputSpacing;
 103   PointType outputOrigin;
 104   for(unsigned int i=0; i<ImageType::ImageDimension; i++) {
 105     outputSize[i] = mOutputSize[i];
 106     outputSpacing[i] = mOutputSpacing[i];
 107     outputOrigin[i] = inputImage->GetOrigin()[i];
 108   }
 109
 110   filter->SetSize(outputSize);
 111   filter->SetOutputSpacing(outputSpacing);
 112   filter->SetOutputOrigin(outputOrigin);
 113   filter->SetDefaultPixelValue(static_cast<typename ImageType::PixelType>(mDefaultPixelValue));//DS TODO//JV comme ça?
 114
 115   // Select interpolator
 116   if (mInterpolatorName == "nn") {
 117     typedef itk::NearestNeighborInterpolateImageFunction<ImageType, double> InterpolatorType;
 118     typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
 119     filter->SetInterpolator(interpolator);
 120   }
 121   else {
 122     if (mInterpolatorName == "linear") {
 123       typedef itk::LinearInterpolateImageFunction<ImageType, double> InterpolatorType;
 124       typename InterpolatorType::Pointer interpolator =  InterpolatorType::New();
 125       filter->SetInterpolator(interpolator);
 126     }
 127     else {
 128       if (mInterpolatorName == "bspline") {
 129         typedef itk::BSplineInterpolateImageFunction<ImageType, double> InterpolatorType;
 130         typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
 131         interpolator->SetSplineOrder(mBSplineOrder);
 132         filter->SetInterpolator(interpolator);
 133       }
 134       else {
 135         if (mInterpolatorName == "blut") {
 136           typedef itk::BSplineInterpolateImageFunctionWithLUT<ImageType, double> InterpolatorType;
 137           typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
 138           interpolator->SetSplineOrder(mBSplineOrder);
 139           interpolator->SetLUTSamplingFactor(mSamplingFactors[0]);
 140           filter->SetInterpolator(interpolator);
 141         }
 142         else {
 143           std::cerr << "Sorry, I do not know the interpolator '" << mInterpolatorName
 144                     << "'. Known interpolators are :  nn, linear, bspline, blut" << std::endl;
 145           exit(0);
 146         }
 147       }
 148     }
 149   }
 150
 151   // Build initial Gaussian bluring (if needed)
 152   typedef itk::RecursiveGaussianImageFilter<ImageType, ImageType> GaussianFilterType;
 153   std::vector<typename GaussianFilterType::Pointer> gaussianFilters;
 154   if (mApplyGaussianFilterBefore) {
 155     for(unsigned int i=0; i<ImageType::ImageDimension; i++) {
 156       // Create filter
 157       gaussianFilters.push_back(GaussianFilterType::New());
 158       // Set options
 159       gaussianFilters[i]->SetDirection(i);
 160       gaussianFilters[i]->SetOrder(GaussianFilterType::ZeroOrder);
 161       gaussianFilters[i]->SetNormalizeAcrossScale(false);
 162       gaussianFilters[i]->SetSigma(mSigma[i]); // in millimeter !
 163       // Set input
 164       if (i==0) gaussianFilters[i]->SetInput(inputImage);
 165       else gaussianFilters[i]->SetInput(gaussianFilters[i-1]->GetOutput());
 166     }
 167     filter->SetInput(gaussianFilters[ImageType::ImageDimension-1]->GetOutput());
 168   }
 169   else {
 170     filter->SetInput(inputImage);
 171   }
 172
 173   // Go !
 174   try {
 175     filter->Update();
 176   }
 177   catch( itk::ExceptionObject & err ) {
 178     std::cerr << "Error while filtering " << mInputFilenames[0].c_str()
 179               << " " << err << std::endl;
 180     exit(0);
 181   }
 182
 183   // Return result
 184   return filter->GetOutput();
 185
 186 }
 187 //--------------------------------------------------------------------
 188
 189
 190
 191 #endif /* end #define CLITKIMAGERESAMPLEGENERICFILTER_TXX */
 192