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 clitkPointTrajectoryGenericFilter_cxx
19 #define clitkPointTrajectoryGenericFilter_cxx
21 /* =================================================
22 * @file clitkPointTrajectoryGenericFilter.cxx
28 ===================================================*/
30 #include "clitkPointTrajectoryGenericFilter.h"
37 //-----------------------------------------------------------
39 //-----------------------------------------------------------
40 PointTrajectoryGenericFilter::PointTrajectoryGenericFilter()
43 // m_InputFileName="";
47 //-----------------------------------------------------------
49 //-----------------------------------------------------------
50 void PointTrajectoryGenericFilter::Update()
53 const unsigned int ImageDimension=4;
54 const unsigned int SpaceDimension=3;
55 typedef itk::Vector<double, SpaceDimension> CoefficientPixelType;
56 typedef itk::Vector<float, SpaceDimension> VectorPixelType;
57 typedef itk::Image<CoefficientPixelType, ImageDimension> CoefficientImageType;
58 typedef itk::Image<VectorPixelType, ImageDimension> VectorFieldType;
61 // -----------------------------------------------
62 // Reference Point List 3D
63 // -----------------------------------------------
64 typedef itk::Point<double, SpaceDimension> SpacePointType;
65 typedef clitk::List<SpacePointType> PointListType;
66 PointListType referencePointList;
67 if (m_Verbose) std::cout<<"Reference point list:"<<std::endl;
68 referencePointList=PointListType(m_ArgsInfo.ref_arg, m_Verbose);
71 // -----------------------------------------------
72 // Transform: based on points, 4DVF, spatio-Temporal transform
73 // -----------------------------------------------
74 typedef itk::Transform<double, 4, 4> TransformType;
75 TransformType::Pointer transform;
76 switch (m_ArgsInfo.transform_arg)
78 // ==========================
80 // ==========================
83 //-----------------------------
85 //-----------------------------
86 typedef itk::Point<double, SpaceDimension> PointType;
87 typedef clitk::List<PointType> PointListType;
88 typedef clitk::Lists<PointType> PointListsType;
89 PointListsType inputPointLists, sortedPointLists;
92 for (unsigned int i=0; i<m_ArgsInfo.points_given; i++)
93 inputPointLists.push_back(PointListType(m_ArgsInfo.points_arg[i], m_Verbose) );
95 // Convert/sort the lists
96 sortedPointLists.resize(inputPointLists[0].size());
97 for (unsigned int i=0; i<inputPointLists[0].size(); i++)
99 sortedPointLists[i].push_back(referencePointList[i]);
100 for (unsigned int j=0; j<inputPointLists.size(); j++)
102 sortedPointLists[i].push_back(inputPointLists[j][i]);
106 // Point list Transform
107 typedef clitk::PointListTransform<double, ImageDimension,ImageDimension> PointListTransformType;
108 PointListTransformType::Pointer pointListTransform=PointListTransformType::New();
109 pointListTransform->SetPointLists(sortedPointLists);
111 // Vector Interpolator
112 typedef PointListTransformType::PointListImageType PointListImageType;
113 typedef clitk::GenericVectorInterpolator<args_info_clitkPointTrajectory,PointListImageType, double> GenericVectorInterpolatorType;
114 GenericVectorInterpolatorType::Pointer genericInterpolator=GenericVectorInterpolatorType::New();
115 genericInterpolator->SetArgsInfo(m_ArgsInfo);
116 typedef itk::VectorInterpolateImageFunction<PointListImageType, double> InterpolatorType;
117 InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer();
118 pointListTransform->SetInterpolator(interpolator);
119 transform=pointListTransform;
124 // ==========================
126 // ==========================
129 // Deformation field transform
130 typedef clitk::DeformationFieldTransform<double, ImageDimension,ImageDimension, SpaceDimension> DeformationFieldTransformType;
131 DeformationFieldTransformType::Pointer deformationFieldTransform=DeformationFieldTransformType::New();
133 // The deformation field
134 typedef DeformationFieldTransformType::DeformationFieldType DeformationFieldType;
135 typedef itk::ImageFileReader<DeformationFieldType> InputReaderType;
136 InputReaderType::Pointer reader = InputReaderType::New();
137 reader->SetFileName( m_ArgsInfo.input_arg);
139 DeformationFieldType::Pointer input= reader->GetOutput();
140 deformationFieldTransform->SetDeformationField(input);
142 // Vector Interpolator
143 typedef clitk::GenericVectorInterpolator<args_info_clitkPointTrajectory,DeformationFieldType, double> GenericVectorInterpolatorType;
144 GenericVectorInterpolatorType::Pointer genericInterpolator=GenericVectorInterpolatorType::New();
145 genericInterpolator->SetArgsInfo(m_ArgsInfo);
146 typedef itk::VectorInterpolateImageFunction<DeformationFieldType, double> InterpolatorType;
147 InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer();
148 deformationFieldTransform->SetInterpolator(interpolator);
149 transform=deformationFieldTransform;
154 // ==========================
155 // Spatio-Temporal transform
156 // ==========================
160 typedef clitk::ShapedBLUTSpatioTemporalDeformableTransform< double, ImageDimension, ImageDimension > TransformType;
161 TransformType::Pointer spatioTemporalTransform = TransformType::New();
164 // Spline orders: Default is cubic splines
165 CoefficientImageType::RegionType::SizeType splineOrders ;
166 splineOrders.Fill(3);
167 if (m_ArgsInfo.order_given)
168 for(unsigned int i=0; i<ImageDimension;i++)
169 splineOrders[i]=m_ArgsInfo.order_arg[i];
170 if (m_Verbose) std::cout<<"Setting the spline orders to "<<splineOrders<<"..."<<std::endl;
173 typedef itk::ImageFileReader<CoefficientImageType> InputReaderType;
174 InputReaderType::Pointer reader = InputReaderType::New();
175 reader->SetFileName( m_ArgsInfo.input_arg);
177 CoefficientImageType::Pointer input= reader->GetOutput();
178 // itk::Vector<double,3> vector;
181 // input->FillBuffer(vector);
184 typedef itk::ImageMaskSpatialObject< ImageDimension > MaskType;
185 MaskType::Pointer spatialObjectMask=NULL;
186 if (m_ArgsInfo.mask_given)
188 typedef itk::Image< unsigned char, ImageDimension > ImageMaskType;
189 typedef itk::ImageFileReader< ImageMaskType > MaskReaderType;
190 MaskReaderType::Pointer maskReader = MaskReaderType::New();
191 maskReader->SetFileName(m_ArgsInfo.mask_arg);
195 maskReader->Update();
197 catch( itk::ExceptionObject & err )
199 std::cerr << "ExceptionObject caught while reading mask !" << std::endl;
200 std::cerr << err << std::endl;
203 if (m_Verbose)std::cout <<"Mask was read..." <<std::endl;
205 // Set the image to the spatialObject
206 spatialObjectMask = MaskType::New();
207 spatialObjectMask->SetImage( maskReader->GetOutput() );
211 CoefficientImageType::SizeType samplingFactors;
212 for (unsigned int i=0; i< ImageDimension-1; i++)
214 samplingFactors[i]= (int) ( input->GetSpacing()[i]/ m_ArgsInfo.spacing_arg);
215 if (m_Verbose) std::cout<<"Setting sampling factor "<<i<<" to "<<samplingFactors[i]<<"..."<<std::endl;
217 samplingFactors[ImageDimension-1]= (int) ( input->GetSpacing()[ImageDimension-1]/ m_ArgsInfo.phaseIncrement_arg);
218 if (m_Verbose) std::cout<<"Setting sampling factor "<<ImageDimension-1<<" to "<<samplingFactors[ImageDimension-1]<<"..."<<std::endl;
219 if( (m_ArgsInfo.shape_arg==3) |
220 (m_ArgsInfo.shape_arg==4) |
221 (m_ArgsInfo.shape_arg==6) |
222 (m_ArgsInfo.shape_arg==8)
223 ) samplingFactors[ImageDimension-1]*=2.5;
226 spatioTemporalTransform->SetTransformShape(m_ArgsInfo.shape_arg);
227 spatioTemporalTransform->SetSplineOrders(splineOrders);
228 spatioTemporalTransform->SetMask(spatialObjectMask);
229 spatioTemporalTransform->SetLUTSamplingFactors(samplingFactors);
230 spatioTemporalTransform->SetCoefficientImage(input);
231 transform=spatioTemporalTransform;
236 // // ==========================
237 // // Spatio-Temporal transform
238 // // ==========================
242 // typedef clitk::BSplineSpatioTemporalDeformableTransform< double, ImageDimension, ImageDimension > TransformType;
243 // TransformType::Pointer spatioTemporalTransform = TransformType::New();
246 // // Spline orders: Default is cubic splines
247 // CoefficientImageType::RegionType::SizeType splineOrders ;
248 // splineOrders.Fill(3);
249 // if (m_ArgsInfo.order_given)
250 // for(unsigned int i=0; i<ImageDimension;i++)
251 // splineOrders[i]=m_ArgsInfo.order_arg[i];
252 // if (m_Verbose) std::cout<<"Setting the spline orders to "<<splineOrders<<"..."<<std::endl;
254 // // Coefficient image
255 // typedef itk::ImageFileReader<CoefficientImageType> InputReaderType;
256 // InputReaderType::Pointer reader = InputReaderType::New();
257 // reader->SetFileName( m_ArgsInfo.input_arg);
259 // CoefficientImageType::Pointer input= reader->GetOutput();
260 // // itk::Vector<double,3> vector;
261 // // vector.Fill(0.);
263 // // input->FillBuffer(vector);
266 // typedef itk::ImageMaskSpatialObject< ImageDimension > MaskType;
267 // MaskType::Pointer spatialObjectMask=NULL;
268 // if (m_ArgsInfo.mask_given)
270 // typedef itk::Image< unsigned char, ImageDimension > ImageMaskType;
271 // typedef itk::ImageFileReader< ImageMaskType > MaskReaderType;
272 // MaskReaderType::Pointer maskReader = MaskReaderType::New();
273 // maskReader->SetFileName(m_ArgsInfo.mask_arg);
277 // maskReader->Update();
279 // catch( itk::ExceptionObject & err )
281 // std::cerr << "ExceptionObject caught while reading mask !" << std::endl;
282 // std::cerr << err << std::endl;
285 // if (m_Verbose)std::cout <<"Mask was read..." <<std::endl;
287 // // Set the image to the spatialObject
288 // spatialObjectMask = MaskType::New();
289 // spatialObjectMask->SetImage( maskReader->GetOutput() );
292 // // Samplingfactors
293 // CoefficientImageType::SizeType samplingFactors;
294 // for (unsigned int i=0; i< ImageDimension-1; i++)
296 // samplingFactors[i]= (int) ( input->GetSpacing()[i]/ m_ArgsInfo.spacing_arg);
297 // if (m_Verbose) std::cout<<"Setting sampling factor "<<i<<" to "<<samplingFactors[i]<<"..."<<std::endl;
299 // samplingFactors[ImageDimension-1]= (int) ( input->GetSpacing()[ImageDimension-1]/ m_ArgsInfo.phaseIncrement_arg);
300 // if (m_Verbose) std::cout<<"Setting sampling factor "<<ImageDimension-1<<" to "<<samplingFactors[ImageDimension-1]<<"..."<<std::endl;
303 // //spatioTemporalTransform->SetTransformShape(m_ArgsInfo.shape_arg);
304 // spatioTemporalTransform->SetSplineOrders(splineOrders);
305 // spatioTemporalTransform->SetMask(spatialObjectMask);
306 // spatioTemporalTransform->SetLUTSamplingFactors(samplingFactors);
307 // spatioTemporalTransform->SetCoefficientImage(input);
308 // transform=spatioTemporalTransform;
315 // -----------------------------------------------
316 // Construct Spatio-Temporal Point lists 4D
317 // -----------------------------------------------
318 typedef itk::Point<double, ImageDimension> SpaceTimePointType;
319 typedef clitk::Lists<SpaceTimePointType> PointListsType;
320 PointListsType pointLists(referencePointList.size());
321 SpaceTimePointType spaceTimePoint;
323 for (unsigned int i=0; i<referencePointList.size(); i++)
325 for (unsigned int d=0; d<SpaceDimension; d++)
327 spaceTimePoint[d]=referencePointList[i][d];
332 spaceTimePoint[ImageDimension-1]=phase;
333 pointLists[i].push_back(spaceTimePoint);
334 phase+=m_ArgsInfo.phaseIncrement_arg;
339 // -----------------------------------------------
341 // -----------------------------------------------
342 typedef itk::Vector<double, ImageDimension> VectorType;
343 typedef clitk::List<VectorType> VectorListType;
344 typedef clitk::Lists<VectorType> VectorListsType;
345 VectorListsType displacementLists(pointLists.size());
346 VectorType displacement;
347 for (unsigned int i=0; i<pointLists.size(); i++)
349 if (m_Verbose) std::cout<<"Transforming point "<<pointLists[i][0]<<"..."<<std::endl;
350 for (unsigned int j=0; j<pointLists[i].size(); j++)
352 spaceTimePoint= transform->TransformPoint(pointLists[i][j]);
353 if (m_Verbose) std::cout<<"Transformed point "<<spaceTimePoint<<"..."<<std::endl;
354 displacement=spaceTimePoint-pointLists[i][j];
355 displacementLists[i].push_back(displacement);
360 // -----------------------------------------------
361 // Write displacements
362 // -----------------------------------------------
363 std::vector<std::string> filenames;
364 for (unsigned int i=0;i<displacementLists.size();i++)
366 std::ostringstream number_ostr;
368 std::string number_str;
369 if (i<10) number_str= "0"+number_ostr.str();
370 else number_str = number_ostr.str();
371 filenames.push_back(m_ArgsInfo.trajectory_arg+number_str);
373 displacementLists.Write(filenames, m_Verbose );
380 #endif //#define clitkPointTrajectoryGenericFilter_cxx