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 CLITKVFRESAMPLEGENERICFILTER_CXX
19 #define CLITKVFRESAMPLEGENERICFILTER_CXX
21 -------------------------------------------------------------------
22 * @file clitkVFResampleGenericFilter.cxx
23 * @author David Sarrut <David.Sarrut@creatis.insa-lyon.fr>
24 * @date 23 Feb 2008 08:37:53
28 -------------------------------------------------------------------*/
30 #include "clitkVFResampleGenericFilter.h"
32 //--------------------------------------------------------------------
33 clitk::VFResampleGenericFilter::VFResampleGenericFilter():
34 clitk::ImageToImageGenericFilter<Self>("VFResample") {
35 InitializeImageType<2>();
36 InitializeImageType<3>();
37 InitializeImageType<4>();
38 mApplyGaussianFilterBefore = false;
39 mDefaultPixelValue = 0.0;
40 mInterpolatorName = "NN";
43 //--------------------------------------------------------------------
46 //--------------------------------------------------------------------
47 template<unsigned int Dim>
48 void clitk::VFResampleGenericFilter::InitializeImageType() {
49 ADD_IMAGE_TYPE(Dim, float);
51 //--------------------------------------------------------------------
54 //--------------------------------------------------------------------
55 template<class ImageType>
56 void clitk::VFResampleGenericFilter::UpdateWithInputImageType() {
58 if (mNbOfComponents == 1) {
59 std::cerr << "Error, only one components ? Use clitkImageResample instead." << std::endl;
62 typedef typename ImageType::PixelType PixelType;
63 if (mNbOfComponents == 2) Update_WithDimAndPixelTypeAndComponent<ImageType::ImageDimension,PixelType,2>();
64 if (mNbOfComponents == 3) Update_WithDimAndPixelTypeAndComponent<ImageType::ImageDimension,PixelType,3>();
65 if (mNbOfComponents == 4) Update_WithDimAndPixelTypeAndComponent<ImageType::ImageDimension,PixelType,4>();
67 //--------------------------------------------------------------------
70 //--------------------------------------------------------------------
71 template<unsigned int Dim, class PixelType, unsigned int DimCompo>
72 void clitk::VFResampleGenericFilter::Update_WithDimAndPixelTypeAndComponent() {
74 typedef itk::Vector<PixelType, DimCompo> DisplacementType;
75 typedef itk::Image< DisplacementType, Dim > ImageType;
77 typename ImageType::Pointer input = clitk::readImage<ImageType>(mInputFilenames, mIOVerbose);
80 typename ImageType::Pointer outputImage = ComputeImage<ImageType>(input);
83 SetNextOutput<ImageType>(outputImage);
85 //--------------------------------------------------------------------
87 //--------------------------------------------------------------------
88 template<class ImageType>
89 typename ImageType::Pointer
90 clitk::VFResampleGenericFilter::ComputeImage(typename ImageType::Pointer inputImage) {
93 static unsigned int dim = ImageType::ImageDimension;
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;
110 typedef typename ImageType::SizeType SizeType;
111 typedef typename ImageType::SpacingType SpacingType;
112 typedef typename ImageType::PointType PointType;
114 // Create Image Filter
115 typedef itk::VectorResampleImageFilter<ImageType,ImageType> FilterType;
116 typename FilterType::Pointer filter = FilterType::New();
118 // Instance of the transform object to be passed to the resample
119 // filter. By default, identity transform is applied
120 typedef itk::AffineTransform<double, ImageType::ImageDimension> TransformType;
121 typename TransformType::Pointer transform = TransformType::New();
122 filter->SetTransform(transform);
124 // Set filter's parameters
126 SpacingType outputSpacing;
127 PointType outputOrigin;
128 for(unsigned int i=0; i<ImageType::ImageDimension; i++) {
129 outputSize[i] = mOutputSize[i];
130 outputSpacing[i] = mOutputSpacing[i];
131 outputOrigin[i] = inputImage->GetOrigin()[i];
134 filter->SetSize(outputSize);
135 filter->SetOutputSpacing(outputSpacing);
136 filter->SetOutputOrigin(outputOrigin);
137 filter->SetDefaultPixelValue(static_cast<typename ImageType::PixelType>(mDefaultPixelValue));
139 // Select interpolator
140 if (mInterpolatorName == "nn") {
141 typedef itk::VectorNearestNeighborInterpolateImageFunction<ImageType, double> InterpolatorType;
142 typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
143 filter->SetInterpolator(interpolator);
146 if (mInterpolatorName == "linear") {
147 typedef itk::VectorLinearInterpolateImageFunction<ImageType, double> InterpolatorType;
148 typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
149 filter->SetInterpolator(interpolator);
152 std::cerr << "Sorry, I do not know the interpolator (for vector field) '" << mInterpolatorName
153 << "'. Known interpolators are : nn, linear" << std::endl;
158 // Build initial Gaussian bluring (if needed)
159 typedef itk::RecursiveGaussianImageFilter<ImageType, ImageType> GaussianFilterType;
160 std::vector<typename GaussianFilterType::Pointer> gaussianFilters;
161 if (mApplyGaussianFilterBefore) {
162 for(unsigned int i=0; i<ImageType::ImageDimension; i++) {
164 gaussianFilters.push_back(GaussianFilterType::New());
166 gaussianFilters[i]->SetDirection(i);
167 gaussianFilters[i]->SetOrder(GaussianFilterType::ZeroOrder);
168 gaussianFilters[i]->SetNormalizeAcrossScale(false);
169 gaussianFilters[i]->SetSigma(mSigma[i]); // in millimeter !
171 if (i==0) gaussianFilters[i]->SetInput(inputImage);
172 else gaussianFilters[i]->SetInput(gaussianFilters[i-1]->GetOutput());
174 filter->SetInput(gaussianFilters[ImageType::ImageDimension-1]->GetOutput());
177 filter->SetInput(inputImage);
184 catch( itk::ExceptionObject & err ) {
185 std::cerr << "Error while filtering " << mInputFilenames[0].c_str()
186 << " " << err << std::endl;
191 return filter->GetOutput();
194 //--------------------------------------------------------------------
197 //--------------------------------------------------------------------
198 void clitk::VFResampleGenericFilter::SetOutputSize(const std::vector<int> & size) {
199 mOutputSize.resize(size.size());
200 std::copy(size.begin(), size.end(), mOutputSize.begin());
202 //--------------------------------------------------------------------
204 //--------------------------------------------------------------------
205 void clitk::VFResampleGenericFilter::SetOutputSpacing(const std::vector<double> & spacing) {
206 mOutputSpacing.resize(spacing.size());
207 std::copy(spacing.begin(), spacing.end(), mOutputSpacing.begin());
209 //--------------------------------------------------------------------
211 //--------------------------------------------------------------------
212 void clitk::VFResampleGenericFilter::SetInterpolationName(const std::string & inter) {
213 mInterpolatorName = inter;
215 //--------------------------------------------------------------------
217 //--------------------------------------------------------------------
218 void clitk::VFResampleGenericFilter::SetGaussianSigma(const std::vector<double> & sigma) {
219 mApplyGaussianFilterBefore = true;
220 mSigma.resize(sigma.size());
221 std::copy(sigma.begin(), sigma.end(), mSigma.begin());
223 //--------------------------------------------------------------------