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 clitkCalculateTREGenericFilter_txx
19 #define clitkCalculateTREGenericFilter_txx
22 #include "clitkBSplineDeformableTransform.h"
25 /* =================================================
26 * @file clitkCalculateTREGenericFilter.txx
32 ===================================================*/
35 //=====================================================================================
37 // Rômulo Pinho - 25/02/2011
38 // The macros below avoid the duplication of code after introducing the function
39 // UpdateCoeffsWithDim(). Although the code generated by these functions will be
40 // duplicated anyway (since they're template functions) the macros make their sources
42 // Ideally, these macros should be made member functions, with the typedefs inside the
43 // the filter class, but that would change the class' interface, since it would need
44 // to become a template class (with parameter Dimension). To avoid breaking legacy code,
45 // I opted for the macros, although I don't quite like them (especially long ones like
46 // these). In any case, it's advisable to think of a better solution, even if it breaks
47 // the class' interface.
48 //=====================================================================================
50 //-----------------------------
52 //-----------------------------
53 #define TypedefsMacro \
55 typedef itk::Point<ValueType, Dimension> PointType; \
56 typedef clitk::List<PointType> PointListType; \
57 typedef clitk::Lists<PointType> PointListsType; \
59 typedef itk::Vector<ValueType, Dimension> VectorType; \
60 typedef clitk::List<VectorType> VectorListType; \
61 typedef clitk::Lists<VectorType> VectorListsType; \
63 typedef itk::Vector<ValueType, Dimension> DeformationVectorType; \
64 typedef itk::Image<DeformationVectorType, Dimension> DeformationFieldType;
66 //-------------------------------------------------------------------
67 // Build the point lists necessary for the computation of distances.
68 //-------------------------------------------------------------------
69 #define BuildPointListsMacro(filenames) \
71 m_NumberOfFields=filenames.size(); \
72 m_NumberOfLists=m_ArgsInfo.input_given; \
73 if ( (m_NumberOfLists!=m_NumberOfFields) && (m_NumberOfLists!=1) ) \
75 std::cerr<<"Error: Number of lists (="<<m_NumberOfLists<<") different from number of DVF's (="<<m_NumberOfFields<<") !"<<std::endl; \
78 else if (m_NumberOfLists==1) \
80 if (m_Verbose) std::cout<<m_NumberOfFields<<" DVFs and 1 point list given..."<<std::endl; \
84 if (m_Verbose) std::cout<<m_NumberOfLists<<" point lists and DVFs given..."<<std::endl; \
88 PointListsType pointLists; \
90 for (unsigned int i=0; i<m_NumberOfFields; i++) \
93 if (m_NumberOfLists==1) \
94 pointLists.push_back(PointListType(m_ArgsInfo.input_arg[0], m_Verbose) ); \
96 pointLists.push_back(PointListType(m_ArgsInfo.input_arg[i], m_Verbose) ); \
99 /* Verify the number of points */\
100 if (i==0) m_NumberOfPoints=pointLists[i].size(); \
103 if (m_NumberOfPoints!=pointLists[i].size()) \
105 std::cerr<<"Size of first list ("<<m_NumberOfPoints \
106 <<") is different from size of list "<<i \
107 <<" ("<<pointLists[i].size()<<")..."<<std::endl; \
114 PointListType referencePointList; \
115 if (m_Verbose) std::cout<<"Reference point list:"<<std::endl; \
116 referencePointList=PointListType(m_ArgsInfo.ref_arg, m_Verbose); \
117 if (m_NumberOfPoints!=referencePointList.size()) \
119 std::cerr<<"Size of the first list ("<<m_NumberOfPoints \
120 <<") is different from size of the reference list (" \
121 << referencePointList.size() <<")..."<<std::endl; \
125 /* //===================================================================================== */\
126 /* // Original distance between points */\
127 /* //===================================================================================== */\
128 #define OriginalDistancesMacro \
129 VectorListsType originalDistanceLists(m_NumberOfFields); \
131 /* //----------------------------- */\
132 /* // Calculate original distances */\
133 /* //----------------------------- */\
134 PointType referencePoint; \
135 VectorType distance; \
136 for (unsigned int phaseIndex=0; phaseIndex<m_NumberOfFields; phaseIndex++) \
138 for (unsigned int pointIndex=0; pointIndex<m_NumberOfPoints; pointIndex++) \
140 referencePoint=referencePointList[pointIndex]; \
141 distance=pointLists[phaseIndex][pointIndex]-referencePoint; \
142 originalDistanceLists[phaseIndex].push_back(distance); \
147 /* //----------------------------- */\
148 /* // Statistics Original */\
149 /* //----------------------------- */\
150 typedef DeformationListStatisticsFilter<VectorType> StatisticsFilterType; \
151 typename StatisticsFilterType::Pointer statisticsFilter=StatisticsFilterType::New(); \
153 /* // Statistics (magnitude) */ \
154 MeasureListType oMeanList, oStdList, oMaxList; \
155 ValueType oMean, oStd, oMax; \
156 statisticsFilter->GetStatistics(originalDistanceLists, oMean, oStd, oMax, oMeanList, oStdList, oMaxList); \
158 /* // Statistics (per component) */\
159 VectorListType oMeanXYZList, oStdXYZList,oMaxXYZList; \
160 VectorType oMeanXYZ, oStdXYZ, oMaxXYZ; \
161 statisticsFilter->GetStatistics(originalDistanceLists, oMeanXYZ, oStdXYZ, oMaxXYZ, oMeanXYZList, oStdXYZList, oMaxXYZList); \
164 /* //----------------------------- */\
166 /* //----------------------------- */\
167 std::vector<std::string> labels; \
168 labels.push_back("MeanX\tSDX"); \
169 labels.push_back("MeanY\tSDY"); \
170 labels.push_back("MeanZ\tSDZ"); \
171 labels.push_back("MeanT\tSDT"); \
172 labels.push_back("MaxX"); \
173 labels.push_back("MaxY"); \
174 labels.push_back("MaxZ"); \
175 labels.push_back("MaxT"); \
178 /* // Output to screen */\
182 /* // Numbers of DVF */\
183 std::cout<<"# Number\tDVF"<<std::endl; \
184 for (unsigned int phaseIndex=0; phaseIndex<m_NumberOfFields; phaseIndex++) \
185 std::cout<<phaseIndex<<"\t"<<filenames[phaseIndex]<<std::endl; \
188 std::cout<<std::endl; \
189 std::cout<<"=================================================="<<std::endl; \
190 std::cout<<"|| Original distance between points ||"<<std::endl; \
191 std::cout<<"=================================================="<<std::endl; \
192 std::cout<<std::endl; \
194 std::cout<<std::setprecision(3); \
198 /* // Labels of the columns */\
199 std::cout<<"#DVF\tMean\tSD\t"<<labels[0]; \
200 for(unsigned int dim=1;dim<Dimension;dim++) std::cout<<"\t"<< labels[dim]; \
201 std::cout<<"\tMax \t"<<labels[4]; \
202 for(unsigned int dim=1;dim<Dimension;dim++) std::cout<<"\t"<< labels[dim+4]; \
203 std::cout<<std::endl; \
204 for (unsigned int phaseIndex=0; phaseIndex<m_NumberOfFields; phaseIndex++) \
206 std::cout<<phaseIndex; \
207 std::cout<<"\t"<<oMeanList[phaseIndex]<<"\t"<<oStdList[phaseIndex]; \
208 std::cout<<"\t"<<oMeanXYZList[phaseIndex][0]<<"\t"<<oStdXYZList[phaseIndex][0]; \
209 for(unsigned int dim=1;dim<Dimension;dim++) \
210 std::cout<<"\t"<<oMeanXYZList[phaseIndex][dim]<<"\t"<<oStdXYZList[phaseIndex][dim]; \
211 std::cout<<"\t"<<oMaxList[phaseIndex]; \
212 std::cout<<"\t"<<oMaxXYZList[phaseIndex][0]; \
213 for(unsigned int dim=1;dim<Dimension;dim++) \
214 std::cout<<"\t"<<oMaxXYZList[phaseIndex][dim]; \
215 std::cout<<std::endl; \
219 std::cout<<std::endl; \
221 std::cout<<"\t"<<oMean<<"\t"<<oStd; \
222 std::cout<<"\t"<<oMeanXYZ[0]<<"\t"<<oStdXYZ[0]; \
223 for(unsigned int dim=1;dim<Dimension;dim++) \
224 std::cout<<"\t"<<oMeanXYZ[dim]<<"\t"<<oStdXYZ[dim]; \
225 std::cout<<"\t"<<oMax; \
226 std::cout<<"\t"<<oMaxXYZ[0]; \
227 for(unsigned int dim=1;dim<Dimension;dim++) \
228 std::cout<<"\t"<<oMaxXYZ[dim]; \
229 std::cout<<std::endl; \
232 /* // Output to file */ \
233 if( m_ArgsInfo.general_given) \
235 /* // Add to the file */\
236 std::ofstream general(m_ArgsInfo.general_arg); \
239 /* // Numbers of DVF */\
240 for (unsigned int phaseIndex=0; phaseIndex<m_NumberOfFields; phaseIndex++) \
241 general<<phaseIndex<<"\t"<< filenames[phaseIndex]<<std::endl; \
243 general<<std::endl; \
244 general<<"=================================================="<<std::endl; \
245 general<<"|| Original distance between points ||"<<std::endl; \
246 general<<"=================================================="<<std::endl; \
247 general<<std::endl; \
249 general<<std::setprecision(3); \
251 /* // Labels of the columns */\
252 general<<"#DVF\tMean\tSD\t"<<labels[0]; \
253 for(unsigned int dim=1;dim<Dimension;dim++) general<<"\t"<< labels[dim]; \
254 general<<"\tMax \t"<<labels[4]; \
255 for(unsigned int dim=1;dim<Dimension;dim++) general<<"\t"<< labels[dim+4]; \
256 general<<std::endl; \
257 for (unsigned int phaseIndex=0; phaseIndex<m_NumberOfFields; phaseIndex++) \
259 general<<phaseIndex; \
260 general<<"\t"<<oMeanList[phaseIndex]<<"\t"<<oStdList[phaseIndex]; \
261 general<<"\t"<<oMeanXYZList[phaseIndex][0]<<"\t"<<oStdXYZList[phaseIndex][0]; \
262 for(unsigned int dim=1;dim<Dimension;dim++) \
263 general<<"\t"<<oMeanXYZList[phaseIndex][dim]<<"\t"<<oStdXYZList[phaseIndex][dim]; \
264 general<<"\t"<<oMaxList[phaseIndex]; \
265 general<<"\t"<<oMaxXYZList[phaseIndex][0]; \
266 for(unsigned int dim=1;dim<Dimension;dim++) \
267 general<<"\t"<<oMaxXYZList[phaseIndex][dim]; \
268 general<<std::endl; \
272 general<<std::endl; \
274 general<<"\t"<<oMean<<"\t"<<oStd; \
275 general<<"\t"<<oMeanXYZ[0]<<"\t"<<oStdXYZ[0]; \
276 for(unsigned int dim=1;dim<Dimension;dim++) \
277 general<<"\t"<<oMeanXYZ[dim]<<"\t"<<oStdXYZ[dim]; \
278 general<<"\t"<<oMax; \
279 general<<"\t"<<oMaxXYZ[0]; \
280 for(unsigned int dim=1;dim<Dimension;dim++) \
281 general<<"\t"<<oMaxXYZ[dim]; \
282 general<<std::endl; \
287 #define StatisticsDisplacementsMacro \
288 /* //----------------------------- */\
289 /* // Statistics displacements */\
290 /* //----------------------------- */\
292 /* // Statistics (magnitude) */\
293 MeasureListType dmeanList, dstdList, dmaxList; \
294 ValueType dmean, dstd, dmax; \
295 statisticsFilter->GetStatistics(displacementLists, dmean, dstd, dmax, dmeanList, dstdList, dmaxList); \
297 /* // Statistics (per component) */\
298 VectorListType dmeanXYZList, dstdXYZList,dmaxXYZList; \
299 VectorType dmeanXYZ, dstdXYZ, dmaxXYZ; \
300 statisticsFilter->GetStatistics(displacementLists, dmeanXYZ, dstdXYZ, dmaxXYZ,dmeanXYZList, dstdXYZList, dmaxXYZList); \
303 /* //----------------------------- */\
304 /* // Statistics TRE */\
305 /* //----------------------------- */\
307 /* // Statistics (magnitude) */\
308 MeasureListType meanList, stdList, maxList; \
309 ValueType mean, std, max; \
310 statisticsFilter->GetStatistics(treLists, mean, std, max, meanList, stdList, maxList); \
312 /* // Statistics (per component) */\
313 VectorListType meanXYZList, stdXYZList,maxXYZList; \
314 VectorType meanXYZ, stdXYZ, maxXYZ; \
315 statisticsFilter->GetStatistics(treLists, meanXYZ, stdXYZ, maxXYZ, meanXYZList, stdXYZList, maxXYZList); \
318 /* // Output to screen */\
322 std::cout<<std::endl; \
323 std::cout<<"=================================================="<<std::endl; \
324 std::cout<<"|| Residual distance between points ||"<<std::endl; \
325 std::cout<<"=================================================="<<std::endl; \
326 std::cout<<std::endl; \
328 std::cout<<std::setprecision(3); \
330 /* // Labels of the columns */\
331 std::cout<<"#DVF\tMean\tSD\t"<<labels[0]; \
332 for(unsigned int dim=1;dim<Dimension;dim++) std::cout<<"\t"<< labels[dim]; \
333 std::cout<<"\tMax \t"<<labels[4]; \
334 for(unsigned int dim=1;dim<Dimension;dim++) std::cout<<"\t"<< labels[dim+4]; \
335 std::cout<<std::endl; \
336 for (unsigned int phaseIndex=0; phaseIndex<m_NumberOfFields; phaseIndex++) \
338 std::cout<<phaseIndex; \
339 std::cout<<"\t"<<meanList[phaseIndex]<<"\t"<<stdList[phaseIndex]; \
340 std::cout<<"\t"<<meanXYZList[phaseIndex][0]<<"\t"<<stdXYZList[phaseIndex][0]; \
341 for(unsigned int dim=1;dim<Dimension;dim++) \
342 std::cout<<"\t"<<meanXYZList[phaseIndex][dim]<<"\t"<<stdXYZList[phaseIndex][dim]; \
343 std::cout<<"\t"<<maxList[phaseIndex]; \
344 std::cout<<"\t"<<maxXYZList[phaseIndex][0]; \
345 for(unsigned int dim=1;dim<Dimension;dim++) \
346 std::cout<<"\t"<<maxXYZList[phaseIndex][dim]; \
347 std::cout<<std::endl; \
351 std::cout<<std::endl; \
353 std::cout<<"\t"<<mean<<"\t"<<std; \
354 std::cout<<"\t"<<meanXYZ[0]<<"\t"<<stdXYZ[0]; \
355 for(unsigned int dim=1;dim<Dimension;dim++) \
356 std::cout<<"\t"<<meanXYZ[dim]<<"\t"<<stdXYZ[dim]; \
357 std::cout<<"\t"<<max; \
358 std::cout<<"\t"<<maxXYZ[0]; \
359 for(unsigned int dim=1;dim<Dimension;dim++) \
360 std::cout<<"\t"<<maxXYZ[dim]; \
361 std::cout<<std::endl; \
364 /* // Output to file */\
365 if( m_ArgsInfo.general_given) \
367 /* // Add to the file */\
368 std::ofstream general(m_ArgsInfo.general_arg, ios_base::app); \
370 general<<std::endl; \
371 general<<"=================================================="<<std::endl; \
372 general<<"|| Residual distance between points ||"<<std::endl; \
373 general<<"=================================================="<<std::endl; \
374 general<<std::endl; \
376 general<<std::setprecision(3); \
378 /* // Labels of the columns */\
379 general<<"#DVF\tMean\tSD\t"<<labels[0]; \
380 for(unsigned int dim=1;dim<Dimension;dim++) general<<"\t"<< labels[dim]; \
381 general<<"\tMax \t"<<labels[4]; \
382 for(unsigned int dim=1;dim<Dimension;dim++) general<<"\t"<< labels[dim+4]; \
383 general<<std::endl; \
384 for (unsigned int phaseIndex=0; phaseIndex<m_NumberOfFields; phaseIndex++) \
386 general<<phaseIndex; \
387 general<<"\t"<<meanList[phaseIndex]<<"\t"<<stdList[phaseIndex]; \
388 general<<"\t"<<meanXYZList[phaseIndex][0]<<"\t"<<stdXYZList[phaseIndex][0]; \
389 for(unsigned int dim=1;dim<Dimension;dim++) \
390 general<<"\t"<<meanXYZList[phaseIndex][dim]<<"\t"<<stdXYZList[phaseIndex][dim]; \
391 general<<"\t"<<maxList[phaseIndex]; \
392 general<<"\t"<<maxXYZList[phaseIndex][0]; \
393 for(unsigned int dim=1;dim<Dimension;dim++) \
394 general<<"\t"<<maxXYZList[phaseIndex][dim]; \
395 general<<std::endl; \
399 general<<std::endl; \
401 general<<"\t"<<mean<<"\t"<<std; \
402 general<<"\t"<<meanXYZ[0]<<"\t"<<stdXYZ[0]; \
403 for(unsigned int dim=1;dim<Dimension;dim++) \
404 general<<"\t"<<meanXYZ[dim]<<"\t"<<stdXYZ[dim]; \
405 general<<"\t"<<max; \
406 general<<"\t"<<maxXYZ[0]; \
407 for(unsigned int dim=1;dim<Dimension;dim++) \
408 general<<"\t"<<maxXYZ[dim]; \
409 general<<std::endl; \
413 /* // Output original points */\
414 if( m_ArgsInfo.original_given) \
417 std::vector<std::string> filenames; \
418 for (unsigned int i=0;i<m_NumberOfFields;i++) \
420 std::ostringstream number_ostr; \
422 std::string number_str = number_ostr.str(); \
423 filenames.push_back(m_ArgsInfo.original_arg+number_str); \
425 originalDistanceLists.Write(filenames, m_Verbose ); \
429 /* // Output original magnitude points */\
430 if( m_ArgsInfo.originalMag_given) \
432 clitk::Lists<itk::FixedArray<ValueType,1> > originalDistanceListsMag=originalDistanceLists.Norm(); \
433 std::vector<std::string> filenames; \
434 for (unsigned int i=0;i<m_NumberOfFields;i++) \
436 std::ostringstream number_ostr; \
438 std::string number_str = number_ostr.str(); \
439 filenames.push_back(m_ArgsInfo.originalMag_arg+number_str); \
441 originalDistanceListsMag.Write(filenames, m_Verbose ); \
444 /* // Output displacement */\
445 if( m_ArgsInfo.displacement_given) \
448 std::vector<std::string> filenames; \
449 for (unsigned int i=0;i<m_NumberOfFields;i++) \
451 std::ostringstream number_ostr; \
453 std::string number_str = number_ostr.str(); \
454 filenames.push_back(m_ArgsInfo.displacement_arg+number_str); \
456 displacementLists.Write(filenames, m_Verbose ); \
460 /* // Output displacement magnitude */\
461 if( m_ArgsInfo.displacementMag_given) \
463 clitk::Lists<itk::FixedArray<ValueType,1> > displacementListsMag=displacementLists.Norm(); \
464 std::vector<std::string> filenames; \
465 for (unsigned int i=0;i<m_NumberOfFields;i++) \
467 std::ostringstream number_ostr; \
469 std::string number_str = number_ostr.str(); \
470 filenames.push_back(m_ArgsInfo.displacementMag_arg+number_str); \
472 displacementListsMag.Write(filenames, m_Verbose ); \
475 /* // Output warped points */\
476 if( m_ArgsInfo.warp_given) \
479 std::vector<std::string> filenames; \
480 for (unsigned int i=0;i<m_ArgsInfo.warp_given;i++) \
482 std::ostringstream number_ostr; \
484 std::string number_str = number_ostr.str(); \
485 filenames.push_back(m_ArgsInfo.warp_arg+number_str); \
487 warpedPointLists.Write(filenames, m_Verbose ); \
491 if( m_ArgsInfo.tre_given) \
494 std::vector<std::string> filenames; \
495 for (unsigned int i=0;i<m_NumberOfFields;i++) \
497 std::ostringstream number_ostr; \
499 std::string number_str = number_ostr.str(); \
500 filenames.push_back(m_ArgsInfo.tre_arg+number_str); \
502 treLists.Write(filenames, m_Verbose ); \
505 /* // Output tre mag */\
506 if( m_ArgsInfo.treMag_given) \
508 clitk::Lists<itk::FixedArray<ValueType,1> > treMagLists=treLists.Norm(); \
510 std::vector<std::string> filenames; \
511 for (unsigned int i=0;i<m_NumberOfFields;i++) \
513 std::ostringstream number_ostr; \
515 std::string number_str = number_ostr.str(); \
516 filenames.push_back(m_ArgsInfo.treMag_arg+number_str); \
518 treMagLists.Write(filenames, m_Verbose ); \
527 //-----------------------------
529 //-----------------------------
530 template<unsigned int Dimension, unsigned int Components >
532 CalculateTREGenericFilter::ReadVectorFields(void)
535 typedef itk::Vector<ValueType, Components> VectorType;
536 typedef itk::Image<VectorType, Dimension> DeformationFieldType;
538 typedef itk::ImageFileReader<DeformationFieldType> DvfReaderType;
539 std::vector<typename DeformationFieldType::Pointer> dvfs;
541 for (unsigned int i=0; i< m_ArgsInfo.vf_given; i++)
543 typename DvfReaderType::Pointer reader = DvfReaderType::New();
544 reader->SetFileName( m_ArgsInfo.vf_arg[i]);
545 if (m_Verbose) std::cout<<"Reading vector field "<< i <<"..."<<std::endl;
547 dvfs.push_back( reader->GetOutput() );
550 ProcessVectorFields<Dimension, Components>(dvfs, m_ArgsInfo.vf_arg);
554 //-----------------------------
556 //-----------------------------
557 template<unsigned int Dimension, unsigned int Components >
559 CalculateTREGenericFilter::ProcessVectorFields(std::vector<typename itk::Image<itk::Vector<CalculateTREGenericFilter::ValueType, Components>, Dimension>::Pointer > dvfs, char** filenames )
561 std::vector<std::string> new_filenames;
562 for (unsigned int i=0;i<m_ArgsInfo.vf_given;i++)
563 new_filenames.push_back(filenames[i]);
564 UpdateDVFWithDim<Dimension>(dvfs, new_filenames);
568 //-----------------------------
569 // Read Coefficient Images
570 //-----------------------------
571 template<unsigned int Dimension, unsigned int Components >
573 CalculateTREGenericFilter::ReadCoefficientImages(void)
575 // if using coefficient images
576 typedef itk::Vector<ValueType, Components> VectorType;
577 typedef itk::Image<VectorType, Dimension> CoefficientImageType;
579 typedef itk::ImageFileReader<CoefficientImageType> CoeffReaderType;
580 std::vector<typename CoefficientImageType::Pointer> coeffs;
582 for (unsigned int i=0; i< m_ArgsInfo.coeff_given; i++)
584 typename CoeffReaderType::Pointer reader = CoeffReaderType::New();
585 reader->SetFileName( m_ArgsInfo.coeff_arg[i]);
586 if (m_Verbose) std::cout<<"Reading coefficient image "<< i <<"..."<<std::endl;
588 coeffs.push_back( reader->GetOutput() );
591 ProcessCoefficientImages<Dimension, Components>(coeffs, m_ArgsInfo.coeff_arg);
594 //-----------------------------
595 // Process Coefficient Images
596 //-----------------------------
597 template<unsigned int Dimension, unsigned int Components >
599 CalculateTREGenericFilter::ProcessCoefficientImages(std::vector<typename itk::Image<itk::Vector<CalculateTREGenericFilter::ValueType, Components>, Dimension>::Pointer > coeffs, char** filenames )
601 std::vector<std::string> new_filenames;
602 for (unsigned int i=0;i<m_ArgsInfo.coeff_given;i++)
603 new_filenames.push_back(filenames[i]);
604 UpdateCoeffsWithDim<Dimension>(coeffs, new_filenames);
608 //-------------------------------------------------------------------
609 // Update DVF with the number of dimensions
610 //-------------------------------------------------------------------
611 template<unsigned int Dimension>
613 CalculateTREGenericFilter::UpdateDVFWithDim(std::vector<typename itk::Image<itk::Vector<CalculateTREGenericFilter::ValueType, Dimension>, Dimension>::Pointer > dvfs, std::vector<std::string> filenames)
615 if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D..."<<std::endl;
618 BuildPointListsMacro(filenames);
619 OriginalDistancesMacro;
621 //=====================================================================================
622 // Distance between points after warp
623 //=====================================================================================
625 //-----------------------------
627 //-----------------------------
628 typedef clitk::GenericVectorInterpolator<args_info_clitkCalculateTRE,DeformationFieldType, ValueType> GenericVectorInterpolatorType;
629 typename GenericVectorInterpolatorType::Pointer genericInterpolator=GenericVectorInterpolatorType::New();
630 genericInterpolator->SetArgsInfo(m_ArgsInfo);
631 typedef itk::VectorInterpolateImageFunction<DeformationFieldType, ValueType> InterpolatorType;
632 typename InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer();
635 PointListsType warpedPointLists(m_NumberOfFields);
636 VectorListsType treLists(m_NumberOfFields), displacementLists(m_NumberOfFields);
638 //-----------------------------
639 // Calculate residual distances
640 //-----------------------------
641 VectorType displacement, tre;
642 PointType warpedPoint;
644 for (unsigned int phaseIndex=0; phaseIndex<m_NumberOfFields; phaseIndex++)
646 interpolator->SetInputImage( dvfs[phaseIndex] );
647 for (unsigned int pointIndex=0; pointIndex<m_NumberOfPoints; pointIndex++)
650 referencePoint=referencePointList[pointIndex];
653 if ( interpolator->IsInsideBuffer(referencePoint) )
654 displacement=interpolator->Evaluate(referencePoint);
656 displacement.Fill(0.0);
659 warpedPoint=referencePoint+displacement;
660 displacementLists[phaseIndex].push_back(displacement);
661 warpedPointLists[phaseIndex].push_back(warpedPoint);
662 tre=pointLists[phaseIndex][pointIndex]-warpedPoint;
663 treLists[phaseIndex].push_back(tre);
667 StatisticsDisplacementsMacro;
670 //-------------------------------------------------------------------
671 // Update coefficients with the number of dimensions
672 //-------------------------------------------------------------------
673 template<unsigned int Dimension>
675 CalculateTREGenericFilter::UpdateCoeffsWithDim(std::vector<typename itk::Image<itk::Vector<CalculateTREGenericFilter::ValueType, Dimension>, Dimension>::Pointer > coeffs, std::vector<std::string> filenames)
679 BuildPointListsMacro(filenames);
680 OriginalDistancesMacro;
682 //=====================================================================================
683 // Distance between points after warp
684 //=====================================================================================
686 //-----------------------------
687 // Transform to get the transformed points from the coeff images
688 //-----------------------------
689 typedef clitk::BSplineDeformableTransform<ValueType, Dimension, Dimension> TransformType;
690 typename TransformType::Pointer transform = TransformType::New();
692 PointListsType warpedPointLists(m_NumberOfFields);
693 VectorListsType treLists(m_NumberOfFields), displacementLists(m_NumberOfFields);
695 //-----------------------------
696 // Calculate residual distances
697 //-----------------------------
698 VectorType displacement, tre;
699 PointType warpedPoint;
701 for (unsigned int phaseIndex=0; phaseIndex<m_NumberOfFields; phaseIndex++)
703 // set the initial parameters (b-spline coefficients from cmd-line
704 transform->SetCoefficientImage(coeffs[phaseIndex]);
706 // transform the input points
707 for (unsigned int pointIndex=0; pointIndex<m_NumberOfPoints; pointIndex++)
710 referencePoint=referencePointList[pointIndex];
713 warpedPoint=transform->TransformPoint(referencePoint);
716 displacement=warpedPoint-referencePoint;
717 displacementLists[phaseIndex].push_back(displacement);
718 warpedPointLists[phaseIndex].push_back(warpedPoint);
719 tre=pointLists[phaseIndex][pointIndex]-warpedPoint;
720 treLists[phaseIndex].push_back(tre);
724 StatisticsDisplacementsMacro;
728 #endif //#define clitkCalculateTREGenericFilter_txx