+//=====================================================================================
+// Macros
+// Rômulo Pinho - 25/02/2011
+// The macros below avoid the duplication of code after introducing the function
+// UpdateCoeffsWithDim(). Although the code generated by these functions will be
+// duplicated anyway (since they're template functions) the macros make their sources
+// more intelligible.
+// Ideally, these macros should be made member functions, with the typedefs inside the
+// the filter class, but that would change the class' interface, since it would need
+// to become a template class (with parameter Dimension). To avoid breaking legacy code,
+// I opted for the macros, although I don't quite like them (especially long ones like
+// these). In any case, it's advisable to think of a better solution, even if it breaks
+// the class' interface.
+//=====================================================================================
+
+//-----------------------------
+// Typedefs
+//-----------------------------
+#define TypedefsMacro \
+\
+ typedef itk::Point<ValueType, Dimension> PointType; \
+ typedef clitk::List<PointType> PointListType; \
+ typedef clitk::Lists<PointType> PointListsType; \
+\
+ typedef itk::Vector<ValueType, Dimension> VectorType; \
+ typedef clitk::List<VectorType> VectorListType; \
+ typedef clitk::Lists<VectorType> VectorListsType; \
+\
+ typedef itk::Vector<ValueType, Dimension> DeformationVectorType; \
+ typedef itk::Image<DeformationVectorType, Dimension> DeformationFieldType;
+
+//-------------------------------------------------------------------
+// Build the point lists necessary for the computation of distances.
+//-------------------------------------------------------------------
+#define BuildPointListsMacro(filenames) \
+\
+ m_NumberOfFields=filenames.size(); \
+ m_NumberOfLists=m_ArgsInfo.input_given; \
+ if ( (m_NumberOfLists!=m_NumberOfFields) && (m_NumberOfLists!=1) ) \
+ { \
+ std::cerr<<"Error: Number of lists (="<<m_NumberOfLists<<") different from number of DVF's (="<<m_NumberOfFields<<") !"<<std::endl; \
+ return; \
+ } \
+ else if (m_NumberOfLists==1) \
+ { \
+ if (m_Verbose) std::cout<<m_NumberOfFields<<" DVFs and 1 point list given..."<<std::endl; \
+ } \
+ else \
+ { \
+ if (m_Verbose) std::cout<<m_NumberOfLists<<" point lists and DVFs given..."<<std::endl; \
+ } \
+\
+\
+ PointListsType pointLists; \
+ m_NumberOfPoints=0; \
+ for (unsigned int i=0; i<m_NumberOfFields; i++) \
+ { \
+ /* Read the lists */\
+ if (m_NumberOfLists==1) \
+ pointLists.push_back(PointListType(m_ArgsInfo.input_arg[0], m_Verbose) ); \
+ else \
+ pointLists.push_back(PointListType(m_ArgsInfo.input_arg[i], m_Verbose) ); \
+\
+\
+ /* Verify the number of points */\
+ if (i==0) m_NumberOfPoints=pointLists[i].size(); \
+ else \
+ { \
+ if (m_NumberOfPoints!=pointLists[i].size()) \
+ { \
+ std::cerr<<"Size of first list ("<<m_NumberOfPoints \
+ <<") is different from size of list "<<i \
+ <<" ("<<pointLists[i].size()<<")..."<<std::endl; \
+ return; \
+ } \
+ } \
+ } \
+\
+\
+ PointListType referencePointList; \
+ if (m_Verbose) std::cout<<"Reference point list:"<<std::endl; \
+ referencePointList=PointListType(m_ArgsInfo.ref_arg, m_Verbose); \
+ if (m_NumberOfPoints!=referencePointList.size()) \
+ { \
+ std::cerr<<"Size of the first list ("<<m_NumberOfPoints \
+ <<") is different from size of the reference list (" \
+ << referencePointList.size() <<")..."<<std::endl; \
+ return; \
+ }
+
+/* //===================================================================================== */\
+/* // Original distance between points */\
+/* //===================================================================================== */\
+#define OriginalDistancesMacro \
+ VectorListsType originalDistanceLists(m_NumberOfFields); \
+\
+ /* //----------------------------- */\
+ /* // Calculate original distances */\
+ /* //----------------------------- */\
+ PointType referencePoint; \
+ VectorType distance; \
+ for (unsigned int phaseIndex=0; phaseIndex<m_NumberOfFields; phaseIndex++) \
+ { \
+ for (unsigned int pointIndex=0; pointIndex<m_NumberOfPoints; pointIndex++) \
+ { \
+ referencePoint=referencePointList[pointIndex]; \
+ distance=pointLists[phaseIndex][pointIndex]-referencePoint; \
+ originalDistanceLists[phaseIndex].push_back(distance); \
+ } \
+ } \
+\
+\
+ /* //----------------------------- */\
+ /* // Statistics Original */\
+ /* //----------------------------- */\
+ typedef DeformationListStatisticsFilter<VectorType> StatisticsFilterType; \
+ typename StatisticsFilterType::Pointer statisticsFilter=StatisticsFilterType::New(); \
+\
+ /* // Statistics (magnitude) */ \
+ MeasureListType oMeanList, oStdList, oMaxList; \
+ ValueType oMean, oStd, oMax; \
+ statisticsFilter->GetStatistics(originalDistanceLists, oMean, oStd, oMax, oMeanList, oStdList, oMaxList); \
+\
+ /* // Statistics (per component) */\
+ VectorListType oMeanXYZList, oStdXYZList,oMaxXYZList; \
+ VectorType oMeanXYZ, oStdXYZ, oMaxXYZ; \
+ statisticsFilter->GetStatistics(originalDistanceLists, oMeanXYZ, oStdXYZ, oMaxXYZ, oMeanXYZList, oStdXYZList, oMaxXYZList); \
+\
+\
+ /* //----------------------------- */\
+ /* // Output */\
+ /* //----------------------------- */\
+ std::vector<std::string> labels; \
+ labels.push_back("MeanX\tSDX"); \
+ labels.push_back("MeanY\tSDY"); \
+ labels.push_back("MeanZ\tSDZ"); \
+ labels.push_back("MeanT\tSDT"); \
+ labels.push_back("MaxX"); \
+ labels.push_back("MaxY"); \
+ labels.push_back("MaxZ"); \
+ labels.push_back("MaxT"); \
+\
+\
+ /* // Output to screen */\
+ if(m_Verbose) \
+ { \
+\
+ /* // Numbers of DVF */\
+ std::cout<<"# Number\tDVF"<<std::endl; \
+ for (unsigned int phaseIndex=0; phaseIndex<m_NumberOfFields; phaseIndex++) \
+ std::cout<<phaseIndex<<"\t"<<filenames[phaseIndex]<<std::endl; \
+\
+\
+ std::cout<<std::endl; \
+ std::cout<<"=================================================="<<std::endl; \
+ std::cout<<"|| Original distance between points ||"<<std::endl; \
+ std::cout<<"=================================================="<<std::endl; \
+ std::cout<<std::endl; \
+\
+ std::cout<<std::setprecision(3); \
+\
+\
+\
+ /* // Labels of the columns */\
+ std::cout<<"#DVF\tMean\tSD\t"<<labels[0]; \
+ for(unsigned int dim=1;dim<Dimension;dim++) std::cout<<"\t"<< labels[dim]; \
+ std::cout<<"\tMax \t"<<labels[4]; \
+ for(unsigned int dim=1;dim<Dimension;dim++) std::cout<<"\t"<< labels[dim+4]; \
+ std::cout<<std::endl; \
+ for (unsigned int phaseIndex=0; phaseIndex<m_NumberOfFields; phaseIndex++) \
+ { \
+ std::cout<<phaseIndex; \
+ std::cout<<"\t"<<oMeanList[phaseIndex]<<"\t"<<oStdList[phaseIndex]; \
+ std::cout<<"\t"<<oMeanXYZList[phaseIndex][0]<<"\t"<<oStdXYZList[phaseIndex][0]; \
+ for(unsigned int dim=1;dim<Dimension;dim++) \
+ std::cout<<"\t"<<oMeanXYZList[phaseIndex][dim]<<"\t"<<oStdXYZList[phaseIndex][dim]; \
+ std::cout<<"\t"<<oMaxList[phaseIndex]; \
+ std::cout<<"\t"<<oMaxXYZList[phaseIndex][0]; \
+ for(unsigned int dim=1;dim<Dimension;dim++) \
+ std::cout<<"\t"<<oMaxXYZList[phaseIndex][dim]; \
+ std::cout<<std::endl; \
+ } \
+\
+ /* // General */\
+ std::cout<<std::endl; \
+ std::cout<<"@"; \
+ std::cout<<"\t"<<oMean<<"\t"<<oStd; \
+ std::cout<<"\t"<<oMeanXYZ[0]<<"\t"<<oStdXYZ[0]; \
+ for(unsigned int dim=1;dim<Dimension;dim++) \
+ std::cout<<"\t"<<oMeanXYZ[dim]<<"\t"<<oStdXYZ[dim]; \
+ std::cout<<"\t"<<oMax; \
+ std::cout<<"\t"<<oMaxXYZ[0]; \
+ for(unsigned int dim=1;dim<Dimension;dim++) \
+ std::cout<<"\t"<<oMaxXYZ[dim]; \
+ std::cout<<std::endl; \
+ } \
+\
+ /* // Output to file */ \
+ if( m_ArgsInfo.general_given) \
+ { \
+ /* // Add to the file */\
+ std::ofstream general(m_ArgsInfo.general_arg); \
+\
+\
+ /* // Numbers of DVF */\
+ for (unsigned int phaseIndex=0; phaseIndex<m_NumberOfFields; phaseIndex++) \
+ general<<phaseIndex<<"\t"<< filenames[phaseIndex]<<std::endl; \
+\
+ general<<std::endl; \
+ general<<"=================================================="<<std::endl; \
+ general<<"|| Original distance between points ||"<<std::endl; \
+ general<<"=================================================="<<std::endl; \
+ general<<std::endl; \
+\
+ general<<std::setprecision(3); \
+\
+ /* // Labels of the columns */\
+ general<<"#DVF\tMean\tSD\t"<<labels[0]; \
+ for(unsigned int dim=1;dim<Dimension;dim++) general<<"\t"<< labels[dim]; \
+ general<<"\tMax \t"<<labels[4]; \
+ for(unsigned int dim=1;dim<Dimension;dim++) general<<"\t"<< labels[dim+4]; \
+ general<<std::endl; \
+ for (unsigned int phaseIndex=0; phaseIndex<m_NumberOfFields; phaseIndex++) \
+ { \
+ general<<phaseIndex; \
+ general<<"\t"<<oMeanList[phaseIndex]<<"\t"<<oStdList[phaseIndex]; \
+ general<<"\t"<<oMeanXYZList[phaseIndex][0]<<"\t"<<oStdXYZList[phaseIndex][0]; \
+ for(unsigned int dim=1;dim<Dimension;dim++) \
+ general<<"\t"<<oMeanXYZList[phaseIndex][dim]<<"\t"<<oStdXYZList[phaseIndex][dim]; \
+ general<<"\t"<<oMaxList[phaseIndex]; \
+ general<<"\t"<<oMaxXYZList[phaseIndex][0]; \
+ for(unsigned int dim=1;dim<Dimension;dim++) \
+ general<<"\t"<<oMaxXYZList[phaseIndex][dim]; \
+ general<<std::endl; \
+ } \
+\
+ /* // General */\
+ general<<std::endl; \
+ general<<"@"; \
+ general<<"\t"<<oMean<<"\t"<<oStd; \
+ general<<"\t"<<oMeanXYZ[0]<<"\t"<<oStdXYZ[0]; \
+ for(unsigned int dim=1;dim<Dimension;dim++) \
+ general<<"\t"<<oMeanXYZ[dim]<<"\t"<<oStdXYZ[dim]; \
+ general<<"\t"<<oMax; \
+ general<<"\t"<<oMaxXYZ[0]; \
+ for(unsigned int dim=1;dim<Dimension;dim++) \
+ general<<"\t"<<oMaxXYZ[dim]; \
+ general<<std::endl; \
+ general.close(); \
+ }
+
+
+#define StatisticsDisplacementsMacro \
+ /* //----------------------------- */\
+ /* // Statistics displacements */\
+ /* //----------------------------- */\
+\
+ /* // Statistics (magnitude) */\
+ MeasureListType dmeanList, dstdList, dmaxList; \
+ ValueType dmean, dstd, dmax; \
+ statisticsFilter->GetStatistics(displacementLists, dmean, dstd, dmax, dmeanList, dstdList, dmaxList); \
+\
+ /* // Statistics (per component) */\
+ VectorListType dmeanXYZList, dstdXYZList,dmaxXYZList; \
+ VectorType dmeanXYZ, dstdXYZ, dmaxXYZ; \
+ statisticsFilter->GetStatistics(displacementLists, dmeanXYZ, dstdXYZ, dmaxXYZ,dmeanXYZList, dstdXYZList, dmaxXYZList); \
+\
+\
+ /* //----------------------------- */\
+ /* // Statistics TRE */\
+ /* //----------------------------- */\
+\
+ /* // Statistics (magnitude) */\
+ MeasureListType meanList, stdList, maxList; \
+ ValueType mean, std, max; \
+ statisticsFilter->GetStatistics(treLists, mean, std, max, meanList, stdList, maxList); \
+\
+ /* // Statistics (per component) */\
+ VectorListType meanXYZList, stdXYZList,maxXYZList; \
+ VectorType meanXYZ, stdXYZ, maxXYZ; \
+ statisticsFilter->GetStatistics(treLists, meanXYZ, stdXYZ, maxXYZ, meanXYZList, stdXYZList, maxXYZList); \
+\
+\
+ /* // Output to screen */\
+ if(m_Verbose) \
+ { \
+\
+ std::cout<<std::endl; \
+ std::cout<<"=================================================="<<std::endl; \
+ std::cout<<"|| Residual distance between points ||"<<std::endl; \
+ std::cout<<"=================================================="<<std::endl; \
+ std::cout<<std::endl; \
+\
+ std::cout<<std::setprecision(3); \
+\
+ /* // Labels of the columns */\
+ std::cout<<"#DVF\tMean\tSD\t"<<labels[0]; \
+ for(unsigned int dim=1;dim<Dimension;dim++) std::cout<<"\t"<< labels[dim]; \
+ std::cout<<"\tMax \t"<<labels[4]; \
+ for(unsigned int dim=1;dim<Dimension;dim++) std::cout<<"\t"<< labels[dim+4]; \
+ std::cout<<std::endl; \
+ for (unsigned int phaseIndex=0; phaseIndex<m_NumberOfFields; phaseIndex++) \
+ { \
+ std::cout<<phaseIndex; \
+ std::cout<<"\t"<<meanList[phaseIndex]<<"\t"<<stdList[phaseIndex]; \
+ std::cout<<"\t"<<meanXYZList[phaseIndex][0]<<"\t"<<stdXYZList[phaseIndex][0]; \
+ for(unsigned int dim=1;dim<Dimension;dim++) \
+ std::cout<<"\t"<<meanXYZList[phaseIndex][dim]<<"\t"<<stdXYZList[phaseIndex][dim]; \
+ std::cout<<"\t"<<maxList[phaseIndex]; \
+ std::cout<<"\t"<<maxXYZList[phaseIndex][0]; \
+ for(unsigned int dim=1;dim<Dimension;dim++) \
+ std::cout<<"\t"<<maxXYZList[phaseIndex][dim]; \
+ std::cout<<std::endl; \
+ } \
+\
+ /* // General */\
+ std::cout<<std::endl; \
+ std::cout<<"@"; \
+ std::cout<<"\t"<<mean<<"\t"<<std; \
+ std::cout<<"\t"<<meanXYZ[0]<<"\t"<<stdXYZ[0]; \
+ for(unsigned int dim=1;dim<Dimension;dim++) \
+ std::cout<<"\t"<<meanXYZ[dim]<<"\t"<<stdXYZ[dim]; \
+ std::cout<<"\t"<<max; \
+ std::cout<<"\t"<<maxXYZ[0]; \
+ for(unsigned int dim=1;dim<Dimension;dim++) \
+ std::cout<<"\t"<<maxXYZ[dim]; \
+ std::cout<<std::endl; \
+ } \
+\
+ /* // Output to file */\
+ if( m_ArgsInfo.general_given) \
+ { \
+ /* // Add to the file */\
+ std::ofstream general(m_ArgsInfo.general_arg, ios_base::app); \
+\
+ general<<std::endl; \
+ general<<"=================================================="<<std::endl; \
+ general<<"|| Residual distance between points ||"<<std::endl; \
+ general<<"=================================================="<<std::endl; \
+ general<<std::endl; \
+\
+ general<<std::setprecision(3); \
+\
+ /* // Labels of the columns */\
+ general<<"#DVF\tMean\tSD\t"<<labels[0]; \
+ for(unsigned int dim=1;dim<Dimension;dim++) general<<"\t"<< labels[dim]; \
+ general<<"\tMax \t"<<labels[4]; \
+ for(unsigned int dim=1;dim<Dimension;dim++) general<<"\t"<< labels[dim+4]; \
+ general<<std::endl; \
+ for (unsigned int phaseIndex=0; phaseIndex<m_NumberOfFields; phaseIndex++) \
+ { \
+ general<<phaseIndex; \
+ general<<"\t"<<meanList[phaseIndex]<<"\t"<<stdList[phaseIndex]; \
+ general<<"\t"<<meanXYZList[phaseIndex][0]<<"\t"<<stdXYZList[phaseIndex][0]; \
+ for(unsigned int dim=1;dim<Dimension;dim++) \
+ general<<"\t"<<meanXYZList[phaseIndex][dim]<<"\t"<<stdXYZList[phaseIndex][dim]; \
+ general<<"\t"<<maxList[phaseIndex]; \
+ general<<"\t"<<maxXYZList[phaseIndex][0]; \
+ for(unsigned int dim=1;dim<Dimension;dim++) \
+ general<<"\t"<<maxXYZList[phaseIndex][dim]; \
+ general<<std::endl; \
+ } \
+\
+ /* // General */\
+ general<<std::endl; \
+ general<<"@"; \
+ general<<"\t"<<mean<<"\t"<<std; \
+ general<<"\t"<<meanXYZ[0]<<"\t"<<stdXYZ[0]; \
+ for(unsigned int dim=1;dim<Dimension;dim++) \
+ general<<"\t"<<meanXYZ[dim]<<"\t"<<stdXYZ[dim]; \
+ general<<"\t"<<max; \
+ general<<"\t"<<maxXYZ[0]; \
+ for(unsigned int dim=1;dim<Dimension;dim++) \
+ general<<"\t"<<maxXYZ[dim]; \
+ general<<std::endl; \
+ general.close(); \
+ } \
+\
+ /* // Output original points */\
+ if( m_ArgsInfo.original_given) \
+ { \
+\
+ std::vector<std::string> filenames; \
+ for (unsigned int i=0;i<m_NumberOfFields;i++) \
+ { \
+ std::ostringstream number_ostr; \
+ number_ostr << i; \
+ std::string number_str = number_ostr.str(); \
+ filenames.push_back(m_ArgsInfo.original_arg+number_str); \
+ } \
+ originalDistanceLists.Write(filenames, m_Verbose ); \
+ } \
+\
+\
+ /* // Output original magnitude points */\
+ if( m_ArgsInfo.originalMag_given) \
+ { \
+ clitk::Lists<itk::FixedArray<ValueType,1> > originalDistanceListsMag=originalDistanceLists.Norm(); \
+ std::vector<std::string> filenames; \
+ for (unsigned int i=0;i<m_NumberOfFields;i++) \
+ { \
+ std::ostringstream number_ostr; \
+ number_ostr << i; \
+ std::string number_str = number_ostr.str(); \
+ filenames.push_back(m_ArgsInfo.originalMag_arg+number_str); \
+ } \
+ originalDistanceListsMag.Write(filenames, m_Verbose ); \
+ } \
+\
+ /* // Output displacement */\
+ if( m_ArgsInfo.displacement_given) \
+ { \
+\
+ std::vector<std::string> filenames; \
+ for (unsigned int i=0;i<m_NumberOfFields;i++) \
+ { \
+ std::ostringstream number_ostr; \
+ number_ostr << i; \
+ std::string number_str = number_ostr.str(); \
+ filenames.push_back(m_ArgsInfo.displacement_arg+number_str); \
+ } \
+ displacementLists.Write(filenames, m_Verbose ); \
+ } \
+\
+\
+ /* // Output displacement magnitude */\
+ if( m_ArgsInfo.displacementMag_given) \
+ { \
+ clitk::Lists<itk::FixedArray<ValueType,1> > displacementListsMag=displacementLists.Norm(); \
+ std::vector<std::string> filenames; \
+ for (unsigned int i=0;i<m_NumberOfFields;i++) \
+ { \
+ std::ostringstream number_ostr; \
+ number_ostr << i; \
+ std::string number_str = number_ostr.str(); \
+ filenames.push_back(m_ArgsInfo.displacementMag_arg+number_str); \
+ } \
+ displacementListsMag.Write(filenames, m_Verbose ); \
+ } \
+\
+ /* // Output warped points */\
+ if( m_ArgsInfo.warp_given) \
+ { \
+\
+ std::vector<std::string> filenames; \
+ for (unsigned int i=0;i<m_ArgsInfo.warp_given;i++) \
+ { \
+ std::ostringstream number_ostr; \
+ number_ostr << i; \
+ std::string number_str = number_ostr.str(); \
+ filenames.push_back(m_ArgsInfo.warp_arg+number_str); \
+ } \
+ warpedPointLists.Write(filenames, m_Verbose ); \
+ } \
+\
+ /* // Output tre */\
+ if( m_ArgsInfo.tre_given) \
+ { \
+\
+ std::vector<std::string> filenames; \
+ for (unsigned int i=0;i<m_NumberOfFields;i++) \
+ { \
+ std::ostringstream number_ostr; \
+ number_ostr << i; \
+ std::string number_str = number_ostr.str(); \
+ filenames.push_back(m_ArgsInfo.tre_arg+number_str); \
+ } \
+ treLists.Write(filenames, m_Verbose ); \
+ } \
+\
+ /* // Output tre mag */\
+ if( m_ArgsInfo.treMag_given) \
+ { \
+ clitk::Lists<itk::FixedArray<ValueType,1> > treMagLists=treLists.Norm(); \
+\
+ std::vector<std::string> filenames; \
+ for (unsigned int i=0;i<m_NumberOfFields;i++) \
+ { \
+ std::ostringstream number_ostr; \
+ number_ostr << i; \
+ std::string number_str = number_ostr.str(); \
+ filenames.push_back(m_ArgsInfo.treMag_arg+number_str); \
+ } \
+ treMagLists.Write(filenames, m_Verbose ); \
+ }
+
+
+
+