Authors belong to:
- University of LYON http://www.universite-lyon.fr/
- - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr
+ - Léon Bérard cancer center http://www.centreleonberard.fr
- CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
This software is distributed WITHOUT ANY WARRANTY; without even
- BSD See included LICENSE.txt file
- CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
-======================================================================-====*/
+===========================================================================**/
#ifndef clitkCalculateTREGenericFilter_txx
#define clitkCalculateTREGenericFilter_txx
+
+#include "clitkBSplineDeformableTransform.h"
+
+
/* =================================================
* @file clitkCalculateTREGenericFilter.txx
* @author
===================================================*/
+//=====================================================================================
+// 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 ); \
+ }
+
+
+
+
namespace clitk
{
void
CalculateTREGenericFilter::ReadVectorFields(void)
{
-
- typedef itk::Vector<float, Components> VectorType;
+ // if using dvfs
+ typedef itk::Vector<ValueType, Components> VectorType;
typedef itk::Image<VectorType, Dimension> DeformationFieldType;
- typedef itk::ImageFileReader<DeformationFieldType> InputReaderType;
+ typedef itk::ImageFileReader<DeformationFieldType> DvfReaderType;
std::vector<typename DeformationFieldType::Pointer> dvfs;
+
for (unsigned int i=0; i< m_ArgsInfo.vf_given; i++)
{
- typename InputReaderType::Pointer reader = InputReaderType::New();
- reader->SetFileName( m_ArgsInfo.vf_arg[i]);
- if (m_Verbose) std::cout<<"Reading vector field "<< i <<"..."<<std::endl;
- reader->Update();
- dvfs.push_back( reader->GetOutput() );
+ typename DvfReaderType::Pointer reader = DvfReaderType::New();
+ reader->SetFileName( m_ArgsInfo.vf_arg[i]);
+ if (m_Verbose) std::cout<<"Reading vector field "<< i <<"..."<<std::endl;
+ reader->Update();
+ dvfs.push_back( reader->GetOutput() );
}
-
+
ProcessVectorFields<Dimension, Components>(dvfs, m_ArgsInfo.vf_arg);
+
}
//-----------------------------
//-----------------------------
template<unsigned int Dimension, unsigned int Components >
void
- CalculateTREGenericFilter::ProcessVectorFields(std::vector<typename itk::Image<itk::Vector<float, Components>, Dimension>::Pointer > dvfs, char** filenames )
+ CalculateTREGenericFilter::ProcessVectorFields(std::vector<typename itk::Image<itk::Vector<CalculateTREGenericFilter::ValueType, Components>, Dimension>::Pointer > dvfs, char** filenames )
{
-
std::vector<std::string> new_filenames;
for (unsigned int i=0;i<m_ArgsInfo.vf_given;i++)
new_filenames.push_back(filenames[i]);
- UpdateWithDim<Dimension>(dvfs, new_filenames);
+ UpdateDVFWithDim<Dimension>(dvfs, new_filenames);
}
- //-------------------------------------------------------------------
- // Update with the number of dimensions
- //-------------------------------------------------------------------
- template<unsigned int Dimension>
+ //-----------------------------
+ // Read Coefficient Images
+ //-----------------------------
+ template<unsigned int Dimension, unsigned int Components >
void
- CalculateTREGenericFilter::UpdateWithDim(std::vector<typename itk::Image<itk::Vector<float, Dimension>, Dimension>::Pointer > dvfs, std::vector<std::string> filenames)
+ CalculateTREGenericFilter::ReadCoefficientImages(void)
{
- if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D..."<<std::endl;
-
- //-----------------------------
- // Typedefs
- //-----------------------------
- typedef double ValueType;
- typedef std::vector<ValueType> MeasureListType;
-
- typedef itk::Point<double, Dimension> PointType;
- typedef clitk::List<PointType> PointListType;
- typedef clitk::Lists<PointType> PointListsType;
-
- typedef itk::Vector<double, Dimension> VectorType;
- typedef clitk::List<VectorType> VectorListType;
- typedef clitk::Lists<VectorType> VectorListsType;
+ // if using coefficient images
+ typedef itk::Vector<ValueType, Components> VectorType;
+ typedef itk::Image<VectorType, Dimension> CoefficientImageType;
- typedef itk::Vector<float, Dimension> DeformationVectorType;
- typedef itk::Image<DeformationVectorType, Dimension> DeformationFieldType;
+ typedef itk::ImageFileReader<CoefficientImageType> CoeffReaderType;
+ std::vector<typename CoefficientImageType::Pointer> coeffs;
- //-----------------------------
- // Number of inputs
- //-----------------------------
- unsigned int numberOfFields=filenames.size();
- unsigned int numberOfLists=m_ArgsInfo.input_given;
- if ( (numberOfLists!=numberOfFields) && (numberOfLists!=1) )
- {
- std::cerr<<"Error: Number of lists (="<<numberOfLists<<") different from number of DVF's (="<<numberOfFields<<") !"<<std::endl;
- return;
- }
- else if (numberOfLists==1)
- {
- if (m_Verbose) std::cout<<numberOfFields<<" DVFs and 1 point list given..."<<std::endl;
- }
- else
+ for (unsigned int i=0; i< m_ArgsInfo.coeff_given; i++)
{
- if (m_Verbose) std::cout<<numberOfLists<<" point lists and DVFs given..."<<std::endl;
+ typename CoeffReaderType::Pointer reader = CoeffReaderType::New();
+ reader->SetFileName( m_ArgsInfo.coeff_arg[i]);
+ if (m_Verbose) std::cout<<"Reading coefficient image "<< i <<"..."<<std::endl;
+ reader->Update();
+ coeffs.push_back( reader->GetOutput() );
}
+
+ ProcessCoefficientImages<Dimension, Components>(coeffs, m_ArgsInfo.coeff_arg);
+ }
-
- //-----------------------------
- // Input point lists
- //-----------------------------
- PointListsType pointLists;
- unsigned int numberOfPoints=0;
- for (unsigned int i=0; i<numberOfFields; i++)
- {
- // Read the lists
- if (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) );
-
+ //-----------------------------
+ // Process Coefficient Images
+ //-----------------------------
+ template<unsigned int Dimension, unsigned int Components >
+ void
+ CalculateTREGenericFilter::ProcessCoefficientImages(std::vector<typename itk::Image<itk::Vector<CalculateTREGenericFilter::ValueType, Components>, Dimension>::Pointer > coeffs, char** filenames )
+ {
+ std::vector<std::string> new_filenames;
+ for (unsigned int i=0;i<m_ArgsInfo.coeff_given;i++)
+ new_filenames.push_back(filenames[i]);
+ UpdateCoeffsWithDim<Dimension>(coeffs, new_filenames);
+ }
+
- // Verify the number of points
- if (i==0) numberOfPoints=pointLists[i].size();
- else
- {
- if (numberOfPoints!=pointLists[i].size())
- {
- std::cerr<<"Size of first list ("<<numberOfPoints
- <<") is different from size of list "<<i
- <<" ("<<pointLists[i].size()<<")..."<<std::endl;
- return;
- }
- }
- }
-
+ //-------------------------------------------------------------------
+ // Update DVF with the number of dimensions
+ //-------------------------------------------------------------------
+ template<unsigned int Dimension>
+ void
+ CalculateTREGenericFilter::UpdateDVFWithDim(std::vector<typename itk::Image<itk::Vector<CalculateTREGenericFilter::ValueType, Dimension>, Dimension>::Pointer > dvfs, std::vector<std::string> filenames)
+ {
+ if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D..."<<std::endl;
- PointListType referencePointList;
- if (m_Verbose) std::cout<<"Reference point list:"<<std::endl;
- referencePointList=PointListType(m_ArgsInfo.ref_arg, m_Verbose);
- if (numberOfPoints!=referencePointList.size())
- {
- std::cerr<<"Size of the first list ("<<numberOfPoints
- <<") is different from size of the reference list ("
- << referencePointList.size() <<")..."<<std::endl;
- return;
- }
+ TypedefsMacro;
+ BuildPointListsMacro(filenames);
+ OriginalDistancesMacro;
+ //=====================================================================================
+ // Distance between points after warp
+ //=====================================================================================
//-----------------------------
// Interpolator
//-----------------------------
- typedef clitk::GenericVectorInterpolator<args_info_clitkCalculateTRE,DeformationFieldType, double> GenericVectorInterpolatorType;
+ typedef clitk::GenericVectorInterpolator<args_info_clitkCalculateTRE,DeformationFieldType, ValueType> GenericVectorInterpolatorType;
typename GenericVectorInterpolatorType::Pointer genericInterpolator=GenericVectorInterpolatorType::New();
genericInterpolator->SetArgsInfo(m_ArgsInfo);
- typedef itk::VectorInterpolateImageFunction<DeformationFieldType, double> InterpolatorType;
+ typedef itk::VectorInterpolateImageFunction<DeformationFieldType, ValueType> InterpolatorType;
typename InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer();
- //=====================================================================================
- // Original distance between points
- //=====================================================================================
- VectorListsType originalDistanceLists(numberOfFields);
-
- //-----------------------------
- // Calculate original distances
- //-----------------------------
- PointType referencePoint;
- VectorType distance;
- for (unsigned int phaseIndex=0; phaseIndex<numberOfFields; phaseIndex++)
- {
- for (unsigned int pointIndex=0; pointIndex<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<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<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<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<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();
- }
-
-
- //=====================================================================================
- // Distance between points after warp
- //=====================================================================================
- PointListsType warpedPointLists(numberOfFields);
- VectorListsType treLists(numberOfFields), displacementLists(numberOfFields);
+ PointListsType warpedPointLists(m_NumberOfFields);
+ VectorListsType treLists(m_NumberOfFields), displacementLists(m_NumberOfFields);
//-----------------------------
// Calculate residual distances
VectorType displacement, tre;
PointType warpedPoint;
- for (unsigned int phaseIndex=0; phaseIndex<numberOfFields; phaseIndex++)
+ for (unsigned int phaseIndex=0; phaseIndex<m_NumberOfFields; phaseIndex++)
{
- interpolator->SetInputImage( dvfs[phaseIndex]);
- for (unsigned int pointIndex=0; pointIndex<numberOfPoints; pointIndex++)
+ interpolator->SetInputImage( dvfs[phaseIndex] );
+ for (unsigned int pointIndex=0; pointIndex<m_NumberOfPoints; pointIndex++)
{
// Reference
referencePoint=referencePointList[pointIndex];
}
}
+ StatisticsDisplacementsMacro;
+ }
+ //-------------------------------------------------------------------
+ // Update coefficients with the number of dimensions
+ //-------------------------------------------------------------------
+ template<unsigned int Dimension>
+ void
+ CalculateTREGenericFilter::UpdateCoeffsWithDim(std::vector<typename itk::Image<itk::Vector<CalculateTREGenericFilter::ValueType, Dimension>, Dimension>::Pointer > coeffs, std::vector<std::string> filenames)
+ {
+ // call macros
+ TypedefsMacro;
+ BuildPointListsMacro(filenames);
+ OriginalDistancesMacro;
+
+ //=====================================================================================
+ // Distance between points after warp
+ //=====================================================================================
+
//-----------------------------
- // Statistics displacements
+ // Transform to get the transformed points from the coeff images
//-----------------------------
+ typedef clitk::BSplineDeformableTransform<ValueType, Dimension, Dimension> TransformType;
+ typename TransformType::Pointer transform = TransformType::New();
+
+ PointListsType warpedPointLists(m_NumberOfFields);
+ VectorListsType treLists(m_NumberOfFields), displacementLists(m_NumberOfFields);
- // 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
+ // Calculate residual distances
//-----------------------------
-
- // 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<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<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<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<double,1> > originalDistanceListsMag=originalDistanceLists.Norm();
- std::vector<std::string> filenames;
- for (unsigned int i=0;i<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<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<double,1> > displacementListsMag=displacementLists.Norm();
- std::vector<std::string> filenames;
- for (unsigned int i=0;i<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)
- {
+ VectorType displacement, tre;
+ PointType warpedPoint;
- std::vector<std::string> filenames;
- for (unsigned int i=0;i<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 );
+ for (unsigned int phaseIndex=0; phaseIndex<m_NumberOfFields; phaseIndex++)
+ {
+ // set the initial parameters (b-spline coefficients from cmd-line
+ transform->SetCoefficientImage(coeffs[phaseIndex]);
+
+ // transform the input points
+ for (unsigned int pointIndex=0; pointIndex<m_NumberOfPoints; pointIndex++)
+ {
+ // Reference
+ referencePoint=referencePointList[pointIndex];
+
+ // Transform
+ warpedPoint=transform->TransformPoint(referencePoint);
+
+ // Warp
+ displacement=warpedPoint-referencePoint;
+ displacementLists[phaseIndex].push_back(displacement);
+ warpedPointLists[phaseIndex].push_back(warpedPoint);
+ tre=pointLists[phaseIndex][pointIndex]-warpedPoint;
+ treLists[phaseIndex].push_back(tre);
+ }
}
- // Output tre mag
- if( m_ArgsInfo.treMag_given)
- {
- clitk::Lists<itk::FixedArray<double,1> > treMagLists=treLists.Norm();
-
- std::vector<std::string> filenames;
- for (unsigned int i=0;i<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 );
- }
-
+ StatisticsDisplacementsMacro;
}
-
-
+
}//end clitk
-
#endif //#define clitkCalculateTREGenericFilter_txx
git://git.creatis.insa-lyon.fr / clitk.git / blobdiff