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://www.centreleonberard.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 #include "clitkVFInterpolateGenericFilter.h"
22 #include "itkInterpolateImageFilter.h"
23 #include "itkLinearInterpolateImageFunction.h"
24 #include "itkNearestNeighborInterpolateImageFunction.h"
25 #include "itkNthElementImageAdaptor.h"
27 //--------------------------------------------------------------------
28 clitk::VFInterpolateGenericFilter::VFInterpolateGenericFilter():
29 clitk::ImageToImageGenericFilter<Self>("VFInterpolate")
31 //InitializeImageType<2>();
32 InitializeImageType<3>();
33 // InitializeImageType<4>();
34 mInterpolatorName = "nn";
38 //--------------------------------------------------------------------
41 //--------------------------------------------------------------------
42 template<unsigned int Dim>
43 void clitk::VFInterpolateGenericFilter::InitializeImageType()
45 //typedef itk::Vector<float,Dim> v3f;
46 //ADD_IMAGE_TYPE(Dim, v3f);
47 // ADD_VEC_IMAGE_TYPE(Dim,2,double)
48 ADD_VEC_IMAGE_TYPE(Dim,3,float)
49 ADD_VEC_IMAGE_TYPE(Dim,3,double)
52 //--------------------------------------------------------------------
55 //--------------------------------------------------------------------
56 template<class ImageType>
57 void clitk::VFInterpolateGenericFilter::UpdateWithInputImageType()
60 if (m_NbOfComponents == 1) {
61 std::cerr << "Error, only one components ? Use clitkImageInterpolate instead." << std::endl;
64 typedef typename ImageType::PixelType PixelType;
65 // if (m_NbOfComponents == 2) Update_WithDimAndPixelTypeAndComponent<ImageType::ImageDimension,PixelType,2>();
66 if (m_NbOfComponents == 3) Update_WithDimAndPixelTypeAndComponent<ImageType::ImageDimension,PixelType,3>();
68 //--------------------------------------------------------------------
71 //--------------------------------------------------------------------
72 template<unsigned int Dim, class PixelType, unsigned int DimCompo>
73 void clitk::VFInterpolateGenericFilter::Update_WithDimAndPixelTypeAndComponent()
76 // typedef itk::Vector<PixelType, DimCompo> DisplacementType;
77 typedef PixelType DisplacementType;
78 typedef itk::Image< DisplacementType, Dim > ImageType;
80 typename ImageType::Pointer input1 = clitk::readImage<ImageType>(m_InputFilenames[0], m_IOVerbose);
81 typename ImageType::Pointer input2 = clitk::readImage<ImageType>(mInputFilename2, m_IOVerbose);
84 typename ImageType::Pointer outputImage = ComputeImage<ImageType>(input1, input2);
87 SetNextOutput<ImageType>(outputImage);
89 //--------------------------------------------------------------------
91 //--------------------------------------------------------------------
92 template<class ImageType>
93 typename ImageType::Pointer
94 clitk::VFInterpolateGenericFilter::ComputeImage(typename ImageType::Pointer inputImage1, typename ImageType::Pointer inputImage2)
98 typedef itk::Image<typename ImageType::PixelType::ValueType, ImageType::ImageDimension> ScalarImageType;
99 typedef itk::Image<typename ImageType::PixelType::ValueType, ImageType::ImageDimension + 1> InterpolationImageType;
100 typedef itk::NthElementImageAdaptor<ImageType, typename ImageType::PixelType::ValueType> ImageAdaptorType;
101 typename ImageAdaptorType::Pointer adaptor1 = ImageAdaptorType::New();
102 typename ImageAdaptorType::Pointer adaptor2 = ImageAdaptorType::New();
103 adaptor1->SetImage(inputImage1);
104 adaptor2->SetImage(inputImage2);
106 // Create Image Filter
107 typedef itk::InterpolateImageFilter<ImageAdaptorType, ScalarImageType> FilterType;
108 typename FilterType::Pointer filter = FilterType::New();
109 filter->SetInput1(adaptor1);
110 filter->SetInput2(adaptor2);
111 filter->SetDistance(mDistance);
113 // Select interpolator
114 if (mInterpolatorName == "nn") {
115 typedef itk::NearestNeighborInterpolateImageFunction<InterpolationImageType> InterpolatorType;
116 typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
117 filter->SetInterpolator(interpolator);
119 if (mInterpolatorName == "linear") {
120 typedef itk::LinearInterpolateImageFunction<InterpolationImageType> InterpolatorType;
121 typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
122 filter->SetInterpolator(interpolator);
124 std::cerr << "Sorry, I do not know the interpolator (for vector field) '" << mInterpolatorName
125 << "'. Known interpolators are : nn, linear" << std::endl;
130 typename ImageType::Pointer output = ImageType::New();
131 typename ImageAdaptorType::Pointer adaptorOutput = ImageAdaptorType::New();
132 output->CopyInformation(inputImage1);
133 output->SetRegions(inputImage1->GetLargestPossibleRegion());
136 typedef itk::ImageRegionIterator<ScalarImageType> IteratorType1;
137 typedef itk::ImageRegionIterator<ImageAdaptorType> IteratorType2;
139 for (unsigned int i = 0; i < ImageType::PixelType::Dimension; i++) {
140 adaptor1->SelectNthElement(i);
141 adaptor2->SelectNthElement(i);
146 } catch( itk::ExceptionObject & err ) {
147 std::cerr << "Error while filtering " << m_InputFilenames[0].c_str()
148 << " " << err << std::endl;
152 adaptorOutput->SelectNthElement(i);
153 adaptorOutput->SetImage(output);
155 IteratorType1 it1(filter->GetOutput(), filter->GetOutput()->GetLargestPossibleRegion());
156 IteratorType2 it2(adaptorOutput, adaptorOutput->GetLargestPossibleRegion());
160 while ( ! it1.IsAtEnd() ) {
171 //--------------------------------------------------------------------
173 //--------------------------------------------------------------------
174 void clitk::VFInterpolateGenericFilter::SetInterpolationName(const std::string & inter)
176 mInterpolatorName = inter;
178 //--------------------------------------------------------------------