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, short);
49 ADD_IMAGE_TYPE(Dim, int);
50 ADD_IMAGE_TYPE(Dim, float);
52 //--------------------------------------------------------------------
55 //--------------------------------------------------------------------
56 template<class InputImageType>
57 void clitk::ImageResampleGenericFilter::UpdateWithInputImageType() {
60 typedef typename InputImageType::SizeType SizeType;
61 typedef typename InputImageType::SpacingType SpacingType;
62 typedef typename InputImageType::PointType PointType;
63 typedef typename InputImageType::PixelType PixelType;
64 static unsigned int dim = InputImageType::ImageDimension;
67 typename InputImageType::Pointer input = this->GetInput<InputImageType>(0);
70 if (!std::numeric_limits<PixelType>::is_signed) {
71 if ((mInterpolatorName == "bspline") || (mInterpolatorName == "blut")) {
72 std::cerr << "Warning : input pixel type is not signed, use bspline interpolation at your own risk ..." << std::endl;
77 if (mOutputSize.size() != dim) {
78 std::cerr << "Please set size with " << dim << " dimensions." << std::endl;
81 if (mOutputSpacing.size() != dim) {
82 std::cerr << "Please set spacing with " << dim << " dimensions." << std::endl;
85 mOutputOrigin.resize(dim);
87 if (mApplyGaussianFilterBefore && mSigma.size() != dim) {
88 std::cerr << "Please set sigma with " << dim << " dimensions." << std::endl;
92 // Create Image Filter
93 typedef itk::ResampleImageFilter<InputImageType,InputImageType> FilterType;
94 typename FilterType::Pointer filter = FilterType::New();
96 // Instance of the transform object to be passed to the resample
97 // filter. By default, identity transform is applied
98 typedef itk::AffineTransform<double, InputImageType::ImageDimension> TransformType;
99 typename TransformType::Pointer transform = TransformType::New();
100 filter->SetTransform(transform);
102 // Set filter's parameters
104 SpacingType outputSpacing;
105 PointType outputOrigin;
106 for(unsigned int i=0; i<InputImageType::ImageDimension; i++) {
107 outputSize[i] = mOutputSize[i];
108 outputSpacing[i] = mOutputSpacing[i];
109 outputOrigin[i] = input->GetOrigin()[i];
112 filter->SetSize(outputSize);
113 filter->SetOutputSpacing(outputSpacing);
114 filter->SetOutputOrigin(outputOrigin);
115 filter->SetDefaultPixelValue(static_cast<PixelType>(mDefaultPixelValue));//DS TODO//JV comme ça?
117 // Select interpolator
118 if (mInterpolatorName == "nn") {
119 typedef itk::NearestNeighborInterpolateImageFunction<InputImageType, double> InterpolatorType;
120 typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
121 filter->SetInterpolator(interpolator);
124 if (mInterpolatorName == "linear") {
125 typedef itk::LinearInterpolateImageFunction<InputImageType, double> InterpolatorType;
126 typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
127 filter->SetInterpolator(interpolator);
130 if (mInterpolatorName == "bspline") {
131 typedef itk::BSplineInterpolateImageFunction<InputImageType, double> InterpolatorType;
132 typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
133 interpolator->SetSplineOrder(mBSplineOrder);
134 filter->SetInterpolator(interpolator);
137 if (mInterpolatorName == "blut") {
138 typedef itk::BSplineInterpolateImageFunctionWithLUT<InputImageType, double> InterpolatorType;
139 typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
140 interpolator->SetSplineOrder(mBSplineOrder);
141 interpolator->SetLUTSamplingFactor(mSamplingFactors[0]);
142 filter->SetInterpolator(interpolator);
145 std::cerr << "Sorry, I do not know the interpolator '" << mInterpolatorName
146 << "'. Known interpolators are : nn, linear, bspline, blut" << std::endl;
153 // Build initial Gaussian bluring (if needed)
154 typedef itk::RecursiveGaussianImageFilter<InputImageType, InputImageType> GaussianFilterType;
155 std::vector<typename GaussianFilterType::Pointer> gaussianFilters;
156 if (mApplyGaussianFilterBefore) {
157 for(unsigned int i=0; i<InputImageType::ImageDimension; i++) {
159 gaussianFilters.push_back(GaussianFilterType::New());
161 gaussianFilters[i]->SetDirection(i);
162 gaussianFilters[i]->SetOrder(GaussianFilterType::ZeroOrder);
163 gaussianFilters[i]->SetNormalizeAcrossScale(false);
164 gaussianFilters[i]->SetSigma(mSigma[i]); // in millimeter !
166 if (i==0) gaussianFilters[i]->SetInput(input);
167 else gaussianFilters[i]->SetInput(gaussianFilters[i-1]->GetOutput());
169 filter->SetInput(gaussianFilters[InputImageType::ImageDimension-1]->GetOutput());
172 filter->SetInput(input);
179 catch( itk::ExceptionObject & err ) {
180 std::cerr << "Error while filtering " << mInputFilenames[0].c_str()
181 << " " << err << std::endl;
186 typename InputImageType::Pointer outputImage = filter->GetOutput();
188 // Write/save results
189 this->SetNextOutput<InputImageType>(outputImage);
192 //--------------------------------------------------------------------
195 //--------------------------------------------------------------------
196 void clitk::ImageResampleGenericFilter::SetOutputSize(const std::vector<int> & size) {
197 mOutputSize.resize(size.size());
198 std::copy(size.begin(), size.end(), mOutputSize.begin());
200 //--------------------------------------------------------------------
202 //--------------------------------------------------------------------
203 void clitk::ImageResampleGenericFilter::SetOutputSpacing(const std::vector<double> & spacing) {
204 mOutputSpacing.resize(spacing.size());
205 std::copy(spacing.begin(), spacing.end(), mOutputSpacing.begin());
207 //--------------------------------------------------------------------
209 //--------------------------------------------------------------------
210 void clitk::ImageResampleGenericFilter::SetInterpolationName(const std::string & inter) {
211 mInterpolatorName = inter;
213 //--------------------------------------------------------------------
215 //--------------------------------------------------------------------
216 void clitk::ImageResampleGenericFilter::SetGaussianSigma(const std::vector<double> & sigma) {
217 mApplyGaussianFilterBefore = true;
218 mSigma.resize(sigma.size());
219 std::copy(sigma.begin(), sigma.end(), mSigma.begin());
221 //--------------------------------------------------------------------