/*=========================================================================
  Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv

  Authors belong to: 
  - University of LYON              http://www.universite-lyon.fr/
  - Léon Bérard cancer center       http://www.centreleonberard.fr
  - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr

  This software is distributed WITHOUT ANY WARRANTY; without even
  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  PURPOSE.  See the copyright notices for more information.

  It is distributed under dual licence

  - BSD        See included LICENSE.txt file
  - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
===========================================================================**/
#ifndef clitkCalculateTREGenericFilter_txx
#define clitkCalculateTREGenericFilter_txx


#include "clitkBSplineDeformableTransform.h"


/* =================================================
 * @file   clitkCalculateTREGenericFilter.txx
 * @author 
 * @date   
 * 
 * @brief 
 * 
 ===================================================*/


//=====================================================================================
// 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
{
  
  //-----------------------------
  // Read DVF
  //-----------------------------
  template<unsigned int Dimension, unsigned int Components >
  void 
  CalculateTREGenericFilter::ReadVectorFields(void)
  {
    // if using dvfs
    typedef itk::Vector<ValueType, Components> VectorType;
    typedef itk::Image<VectorType, Dimension> DeformationFieldType;

    typedef itk::ImageFileReader<DeformationFieldType> DvfReaderType;    
    std::vector<typename DeformationFieldType::Pointer> dvfs;

    for (unsigned int i=0; i< m_ArgsInfo.vf_given; i++)
      {
        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);

  }

  //-----------------------------
  // Process DVF
  //-----------------------------
  template<unsigned int Dimension, unsigned int Components >
  void 
  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]);
    UpdateDVFWithDim<Dimension>(dvfs, new_filenames); 
  }

    
  //-----------------------------
  // Read Coefficient Images
  //-----------------------------
  template<unsigned int Dimension, unsigned int Components >
  void 
  CalculateTREGenericFilter::ReadCoefficientImages(void)
  {
    // if using coefficient images
    typedef itk::Vector<ValueType, Components> VectorType;
    typedef itk::Image<VectorType, Dimension> CoefficientImageType;
    
    typedef itk::ImageFileReader<CoefficientImageType> CoeffReaderType;
    std::vector<typename CoefficientImageType::Pointer> coeffs;

    for (unsigned int i=0; i< m_ArgsInfo.coeff_given; i++)
      {
        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);
  }

  //-----------------------------
  // 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); 
  }
  

  //-------------------------------------------------------------------
  // 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;

    TypedefsMacro;
    BuildPointListsMacro(filenames);
    OriginalDistancesMacro;

    //=====================================================================================
    // Distance between points after warp
    //=====================================================================================
    
    //-----------------------------
    // Interpolator
    //-----------------------------
    typedef clitk::GenericVectorInterpolator<args_info_clitkCalculateTRE,DeformationFieldType, ValueType> GenericVectorInterpolatorType;
    typename GenericVectorInterpolatorType::Pointer genericInterpolator=GenericVectorInterpolatorType::New();
    genericInterpolator->SetArgsInfo(m_ArgsInfo);
    typedef itk::VectorInterpolateImageFunction<DeformationFieldType, ValueType> InterpolatorType; 
    typename InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer();
    

    PointListsType warpedPointLists(m_NumberOfFields);
    VectorListsType treLists(m_NumberOfFields), displacementLists(m_NumberOfFields);
    
    //-----------------------------
    // Calculate residual distances
    //-----------------------------
    VectorType displacement, tre;
    PointType  warpedPoint;

    for (unsigned int phaseIndex=0; phaseIndex<m_NumberOfFields; phaseIndex++)
      { 
	interpolator->SetInputImage( dvfs[phaseIndex] );
	for (unsigned int pointIndex=0; pointIndex<m_NumberOfPoints; pointIndex++)   
	  {
	    // Reference
	    referencePoint=referencePointList[pointIndex];

	    // Inside?
	    if ( interpolator->IsInsideBuffer(referencePoint) )
	      displacement=interpolator->Evaluate(referencePoint); 
	    else
	      displacement.Fill(0.0); 

	    // Warp
	    warpedPoint=referencePoint+displacement;
	    displacementLists[phaseIndex].push_back(displacement);
	    warpedPointLists[phaseIndex].push_back(warpedPoint);
	    tre=pointLists[phaseIndex][pointIndex]-warpedPoint;
	    treLists[phaseIndex].push_back(tre);
	  }
      }

    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
    //=====================================================================================
    
    //-----------------------------
    // 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);
    
    //-----------------------------
    // Calculate residual distances
    //-----------------------------
    VectorType displacement, tre;
    PointType  warpedPoint;

    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);
          }
      }

    StatisticsDisplacementsMacro;
  }
  
}//end clitk
#endif //#define clitkCalculateTREGenericFilter_txx