]> Creatis software - clitk.git/commitdiff
Romulo:
authordelmon <delmon>
Thu, 10 Mar 2011 16:16:15 +0000 (16:16 +0000)
committerdelmon <delmon>
Thu, 10 Mar 2011 16:16:15 +0000 (16:16 +0000)
- 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)

registration/clitkCalculateTRE.ggo
registration/clitkCalculateTREGenericFilter.cxx
registration/clitkCalculateTREGenericFilter.h
registration/clitkCalculateTREGenericFilter.txx

index e6b42f694481561dfbb5bcba2276e22b7b710585..7acbe02a2d6a57268ba7fb7af3da4c59c15c7f03 100755 (executable)
@@ -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"
index 1da773bb2e97ceae7cb572e1e5bdb2f2bc223e78..413823d5c0a966d38b988ac35d99cda0cf1bbf45 100755 (executable)
@@ -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); 
+
+  }
 
 
   
index d9caea1890917ecca206c91d78637bbc35cfe2a2..e504817af7fdfa791aa20fe9aa33e01daf3cfdf7 100755 (executable)
@@ -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"
index 1198ac83cae18e6e835dfc20fe7c14e3f1fbb839..f75a3f29b152894f6b048223e92aa43e49c37bff 100755 (executable)
 #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
 {
   
@@ -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