1 #ifndef CLITKIMAGERESAMPLEGENERICFILTER2_CXX
2 #define CLITKIMAGERESAMPLEGENERICFILTER2_CXX
5 -------------------------------------------------------------------
6 * @file clitkImageResampleGenericFilter.cxx
7 * @author David Sarrut <David.Sarrut@creatis.insa-lyon.fr>
8 * @date 23 Feb 2008 08:37:53
12 -------------------------------------------------------------------*/
14 #include "clitkImageResampleGenericFilter.h"
18 #include "itkImageFileReader.h"
19 #include "itkImageSeriesReader.h"
20 #include "itkImageFileWriter.h"
21 #include "itkRecursiveGaussianImageFilter.h"
22 #include "itkResampleImageFilter.h"
23 #include "itkAffineTransform.h"
24 #include "itkNearestNeighborInterpolateImageFunction.h"
25 #include "itkLinearInterpolateImageFunction.h"
26 #include "itkBSplineInterpolateImageFunction.h"
27 #include "itkBSplineInterpolateImageFunctionWithLUT.h"
28 #include "itkCommand.h"
30 //--------------------------------------------------------------------
31 clitk::ImageResampleGenericFilter::ImageResampleGenericFilter():
32 ImageToImageGenericFilter<Self>("ImageResample") {
33 mApplyGaussianFilterBefore = false;
34 mDefaultPixelValue = 0.0;
35 mInterpolatorName = "NN";
37 InitializeImageTypeWithDim<2>();
38 InitializeImageTypeWithDim<3>();
39 InitializeImageTypeWithDim<4>();
41 //--------------------------------------------------------------------
44 //--------------------------------------------------------------------
45 template<unsigned int Dim>
46 void clitk::ImageResampleGenericFilter::InitializeImageTypeWithDim() {
47 ADD_IMAGE_TYPE(Dim, char);
48 ADD_IMAGE_TYPE(Dim, unsigned char);
49 ADD_IMAGE_TYPE(Dim, short);
50 ADD_IMAGE_TYPE(Dim, int);
51 ADD_IMAGE_TYPE(Dim, float);
53 //--------------------------------------------------------------------
56 //--------------------------------------------------------------------
57 template<class InputImageType>
58 void clitk::ImageResampleGenericFilter::UpdateWithInputImageType() {
61 typedef typename InputImageType::SizeType SizeType;
62 typedef typename InputImageType::SpacingType SpacingType;
63 typedef typename InputImageType::PointType PointType;
64 typedef typename InputImageType::PixelType PixelType;
65 static unsigned int dim = InputImageType::ImageDimension;
68 typename InputImageType::Pointer input = this->GetInput<InputImageType>(0);
71 if (!std::numeric_limits<PixelType>::is_signed) {
72 if ((mInterpolatorName == "bspline") || (mInterpolatorName == "blut")) {
73 std::cerr << "Warning : input pixel type is not signed, use bspline interpolation at your own risk ..." << std::endl;
78 if (mOutputSize.size() != dim) {
79 std::cerr << "Please set size with " << dim << " dimensions." << std::endl;
82 if (mOutputSpacing.size() != dim) {
83 std::cerr << "Please set spacing with " << dim << " dimensions." << std::endl;
86 mOutputOrigin.resize(dim);
88 if (mApplyGaussianFilterBefore && mSigma.size() != dim) {
89 std::cerr << "Please set sigma with " << dim << " dimensions." << std::endl;
93 // Create Image Filter
94 typedef itk::ResampleImageFilter<InputImageType,InputImageType> FilterType;
95 typename FilterType::Pointer filter = FilterType::New();
97 // Instance of the transform object to be passed to the resample
98 // filter. By default, identity transform is applied
99 typedef itk::AffineTransform<double, InputImageType::ImageDimension> TransformType;
100 typename TransformType::Pointer transform = TransformType::New();
101 filter->SetTransform(transform);
103 // Set filter's parameters
105 SpacingType outputSpacing;
106 PointType outputOrigin;
107 for(unsigned int i=0; i<InputImageType::ImageDimension; i++) {
108 outputSize[i] = mOutputSize[i];
109 outputSpacing[i] = mOutputSpacing[i];
110 outputOrigin[i] = input->GetOrigin()[i];
113 filter->SetSize(outputSize);
114 filter->SetOutputSpacing(outputSpacing);
115 filter->SetOutputOrigin(outputOrigin);
116 filter->SetDefaultPixelValue(static_cast<PixelType>(mDefaultPixelValue));//DS TODO//JV comme ça?
118 // Select interpolator
119 if (mInterpolatorName == "nn") {
120 typedef itk::NearestNeighborInterpolateImageFunction<InputImageType, double> InterpolatorType;
121 typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
122 filter->SetInterpolator(interpolator);
125 if (mInterpolatorName == "linear") {
126 typedef itk::LinearInterpolateImageFunction<InputImageType, double> InterpolatorType;
127 typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
128 filter->SetInterpolator(interpolator);
131 if (mInterpolatorName == "bspline") {
132 typedef itk::BSplineInterpolateImageFunction<InputImageType, double> InterpolatorType;
133 typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
134 interpolator->SetSplineOrder(mBSplineOrder);
135 filter->SetInterpolator(interpolator);
138 if (mInterpolatorName == "blut") {
139 typedef itk::BSplineInterpolateImageFunctionWithLUT<InputImageType, double> InterpolatorType;
140 typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
141 interpolator->SetSplineOrder(mBSplineOrder);
142 interpolator->SetLUTSamplingFactor(mSamplingFactors[0]);
143 filter->SetInterpolator(interpolator);
146 std::cerr << "Sorry, I do not know the interpolator '" << mInterpolatorName
147 << "'. Known interpolators are : nn, linear, bspline, blut" << std::endl;
154 // Build initial Gaussian bluring (if needed)
155 typedef itk::RecursiveGaussianImageFilter<InputImageType, InputImageType> GaussianFilterType;
156 std::vector<typename GaussianFilterType::Pointer> gaussianFilters;
157 if (mApplyGaussianFilterBefore) {
158 for(unsigned int i=0; i<InputImageType::ImageDimension; i++) {
160 gaussianFilters.push_back(GaussianFilterType::New());
162 gaussianFilters[i]->SetDirection(i);
163 gaussianFilters[i]->SetOrder(GaussianFilterType::ZeroOrder);
164 gaussianFilters[i]->SetNormalizeAcrossScale(false);
165 gaussianFilters[i]->SetSigma(mSigma[i]); // in millimeter !
167 if (i==0) gaussianFilters[i]->SetInput(input);
168 else gaussianFilters[i]->SetInput(gaussianFilters[i-1]->GetOutput());
170 filter->SetInput(gaussianFilters[InputImageType::ImageDimension-1]->GetOutput());
173 filter->SetInput(input);
180 catch( itk::ExceptionObject & err ) {
181 std::cerr << "Error while filtering " << mInputFilenames[0].c_str()
182 << " " << err << std::endl;
187 typename InputImageType::Pointer outputImage = filter->GetOutput();
189 // Write/save results
190 this->SetNextOutput<InputImageType>(outputImage);
193 //--------------------------------------------------------------------
196 //--------------------------------------------------------------------
197 void clitk::ImageResampleGenericFilter::SetOutputSize(const std::vector<int> & size) {
198 mOutputSize.resize(size.size());
199 std::copy(size.begin(), size.end(), mOutputSize.begin());
201 //--------------------------------------------------------------------
203 //--------------------------------------------------------------------
204 void clitk::ImageResampleGenericFilter::SetOutputSpacing(const std::vector<double> & spacing) {
205 mOutputSpacing.resize(spacing.size());
206 std::copy(spacing.begin(), spacing.end(), mOutputSpacing.begin());
208 //--------------------------------------------------------------------
210 //--------------------------------------------------------------------
211 void clitk::ImageResampleGenericFilter::SetInterpolationName(const std::string & inter) {
212 mInterpolatorName = inter;
214 //--------------------------------------------------------------------
216 //--------------------------------------------------------------------
217 void clitk::ImageResampleGenericFilter::SetGaussianSigma(const std::vector<double> & sigma) {
218 mApplyGaussianFilterBefore = true;
219 mSigma.resize(sigma.size());
220 std::copy(sigma.begin(), sigma.end(), mSigma.begin());
222 //--------------------------------------------------------------------