From: delmon <delmon>
Date: Thu, 10 Mar 2011 16:16:15 +0000 (+0000)
Subject: Romulo:
X-Git-Tag: v1.2.0~183
X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=4af47ddf6c7819cd561e97b454521c29c21ddab8;p=clitk.git

Romulo:
- Added functionality to transform points based on the BSpline coefficients
  instead of the vector field
- Added new command-line option --coeff (mutual exclusive wrt --vf)
---

diff --git a/registration/clitkCalculateTRE.ggo b/registration/clitkCalculateTRE.ggo
index e6b42f6..7acbe02 100755
--- a/registration/clitkCalculateTRE.ggo
+++ b/registration/clitkCalculateTRE.ggo
@@ -11,7 +11,8 @@ section "Input"
 
 option "ref"   		-   	"List of points in reference"	 	string		yes 			
 option "input"		i	"Lists of points in targets"	  	string  	multiple 	yes
-option "vf"    		-   	"Input deformation fields"	  	string 		multiple	yes
+option "vf"    		-   	"Input deformation fields"	  	string 		multiple	no
+option "coeff"          -       "Input coefficient images"              string          multiple        no
 option "skip"		-	"Skip a phase of a 4D DVF"		int		no	
 
 section "Interpolation"
diff --git a/registration/clitkCalculateTREGenericFilter.cxx b/registration/clitkCalculateTREGenericFilter.cxx
index 1da773b..413823d 100755
--- a/registration/clitkCalculateTREGenericFilter.cxx
+++ b/registration/clitkCalculateTREGenericFilter.cxx
@@ -41,6 +41,9 @@ namespace clitk
   {
     m_Verbose=false;
     m_InputFileName="";
+    m_NumberOfFields=0;
+    m_NumberOfLists=0;
+    m_NumberOfPoints=0;
   }
 
 
@@ -50,39 +53,64 @@ namespace clitk
   void CalculateTREGenericFilter::Update()
   {
     // Read the Dimension and PixelType
+    if (m_ArgsInfo.vf_given && m_ArgsInfo.coeff_given)
+      {
+        std::cout << "Error, You have to supply either vector fields or coefficient images. Supplying both is not allowed. See clitkCalculateTRE --help for details."<<std::endl;;
+        return;
+      }
+    else if (!m_ArgsInfo.vf_given && !m_ArgsInfo.coeff_given)
+      {
+        std::cout << "Error, You have to supply either vector fields or coefficient images. Supplying none is not allowed. See clitkCalculateTRE --help for details."<<std::endl;;
+        return;
+      }
+    else if (m_ArgsInfo.vf_given)
+      m_InputFileName=m_ArgsInfo.vf_arg[0];
+    else
+      {
+        m_InputFileName=m_ArgsInfo.coeff_arg[0];
+      }
+
     int Dimension, Components;
     std::string PixelType;
     ReadImageDimensionAndPixelType(m_InputFileName, Dimension, PixelType, Components);
+    if (Dimension < 2 || Dimension > 4)
+      {
+        std::cout<<"Error, Only for 2, 3 or 4  Dimensions!!!"<<std::endl ;
+        return;
+      }
 
-    
     // Call UpdateWithDim
-    if(Dimension==2) ReadVectorFields<2,2>();
-    else if(Dimension==3) ReadVectorFields<3,3>();
-    else if (Dimension==4)ReadVectorFields<4,3>(); 
-    else 
+    if (m_ArgsInfo.vf_given)
       {
-	std::cout<<"Error, Only for 2, 3 or 4  Dimensions!!!"<<std::endl ;
-	return;
+        if(Dimension==2) ReadVectorFields<2,2>();
+        else 
+        if(Dimension==3) ReadVectorFields<3,3>();
+        else if (Dimension==4)ReadVectorFields<4,3>(); 
       }
+    else
+      {
+        if(Dimension==2) ReadCoefficientImages<2,2>();
+        else 
+        if(Dimension==3) ReadCoefficientImages<3,3>();
+        else if (Dimension==4)ReadCoefficientImages<4,3>(); 
+      }        
   }
 
-
-
   //-----------------------------
   // Process DVF
   //-----------------------------
   template< >
   void 
-  CalculateTREGenericFilter::ProcessVectorFields<4,3>(std::vector<itk::Image<itk::Vector<float, 3>, 4>::Pointer > dvfs,  char** filenames )
+  CalculateTREGenericFilter::ProcessVectorFields<4,3>(std::vector<itk::Image<itk::Vector<double, 3>, 4>::Pointer > dvfs,  char** filenames )
   {
     // Typedefs
-    typedef itk::Vector<float,3> PixelType;
+    typedef itk::Vector<double,3> PixelType;
     typedef itk::Image<PixelType, 4> InputImageType;
     typedef itk::Image<PixelType, 3> OutputImageType;
 
     // IO
     InputImageType::Pointer input=dvfs[0];
-    std::vector<itk::Image<itk::Vector<float, 3>, 3>::Pointer > new_dvfs;
+    std::vector<itk::Image<itk::Vector<double, 3>, 3>::Pointer > new_dvfs;
     
     // Split vector field
     typedef itk::ExtractImageFilter<InputImageType,OutputImageType> FilterType;
@@ -124,11 +152,69 @@ namespace clitk
       }
 
     // Update
-    this->UpdateWithDim<3>(new_dvfs, new_filenames); 
+    this->UpdateDVFWithDim<3>(new_dvfs, new_filenames); 
 
   }
 
+  //-----------------------------
+  // Process Coefficient images
+  //-----------------------------
+  template< >
+  void 
+  CalculateTREGenericFilter::ProcessCoefficientImages<4,3>(std::vector<itk::Image<itk::Vector<double, 3>, 4>::Pointer > coeffs,  char** filenames )
+  {
+    // Typedefs
+    typedef itk::Vector<double,3> PixelType;
+    typedef itk::Image<PixelType, 4> InputImageType;
+    typedef itk::Image<PixelType, 3> OutputImageType;
 
+    // IO
+    InputImageType::Pointer input=coeffs[0];
+    std::vector<itk::Image<itk::Vector<double, 3>, 3>::Pointer > new_coeffs;
+    
+    // Split vector field
+    typedef itk::ExtractImageFilter<InputImageType,OutputImageType> FilterType;
+    unsigned int splitDimension=3;
+
+    // Make new file names
+    std::vector<std::string> new_filenames;
+    std::string base = filenames[0];
+
+    // Set the extract region
+    InputImageType::SizeType size=input->GetLargestPossibleRegion().GetSize();
+    size[splitDimension]=0;
+    InputImageType::RegionType extracted_region;
+    extracted_region.SetSize(size);
+    InputImageType::IndexType index=input->GetLargestPossibleRegion().GetIndex();
+    
+  
+    // Loop
+    for (unsigned int i=0;i<input->GetLargestPossibleRegion().GetSize()[splitDimension];i++)
+      {
+        
+        // Skip?
+        if (m_ArgsInfo.skip_given && i==(unsigned int) m_ArgsInfo.skip_arg) continue;
+
+        // extract dvf
+         FilterType::Pointer filter= FilterType::New();
+        filter->SetInput(input);
+        index[splitDimension]=i;
+        extracted_region.SetIndex(index);
+        filter->SetExtractionRegion(extracted_region);
+        filter->Update();
+        new_coeffs.push_back(filter->GetOutput());
+        
+        // make name
+        std::ostringstream number_dvf;
+        number_dvf << i;
+        std::string number =  number_dvf.str();
+        new_filenames.push_back(base+"_"+number);
+      }
+
+    // Update
+    this->UpdateCoeffsWithDim<3>(new_coeffs, new_filenames); 
+
+  }
 
 
   
diff --git a/registration/clitkCalculateTREGenericFilter.h b/registration/clitkCalculateTREGenericFilter.h
index d9caea1..e504817 100755
--- a/registration/clitkCalculateTREGenericFilter.h
+++ b/registration/clitkCalculateTREGenericFilter.h
@@ -67,9 +67,11 @@ namespace clitk
     itkTypeMacro( CalculateTREGenericFilter, LightObject );
 
 
-    //----------------------------------------
+    //-----------------------------
     // Typedefs
-    //----------------------------------------
+    //-----------------------------
+    typedef double ValueType;
+    typedef std::vector<ValueType> MeasureListType;
 
 
     //----------------------------------------
@@ -79,7 +81,6 @@ namespace clitk
     {
       m_ArgsInfo=a;
       m_Verbose=m_ArgsInfo.verbose_flag;
-      m_InputFileName=m_ArgsInfo.vf_arg[0];
     }
     
     
@@ -101,9 +102,12 @@ namespace clitk
     // Templated members
     //----------------------------------------  
     template <unsigned int Dimension,unsigned int Components> void ReadVectorFields(void);
-    template <unsigned int Dimension,unsigned int Components> void ProcessVectorFields(std::vector< typename itk::Image<itk::Vector<float, Components>, Dimension>::Pointer > dvfs, char** filenames);
-    template <unsigned int Dimension> void UpdateWithDim( std::vector<typename itk::Image<itk::Vector<float, Dimension>, Dimension>::Pointer > dvfs,  std::vector<std::string> filenames);
-    
+    template <unsigned int Dimension,unsigned int Components> void ProcessVectorFields(std::vector< typename itk::Image<itk::Vector<double, Components>, Dimension>::Pointer > dvfs, char** filenames);
+    template <unsigned int Dimension> void UpdateDVFWithDim( std::vector<typename itk::Image<itk::Vector<ValueType, Dimension>, Dimension>::Pointer > dvfs,  std::vector<std::string> filenames);
+    template <unsigned int Dimension,unsigned int Components> void ReadCoefficientImages(void);
+    template <unsigned int Dimension,unsigned int Components> void ProcessCoefficientImages(std::vector< typename itk::Image<itk::Vector<double, Components>, Dimension>::Pointer > dvfs, char** filenames);
+    template <unsigned int Dimension> void UpdateCoeffsWithDim( std::vector<typename itk::Image<itk::Vector<ValueType, Dimension>, Dimension>::Pointer > dvfs,  std::vector<std::string> filenames);
+    template<unsigned int Dimension> void BuildPointLists(std::vector<std::string>& filenames);
 
     //----------------------------------------  
     // Data members
@@ -111,11 +115,14 @@ namespace clitk
     args_info_clitkCalculateTRE m_ArgsInfo;
     bool m_Verbose;
     std::string m_InputFileName;
+    unsigned int m_NumberOfFields;
+    unsigned int m_NumberOfLists;
+    unsigned int m_NumberOfPoints;
 
   };
 
 
-} // end namespace clitk
+} // end namespace clitk 
 
 #ifndef ITK_MANUAL_INSTANTIATION
 #include "clitkCalculateTREGenericFilter.txx"
diff --git a/registration/clitkCalculateTREGenericFilter.txx b/registration/clitkCalculateTREGenericFilter.txx
index 1198ac8..f75a3f2 100755
--- a/registration/clitkCalculateTREGenericFilter.txx
+++ b/registration/clitkCalculateTREGenericFilter.txx
@@ -18,6 +18,10 @@
 #ifndef clitkCalculateTREGenericFilter_txx
 #define clitkCalculateTREGenericFilter_txx
 
+
+#include "clitkBSplineDeformableTransform.h"
+
+
 /* =================================================
  * @file   clitkCalculateTREGenericFilter.txx
  * @author 
@@ -28,6 +32,495 @@
  ===================================================*/
 
 
+//=====================================================================================
+// 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
 {
   
@@ -38,22 +531,24 @@ 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);
+
   }
 
   //-----------------------------
@@ -61,279 +556,84 @@ namespace clitk
   //-----------------------------
   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
@@ -341,10 +641,10 @@ namespace clitk
     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];
@@ -364,243 +664,65 @@ namespace clitk
 	  }
       }
 
+    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