filters/clitkVFResampleGenericFilter.txx

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