registration/clitkCalculateTREGenericFilter.txx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   5   - University of LYON              http://www.universite-lyon.fr/
   6   - Léon Bérard cancer center       http://www.centreleonberard.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   9   This software is distributed WITHOUT ANY WARRANTY; without even
  10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11   PURPOSE.  See the copyright notices for more information.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17 ===========================================================================**/
  18 #ifndef clitkCalculateTREGenericFilter_txx
  19 #define clitkCalculateTREGenericFilter_txx
  20
  21
  22 #include "clitkBSplineDeformableTransform.h"
  23
  24
  25 /* =================================================
  26  * @file   clitkCalculateTREGenericFilter.txx
  27  * @author
  28  * @date
  29  *
  30  * @brief
  31  *
  32  ===================================================*/
  33
  34
  35 //=====================================================================================
  36 // Macros
  37 //  Rômulo Pinho - 25/02/2011
  38 //  The macros below avoid the duplication of code after introducing the function
  39 // UpdateCoeffsWithDim(). Although the code generated by these functions will be
  40 // duplicated anyway (since they're template functions) the macros make their sources
  41 // more intelligible.
  42 // Ideally, these macros should be made member functions, with the typedefs inside the
  43 // the filter class, but that would change the class' interface, since it would need
  44 // to become a template class (with parameter Dimension). To avoid breaking legacy code,
  45 // I opted for the macros, although I don't quite like them (especially long ones like
  46 // these). In any case, it's advisable to think of a better solution, even if it breaks
  47 // the class' interface.
  48 //=====================================================================================
  49
  50 //-----------------------------
  51 // Typedefs
  52 //-----------------------------
  53 #define TypedefsMacro \
  54 \
  55     typedef itk::Point<ValueType, Dimension> PointType; \
  56     typedef clitk::List<PointType> PointListType; \
  57     typedef clitk::Lists<PointType> PointListsType; \
  58 \
  59     typedef itk::Vector<ValueType, Dimension> VectorType; \
  60     typedef clitk::List<VectorType> VectorListType; \
  61     typedef clitk::Lists<VectorType> VectorListsType; \
  62 \
  63     typedef itk::Vector<ValueType, Dimension> DeformationVectorType; \
  64     typedef itk::Image<DeformationVectorType, Dimension> DeformationFieldType;
  65
  66 //-------------------------------------------------------------------
  67 // Build the point lists necessary for the computation of distances.
  68 //-------------------------------------------------------------------
  69 #define BuildPointListsMacro(filenames) \
  70 \
  71   m_NumberOfFields=filenames.size(); \
  72   m_NumberOfLists=m_ArgsInfo.input_given; \
  73   if ( (m_NumberOfLists!=m_NumberOfFields)  &&  (m_NumberOfLists!=1) ) \
  74     { \
  75       std::cerr<<"Error: Number of lists (="<<m_NumberOfLists<<") different from number of DVF's (="<<m_NumberOfFields<<") !"<<std::endl; \
  76       return; \
  77     } \
  78   else if (m_NumberOfLists==1) \
  79     { \
  80       if (m_Verbose) std::cout<<m_NumberOfFields<<" DVFs and 1 point list given..."<<std::endl; \
  81     } \
  82   else \
  83     { \
  84       if (m_Verbose) std::cout<<m_NumberOfLists<<" point lists and DVFs given..."<<std::endl; \
  85     } \
  86 \
  87 \
  88   PointListsType pointLists; \
  89   m_NumberOfPoints=0; \
  90   for (unsigned int i=0; i<m_NumberOfFields; i++) \
  91     { \
  92       /* Read the lists */\
  93       if (m_NumberOfLists==1) \
  94         pointLists.push_back(PointListType(m_ArgsInfo.input_arg[0], m_Verbose) ); \
  95       else \
  96         pointLists.push_back(PointListType(m_ArgsInfo.input_arg[i], m_Verbose) ); \
  97 \
  98 \
  99       /* Verify the number of points */\
 100       if (i==0) m_NumberOfPoints=pointLists[i].size(); \
 101       else \
 102         { \
 103           if  (m_NumberOfPoints!=pointLists[i].size()) \
 104             { \
 105               std::cerr<<"Size of first list ("<<m_NumberOfPoints \
 106                         <<") is different from size of list "<<i \
 107                         <<" ("<<pointLists[i].size()<<")..."<<std::endl; \
 108               return; \
 109             } \
 110         } \
 111     }    \
 112 \
 113 \
 114   PointListType referencePointList; \
 115   if (m_Verbose) std::cout<<"Reference point list:"<<std::endl; \
 116   referencePointList=PointListType(m_ArgsInfo.ref_arg, m_Verbose); \
 117   if  (m_NumberOfPoints!=referencePointList.size()) \
 118     { \
 119       std::cerr<<"Size of the first list ("<<m_NumberOfPoints \
 120                 <<") is different from size of the reference list (" \
 121                 << referencePointList.size() <<")..."<<std::endl; \
 122       return; \
 123     }
 124
 125 /* //===================================================================================== */\
 126 /* // Original distance between points */\
 127 /* //===================================================================================== */\
 128 #define OriginalDistancesMacro \
 129     VectorListsType originalDistanceLists(m_NumberOfFields); \
 130 \
 131     /* //----------------------------- */\
 132     /* // Calculate original distances */\
 133     /* //----------------------------- */\
 134     PointType referencePoint; \
 135     VectorType distance; \
 136     for (unsigned int phaseIndex=0; phaseIndex<m_NumberOfFields; phaseIndex++) \
 137       { \
 138         for (unsigned int pointIndex=0; pointIndex<m_NumberOfPoints; pointIndex++)    \
 139           { \
 140             referencePoint=referencePointList[pointIndex]; \
 141             distance=pointLists[phaseIndex][pointIndex]-referencePoint;  \
 142             originalDistanceLists[phaseIndex].push_back(distance);  \
 143           } \
 144       } \
 145 \
 146 \
 147     /* //----------------------------- */\
 148     /* // Statistics Original */\
 149     /* //----------------------------- */\
 150     typedef DeformationListStatisticsFilter<VectorType> StatisticsFilterType; \
 151     typename StatisticsFilterType::Pointer statisticsFilter=StatisticsFilterType::New(); \
 152 \
 153     /* // Statistics  (magnitude) */ \
 154     MeasureListType oMeanList, oStdList, oMaxList; \
 155     ValueType oMean, oStd, oMax; \
 156     statisticsFilter->GetStatistics(originalDistanceLists, oMean, oStd, oMax, oMeanList, oStdList, oMaxList); \
 157 \
 158     /* // Statistics  (per component) */\
 159     VectorListType oMeanXYZList, oStdXYZList,oMaxXYZList; \
 160     VectorType oMeanXYZ, oStdXYZ, oMaxXYZ; \
 161     statisticsFilter->GetStatistics(originalDistanceLists, oMeanXYZ, oStdXYZ, oMaxXYZ, oMeanXYZList, oStdXYZList, oMaxXYZList); \
 162 \
 163 \
 164     /* //----------------------------- */\
 165     /* // Output */\
 166     /* //----------------------------- */\
 167     std::vector<std::string> labels; \
 168     labels.push_back("MeanX\tSDX"); \
 169     labels.push_back("MeanY\tSDY"); \
 170     labels.push_back("MeanZ\tSDZ"); \
 171     labels.push_back("MeanT\tSDT"); \
 172     labels.push_back("MaxX"); \
 173     labels.push_back("MaxY"); \
 174     labels.push_back("MaxZ"); \
 175     labels.push_back("MaxT"); \
 176 \
 177 \
 178     /* // Output to screen */\
 179     if(m_Verbose) \
 180       { \
 181 \
 182         /* // Numbers of DVF */\
 183         std::cout<<"# Number\tDVF"<<std::endl; \
 184         for (unsigned int phaseIndex=0; phaseIndex<m_NumberOfFields; phaseIndex++) \
 185           std::cout<<phaseIndex<<"\t"<<filenames[phaseIndex]<<std::endl; \
 186 \
 187 \
 188         std::cout<<std::endl; \
 189         std::cout<<"=================================================="<<std::endl; \
 190         std::cout<<"||       Original distance between points       ||"<<std::endl; \
 191         std::cout<<"=================================================="<<std::endl; \
 192         std::cout<<std::endl; \
 193 \
 194         std::cout<<std::setprecision(3); \
 195 \
 196 \
 197 \
 198         /* // Labels of the columns */\
 199         std::cout<<"#DVF\tMean\tSD\t"<<labels[0]; \
 200         for(unsigned int dim=1;dim<Dimension;dim++)     std::cout<<"\t"<< labels[dim]; \
 201         std::cout<<"\tMax \t"<<labels[4]; \
 202         for(unsigned int dim=1;dim<Dimension;dim++)     std::cout<<"\t"<< labels[dim+4]; \
 203         std::cout<<std::endl; \
 204         for (unsigned int phaseIndex=0; phaseIndex<m_NumberOfFields; phaseIndex++) \
 205           { \
 206             std::cout<<phaseIndex; \
 207             std::cout<<"\t"<<oMeanList[phaseIndex]<<"\t"<<oStdList[phaseIndex]; \
 208             std::cout<<"\t"<<oMeanXYZList[phaseIndex][0]<<"\t"<<oStdXYZList[phaseIndex][0]; \
 209             for(unsigned int dim=1;dim<Dimension;dim++) \
 210               std::cout<<"\t"<<oMeanXYZList[phaseIndex][dim]<<"\t"<<oStdXYZList[phaseIndex][dim]; \
 211             std::cout<<"\t"<<oMaxList[phaseIndex]; \
 212             std::cout<<"\t"<<oMaxXYZList[phaseIndex][0]; \
 213             for(unsigned int dim=1;dim<Dimension;dim++) \
 214               std::cout<<"\t"<<oMaxXYZList[phaseIndex][dim]; \
 215             std::cout<<std::endl; \
 216           } \
 217 \
 218         /* // General */\
 219         std::cout<<std::endl; \
 220         std::cout<<"@"; \
 221         std::cout<<"\t"<<oMean<<"\t"<<oStd; \
 222         std::cout<<"\t"<<oMeanXYZ[0]<<"\t"<<oStdXYZ[0]; \
 223         for(unsigned int dim=1;dim<Dimension;dim++) \
 224           std::cout<<"\t"<<oMeanXYZ[dim]<<"\t"<<oStdXYZ[dim]; \
 225         std::cout<<"\t"<<oMax; \
 226         std::cout<<"\t"<<oMaxXYZ[0]; \
 227         for(unsigned int dim=1;dim<Dimension;dim++) \
 228           std::cout<<"\t"<<oMaxXYZ[dim]; \
 229         std::cout<<std::endl; \
 230       } \
 231 \
 232     /* // Output to file */ \
 233     if( m_ArgsInfo.general_given) \
 234       { \
 235         /* // Add to the file */\
 236         std::ofstream general(m_ArgsInfo.general_arg); \
 237 \
 238 \
 239         /* // Numbers of DVF */\
 240         for (unsigned int phaseIndex=0; phaseIndex<m_NumberOfFields; phaseIndex++) \
 241           general<<phaseIndex<<"\t"<< filenames[phaseIndex]<<std::endl; \
 242 \
 243         general<<std::endl; \
 244         general<<"=================================================="<<std::endl; \
 245         general<<"||       Original distance between points       ||"<<std::endl; \
 246         general<<"=================================================="<<std::endl; \
 247         general<<std::endl; \
 248 \
 249         general<<std::setprecision(3); \
 250 \
 251         /* // Labels of the columns */\
 252         general<<"#DVF\tMean\tSD\t"<<labels[0]; \
 253         for(unsigned int dim=1;dim<Dimension;dim++)     general<<"\t"<< labels[dim]; \
 254         general<<"\tMax \t"<<labels[4]; \
 255         for(unsigned int dim=1;dim<Dimension;dim++)     general<<"\t"<< labels[dim+4]; \
 256         general<<std::endl; \
 257         for (unsigned int phaseIndex=0; phaseIndex<m_NumberOfFields; phaseIndex++) \
 258           { \
 259             general<<phaseIndex; \
 260             general<<"\t"<<oMeanList[phaseIndex]<<"\t"<<oStdList[phaseIndex]; \
 261             general<<"\t"<<oMeanXYZList[phaseIndex][0]<<"\t"<<oStdXYZList[phaseIndex][0]; \
 262             for(unsigned int dim=1;dim<Dimension;dim++) \
 263               general<<"\t"<<oMeanXYZList[phaseIndex][dim]<<"\t"<<oStdXYZList[phaseIndex][dim]; \
 264             general<<"\t"<<oMaxList[phaseIndex]; \
 265             general<<"\t"<<oMaxXYZList[phaseIndex][0]; \
 266             for(unsigned int dim=1;dim<Dimension;dim++) \
 267               general<<"\t"<<oMaxXYZList[phaseIndex][dim]; \
 268             general<<std::endl; \
 269           } \
 270 \
 271         /* // General */\
 272         general<<std::endl; \
 273         general<<"@"; \
 274         general<<"\t"<<oMean<<"\t"<<oStd; \
 275         general<<"\t"<<oMeanXYZ[0]<<"\t"<<oStdXYZ[0]; \
 276         for(unsigned int dim=1;dim<Dimension;dim++) \
 277           general<<"\t"<<oMeanXYZ[dim]<<"\t"<<oStdXYZ[dim]; \
 278         general<<"\t"<<oMax; \
 279         general<<"\t"<<oMaxXYZ[0]; \
 280         for(unsigned int dim=1;dim<Dimension;dim++) \
 281           general<<"\t"<<oMaxXYZ[dim]; \
 282         general<<std::endl; \
 283         general.close(); \
 284       }
 285
 286
 287 #define StatisticsDisplacementsMacro  \
 288     /* //----------------------------- */\
 289     /* // Statistics displacements */\
 290     /* //----------------------------- */\
 291 \
 292     /* // Statistics  (magnitude) */\
 293     MeasureListType dmeanList, dstdList, dmaxList; \
 294     ValueType dmean, dstd, dmax; \
 295     statisticsFilter->GetStatistics(displacementLists, dmean, dstd, dmax, dmeanList, dstdList, dmaxList); \
 296 \
 297     /* // Statistics  (per component) */\
 298     VectorListType dmeanXYZList, dstdXYZList,dmaxXYZList; \
 299     VectorType dmeanXYZ, dstdXYZ, dmaxXYZ; \
 300     statisticsFilter->GetStatistics(displacementLists, dmeanXYZ, dstdXYZ, dmaxXYZ,dmeanXYZList, dstdXYZList, dmaxXYZList); \
 301 \
 302 \
 303     /* //----------------------------- */\
 304     /* // Statistics TRE */\
 305     /* //----------------------------- */\
 306 \
 307     /* // Statistics  (magnitude) */\
 308     MeasureListType meanList, stdList, maxList;  \
 309     ValueType mean, std, max; \
 310     statisticsFilter->GetStatistics(treLists, mean, std, max, meanList, stdList, maxList); \
 311 \
 312     /* // Statistics  (per component) */\
 313     VectorListType meanXYZList, stdXYZList,maxXYZList; \
 314     VectorType meanXYZ, stdXYZ, maxXYZ; \
 315     statisticsFilter->GetStatistics(treLists, meanXYZ, stdXYZ, maxXYZ, meanXYZList, stdXYZList, maxXYZList); \
 316 \
 317 \
 318     /* // Output to screen */\
 319     if(m_Verbose) \
 320       { \
 321 \
 322         std::cout<<std::endl; \
 323         std::cout<<"=================================================="<<std::endl; \
 324         std::cout<<"||       Residual distance between points       ||"<<std::endl; \
 325         std::cout<<"=================================================="<<std::endl; \
 326         std::cout<<std::endl; \
 327 \
 328         std::cout<<std::setprecision(3); \
 329 \
 330         /* // Labels of the columns */\
 331         std::cout<<"#DVF\tMean\tSD\t"<<labels[0]; \
 332         for(unsigned int dim=1;dim<Dimension;dim++)     std::cout<<"\t"<< labels[dim]; \
 333         std::cout<<"\tMax \t"<<labels[4]; \
 334         for(unsigned int dim=1;dim<Dimension;dim++)     std::cout<<"\t"<< labels[dim+4]; \
 335         std::cout<<std::endl; \
 336         for (unsigned int phaseIndex=0; phaseIndex<m_NumberOfFields; phaseIndex++) \
 337           { \
 338             std::cout<<phaseIndex; \
 339             std::cout<<"\t"<<meanList[phaseIndex]<<"\t"<<stdList[phaseIndex]; \
 340             std::cout<<"\t"<<meanXYZList[phaseIndex][0]<<"\t"<<stdXYZList[phaseIndex][0]; \
 341             for(unsigned int dim=1;dim<Dimension;dim++) \
 342               std::cout<<"\t"<<meanXYZList[phaseIndex][dim]<<"\t"<<stdXYZList[phaseIndex][dim]; \
 343             std::cout<<"\t"<<maxList[phaseIndex]; \
 344             std::cout<<"\t"<<maxXYZList[phaseIndex][0]; \
 345             for(unsigned int dim=1;dim<Dimension;dim++) \
 346               std::cout<<"\t"<<maxXYZList[phaseIndex][dim]; \
 347             std::cout<<std::endl; \
 348           } \
 349 \
 350         /* // General */\
 351         std::cout<<std::endl; \
 352         std::cout<<"@"; \
 353         std::cout<<"\t"<<mean<<"\t"<<std; \
 354         std::cout<<"\t"<<meanXYZ[0]<<"\t"<<stdXYZ[0]; \
 355         for(unsigned int dim=1;dim<Dimension;dim++) \
 356           std::cout<<"\t"<<meanXYZ[dim]<<"\t"<<stdXYZ[dim]; \
 357         std::cout<<"\t"<<max; \
 358         std::cout<<"\t"<<maxXYZ[0]; \
 359         for(unsigned int dim=1;dim<Dimension;dim++) \
 360           std::cout<<"\t"<<maxXYZ[dim]; \
 361         std::cout<<std::endl; \
 362       } \
 363 \
 364     /* // Output to file */\
 365     if( m_ArgsInfo.general_given) \
 366       { \
 367         /*  // Add to the file */\
 368         std::ofstream general(m_ArgsInfo.general_arg, ios_base::app); \
 369 \
 370         general<<std::endl; \
 371         general<<"=================================================="<<std::endl; \
 372         general<<"||       Residual distance between points       ||"<<std::endl; \
 373         general<<"=================================================="<<std::endl; \
 374         general<<std::endl; \
 375 \
 376         general<<std::setprecision(3); \
 377 \
 378         /* // Labels of the columns */\
 379         general<<"#DVF\tMean\tSD\t"<<labels[0]; \
 380         for(unsigned int dim=1;dim<Dimension;dim++)     general<<"\t"<< labels[dim]; \
 381         general<<"\tMax \t"<<labels[4]; \
 382         for(unsigned int dim=1;dim<Dimension;dim++)     general<<"\t"<< labels[dim+4]; \
 383         general<<std::endl; \
 384         for (unsigned int phaseIndex=0; phaseIndex<m_NumberOfFields; phaseIndex++) \
 385           { \
 386             general<<phaseIndex; \
 387             general<<"\t"<<meanList[phaseIndex]<<"\t"<<stdList[phaseIndex]; \
 388             general<<"\t"<<meanXYZList[phaseIndex][0]<<"\t"<<stdXYZList[phaseIndex][0]; \
 389             for(unsigned int dim=1;dim<Dimension;dim++) \
 390               general<<"\t"<<meanXYZList[phaseIndex][dim]<<"\t"<<stdXYZList[phaseIndex][dim]; \
 391             general<<"\t"<<maxList[phaseIndex]; \
 392             general<<"\t"<<maxXYZList[phaseIndex][0]; \
 393             for(unsigned int dim=1;dim<Dimension;dim++) \
 394               general<<"\t"<<maxXYZList[phaseIndex][dim]; \
 395             general<<std::endl; \
 396           } \
 397 \
 398         /* // General */\
 399         general<<std::endl; \
 400         general<<"@"; \
 401         general<<"\t"<<mean<<"\t"<<std; \
 402         general<<"\t"<<meanXYZ[0]<<"\t"<<stdXYZ[0]; \
 403         for(unsigned int dim=1;dim<Dimension;dim++) \
 404           general<<"\t"<<meanXYZ[dim]<<"\t"<<stdXYZ[dim]; \
 405         general<<"\t"<<max; \
 406         general<<"\t"<<maxXYZ[0]; \
 407         for(unsigned int dim=1;dim<Dimension;dim++) \
 408           general<<"\t"<<maxXYZ[dim]; \
 409         general<<std::endl; \
 410         general.close(); \
 411       } \
 412 \
 413     /* // Output original points */\
 414     if( m_ArgsInfo.original_given) \
 415       { \
 416 \
 417         std::vector<std::string> filenames; \
 418         for (unsigned int i=0;i<m_NumberOfFields;i++) \
 419           { \
 420             std::ostringstream number_ostr; \
 421             number_ostr << i; \
 422             std::string number_str =  number_ostr.str(); \
 423             filenames.push_back(m_ArgsInfo.original_arg+number_str); \
 424           } \
 425         originalDistanceLists.Write(filenames, m_Verbose ); \
 426       } \
 427 \
 428 \
 429     /*  // Output original magnitude points */\
 430     if( m_ArgsInfo.originalMag_given) \
 431       { \
 432         clitk::Lists<itk::FixedArray<ValueType,1> > originalDistanceListsMag=originalDistanceLists.Norm(); \
 433         std::vector<std::string> filenames; \
 434         for (unsigned int i=0;i<m_NumberOfFields;i++) \
 435           { \
 436             std::ostringstream number_ostr; \
 437             number_ostr << i; \
 438             std::string number_str =  number_ostr.str(); \
 439             filenames.push_back(m_ArgsInfo.originalMag_arg+number_str); \
 440           } \
 441         originalDistanceListsMag.Write(filenames, m_Verbose ); \
 442       } \
 443 \
 444     /*  // Output displacement */\
 445     if( m_ArgsInfo.displacement_given) \
 446       { \
 447 \
 448         std::vector<std::string> filenames; \
 449         for (unsigned int i=0;i<m_NumberOfFields;i++) \
 450           { \
 451             std::ostringstream number_ostr; \
 452             number_ostr << i; \
 453             std::string number_str =  number_ostr.str(); \
 454             filenames.push_back(m_ArgsInfo.displacement_arg+number_str); \
 455           } \
 456         displacementLists.Write(filenames, m_Verbose ); \
 457       } \
 458 \
 459 \
 460     /*  // Output displacement magnitude  */\
 461     if( m_ArgsInfo.displacementMag_given) \
 462       { \
 463         clitk::Lists<itk::FixedArray<ValueType,1> > displacementListsMag=displacementLists.Norm(); \
 464         std::vector<std::string> filenames; \
 465         for (unsigned int i=0;i<m_NumberOfFields;i++) \
 466           { \
 467             std::ostringstream number_ostr; \
 468             number_ostr << i; \
 469             std::string number_str =  number_ostr.str(); \
 470             filenames.push_back(m_ArgsInfo.displacementMag_arg+number_str); \
 471           } \
 472         displacementListsMag.Write(filenames, m_Verbose ); \
 473       } \
 474 \
 475     /*  // Output warped points */\
 476     if( m_ArgsInfo.warp_given) \
 477       { \
 478 \
 479         std::vector<std::string> filenames; \
 480         for (unsigned int i=0;i<m_ArgsInfo.warp_given;i++) \
 481           { \
 482             std::ostringstream number_ostr; \
 483             number_ostr << i; \
 484             std::string number_str =  number_ostr.str(); \
 485             filenames.push_back(m_ArgsInfo.warp_arg+number_str); \
 486           } \
 487         warpedPointLists.Write(filenames, m_Verbose ); \
 488       } \
 489 \
 490     /*  // Output tre  */\
 491     if( m_ArgsInfo.tre_given) \
 492       { \
 493 \
 494         std::vector<std::string> filenames; \
 495         for (unsigned int i=0;i<m_NumberOfFields;i++) \
 496           { \
 497             std::ostringstream number_ostr; \
 498             number_ostr << i; \
 499             std::string number_str =  number_ostr.str(); \
 500             filenames.push_back(m_ArgsInfo.tre_arg+number_str); \
 501           } \
 502         treLists.Write(filenames, m_Verbose ); \
 503       } \
 504 \
 505     /*  // Output tre mag */\
 506     if( m_ArgsInfo.treMag_given) \
 507       { \
 508         clitk::Lists<itk::FixedArray<ValueType,1> > treMagLists=treLists.Norm(); \
 509 \
 510         std::vector<std::string> filenames; \
 511         for (unsigned int i=0;i<m_NumberOfFields;i++) \
 512           { \
 513             std::ostringstream number_ostr; \
 514             number_ostr << i; \
 515             std::string number_str =  number_ostr.str(); \
 516             filenames.push_back(m_ArgsInfo.treMag_arg+number_str); \
 517           } \
 518         treMagLists.Write(filenames, m_Verbose ); \
 519       }
 520
 521
 522
 523
 524 namespace clitk
 525 {
 526
 527   //-----------------------------
 528   // Read DVF
 529   //-----------------------------
 530   template<unsigned int Dimension, unsigned int Components >
 531   void
 532   CalculateTREGenericFilter::ReadVectorFields(void)
 533   {
 534     // if using dvfs
 535     typedef itk::Vector<ValueType, Components> VectorType;
 536     typedef itk::Image<VectorType, Dimension> DeformationFieldType;
 537
 538     typedef itk::ImageFileReader<DeformationFieldType> DvfReaderType;
 539     std::vector<typename DeformationFieldType::Pointer> dvfs;
 540
 541     for (unsigned int i=0; i< m_ArgsInfo.vf_given; i++)
 542       {
 543         typename DvfReaderType::Pointer reader = DvfReaderType::New();
 544         reader->SetFileName( m_ArgsInfo.vf_arg[i]);
 545         if (m_Verbose) std::cout<<"Reading vector field "<< i <<"..."<<std::endl;
 546         reader->Update();
 547         dvfs.push_back( reader->GetOutput() );
 548       }
 549
 550     ProcessVectorFields<Dimension, Components>(dvfs, m_ArgsInfo.vf_arg);
 551
 552   }
 553
 554   //-----------------------------
 555   // Process DVF
 556   //-----------------------------
 557   template<unsigned int Dimension, unsigned int Components >
 558   void
 559   CalculateTREGenericFilter::ProcessVectorFields(std::vector<typename itk::Image<itk::Vector<CalculateTREGenericFilter::ValueType, Components>, Dimension>::Pointer > dvfs, char** filenames )
 560   {
 561     std::vector<std::string> new_filenames;
 562     for (unsigned int i=0;i<m_ArgsInfo.vf_given;i++)
 563       new_filenames.push_back(filenames[i]);
 564     UpdateDVFWithDim<Dimension>(dvfs, new_filenames);
 565   }
 566
 567
 568   //-----------------------------
 569   // Read Coefficient Images
 570   //-----------------------------
 571   template<unsigned int Dimension, unsigned int Components >
 572   void
 573   CalculateTREGenericFilter::ReadCoefficientImages(void)
 574   {
 575     // if using coefficient images
 576     typedef itk::Vector<ValueType, Components> VectorType;
 577     typedef itk::Image<VectorType, Dimension> CoefficientImageType;
 578
 579     typedef itk::ImageFileReader<CoefficientImageType> CoeffReaderType;
 580     std::vector<typename CoefficientImageType::Pointer> coeffs;
 581
 582     for (unsigned int i=0; i< m_ArgsInfo.coeff_given; i++)
 583       {
 584         typename CoeffReaderType::Pointer reader = CoeffReaderType::New();
 585         reader->SetFileName( m_ArgsInfo.coeff_arg[i]);
 586         if (m_Verbose) std::cout<<"Reading coefficient image "<< i <<"..."<<std::endl;
 587         reader->Update();
 588         coeffs.push_back( reader->GetOutput() );
 589       }
 590
 591     ProcessCoefficientImages<Dimension, Components>(coeffs, m_ArgsInfo.coeff_arg);
 592   }
 593
 594   //-----------------------------
 595   // Process Coefficient Images
 596   //-----------------------------
 597   template<unsigned int Dimension, unsigned int Components >
 598   void
 599   CalculateTREGenericFilter::ProcessCoefficientImages(std::vector<typename itk::Image<itk::Vector<CalculateTREGenericFilter::ValueType, Components>, Dimension>::Pointer > coeffs, char** filenames )
 600   {
 601     std::vector<std::string> new_filenames;
 602     for (unsigned int i=0;i<m_ArgsInfo.coeff_given;i++)
 603       new_filenames.push_back(filenames[i]);
 604     UpdateCoeffsWithDim<Dimension>(coeffs, new_filenames);
 605   }
 606
 607
 608   //-------------------------------------------------------------------
 609   // Update DVF with the number of dimensions
 610   //-------------------------------------------------------------------
 611   template<unsigned int Dimension>
 612   void
 613   CalculateTREGenericFilter::UpdateDVFWithDim(std::vector<typename itk::Image<itk::Vector<CalculateTREGenericFilter::ValueType, Dimension>, Dimension>::Pointer > dvfs,  std::vector<std::string> filenames)
 614   {
 615     if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D..."<<std::endl;
 616
 617     TypedefsMacro;
 618     BuildPointListsMacro(filenames);
 619     OriginalDistancesMacro;
 620
 621     //=====================================================================================
 622     // Distance between points after warp
 623     //=====================================================================================
 624
 625     //-----------------------------
 626     // Interpolator
 627     //-----------------------------
 628     typedef clitk::GenericVectorInterpolator<args_info_clitkCalculateTRE,DeformationFieldType, ValueType> GenericVectorInterpolatorType;
 629     typename GenericVectorInterpolatorType::Pointer genericInterpolator=GenericVectorInterpolatorType::New();
 630     genericInterpolator->SetArgsInfo(m_ArgsInfo);
 631     typedef itk::VectorInterpolateImageFunction<DeformationFieldType, ValueType> InterpolatorType;
 632     typename InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer();
 633
 634
 635     PointListsType warpedPointLists(m_NumberOfFields);
 636     VectorListsType treLists(m_NumberOfFields), displacementLists(m_NumberOfFields);
 637
 638     //-----------------------------
 639     // Calculate residual distances
 640     //-----------------------------
 641     VectorType displacement, tre;
 642     PointType  warpedPoint;
 643
 644     for (unsigned int phaseIndex=0; phaseIndex<m_NumberOfFields; phaseIndex++)
 645       {
 646         interpolator->SetInputImage( dvfs[phaseIndex] );
 647         for (unsigned int pointIndex=0; pointIndex<m_NumberOfPoints; pointIndex++)
 648           {
 649             // Reference
 650             referencePoint=referencePointList[pointIndex];
 651
 652             // Inside?
 653             if ( interpolator->IsInsideBuffer(referencePoint) )
 654               displacement=interpolator->Evaluate(referencePoint);
 655             else
 656               displacement.Fill(0.0);
 657
 658             // Warp
 659             warpedPoint=referencePoint+displacement;
 660             displacementLists[phaseIndex].push_back(displacement);
 661             warpedPointLists[phaseIndex].push_back(warpedPoint);
 662             tre=pointLists[phaseIndex][pointIndex]-warpedPoint;
 663             treLists[phaseIndex].push_back(tre);
 664           }
 665       }
 666
 667     StatisticsDisplacementsMacro;
 668   }
 669
 670   //-------------------------------------------------------------------
 671   // Update coefficients with the number of dimensions
 672   //-------------------------------------------------------------------
 673   template<unsigned int Dimension>
 674   void
 675   CalculateTREGenericFilter::UpdateCoeffsWithDim(std::vector<typename itk::Image<itk::Vector<CalculateTREGenericFilter::ValueType, Dimension>, Dimension>::Pointer >  coeffs,  std::vector<std::string> filenames)
 676   {
 677     // call macros
 678     TypedefsMacro;
 679     BuildPointListsMacro(filenames);
 680     OriginalDistancesMacro;
 681
 682     //=====================================================================================
 683     // Distance between points after warp
 684     //=====================================================================================
 685
 686     //-----------------------------
 687     // Transform to get the transformed points from the coeff images
 688     //-----------------------------
 689     typedef clitk::BSplineDeformableTransform<ValueType, Dimension, Dimension> TransformType;
 690     typename TransformType::Pointer transform = TransformType::New();
 691
 692     PointListsType warpedPointLists(m_NumberOfFields);
 693     VectorListsType treLists(m_NumberOfFields), displacementLists(m_NumberOfFields);
 694
 695     //-----------------------------
 696     // Calculate residual distances
 697     //-----------------------------
 698     VectorType displacement, tre;
 699     PointType  warpedPoint;
 700
 701     for (unsigned int phaseIndex=0; phaseIndex<m_NumberOfFields; phaseIndex++)
 702       {
 703         // set the initial parameters (b-spline coefficients from cmd-line
 704         transform->SetCoefficientImage(coeffs[phaseIndex]);
 705
 706         // transform the input points
 707         for (unsigned int pointIndex=0; pointIndex<m_NumberOfPoints; pointIndex++)
 708           {
 709             // Reference
 710             referencePoint=referencePointList[pointIndex];
 711
 712             // Transform
 713             warpedPoint=transform->TransformPoint(referencePoint);
 714
 715             // Warp
 716             displacement=warpedPoint-referencePoint;
 717             displacementLists[phaseIndex].push_back(displacement);
 718             warpedPointLists[phaseIndex].push_back(warpedPoint);
 719             tre=pointLists[phaseIndex][pointIndex]-warpedPoint;
 720             treLists[phaseIndex].push_back(tre);
 721           }
 722       }
 723
 724     StatisticsDisplacementsMacro;
 725   }
 726
 727 }//end clitk
 728 #endif //#define clitkCalculateTREGenericFilter_txx