1 #ifndef CLITKVFRESAMPLEGENERICFILTER_CXX
2 #define CLITKVFRESAMPLEGENERICFILTER_CXX
4 -------------------------------------------------------------------
5 * @file clitkVFResampleGenericFilter.cxx
6 * @author David Sarrut <David.Sarrut@creatis.insa-lyon.fr>
7 * @date 23 Feb 2008 08:37:53
11 -------------------------------------------------------------------*/
13 #include "clitkVFResampleGenericFilter.h"
15 //--------------------------------------------------------------------
16 clitk::VFResampleGenericFilter::VFResampleGenericFilter():
17 clitk::ImageToImageGenericFilter<Self>("VFResample") {
18 InitializeImageType<2>();
19 InitializeImageType<3>();
20 InitializeImageType<4>();
21 mApplyGaussianFilterBefore = false;
22 mDefaultPixelValue = 0.0;
23 mInterpolatorName = "NN";
26 //--------------------------------------------------------------------
29 //--------------------------------------------------------------------
30 template<unsigned int Dim>
31 void clitk::VFResampleGenericFilter::InitializeImageType() {
32 ADD_IMAGE_TYPE(Dim, float);
34 //--------------------------------------------------------------------
37 //--------------------------------------------------------------------
38 template<class ImageType>
39 void clitk::VFResampleGenericFilter::UpdateWithInputImageType() {
41 if (mNbOfComponents == 1) {
42 std::cerr << "Error, only one components ? Use clitkImageResample instead." << std::endl;
45 typedef typename ImageType::PixelType PixelType;
46 if (mNbOfComponents == 2) Update_WithDimAndPixelTypeAndComponent<ImageType::ImageDimension,PixelType,2>();
47 if (mNbOfComponents == 3) Update_WithDimAndPixelTypeAndComponent<ImageType::ImageDimension,PixelType,3>();
48 if (mNbOfComponents == 4) Update_WithDimAndPixelTypeAndComponent<ImageType::ImageDimension,PixelType,4>();
50 //--------------------------------------------------------------------
53 //--------------------------------------------------------------------
54 template<unsigned int Dim, class PixelType, unsigned int DimCompo>
55 void clitk::VFResampleGenericFilter::Update_WithDimAndPixelTypeAndComponent() {
57 typedef itk::Vector<PixelType, DimCompo> DisplacementType;
58 typedef itk::Image< DisplacementType, Dim > ImageType;
60 typename ImageType::Pointer input = clitk::readImage<ImageType>(mInputFilenames, mIOVerbose);
63 typename ImageType::Pointer outputImage = ComputeImage<ImageType>(input);
66 SetNextOutput<ImageType>(outputImage);
68 //--------------------------------------------------------------------
70 //--------------------------------------------------------------------
71 template<class ImageType>
72 typename ImageType::Pointer
73 clitk::VFResampleGenericFilter::ComputeImage(typename ImageType::Pointer inputImage) {
76 static unsigned int dim = ImageType::ImageDimension;
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;
93 typedef typename ImageType::SizeType SizeType;
94 typedef typename ImageType::SpacingType SpacingType;
95 typedef typename ImageType::PointType PointType;
97 // Create Image Filter
98 typedef itk::VectorResampleImageFilter<ImageType,ImageType> FilterType;
99 typename FilterType::Pointer filter = FilterType::New();
101 // Instance of the transform object to be passed to the resample
102 // filter. By default, identity transform is applied
103 typedef itk::AffineTransform<double, ImageType::ImageDimension> TransformType;
104 typename TransformType::Pointer transform = TransformType::New();
105 filter->SetTransform(transform);
107 // Set filter's parameters
109 SpacingType outputSpacing;
110 PointType outputOrigin;
111 for(unsigned int i=0; i<ImageType::ImageDimension; i++) {
112 outputSize[i] = mOutputSize[i];
113 outputSpacing[i] = mOutputSpacing[i];
114 outputOrigin[i] = inputImage->GetOrigin()[i];
117 filter->SetSize(outputSize);
118 filter->SetOutputSpacing(outputSpacing);
119 filter->SetOutputOrigin(outputOrigin);
120 filter->SetDefaultPixelValue(static_cast<typename ImageType::PixelType>(mDefaultPixelValue));
122 // Select interpolator
123 if (mInterpolatorName == "nn") {
124 typedef itk::VectorNearestNeighborInterpolateImageFunction<ImageType, double> InterpolatorType;
125 typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
126 filter->SetInterpolator(interpolator);
129 if (mInterpolatorName == "linear") {
130 typedef itk::VectorLinearInterpolateImageFunction<ImageType, double> InterpolatorType;
131 typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
132 filter->SetInterpolator(interpolator);
135 std::cerr << "Sorry, I do not know the interpolator (for vector field) '" << mInterpolatorName
136 << "'. Known interpolators are : nn, linear" << std::endl;
141 // Build initial Gaussian bluring (if needed)
142 typedef itk::RecursiveGaussianImageFilter<ImageType, ImageType> GaussianFilterType;
143 std::vector<typename GaussianFilterType::Pointer> gaussianFilters;
144 if (mApplyGaussianFilterBefore) {
145 for(unsigned int i=0; i<ImageType::ImageDimension; i++) {
147 gaussianFilters.push_back(GaussianFilterType::New());
149 gaussianFilters[i]->SetDirection(i);
150 gaussianFilters[i]->SetOrder(GaussianFilterType::ZeroOrder);
151 gaussianFilters[i]->SetNormalizeAcrossScale(false);
152 gaussianFilters[i]->SetSigma(mSigma[i]); // in millimeter !
154 if (i==0) gaussianFilters[i]->SetInput(inputImage);
155 else gaussianFilters[i]->SetInput(gaussianFilters[i-1]->GetOutput());
157 filter->SetInput(gaussianFilters[ImageType::ImageDimension-1]->GetOutput());
160 filter->SetInput(inputImage);
167 catch( itk::ExceptionObject & err ) {
168 std::cerr << "Error while filtering " << mInputFilenames[0].c_str()
169 << " " << err << std::endl;
174 return filter->GetOutput();
177 //--------------------------------------------------------------------
180 //--------------------------------------------------------------------
181 void clitk::VFResampleGenericFilter::SetOutputSize(const std::vector<int> & size) {
182 mOutputSize.resize(size.size());
183 std::copy(size.begin(), size.end(), mOutputSize.begin());
185 //--------------------------------------------------------------------
187 //--------------------------------------------------------------------
188 void clitk::VFResampleGenericFilter::SetOutputSpacing(const std::vector<double> & spacing) {
189 mOutputSpacing.resize(spacing.size());
190 std::copy(spacing.begin(), spacing.end(), mOutputSpacing.begin());
192 //--------------------------------------------------------------------
194 //--------------------------------------------------------------------
195 void clitk::VFResampleGenericFilter::SetInterpolationName(const std::string & inter) {
196 mInterpolatorName = inter;
198 //--------------------------------------------------------------------
200 //--------------------------------------------------------------------
201 void clitk::VFResampleGenericFilter::SetGaussianSigma(const std::vector<double> & sigma) {
202 mApplyGaussianFilterBefore = true;
203 mSigma.resize(sigma.size());
204 std::copy(sigma.begin(), sigma.end(), mSigma.begin());
206 //--------------------------------------------------------------------