1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
5 - University of LYON http://www.universite-lyon.fr/
6 - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr
7 - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
9 This software is distributed WITHOUT ANY WARRANTY; without even
10 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11 PURPOSE. See the copyright notices for more information.
13 It is distributed under dual licence
15 - BSD See included LICENSE.txt file
16 - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ======================================================================-====*/
18 #ifndef CLITKIMAGERESAMPLEGENERICFILTER2_CXX
19 #define CLITKIMAGERESAMPLEGENERICFILTER2_CXX
21 -------------------------------------------------------------------
22 * @file clitkImageResampleGenericFilter.cxx
23 * @author David Sarrut <David.Sarrut@creatis.insa-lyon.fr>
24 * @date 23 Feb 2008 08:37:53
28 -------------------------------------------------------------------*/
30 #include "clitkImageResampleGenericFilter.h"
34 #include "itkImageFileReader.h"
35 #include "itkImageSeriesReader.h"
36 #include "itkImageFileWriter.h"
37 #include "itkRecursiveGaussianImageFilter.h"
38 #include "itkResampleImageFilter.h"
39 #include "itkAffineTransform.h"
40 #include "itkNearestNeighborInterpolateImageFunction.h"
41 #include "itkWindowedSincInterpolateImageFunction.h"
42 #include "itkLinearInterpolateImageFunction.h"
43 #include "itkBSplineInterpolateImageFunction.h"
44 #include "itkBSplineInterpolateImageFunctionWithLUT.h"
45 #include "itkCommand.h"
47 //--------------------------------------------------------------------
48 clitk::ImageResampleGenericFilter::ImageResampleGenericFilter():
49 ImageToImageGenericFilter<Self>("ImageResample")
51 mApplyGaussianFilterBefore = false;
52 mDefaultPixelValue = 0.0;
53 mInterpolatorName = "NN";
55 InitializeImageTypeWithDim<2>();
56 InitializeImageTypeWithDim<3>();
57 InitializeImageTypeWithDim<4>();
59 //--------------------------------------------------------------------
62 //--------------------------------------------------------------------
63 template<unsigned int Dim>
64 void clitk::ImageResampleGenericFilter::InitializeImageTypeWithDim()
66 ADD_DEFAULT_IMAGE_TYPES(Dim);
68 //--------------------------------------------------------------------
71 //--------------------------------------------------------------------
72 template<class InputImageType>
73 void clitk::ImageResampleGenericFilter::UpdateWithInputImageType()
77 typedef typename InputImageType::SizeType SizeType;
78 typedef typename InputImageType::SpacingType SpacingType;
79 typedef typename InputImageType::PointType PointType;
80 typedef typename InputImageType::PixelType PixelType;
81 static unsigned int dim = InputImageType::ImageDimension;
84 typename InputImageType::Pointer input = this->GetInput<InputImageType>(0);
87 if (!std::numeric_limits<PixelType>::is_signed) {
88 if ((mInterpolatorName == "bspline") || (mInterpolatorName == "blut")) {
89 std::cerr << "Warning : input pixel type is not signed, use bspline interpolation at your own risk ..." << std::endl;
94 if (mOutputSize.size() != dim) {
95 std::cerr << "Please set size with " << dim << " dimensions." << std::endl;
98 if (mOutputSpacing.size() != dim) {
99 std::cerr << "Please set spacing with " << dim << " dimensions." << std::endl;
102 mOutputOrigin.resize(dim);
104 if (mApplyGaussianFilterBefore && mSigma.size() != dim) {
105 std::cerr << "Please set sigma with " << dim << " dimensions." << std::endl;
109 // Create Image Filter
110 typedef itk::ResampleImageFilter<InputImageType,InputImageType> FilterType;
111 typename FilterType::Pointer filter = FilterType::New();
113 // Instance of the transform object to be passed to the resample
114 // filter. By default, identity transform is applied
115 typedef itk::AffineTransform<double, InputImageType::ImageDimension> TransformType;
116 typename TransformType::Pointer transform = TransformType::New();
117 filter->SetTransform(transform);
119 // Set filter's parameters
121 SpacingType outputSpacing;
122 PointType outputOrigin;
123 for(unsigned int i=0; i<InputImageType::ImageDimension; i++) {
124 outputSize[i] = mOutputSize[i];
125 outputSpacing[i] = mOutputSpacing[i];
126 outputOrigin[i] = input->GetOrigin()[i];
129 filter->SetSize(outputSize);
130 filter->SetOutputSpacing(outputSpacing);
131 filter->SetOutputOrigin(outputOrigin);
132 filter->SetDefaultPixelValue(static_cast<PixelType>(mDefaultPixelValue));//DS TODO//JV comme ça?
134 // Select interpolator
135 if (mInterpolatorName == "nn") {
136 typedef itk::NearestNeighborInterpolateImageFunction<InputImageType, double> InterpolatorType;
137 typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
138 filter->SetInterpolator(interpolator);
140 if (mInterpolatorName == "linear") {
141 typedef itk::LinearInterpolateImageFunction<InputImageType, double> InterpolatorType;
142 typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
143 filter->SetInterpolator(interpolator);
145 if (mInterpolatorName == "windowed sinc") {
146 typedef itk::WindowedSincInterpolateImageFunction<InputImageType, 4> InterpolatorType;
147 typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
148 filter->SetInterpolator(interpolator);
150 if (mInterpolatorName == "bspline") {
151 typedef itk::BSplineInterpolateImageFunction<InputImageType, double> InterpolatorType;
152 typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
153 interpolator->SetSplineOrder(mBSplineOrder);
154 filter->SetInterpolator(interpolator);
156 if (mInterpolatorName == "blut") {
157 typedef itk::BSplineInterpolateImageFunctionWithLUT<InputImageType, double> InterpolatorType;
158 typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
159 interpolator->SetSplineOrder(mBSplineOrder);
160 interpolator->SetLUTSamplingFactor(mSamplingFactors[0]);
161 filter->SetInterpolator(interpolator);
163 std::cerr << "Sorry, I do not know the interpolator '" << mInterpolatorName
164 << "'. Known interpolators are : nn, linear, bspline, blut" << std::endl;
172 // Build initial Gaussian bluring (if needed)
173 typedef itk::RecursiveGaussianImageFilter<InputImageType, InputImageType> GaussianFilterType;
174 std::vector<typename GaussianFilterType::Pointer> gaussianFilters;
175 if (mApplyGaussianFilterBefore) {
176 for(unsigned int i=0; i<InputImageType::ImageDimension; i++) {
178 gaussianFilters.push_back(GaussianFilterType::New());
180 gaussianFilters[i]->SetDirection(i);
181 gaussianFilters[i]->SetOrder(GaussianFilterType::ZeroOrder);
182 gaussianFilters[i]->SetNormalizeAcrossScale(false);
183 gaussianFilters[i]->SetSigma(mSigma[i]); // in millimeter !
185 if (i==0) gaussianFilters[i]->SetInput(input);
186 else gaussianFilters[i]->SetInput(gaussianFilters[i-1]->GetOutput());
188 filter->SetInput(gaussianFilters[InputImageType::ImageDimension-1]->GetOutput());
190 filter->SetInput(input);
196 } catch( itk::ExceptionObject & err ) {
197 std::cerr << "Error while filtering " << mInputFilenames[0].c_str()
198 << " " << err << std::endl;
203 typename InputImageType::Pointer outputImage = filter->GetOutput();
205 // Write/save results
206 this->SetNextOutput<InputImageType>(outputImage);
209 //--------------------------------------------------------------------
212 //--------------------------------------------------------------------
213 void clitk::ImageResampleGenericFilter::SetOutputSize(const std::vector<int> & size)
215 mOutputSize.resize(size.size());
216 std::copy(size.begin(), size.end(), mOutputSize.begin());
218 //--------------------------------------------------------------------
220 //--------------------------------------------------------------------
221 void clitk::ImageResampleGenericFilter::SetOutputSpacing(const std::vector<double> & spacing)
223 mOutputSpacing.resize(spacing.size());
224 std::copy(spacing.begin(), spacing.end(), mOutputSpacing.begin());
226 //--------------------------------------------------------------------
228 //--------------------------------------------------------------------
229 void clitk::ImageResampleGenericFilter::SetInterpolationName(const std::string & inter)
231 mInterpolatorName = inter;
233 //--------------------------------------------------------------------
235 //--------------------------------------------------------------------
236 void clitk::ImageResampleGenericFilter::SetGaussianSigma(const std::vector<double> & sigma)
238 mApplyGaussianFilterBefore = true;
239 mSigma.resize(sigma.size());
240 std::copy(sigma.begin(), sigma.end(), mSigma.begin());
242 //--------------------------------------------------------------------