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 clitkCalculateTREGenericFilter_txx
19 #define clitkCalculateTREGenericFilter_txx
21 /* =================================================
22 * @file clitkCalculateTREGenericFilter.txx
28 ===================================================*/
34 //-----------------------------
36 //-----------------------------
37 template<unsigned int Dimension, unsigned int Components >
39 CalculateTREGenericFilter::ReadVectorFields(void)
42 typedef itk::Vector<float, Components> VectorType;
43 typedef itk::Image<VectorType, Dimension> DeformationFieldType;
45 typedef itk::ImageFileReader<DeformationFieldType> InputReaderType;
46 std::vector<typename DeformationFieldType::Pointer> dvfs;
47 for (unsigned int i=0; i< m_ArgsInfo.vf_given; i++)
49 typename InputReaderType::Pointer reader = InputReaderType::New();
50 reader->SetFileName( m_ArgsInfo.vf_arg[i]);
51 if (m_Verbose) std::cout<<"Reading vector field "<< i <<"..."<<std::endl;
53 dvfs.push_back( reader->GetOutput() );
56 ProcessVectorFields<Dimension, Components>(dvfs, m_ArgsInfo.vf_arg);
59 //-----------------------------
61 //-----------------------------
62 template<unsigned int Dimension, unsigned int Components >
64 CalculateTREGenericFilter::ProcessVectorFields(std::vector<typename itk::Image<itk::Vector<float, Components>, Dimension>::Pointer > dvfs, char** filenames )
67 std::vector<std::string> new_filenames;
68 for (unsigned int i=0;i<m_ArgsInfo.vf_given;i++)
69 new_filenames.push_back(filenames[i]);
70 UpdateWithDim<Dimension>(dvfs, new_filenames);
74 //-------------------------------------------------------------------
75 // Update with the number of dimensions
76 //-------------------------------------------------------------------
77 template<unsigned int Dimension>
79 CalculateTREGenericFilter::UpdateWithDim(std::vector<typename itk::Image<itk::Vector<float, Dimension>, Dimension>::Pointer > dvfs, std::vector<std::string> filenames)
81 if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D..."<<std::endl;
83 //-----------------------------
85 //-----------------------------
86 typedef double ValueType;
87 typedef std::vector<ValueType> MeasureListType;
89 typedef itk::Point<double, Dimension> PointType;
90 typedef clitk::List<PointType> PointListType;
91 typedef clitk::Lists<PointType> PointListsType;
93 typedef itk::Vector<double, Dimension> VectorType;
94 typedef clitk::List<VectorType> VectorListType;
95 typedef clitk::Lists<VectorType> VectorListsType;
97 typedef itk::Vector<float, Dimension> DeformationVectorType;
98 typedef itk::Image<DeformationVectorType, Dimension> DeformationFieldType;
100 //-----------------------------
102 //-----------------------------
103 unsigned int numberOfFields=filenames.size();
104 unsigned int numberOfLists=m_ArgsInfo.input_given;
105 if ( (numberOfLists!=numberOfFields) && (numberOfLists!=1) )
107 std::cerr<<"Error: Number of lists (="<<numberOfLists<<") different from number of DVF's (="<<numberOfFields<<") !"<<std::endl;
110 else if (numberOfLists==1)
112 if (m_Verbose) std::cout<<numberOfFields<<" DVFs and 1 point list given..."<<std::endl;
116 if (m_Verbose) std::cout<<numberOfLists<<" point lists and DVFs given..."<<std::endl;
120 //-----------------------------
122 //-----------------------------
123 PointListsType pointLists;
124 unsigned int numberOfPoints=0;
125 for (unsigned int i=0; i<numberOfFields; i++)
128 if (numberOfLists==1)
129 pointLists.push_back(PointListType(m_ArgsInfo.input_arg[0], m_Verbose) );
131 pointLists.push_back(PointListType(m_ArgsInfo.input_arg[i], m_Verbose) );
134 // Verify the number of points
135 if (i==0) numberOfPoints=pointLists[i].size();
138 if (numberOfPoints!=pointLists[i].size())
140 std::cerr<<"Size of first list ("<<numberOfPoints
141 <<") is different from size of list "<<i
142 <<" ("<<pointLists[i].size()<<")..."<<std::endl;
149 PointListType referencePointList;
150 if (m_Verbose) std::cout<<"Reference point list:"<<std::endl;
151 referencePointList=PointListType(m_ArgsInfo.ref_arg, m_Verbose);
152 if (numberOfPoints!=referencePointList.size())
154 std::cerr<<"Size of the first list ("<<numberOfPoints
155 <<") is different from size of the reference list ("
156 << referencePointList.size() <<")..."<<std::endl;
161 //-----------------------------
163 //-----------------------------
164 typedef clitk::GenericVectorInterpolator<args_info_clitkCalculateTRE,DeformationFieldType, double> GenericVectorInterpolatorType;
165 typename GenericVectorInterpolatorType::Pointer genericInterpolator=GenericVectorInterpolatorType::New();
166 genericInterpolator->SetArgsInfo(m_ArgsInfo);
167 typedef itk::VectorInterpolateImageFunction<DeformationFieldType, double> InterpolatorType;
168 typename InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer();
171 //=====================================================================================
172 // Original distance between points
173 //=====================================================================================
174 VectorListsType originalDistanceLists(numberOfFields);
176 //-----------------------------
177 // Calculate original distances
178 //-----------------------------
179 PointType referencePoint;
181 for (unsigned int phaseIndex=0; phaseIndex<numberOfFields; phaseIndex++)
183 for (unsigned int pointIndex=0; pointIndex<numberOfPoints; pointIndex++)
185 referencePoint=referencePointList[pointIndex];
186 distance=pointLists[phaseIndex][pointIndex]-referencePoint;
187 originalDistanceLists[phaseIndex].push_back(distance);
192 //-----------------------------
193 // Statistics Original
194 //-----------------------------
195 typedef DeformationListStatisticsFilter<VectorType> StatisticsFilterType;
196 typename StatisticsFilterType::Pointer statisticsFilter=StatisticsFilterType::New();
198 // Statistics (magnitude)
199 MeasureListType oMeanList, oStdList, oMaxList;
200 ValueType oMean, oStd, oMax;
201 statisticsFilter->GetStatistics(originalDistanceLists, oMean, oStd, oMax, oMeanList, oStdList, oMaxList);
203 // Statistics (per component)
204 VectorListType oMeanXYZList, oStdXYZList,oMaxXYZList;
205 VectorType oMeanXYZ, oStdXYZ, oMaxXYZ;
206 statisticsFilter->GetStatistics(originalDistanceLists, oMeanXYZ, oStdXYZ, oMaxXYZ, oMeanXYZList, oStdXYZList, oMaxXYZList);
209 //-----------------------------
211 //-----------------------------
212 std::vector<std::string> labels;
213 labels.push_back("MeanX\tSDX");
214 labels.push_back("MeanY\tSDY");
215 labels.push_back("MeanZ\tSDZ");
216 labels.push_back("MeanT\tSDT");
217 labels.push_back("MaxX");
218 labels.push_back("MaxY");
219 labels.push_back("MaxZ");
220 labels.push_back("MaxT");
228 std::cout<<"# Number\tDVF"<<std::endl;
229 for (unsigned int phaseIndex=0; phaseIndex<numberOfFields; phaseIndex++)
230 std::cout<<phaseIndex<<"\t"<<filenames[phaseIndex]<<std::endl;
233 std::cout<<std::endl;
234 std::cout<<"=================================================="<<std::endl;
235 std::cout<<"|| Original distance between points ||"<<std::endl;
236 std::cout<<"=================================================="<<std::endl;
237 std::cout<<std::endl;
239 std::cout<<std::setprecision(3);
243 // Labels of the columns
244 std::cout<<"#DVF\tMean\tSD\t"<<labels[0];
245 for(unsigned int dim=1;dim<Dimension;dim++) std::cout<<"\t"<< labels[dim];
246 std::cout<<"\tMax \t"<<labels[4];
247 for(unsigned int dim=1;dim<Dimension;dim++) std::cout<<"\t"<< labels[dim+4];
248 std::cout<<std::endl;
249 for (unsigned int phaseIndex=0; phaseIndex<numberOfFields; phaseIndex++)
251 std::cout<<phaseIndex;
252 std::cout<<"\t"<<oMeanList[phaseIndex]<<"\t"<<oStdList[phaseIndex];
253 std::cout<<"\t"<<oMeanXYZList[phaseIndex][0]<<"\t"<<oStdXYZList[phaseIndex][0];
254 for(unsigned int dim=1;dim<Dimension;dim++)
255 std::cout<<"\t"<<oMeanXYZList[phaseIndex][dim]<<"\t"<<oStdXYZList[phaseIndex][dim];
256 std::cout<<"\t"<<oMaxList[phaseIndex];
257 std::cout<<"\t"<<oMaxXYZList[phaseIndex][0];
258 for(unsigned int dim=1;dim<Dimension;dim++)
259 std::cout<<"\t"<<oMaxXYZList[phaseIndex][dim];
260 std::cout<<std::endl;
264 std::cout<<std::endl;
266 std::cout<<"\t"<<oMean<<"\t"<<oStd;
267 std::cout<<"\t"<<oMeanXYZ[0]<<"\t"<<oStdXYZ[0];
268 for(unsigned int dim=1;dim<Dimension;dim++)
269 std::cout<<"\t"<<oMeanXYZ[dim]<<"\t"<<oStdXYZ[dim];
270 std::cout<<"\t"<<oMax;
271 std::cout<<"\t"<<oMaxXYZ[0];
272 for(unsigned int dim=1;dim<Dimension;dim++)
273 std::cout<<"\t"<<oMaxXYZ[dim];
274 std::cout<<std::endl;
278 if( m_ArgsInfo.general_given)
281 std::ofstream general(m_ArgsInfo.general_arg);
285 for (unsigned int phaseIndex=0; phaseIndex<numberOfFields; phaseIndex++)
286 general<<phaseIndex<<"\t"<< filenames[phaseIndex]<<std::endl;
289 general<<"=================================================="<<std::endl;
290 general<<"|| Original distance between points ||"<<std::endl;
291 general<<"=================================================="<<std::endl;
294 general<<std::setprecision(3);
296 // Labels of the columns
297 general<<"#DVF\tMean\tSD\t"<<labels[0];
298 for(unsigned int dim=1;dim<Dimension;dim++) general<<"\t"<< labels[dim];
299 general<<"\tMax \t"<<labels[4];
300 for(unsigned int dim=1;dim<Dimension;dim++) general<<"\t"<< labels[dim+4];
302 for (unsigned int phaseIndex=0; phaseIndex<numberOfFields; phaseIndex++)
305 general<<"\t"<<oMeanList[phaseIndex]<<"\t"<<oStdList[phaseIndex];
306 general<<"\t"<<oMeanXYZList[phaseIndex][0]<<"\t"<<oStdXYZList[phaseIndex][0];
307 for(unsigned int dim=1;dim<Dimension;dim++)
308 general<<"\t"<<oMeanXYZList[phaseIndex][dim]<<"\t"<<oStdXYZList[phaseIndex][dim];
309 general<<"\t"<<oMaxList[phaseIndex];
310 general<<"\t"<<oMaxXYZList[phaseIndex][0];
311 for(unsigned int dim=1;dim<Dimension;dim++)
312 general<<"\t"<<oMaxXYZList[phaseIndex][dim];
319 general<<"\t"<<oMean<<"\t"<<oStd;
320 general<<"\t"<<oMeanXYZ[0]<<"\t"<<oStdXYZ[0];
321 for(unsigned int dim=1;dim<Dimension;dim++)
322 general<<"\t"<<oMeanXYZ[dim]<<"\t"<<oStdXYZ[dim];
324 general<<"\t"<<oMaxXYZ[0];
325 for(unsigned int dim=1;dim<Dimension;dim++)
326 general<<"\t"<<oMaxXYZ[dim];
332 //=====================================================================================
333 // Distance between points after warp
334 //=====================================================================================
335 PointListsType warpedPointLists(numberOfFields);
336 VectorListsType treLists(numberOfFields), displacementLists(numberOfFields);
338 //-----------------------------
339 // Calculate residual distances
340 //-----------------------------
341 VectorType displacement, tre;
342 PointType warpedPoint;
344 for (unsigned int phaseIndex=0; phaseIndex<numberOfFields; phaseIndex++)
346 interpolator->SetInputImage( dvfs[phaseIndex]);
347 for (unsigned int pointIndex=0; pointIndex<numberOfPoints; pointIndex++)
350 referencePoint=referencePointList[pointIndex];
353 if ( interpolator->IsInsideBuffer(referencePoint) )
354 displacement=interpolator->Evaluate(referencePoint);
356 displacement.Fill(0.0);
359 warpedPoint=referencePoint+displacement;
360 displacementLists[phaseIndex].push_back(displacement);
361 warpedPointLists[phaseIndex].push_back(warpedPoint);
362 tre=pointLists[phaseIndex][pointIndex]-warpedPoint;
363 treLists[phaseIndex].push_back(tre);
368 //-----------------------------
369 // Statistics displacements
370 //-----------------------------
372 // Statistics (magnitude)
373 MeasureListType dmeanList, dstdList, dmaxList;
374 ValueType dmean, dstd, dmax;
375 statisticsFilter->GetStatistics(displacementLists, dmean, dstd, dmax, dmeanList, dstdList, dmaxList);
377 // Statistics (per component)
378 VectorListType dmeanXYZList, dstdXYZList,dmaxXYZList;
379 VectorType dmeanXYZ, dstdXYZ, dmaxXYZ;
380 statisticsFilter->GetStatistics(displacementLists, dmeanXYZ, dstdXYZ, dmaxXYZ,dmeanXYZList, dstdXYZList, dmaxXYZList);
383 //-----------------------------
385 //-----------------------------
387 // Statistics (magnitude)
388 MeasureListType meanList, stdList, maxList;
389 ValueType mean, std, max;
390 statisticsFilter->GetStatistics(treLists, mean, std, max, meanList, stdList, maxList);
392 // Statistics (per component)
393 VectorListType meanXYZList, stdXYZList,maxXYZList;
394 VectorType meanXYZ, stdXYZ, maxXYZ;
395 statisticsFilter->GetStatistics(treLists, meanXYZ, stdXYZ, maxXYZ, meanXYZList, stdXYZList, maxXYZList);
402 std::cout<<std::endl;
403 std::cout<<"=================================================="<<std::endl;
404 std::cout<<"|| Residual distance between points ||"<<std::endl;
405 std::cout<<"=================================================="<<std::endl;
406 std::cout<<std::endl;
408 std::cout<<std::setprecision(3);
410 // Labels of the columns
411 std::cout<<"#DVF\tMean\tSD\t"<<labels[0];
412 for(unsigned int dim=1;dim<Dimension;dim++) std::cout<<"\t"<< labels[dim];
413 std::cout<<"\tMax \t"<<labels[4];
414 for(unsigned int dim=1;dim<Dimension;dim++) std::cout<<"\t"<< labels[dim+4];
415 std::cout<<std::endl;
416 for (unsigned int phaseIndex=0; phaseIndex<numberOfFields; phaseIndex++)
418 std::cout<<phaseIndex;
419 std::cout<<"\t"<<meanList[phaseIndex]<<"\t"<<stdList[phaseIndex];
420 std::cout<<"\t"<<meanXYZList[phaseIndex][0]<<"\t"<<stdXYZList[phaseIndex][0];
421 for(unsigned int dim=1;dim<Dimension;dim++)
422 std::cout<<"\t"<<meanXYZList[phaseIndex][dim]<<"\t"<<stdXYZList[phaseIndex][dim];
423 std::cout<<"\t"<<maxList[phaseIndex];
424 std::cout<<"\t"<<maxXYZList[phaseIndex][0];
425 for(unsigned int dim=1;dim<Dimension;dim++)
426 std::cout<<"\t"<<maxXYZList[phaseIndex][dim];
427 std::cout<<std::endl;
431 std::cout<<std::endl;
433 std::cout<<"\t"<<mean<<"\t"<<std;
434 std::cout<<"\t"<<meanXYZ[0]<<"\t"<<stdXYZ[0];
435 for(unsigned int dim=1;dim<Dimension;dim++)
436 std::cout<<"\t"<<meanXYZ[dim]<<"\t"<<stdXYZ[dim];
437 std::cout<<"\t"<<max;
438 std::cout<<"\t"<<maxXYZ[0];
439 for(unsigned int dim=1;dim<Dimension;dim++)
440 std::cout<<"\t"<<maxXYZ[dim];
441 std::cout<<std::endl;
445 if( m_ArgsInfo.general_given)
448 std::ofstream general(m_ArgsInfo.general_arg, ios_base::app);
451 general<<"=================================================="<<std::endl;
452 general<<"|| Residual distance between points ||"<<std::endl;
453 general<<"=================================================="<<std::endl;
456 general<<std::setprecision(3);
458 // Labels of the columns
459 general<<"#DVF\tMean\tSD\t"<<labels[0];
460 for(unsigned int dim=1;dim<Dimension;dim++) general<<"\t"<< labels[dim];
461 general<<"\tMax \t"<<labels[4];
462 for(unsigned int dim=1;dim<Dimension;dim++) general<<"\t"<< labels[dim+4];
464 for (unsigned int phaseIndex=0; phaseIndex<numberOfFields; phaseIndex++)
467 general<<"\t"<<meanList[phaseIndex]<<"\t"<<stdList[phaseIndex];
468 general<<"\t"<<meanXYZList[phaseIndex][0]<<"\t"<<stdXYZList[phaseIndex][0];
469 for(unsigned int dim=1;dim<Dimension;dim++)
470 general<<"\t"<<meanXYZList[phaseIndex][dim]<<"\t"<<stdXYZList[phaseIndex][dim];
471 general<<"\t"<<maxList[phaseIndex];
472 general<<"\t"<<maxXYZList[phaseIndex][0];
473 for(unsigned int dim=1;dim<Dimension;dim++)
474 general<<"\t"<<maxXYZList[phaseIndex][dim];
481 general<<"\t"<<mean<<"\t"<<std;
482 general<<"\t"<<meanXYZ[0]<<"\t"<<stdXYZ[0];
483 for(unsigned int dim=1;dim<Dimension;dim++)
484 general<<"\t"<<meanXYZ[dim]<<"\t"<<stdXYZ[dim];
486 general<<"\t"<<maxXYZ[0];
487 for(unsigned int dim=1;dim<Dimension;dim++)
488 general<<"\t"<<maxXYZ[dim];
493 // Output original points
494 if( m_ArgsInfo.original_given)
497 std::vector<std::string> filenames;
498 for (unsigned int i=0;i<numberOfFields;i++)
500 std::ostringstream number_ostr;
502 std::string number_str = number_ostr.str();
503 filenames.push_back(m_ArgsInfo.original_arg+number_str);
505 originalDistanceLists.Write(filenames, m_Verbose );
509 // Output original magnitude points
510 if( m_ArgsInfo.originalMag_given)
512 clitk::Lists<itk::FixedArray<double,1> > originalDistanceListsMag=originalDistanceLists.Norm();
513 std::vector<std::string> filenames;
514 for (unsigned int i=0;i<numberOfFields;i++)
516 std::ostringstream number_ostr;
518 std::string number_str = number_ostr.str();
519 filenames.push_back(m_ArgsInfo.originalMag_arg+number_str);
521 originalDistanceListsMag.Write(filenames, m_Verbose );
524 // Output displacement
525 if( m_ArgsInfo.displacement_given)
528 std::vector<std::string> filenames;
529 for (unsigned int i=0;i<numberOfFields;i++)
531 std::ostringstream number_ostr;
533 std::string number_str = number_ostr.str();
534 filenames.push_back(m_ArgsInfo.displacement_arg+number_str);
536 displacementLists.Write(filenames, m_Verbose );
540 // Output displacement magnitude
541 if( m_ArgsInfo.displacementMag_given)
543 clitk::Lists<itk::FixedArray<double,1> > displacementListsMag=displacementLists.Norm();
544 std::vector<std::string> filenames;
545 for (unsigned int i=0;i<numberOfFields;i++)
547 std::ostringstream number_ostr;
549 std::string number_str = number_ostr.str();
550 filenames.push_back(m_ArgsInfo.displacementMag_arg+number_str);
552 displacementListsMag.Write(filenames, m_Verbose );
555 // Output warped points
556 if( m_ArgsInfo.warp_given)
559 std::vector<std::string> filenames;
560 for (unsigned int i=0;i<m_ArgsInfo.warp_given;i++)
562 std::ostringstream number_ostr;
564 std::string number_str = number_ostr.str();
565 filenames.push_back(m_ArgsInfo.warp_arg+number_str);
567 warpedPointLists.Write(filenames, m_Verbose );
571 if( m_ArgsInfo.tre_given)
574 std::vector<std::string> filenames;
575 for (unsigned int i=0;i<numberOfFields;i++)
577 std::ostringstream number_ostr;
579 std::string number_str = number_ostr.str();
580 filenames.push_back(m_ArgsInfo.tre_arg+number_str);
582 treLists.Write(filenames, m_Verbose );
586 if( m_ArgsInfo.treMag_given)
588 clitk::Lists<itk::FixedArray<double,1> > treMagLists=treLists.Norm();
590 std::vector<std::string> filenames;
591 for (unsigned int i=0;i<numberOfFields;i++)
593 std::ostringstream number_ostr;
595 std::string number_str = number_ostr.str();
596 filenames.push_back(m_ArgsInfo.treMag_arg+number_str);
598 treMagLists.Write(filenames, m_Verbose );
606 #endif //#define clitkCalculateTREGenericFilter_txx