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://oncora1.lyon.fnclcc.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  * @file   clitkCalculateTREGenericFilter.txx
  23  * @author
  24  * @date
  25  *
  26  * @brief
  27  *
  28  ===================================================*/
  29
  30
  31 namespace clitk
  32 {
  33
  34   //-----------------------------
  35   // Read DVF
  36   //-----------------------------
  37   template<unsigned int Dimension, unsigned int Components >
  38   void
  39   CalculateTREGenericFilter::ReadVectorFields(void)
  40   {
  41
  42     typedef itk::Vector<float, Components> VectorType;
  43     typedef itk::Image<VectorType, Dimension> DeformationFieldType;
  44
  45     typedef itk::ImageFileReader<DeformationFieldType> InputReaderType;
  46     std::vector<typename DeformationFieldType::Pointer> dvfs;
  47     for (unsigned int i=0; i< m_ArgsInfo.vf_given; i++)
  48       {
  49         typename InputReaderType::Pointer reader = InputReaderType::New();
  50         reader->SetFileName( m_ArgsInfo.vf_arg[i]);
  51         if (m_Verbose) std::cout<<"Reading vector field "<< i <<"..."<<std::endl;
  52         reader->Update();
  53         dvfs.push_back( reader->GetOutput() );
  54       }
  55
  56     ProcessVectorFields<Dimension, Components>(dvfs, m_ArgsInfo.vf_arg);
  57   }
  58
  59   //-----------------------------
  60   // Process DVF
  61   //-----------------------------
  62   template<unsigned int Dimension, unsigned int Components >
  63   void
  64   CalculateTREGenericFilter::ProcessVectorFields(std::vector<typename itk::Image<itk::Vector<float, Components>, Dimension>::Pointer > dvfs, char** filenames )
  65   {
  66
  67     std::vector<std::string> new_filenames;
  68     for (unsigned int i=0;i<m_ArgsInfo.vf_given;i++)
  69       new_filenames.push_back(filenames[i]);
  70     UpdateWithDim<Dimension>(dvfs, new_filenames);
  71   }
  72
  73
  74   //-------------------------------------------------------------------
  75   // Update with the number of dimensions
  76   //-------------------------------------------------------------------
  77   template<unsigned int Dimension>
  78   void
  79   CalculateTREGenericFilter::UpdateWithDim(std::vector<typename itk::Image<itk::Vector<float, Dimension>, Dimension>::Pointer > dvfs,  std::vector<std::string> filenames)
  80   {
  81     if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D..."<<std::endl;
  82
  83     //-----------------------------
  84     // Typedefs
  85     //-----------------------------
  86     typedef double ValueType;
  87     typedef std::vector<ValueType> MeasureListType;
  88
  89     typedef itk::Point<double, Dimension> PointType;
  90     typedef clitk::List<PointType> PointListType;
  91     typedef clitk::Lists<PointType> PointListsType;
  92
  93     typedef itk::Vector<double, Dimension> VectorType;
  94     typedef clitk::List<VectorType> VectorListType;
  95     typedef clitk::Lists<VectorType> VectorListsType;
  96
  97     typedef itk::Vector<float, Dimension> DeformationVectorType;
  98     typedef itk::Image<DeformationVectorType, Dimension> DeformationFieldType;
  99
 100     //-----------------------------
 101     // Number of inputs
 102     //-----------------------------
 103     unsigned int numberOfFields=filenames.size();
 104     unsigned int numberOfLists=m_ArgsInfo.input_given;
 105     if ( (numberOfLists!=numberOfFields)  &&  (numberOfLists!=1) )
 106       {
 107         std::cerr<<"Error: Number of lists (="<<numberOfLists<<") different from number of DVF's (="<<numberOfFields<<") !"<<std::endl;
 108         return;
 109       }
 110     else if (numberOfLists==1)
 111       {
 112         if (m_Verbose) std::cout<<numberOfFields<<" DVFs and 1 point list given..."<<std::endl;
 113       }
 114     else
 115       {
 116         if (m_Verbose) std::cout<<numberOfLists<<" point lists and DVFs given..."<<std::endl;
 117       }
 118
 119
 120     //-----------------------------
 121     // Input point lists
 122     //-----------------------------
 123     PointListsType pointLists;
 124     unsigned int numberOfPoints=0;
 125     for (unsigned int i=0; i<numberOfFields; i++)
 126       {
 127         // Read the lists
 128         if (numberOfLists==1)
 129           pointLists.push_back(PointListType(m_ArgsInfo.input_arg[0], m_Verbose) );
 130         else
 131           pointLists.push_back(PointListType(m_ArgsInfo.input_arg[i], m_Verbose) );
 132
 133
 134         // Verify the number of points
 135         if (i==0) numberOfPoints=pointLists[i].size();
 136         else
 137           {
 138             if  (numberOfPoints!=pointLists[i].size())
 139               {
 140                 std::cerr<<"Size of first list ("<<numberOfPoints
 141                          <<") is different from size of list "<<i
 142                          <<" ("<<pointLists[i].size()<<")..."<<std::endl;
 143                 return;
 144               }
 145           }
 146       }
 147
 148
 149     PointListType referencePointList;
 150     if (m_Verbose) std::cout<<"Reference point list:"<<std::endl;
 151     referencePointList=PointListType(m_ArgsInfo.ref_arg, m_Verbose);
 152     if  (numberOfPoints!=referencePointList.size())
 153       {
 154         std::cerr<<"Size of the first list ("<<numberOfPoints
 155                  <<") is different from size of the reference list ("
 156                  << referencePointList.size() <<")..."<<std::endl;
 157         return;
 158       }
 159
 160
 161     //-----------------------------
 162     // Interpolator
 163     //-----------------------------
 164     typedef clitk::GenericVectorInterpolator<args_info_clitkCalculateTRE,DeformationFieldType, double> GenericVectorInterpolatorType;
 165     typename GenericVectorInterpolatorType::Pointer genericInterpolator=GenericVectorInterpolatorType::New();
 166     genericInterpolator->SetArgsInfo(m_ArgsInfo);
 167     typedef itk::VectorInterpolateImageFunction<DeformationFieldType, double> InterpolatorType;
 168     typename InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer();
 169
 170
 171     //=====================================================================================
 172     // Original distance between points
 173     //=====================================================================================
 174     VectorListsType originalDistanceLists(numberOfFields);
 175
 176     //-----------------------------
 177     // Calculate original distances
 178     //-----------------------------
 179     PointType referencePoint;
 180     VectorType distance;
 181     for (unsigned int phaseIndex=0; phaseIndex<numberOfFields; phaseIndex++)
 182       {
 183         for (unsigned int pointIndex=0; pointIndex<numberOfPoints; pointIndex++)
 184           {
 185             referencePoint=referencePointList[pointIndex];
 186             distance=pointLists[phaseIndex][pointIndex]-referencePoint;
 187             originalDistanceLists[phaseIndex].push_back(distance);
 188           }
 189       }
 190
 191
 192     //-----------------------------
 193     // Statistics Original
 194     //-----------------------------
 195     typedef DeformationListStatisticsFilter<VectorType> StatisticsFilterType;
 196     typename StatisticsFilterType::Pointer statisticsFilter=StatisticsFilterType::New();
 197
 198     // Statistics  (magnitude)
 199     MeasureListType oMeanList, oStdList, oMaxList;
 200     ValueType oMean, oStd, oMax;
 201     statisticsFilter->GetStatistics(originalDistanceLists, oMean, oStd, oMax, oMeanList, oStdList, oMaxList);
 202
 203     // Statistics  (per component)
 204     VectorListType oMeanXYZList, oStdXYZList,oMaxXYZList;
 205     VectorType oMeanXYZ, oStdXYZ, oMaxXYZ;
 206     statisticsFilter->GetStatistics(originalDistanceLists, oMeanXYZ, oStdXYZ, oMaxXYZ, oMeanXYZList, oStdXYZList, oMaxXYZList);
 207
 208
 209     //-----------------------------
 210     // Output
 211     //-----------------------------
 212     std::vector<std::string> labels;
 213     labels.push_back("MeanX\tSDX");
 214     labels.push_back("MeanY\tSDY");
 215     labels.push_back("MeanZ\tSDZ");
 216     labels.push_back("MeanT\tSDT");
 217     labels.push_back("MaxX");
 218     labels.push_back("MaxY");
 219     labels.push_back("MaxZ");
 220     labels.push_back("MaxT");
 221
 222
 223     // Output to screen
 224     if(m_Verbose)
 225       {
 226
 227         // Numbers of DVF
 228         std::cout<<"# Number\tDVF"<<std::endl;
 229         for (unsigned int phaseIndex=0; phaseIndex<numberOfFields; phaseIndex++)
 230           std::cout<<phaseIndex<<"\t"<<filenames[phaseIndex]<<std::endl;
 231
 232
 233         std::cout<<std::endl;
 234         std::cout<<"=================================================="<<std::endl;
 235         std::cout<<"||       Original distance between points       ||"<<std::endl;
 236         std::cout<<"=================================================="<<std::endl;
 237         std::cout<<std::endl;
 238
 239         std::cout<<std::setprecision(3);
 240
 241
 242
 243         // Labels of the columns
 244         std::cout<<"#DVF\tMean\tSD\t"<<labels[0];
 245         for(unsigned int dim=1;dim<Dimension;dim++)     std::cout<<"\t"<< labels[dim];
 246         std::cout<<"\tMax \t"<<labels[4];
 247         for(unsigned int dim=1;dim<Dimension;dim++)     std::cout<<"\t"<< labels[dim+4];
 248         std::cout<<std::endl;
 249         for (unsigned int phaseIndex=0; phaseIndex<numberOfFields; phaseIndex++)
 250           {
 251             std::cout<<phaseIndex;
 252             std::cout<<"\t"<<oMeanList[phaseIndex]<<"\t"<<oStdList[phaseIndex];
 253             std::cout<<"\t"<<oMeanXYZList[phaseIndex][0]<<"\t"<<oStdXYZList[phaseIndex][0];
 254             for(unsigned int dim=1;dim<Dimension;dim++)
 255               std::cout<<"\t"<<oMeanXYZList[phaseIndex][dim]<<"\t"<<oStdXYZList[phaseIndex][dim];
 256             std::cout<<"\t"<<oMaxList[phaseIndex];
 257             std::cout<<"\t"<<oMaxXYZList[phaseIndex][0];
 258             for(unsigned int dim=1;dim<Dimension;dim++)
 259               std::cout<<"\t"<<oMaxXYZList[phaseIndex][dim];
 260             std::cout<<std::endl;
 261           }
 262
 263         // General
 264         std::cout<<std::endl;
 265         std::cout<<"@";
 266         std::cout<<"\t"<<oMean<<"\t"<<oStd;
 267         std::cout<<"\t"<<oMeanXYZ[0]<<"\t"<<oStdXYZ[0];
 268         for(unsigned int dim=1;dim<Dimension;dim++)
 269           std::cout<<"\t"<<oMeanXYZ[dim]<<"\t"<<oStdXYZ[dim];
 270         std::cout<<"\t"<<oMax;
 271         std::cout<<"\t"<<oMaxXYZ[0];
 272         for(unsigned int dim=1;dim<Dimension;dim++)
 273           std::cout<<"\t"<<oMaxXYZ[dim];
 274         std::cout<<std::endl;
 275       }
 276
 277     // Output to file
 278     if( m_ArgsInfo.general_given)
 279       {
 280         // Add to the file
 281         std::ofstream general(m_ArgsInfo.general_arg);
 282
 283
 284         // Numbers of DVF
 285         for (unsigned int phaseIndex=0; phaseIndex<numberOfFields; phaseIndex++)
 286           general<<phaseIndex<<"\t"<< filenames[phaseIndex]<<std::endl;
 287
 288         general<<std::endl;
 289         general<<"=================================================="<<std::endl;
 290         general<<"||       Original distance between points       ||"<<std::endl;
 291         general<<"=================================================="<<std::endl;
 292         general<<std::endl;
 293
 294         general<<std::setprecision(3);
 295
 296         // Labels of the columns
 297         general<<"#DVF\tMean\tSD\t"<<labels[0];
 298         for(unsigned int dim=1;dim<Dimension;dim++)     general<<"\t"<< labels[dim];
 299         general<<"\tMax \t"<<labels[4];
 300         for(unsigned int dim=1;dim<Dimension;dim++)     general<<"\t"<< labels[dim+4];
 301         general<<std::endl;
 302         for (unsigned int phaseIndex=0; phaseIndex<numberOfFields; phaseIndex++)
 303           {
 304             general<<phaseIndex;
 305             general<<"\t"<<oMeanList[phaseIndex]<<"\t"<<oStdList[phaseIndex];
 306             general<<"\t"<<oMeanXYZList[phaseIndex][0]<<"\t"<<oStdXYZList[phaseIndex][0];
 307             for(unsigned int dim=1;dim<Dimension;dim++)
 308               general<<"\t"<<oMeanXYZList[phaseIndex][dim]<<"\t"<<oStdXYZList[phaseIndex][dim];
 309             general<<"\t"<<oMaxList[phaseIndex];
 310             general<<"\t"<<oMaxXYZList[phaseIndex][0];
 311             for(unsigned int dim=1;dim<Dimension;dim++)
 312               general<<"\t"<<oMaxXYZList[phaseIndex][dim];
 313             general<<std::endl;
 314           }
 315
 316         // General
 317         general<<std::endl;
 318         general<<"@";
 319         general<<"\t"<<oMean<<"\t"<<oStd;
 320         general<<"\t"<<oMeanXYZ[0]<<"\t"<<oStdXYZ[0];
 321         for(unsigned int dim=1;dim<Dimension;dim++)
 322           general<<"\t"<<oMeanXYZ[dim]<<"\t"<<oStdXYZ[dim];
 323         general<<"\t"<<oMax;
 324         general<<"\t"<<oMaxXYZ[0];
 325         for(unsigned int dim=1;dim<Dimension;dim++)
 326           general<<"\t"<<oMaxXYZ[dim];
 327         general<<std::endl;
 328         general.close();
 329       }
 330
 331
 332     //=====================================================================================
 333     // Distance between points after warp
 334     //=====================================================================================
 335     PointListsType warpedPointLists(numberOfFields);
 336     VectorListsType treLists(numberOfFields), displacementLists(numberOfFields);
 337
 338     //-----------------------------
 339     // Calculate residual distances
 340     //-----------------------------
 341     VectorType displacement, tre;
 342     PointType  warpedPoint;
 343
 344     for (unsigned int phaseIndex=0; phaseIndex<numberOfFields; phaseIndex++)
 345       {
 346         interpolator->SetInputImage( dvfs[phaseIndex]);
 347         for (unsigned int pointIndex=0; pointIndex<numberOfPoints; pointIndex++)
 348           {
 349             // Reference
 350             referencePoint=referencePointList[pointIndex];
 351
 352             // Inside?
 353             if ( interpolator->IsInsideBuffer(referencePoint) )
 354               displacement=interpolator->Evaluate(referencePoint);
 355             else
 356               displacement.Fill(0.0);
 357
 358             // Warp
 359             warpedPoint=referencePoint+displacement;
 360             displacementLists[phaseIndex].push_back(displacement);
 361             warpedPointLists[phaseIndex].push_back(warpedPoint);
 362             tre=pointLists[phaseIndex][pointIndex]-warpedPoint;
 363             treLists[phaseIndex].push_back(tre);
 364           }
 365       }
 366
 367
 368     //-----------------------------
 369     // Statistics displacements
 370     //-----------------------------
 371
 372     // Statistics  (magnitude)
 373     MeasureListType dmeanList, dstdList, dmaxList;
 374     ValueType dmean, dstd, dmax;
 375     statisticsFilter->GetStatistics(displacementLists, dmean, dstd, dmax, dmeanList, dstdList, dmaxList);
 376
 377     // Statistics  (per component)
 378     VectorListType dmeanXYZList, dstdXYZList,dmaxXYZList;
 379     VectorType dmeanXYZ, dstdXYZ, dmaxXYZ;
 380     statisticsFilter->GetStatistics(displacementLists, dmeanXYZ, dstdXYZ, dmaxXYZ,dmeanXYZList, dstdXYZList, dmaxXYZList);
 381
 382
 383     //-----------------------------
 384     // Statistics TRE
 385     //-----------------------------
 386
 387     // Statistics  (magnitude)
 388     MeasureListType meanList, stdList, maxList;
 389     ValueType mean, std, max;
 390     statisticsFilter->GetStatistics(treLists, mean, std, max, meanList, stdList, maxList);
 391
 392     // Statistics  (per component)
 393     VectorListType meanXYZList, stdXYZList,maxXYZList;
 394     VectorType meanXYZ, stdXYZ, maxXYZ;
 395     statisticsFilter->GetStatistics(treLists, meanXYZ, stdXYZ, maxXYZ, meanXYZList, stdXYZList, maxXYZList);
 396
 397
 398     // Output to screen
 399     if(m_Verbose)
 400       {
 401
 402         std::cout<<std::endl;
 403         std::cout<<"=================================================="<<std::endl;
 404         std::cout<<"||       Residual distance between points       ||"<<std::endl;
 405         std::cout<<"=================================================="<<std::endl;
 406         std::cout<<std::endl;
 407
 408         std::cout<<std::setprecision(3);
 409
 410         // Labels of the columns
 411         std::cout<<"#DVF\tMean\tSD\t"<<labels[0];
 412         for(unsigned int dim=1;dim<Dimension;dim++)     std::cout<<"\t"<< labels[dim];
 413         std::cout<<"\tMax \t"<<labels[4];
 414         for(unsigned int dim=1;dim<Dimension;dim++)     std::cout<<"\t"<< labels[dim+4];
 415         std::cout<<std::endl;
 416         for (unsigned int phaseIndex=0; phaseIndex<numberOfFields; phaseIndex++)
 417           {
 418             std::cout<<phaseIndex;
 419             std::cout<<"\t"<<meanList[phaseIndex]<<"\t"<<stdList[phaseIndex];
 420             std::cout<<"\t"<<meanXYZList[phaseIndex][0]<<"\t"<<stdXYZList[phaseIndex][0];
 421             for(unsigned int dim=1;dim<Dimension;dim++)
 422               std::cout<<"\t"<<meanXYZList[phaseIndex][dim]<<"\t"<<stdXYZList[phaseIndex][dim];
 423             std::cout<<"\t"<<maxList[phaseIndex];
 424             std::cout<<"\t"<<maxXYZList[phaseIndex][0];
 425             for(unsigned int dim=1;dim<Dimension;dim++)
 426               std::cout<<"\t"<<maxXYZList[phaseIndex][dim];
 427             std::cout<<std::endl;
 428           }
 429
 430         // General
 431         std::cout<<std::endl;
 432         std::cout<<"@";
 433         std::cout<<"\t"<<mean<<"\t"<<std;
 434         std::cout<<"\t"<<meanXYZ[0]<<"\t"<<stdXYZ[0];
 435         for(unsigned int dim=1;dim<Dimension;dim++)
 436           std::cout<<"\t"<<meanXYZ[dim]<<"\t"<<stdXYZ[dim];
 437         std::cout<<"\t"<<max;
 438         std::cout<<"\t"<<maxXYZ[0];
 439         for(unsigned int dim=1;dim<Dimension;dim++)
 440           std::cout<<"\t"<<maxXYZ[dim];
 441         std::cout<<std::endl;
 442       }
 443
 444     // Output to file
 445     if( m_ArgsInfo.general_given)
 446       {
 447         // Add to the file
 448         std::ofstream general(m_ArgsInfo.general_arg, ios_base::app);
 449
 450         general<<std::endl;
 451         general<<"=================================================="<<std::endl;
 452         general<<"||       Residual distance between points       ||"<<std::endl;
 453         general<<"=================================================="<<std::endl;
 454         general<<std::endl;
 455
 456         general<<std::setprecision(3);
 457
 458         // Labels of the columns
 459         general<<"#DVF\tMean\tSD\t"<<labels[0];
 460         for(unsigned int dim=1;dim<Dimension;dim++)     general<<"\t"<< labels[dim];
 461         general<<"\tMax \t"<<labels[4];
 462         for(unsigned int dim=1;dim<Dimension;dim++)     general<<"\t"<< labels[dim+4];
 463         general<<std::endl;
 464         for (unsigned int phaseIndex=0; phaseIndex<numberOfFields; phaseIndex++)
 465           {
 466             general<<phaseIndex;
 467             general<<"\t"<<meanList[phaseIndex]<<"\t"<<stdList[phaseIndex];
 468             general<<"\t"<<meanXYZList[phaseIndex][0]<<"\t"<<stdXYZList[phaseIndex][0];
 469             for(unsigned int dim=1;dim<Dimension;dim++)
 470               general<<"\t"<<meanXYZList[phaseIndex][dim]<<"\t"<<stdXYZList[phaseIndex][dim];
 471             general<<"\t"<<maxList[phaseIndex];
 472             general<<"\t"<<maxXYZList[phaseIndex][0];
 473             for(unsigned int dim=1;dim<Dimension;dim++)
 474               general<<"\t"<<maxXYZList[phaseIndex][dim];
 475             general<<std::endl;
 476           }
 477
 478         // General
 479         general<<std::endl;
 480         general<<"@";
 481         general<<"\t"<<mean<<"\t"<<std;
 482         general<<"\t"<<meanXYZ[0]<<"\t"<<stdXYZ[0];
 483         for(unsigned int dim=1;dim<Dimension;dim++)
 484           general<<"\t"<<meanXYZ[dim]<<"\t"<<stdXYZ[dim];
 485         general<<"\t"<<max;
 486         general<<"\t"<<maxXYZ[0];
 487         for(unsigned int dim=1;dim<Dimension;dim++)
 488           general<<"\t"<<maxXYZ[dim];
 489         general<<std::endl;
 490         general.close();
 491       }
 492
 493     // Output original points
 494     if( m_ArgsInfo.original_given)
 495       {
 496
 497         std::vector<std::string> filenames;
 498         for (unsigned int i=0;i<numberOfFields;i++)
 499           {
 500             std::ostringstream number_ostr;
 501             number_ostr << i;
 502             std::string number_str =  number_ostr.str();
 503             filenames.push_back(m_ArgsInfo.original_arg+number_str);
 504           }
 505         originalDistanceLists.Write(filenames, m_Verbose );
 506       }
 507
 508
 509     // Output original magnitude points
 510     if( m_ArgsInfo.originalMag_given)
 511       {
 512         clitk::Lists<itk::FixedArray<double,1> > originalDistanceListsMag=originalDistanceLists.Norm();
 513         std::vector<std::string> filenames;
 514         for (unsigned int i=0;i<numberOfFields;i++)
 515           {
 516             std::ostringstream number_ostr;
 517             number_ostr << i;
 518             std::string number_str =  number_ostr.str();
 519             filenames.push_back(m_ArgsInfo.originalMag_arg+number_str);
 520           }
 521         originalDistanceListsMag.Write(filenames, m_Verbose );
 522       }
 523
 524     // Output displacement
 525     if( m_ArgsInfo.displacement_given)
 526       {
 527
 528         std::vector<std::string> filenames;
 529         for (unsigned int i=0;i<numberOfFields;i++)
 530           {
 531             std::ostringstream number_ostr;
 532             number_ostr << i;
 533             std::string number_str =  number_ostr.str();
 534             filenames.push_back(m_ArgsInfo.displacement_arg+number_str);
 535           }
 536         displacementLists.Write(filenames, m_Verbose );
 537       }
 538
 539
 540     // Output displacement magnitude
 541     if( m_ArgsInfo.displacementMag_given)
 542       {
 543         clitk::Lists<itk::FixedArray<double,1> > displacementListsMag=displacementLists.Norm();
 544         std::vector<std::string> filenames;
 545         for (unsigned int i=0;i<numberOfFields;i++)
 546           {
 547             std::ostringstream number_ostr;
 548             number_ostr << i;
 549             std::string number_str =  number_ostr.str();
 550             filenames.push_back(m_ArgsInfo.displacementMag_arg+number_str);
 551           }
 552         displacementListsMag.Write(filenames, m_Verbose );
 553       }
 554
 555     // Output warped points
 556     if( m_ArgsInfo.warp_given)
 557       {
 558
 559         std::vector<std::string> filenames;
 560         for (unsigned int i=0;i<m_ArgsInfo.warp_given;i++)
 561           {
 562             std::ostringstream number_ostr;
 563             number_ostr << i;
 564             std::string number_str =  number_ostr.str();
 565             filenames.push_back(m_ArgsInfo.warp_arg+number_str);
 566           }
 567         warpedPointLists.Write(filenames, m_Verbose );
 568       }
 569
 570     // Output tre
 571     if( m_ArgsInfo.tre_given)
 572       {
 573
 574         std::vector<std::string> filenames;
 575         for (unsigned int i=0;i<numberOfFields;i++)
 576           {
 577             std::ostringstream number_ostr;
 578             number_ostr << i;
 579             std::string number_str =  number_ostr.str();
 580             filenames.push_back(m_ArgsInfo.tre_arg+number_str);
 581           }
 582         treLists.Write(filenames, m_Verbose );
 583       }
 584
 585     // Output tre mag
 586     if( m_ArgsInfo.treMag_given)
 587       {
 588         clitk::Lists<itk::FixedArray<double,1> > treMagLists=treLists.Norm();
 589
 590         std::vector<std::string> filenames;
 591         for (unsigned int i=0;i<numberOfFields;i++)
 592           {
 593             std::ostringstream number_ostr;
 594             number_ostr << i;
 595             std::string number_str =  number_ostr.str();
 596             filenames.push_back(m_ArgsInfo.treMag_arg+number_str);
 597           }
 598         treMagLists.Write(filenames, m_Verbose );
 599       }
 600
 601   }
 602
 603
 604 }//end clitk
 605
 606 #endif //#define clitkCalculateTREGenericFilter_txx