vtk/vtkGdcm4DSplitter.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: vtkGdcm4DSplitter.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2011/04/15 15:14:40 $
   7   Version:   $Revision: 1.11 $
   8
   9   Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
  10   l'Image). All rights reserved. See Doc/License.txt or
  11   http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
  12
  13      This software is distributed WITHOUT ANY WARRANTY; without even
  14      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  15      PURPOSE.  See the above copyright notices for more information.
  16
  17 =========================================================================*/
  18
  19 /* Raisons ne pas utiliser itkImageSeriesReader:
  20
  21 On Wed, Feb 16, 2011 at 11:51 AM, Roger Bramon Feixas <rogerbramon@gmail.com>
  22     Hi,
  23     I'm developing with ITK 3.20 + GDCM 2.0.17 + VTK 5.6 and I've noticed
  24     itkImageSeriesReader is ~2x slower than vtkGDCMImageReader (from GDCM2).
  25     I compared both codes and I think the difference is the extra copy which
  26     itkImageSeriesReader makes from ImageFileReader's output to its own output
  27     (ImageSeriesReader::GenerateData() line 393).
  28 */
  29
  30
  31 /* ====================================================================
  32 vtkGdcm4DSplitter
  33
  34 3D, 2D+T, 3D+T, n*2D+T, 4D images are not always stored the same way :
  35         a single 'Dicom Serie',
  36         several 'Dicom series' within a single directory
  37         several 'Dicom series' within several directories
  38 A 'Dicom Serie' doesn't mean always the same thing :
  39         a given Slice along the time
  40         a given Volume at a given time
  41 Sometimes, an image within a serie is so artefacted than user decides to replace it
  42            by an other image.
  43
  44 User needs to be aware, *only him* knows want he wants to do.
  45 vtkGdcm4DSplitter class does the job for him
  46 (despite its name, it works on 3D or 2D+T images too)
  47
  48 ==> To (try to) understand how the file store is organised, user is deeply encourage to use something like:
  49
  50 PrintFile dirin=Root_Directory_holding_the_images rec > 1.txt
  51
  52 open it with any test editor, and/or (Linux users)
  53 grep a_supposed_to_be_string_of_interest 1.txt
  54
  55
  56 Aware user will have to specify some points :
  57
  58 . Choose input data
  59 -------------------
  60
  61 - a single directory
  62        bool setDirName(std::string &dirName);
  63 - a list of directories
  64        bool setVectDirName(std::vector<std::string> &vectDirName);
  65 - a list of files
  66        bool setVectFileName(std::vector<std::string> &vectFileName);
  67 - a list of gdcm::File*
  68        bool setVectGdcmFile(std::vector<GDCM_NAME_SPACE::File *> &vectGdcmFile);
  69
  70 - Recursive directory exploration
  71        void setRecursive(bool recursive);
  72
  73 . Choose 'split' criterion :
  74 ---------------------------
  75
  76  - ImagePositionPatient
  77         void setSplitOnPosition();
  78  - ImageOrientationPatient
  79         void setSplitOnOrientation();
  80  - User choosen tag
  81         void setSplitOnTag(unsigned short splitGroup, unsigned short splitElem);
  82         void setSplitConvertToFloat(bool conv);
  83  - UserDefined Function
  84         void setSortOnUserFunction (FoncComp f);
  85
  86  for 'true' 3D image sets :
  87    - if you want to get a single 3D vtkImageData, use SplitOnOrientation -i.e. no split-
  88    - if you want to get a vector of 2D vtkImageData, use SplitOnPosition  -i.e. one slice in each 'XCoherent fileset'-
  89
  90  for 'true' 4D multi-orientation image sets (i.e. a stack of n axial + m sagital + p coronal images, at different instants ...)
  91    --> this is 5D, right?
  92    (almost) nothing done, yet :
  93    . use setSplitOnly()
  94    . Use a first time vtkGdcm4DSplitter with setSplitOnOrientation();
  95    . Get the VectGdcmFileLists (a std::vector of 'XCoherent fileset')
  96    . use vtkGdcm4DSplitter, with as many setVectGdcmFile(std::vector<GDCM_NAME_SPACE::File *> &vectGdcmFile) you need
  97         one per element of std::vector<GDCM_NAME_SPACE::File *>
  98         think on 'spliting' and 'sorting' it, according to your needs.
  99
 100 . Choose 'sort' criterion :
 101 --------------------------
 102
 103  - ImagePositionPatient
 104         void setSortOnPosition();
 105  - ImageOrientationPatient
 106        ==> Only in your dreams!
 107        ==> or, please, write a IOP sorter ...
 108  - UserDefined Function
 109         void setSortOnUserFunction (FoncComp f);
 110  - File name
 111         void setSortOnFileName()
 112
 113 - Choose 'flip' option
 114 ----------------------
 115 VTK reader 'flop upside down' the images to display them correctly ?!?
 116 When a 'math guy' gets the image, the first thing he does is flipping it again.
 117 Avoid loosing CPU time using :
 118       bool SetFlipY(bool); // Wrong = no flip
 119
 120 . Execute :
 121 -----------
 122         bool Go();
 123
 124 . Get the result
 125 ----------------
 126
 127 // -a single vtkImageData:
 128 //        vtkImageData *GetImageData();
 129 - a vector of vtkImageData
 130         std::vector<vtkImageData*> *GetImageDataVector();
 131
 132   ===================================================================== */
 133
 134 #include "vtkGdcmReader.h"
 135 #include "vtkGdcm4DSplitter.h"
 136 #include <algorithm>
 137 #include "gdcmSerieHelper.h" // for ImagePositionPatientOrdering()
 138 #include <stdlib.h> // for atof
 139
 140 // Constructor / Destructor
 141 /**
 142  * \brief   Constructor from a given vtkGdcm4DSplitter
 143  */
 144  vtkGdcm4DSplitter::vtkGdcm4DSplitter() :
 145                  SplitOnPosition(false), SplitOnOrientation(false), SplitOnTag(false),
 146                  SplitGroup(0), SplitElem(0), SplitConvertToFloat(false),
 147
 148                  SortOnPosition(false),  SortOnOrientation(false),  SortOnTag(false), SortOnFileName(false), SortOnUserFunction(false),
 149                  SortGroup(0),  SortElem(0), SortConvertToFloat(false),
 150
 151                  Recursive(false), TypeDir(0),
 152                  verbose(false)
 153  {
 154
 155  }
 156
 157 /**
 158  * \brief   Canonical destructor.
 159  */
 160  vtkGdcm4DSplitter::~vtkGdcm4DSplitter()
 161  {
 162     /// \TODO : delete everything that must be!
 163  }
 164
 165        // Locate Data to process
 166        // ======================
 167 /**
 168  * \brief sets the directories exploration mode
 169  * @param recursive whether we want explore recursively the root Directory
 170  */
 171 void  vtkGdcm4DSplitter::setRecursive(bool recursive)
 172 {
 173    Recursive=recursive;
 174 }
 175
 176 /**
 177  * \brief Sets the root Directory to get the images from
 178  * @param   dirName name of the directory to deal with
 179  */
 180  bool vtkGdcm4DSplitter::setDirName(std::string &dirName)
 181  {
 182     if ( ! GDCM_NAME_SPACE::DirList::IsDirectory(dirName) )
 183     {
 184        std::cout << "[" << dirName << "] is NOT a directory" << std::endl;
 185        return false;
 186     }
 187     DirName = dirName;
 188     TypeDir=1;
 189     return true;
 190  }
 191
 192 /**
 193  * \brief Sets a list of Directories to get the images from
 194  * @param   vectDirName vector of directory names to deal with
 195  */
 196  bool vtkGdcm4DSplitter::setVectDirName(std::vector<std::string> &vectDirName)
 197  {
 198     int nbDir = vectDirName.size();
 199     for (int iDir=0; iDir<nbDir; iDir++)
 200     {
 201        if ( ! GDCM_NAME_SPACE::DirList::IsDirectory(vectDirName[iDir]) )
 202        {
 203           std::cout << "[" << vectDirName[iDir] << "] is NOT a directory" << std::endl;
 204           return false;
 205        }
 206     }
 207
 208     VectDirName = vectDirName;
 209     TypeDir=2;
 210     return true;
 211  }
 212
 213 /**
 214  * \brief Sets a list of files read
 215  * @param   vectFileName vector of file names to deal with
 216  */
 217  bool vtkGdcm4DSplitter::setVectFileName(std::vector<std::string> &vectFileName)
 218  {
 219     if ( vectFileName.size() == 0)
 220     {
 221           std::cout << "[ vectFileName ] : empty list" << std::endl;
 222           return false;
 223     }
 224     VectFileName = vectFileName;
 225     TypeDir=3;
 226     return true;
 227  }
 228
 229 /**
 230  * \brief Sets a vector of gdcm::File *
 231  * @param   vectGdcmFileName vector of gdcm::File *
 232  */
 233
 234  bool vtkGdcm4DSplitter::setVectGdcmFile(GDCM_NAME_SPACE::FileList *vectGdcmFile)
 235  {
 236     if ( vectGdcmFile->size() == 0)
 237     {
 238           std::cout << "[ vectGdcmFile ] : empty list" << std::endl;
 239           return false;
 240     }
 241     TypeDir=4;
 242     VectGdcmFile = vectGdcmFile;
 243  }
 244
 245
 246        // Split
 247        // =====
 248
 249 /**
 250  * \brief   asks for splitting  Filesets according to the Position
 251  */
 252
 253  void  vtkGdcm4DSplitter::setSplitOnPosition()
 254  {
 255     SplitOnPosition=true;
 256     SplitOnOrientation=false;
 257     SplitOnTag=false;
 258  }
 259 /**
 260  * \brief   asks for splitting  Filesets according to the Orientation
 261  */
 262  void  vtkGdcm4DSplitter::setSplitOnOrientation()
 263  {
 264     SplitOnPosition=false;
 265     SplitOnOrientation=true;
 266     SplitOnTag=false;
 267  }
 268 /**
 269  * \brief   asks for splitting  Filesets according to the value of a given Tag
 270  * @param   group  group number of the target Element
 271  * @param   element element number of the target Element
 272  */
 273  void  vtkGdcm4DSplitter::setSplitOnTag(unsigned short int splitGroup, unsigned short int splitElem)
 274  {
 275     SplitOnPosition=false;
 276     SplitOnOrientation=false;
 277     SplitOnTag=true;
 278     SplitGroup=splitGroup;
 279     SplitElem=splitElem;
 280  }
 281 /**
 282  * \brief   asks for converting to 'float' the tag values used as a splitting criteria (lexicographic order may not be suitable)
 283  */
 284  void  vtkGdcm4DSplitter::setSplitConvertToFloat(bool conv) {SplitConvertToFloat=conv;}
 285
 286 /**
 287  * \brief   asks for splitting Filesets according to what was asked for (no sorting, no reading data)
 288  */
 289  void  vtkGdcm4DSplitter::setSplitOnly(bool s)
 290  {
 291     SplitOnly = s;
 292  }
 293        // Sort
 294        // ====
 295
 296  void  vtkGdcm4DSplitter::setSortOnPosition()
 297  {
 298     SortOnPosition=true;
 299     SortOnOrientation=false;
 300     SortOnTag=false;
 301     SortOnFileName=false;
 302     SortOnUserFunction=false;
 303     SortOnPosition=true;
 304  }
 305
 306       // use setSortOnUserFunction, instead!
 307       // void setSortOnTag(unsigned short int sortGroup, unsigned short int sortElem)
 308       // {
 309       //    SortOnPosition=false;
 310       //    SortOnOrientation=false;
 311       //    SortOnTag=true;
 312       //    SortOnFileName=false;
 313       //    SortOnUserFunction=false;
 314       //    SortGroup=sortGroup;  SortElem=sortElem;
 315       // }
 316
 317
 318 /**
 319  * \brief sets a user supplied function (comparison)
 320  * @param f comparison function
 321  */
 322  void  vtkGdcm4DSplitter::setSortOnUserFunction (FoncComp f)
 323  {
 324     UserCompareFunction=f;
 325     SortOnPosition=false;
 326     SortOnOrientation=false;
 327     SortOnTag=false;
 328     SortOnFileName=false;
 329     SortOnUserFunction=true;
 330   }
 331
 332
 333  //  void setSortConvertToFloat(bool conv)
 334  //  {
 335  //     SortConvertToFloat=conv;
 336  //  }
 337
 338 /**
 339  * \brief asks for sorting  the images, according to their File Name
 340  */
 341  void  vtkGdcm4DSplitter::setSortOnFileName()
 342  {
 343     SortOnPosition=false;
 344     SortOnOrientation=false;
 345     SortOnTag=false;
 346     SortOnFileName=true;
 347     SortOnUserFunction=false;
 348  }
 349
 350
 351  std::vector<vtkImageData*> * vtkGdcm4DSplitter::GetImageDataVector()
 352  {
 353     if (SplitOnly)
 354        return NULL;
 355
 356      return ImageDataVector;
 357  }
 358
 359  std::vector<GDCM_NAME_SPACE::FileList *> *vtkGdcm4DSplitter::GetVectGdcmFileLists()
 360  {
 361     if (SplitOnly)
 362         return NULL;
 363
 364     GDCM_NAME_SPACE::XCoherentFileSetmap::iterator it;
 365     for ( it = xcm.begin();
 366           it != xcm.end();
 367         ++it)
 368     {
 369        VectGdcmFileLists.push_back((*it).second);
 370     }
 371     return  &VectGdcmFileLists;
 372  }
 373
 374
 375  vtkImageData *vtkGdcm4DSplitter::GetImageData()
 376  {
 377    if (SplitOnly)
 378       return NULL;
 379    return (*ImageDataVector)[0];
 380  }
 381
 382
 383  bool vtkGdcm4DSplitter::CompareOnSortTagConvertToFloat(GDCM_NAME_SPACE::File *file1, GDCM_NAME_SPACE::File *file2)
 384  {
 385   /* if (verbose) printf ("%04x %04x\n", this->SortGroup,this->SortElem);
 386      if (verbose) std :: cout << file1->GetEntryString(SortGroup,SortElem).c_str() << " : "
 387                             << atof(file1->GetEntryString(SortGroup,SortElem).c_str())
 388                             << std::endl;
 389 */
 390 //   return atof(file1->GetEntryString(vtkGdcm4DSplitter::SortGroup,vtkGdcm4DSplitter::SortElem).c_str()) < atof(file2->GetEntryString(vtkGdcm4DSplitter::SortGroup,vtkGdcm4DSplitter::SortElem).c_str());
 391    return atof(file1->GetEntryString(SortGroup,SortElem).c_str()) < atof(file2->GetEntryString(SortGroup,SortElem).c_str());
 392  }
 393
 394  bool vtkGdcm4DSplitter::CompareOnSortTag(GDCM_NAME_SPACE::File *file1, GDCM_NAME_SPACE::File *file2)
 395  {
 396    return file1->GetEntryString(vtkGdcm4DSplitter::SortGroup,vtkGdcm4DSplitter::SortElem) < file2->GetEntryString(vtkGdcm4DSplitter::SortGroup,vtkGdcm4DSplitter::SortElem);
 397  }
 398
 399
 400  bool vtkGdcm4DSplitter::Go()
 401  {
 402    if (!SplitOnPosition && !SplitOnOrientation && !SplitOnTag)
 403    {
 404        ///\TODO (?) Throw an exception "Choose Splitting mode before!"
 405        std::cout << "Choose Splitting mode before!" << std::endl;
 406        return false;
 407    }
 408
 409    /// How To :
 410    /*
 411    entree nom de directory / Vecteur de noms?
 412    recursif oui/non
 413    recuperer la liste des gdcm::File*
 414    passer a SerieHelper (pas de check du Serie UID)
 415    set critere de split
 416
 417    trier chaque Coherent file set
 418    passer chacun a un vtkGcdmReader
 419    retourner le (vecteur de) vtkImageData
 420    */
 421
 422    GDCM_NAME_SPACE::SerieHelper *s;
 423    s = GDCM_NAME_SPACE::SerieHelper::New();
 424
 425    GDCM_NAME_SPACE::File *f;
 426    GDCM_NAME_SPACE::DirListType fileNames;
 427
 428  //
 429  // Fill  fileNames with the user supplied file names (in any)
 430  // --------------------------------------
 431  //
 432    if (TypeDir == 0 )
 433    {
 434       ///\TODO (?) Throw an exception "Set input Directory name(s) / file names  before!"
 435       std::cout << "Set input Directory name(s) / file names  before!" << std::endl;
 436       return false;
 437    }
 438    else if (TypeDir == 1)
 439    {
 440       GDCM_NAME_SPACE::DirList dirlist(DirName, Recursive); // NO recursive exploration
 441       fileNames = dirlist.GetFilenames(); // all the file names
 442    }
 443
 444    else if (TypeDir == 2)
 445    {
 446       int nbDir = VectDirName.size();
 447       GDCM_NAME_SPACE::DirListType tmpFileNames;
 448       for (int iDir=0; iDir<nbDir; iDir++)
 449       {
 450         GDCM_NAME_SPACE::DirList dirlist(VectDirName[iDir], Recursive);
 451         tmpFileNames = dirlist.GetFilenames();
 452         // Concat two std::vector
 453         //vector1.insert( vector1.end(), vector2.begin(), vector2.end() );
 454         fileNames.insert( fileNames.end(), tmpFileNames.begin(), tmpFileNames.end() );
 455       }
 456    }
 457    else if (TypeDir == 3)
 458    {
 459       fileNames=VectFileName;
 460    }
 461
 462  //
 463  // Fill l with the gdcm::File* corresponding to the files
 464  // --------------------------------------
 465  //
 466
 467    GDCM_NAME_SPACE::FileList *l = new GDCM_NAME_SPACE::FileList; // (set of gdcm::File)
 468
 469    if (TypeDir == 4)
 470    {
 471    // User passed a vector of gdcm::File*
 472       l = VectGdcmFile;
 473    }
 474    else
 475    {
 476       double floatTagvalue;
 477       // Loop on all the gdcm-readable files
 478       for (GDCM_NAME_SPACE::DirListType::iterator it = fileNames.begin();
 479                                                   it != fileNames.end();
 480                                                 ++it)
 481       {
 482          int maxSize  = 0x7fff;         // load Elements of any length
 483          f = GDCM_NAME_SPACE::File::New();
 484          f->SetMaxSizeLoadEntry(maxSize);
 485          f->SetFileName( *it );
 486          if (f->Load())
 487             l->push_back(f);
 488          else
 489             std::cout << " Fail to load [" <<  *it << "]" << std::endl;
 490       }
 491    }
 492
 493
 494 //
 495 // Split the gdcm::File* set, according to user's requierements
 496 // ------------------------------------------------------------
 497 //
 498    if (SplitOnOrientation)
 499    {
 500       s->SetDropDuplicatePositions(false);
 501       xcm = s->SplitOnOrientation(l);
 502    }
 503    else if (SplitOnPosition)
 504    {
 505       s->SetDropDuplicatePositions(true);
 506       xcm = s->SplitOnPosition(l);
 507
 508       // the key of xcm follows lexicographical order
 509       // (that may be different than the 'distance' order)
 510       // we have to reorganize it!
 511
 512       reorgXCoherentFileSetmap(xcm);
 513    }
 514    else if (SplitOnTag)
 515    {
 516          s->SetDropDuplicatePositions(false);
 517
 518          // Crashes if DataElement not found
 519          //std:: cout << GDCM_NAME_SPACE::Global::GetDicts()->GetDefaultPubDict()->GetEntry(groupelem[0],groupelem[1])->GetName() << std::endl;
 520             if ( ! SplitConvertToFloat )
 521                xcm = s->SplitOnTagValue(l, SplitGroup, SplitElem);
 522             else
 523             {
 524                 xcm = s->SplitOnTagValueConvertToFloat(l, SplitGroup, SplitElem);
 525             }
 526    }
 527
 528    if (xcm.size() == 0)
 529    {
 530       if(verbose) std::cout << "Empty XCoherent File Set after 'split' ?!?" << std::endl;
 531       return false;
 532    }
 533 /*
 534    else if (xcm.size() == 1)
 535       TypeResult=1;
 536    else
 537       TypeResult=2;
 538 */
 539
 540    if(SplitOnly)
 541       return true;
 542 //
 543 //
 544 // ------------------------------------------------------------
 545 //
 546    ImageDataVector = new std::vector<vtkImageData*>;
 547    /// \TODO move inside the loop, or be clever using vtk!
 548   // vtkGdcmReader *reader = vtkGdcmReader::New(); // move inside the loop, or be clever using vtk!
 549
 550  // XCoherentFileSetmap map < critère de split, FileList (= std::vector de gdcm::File*) >
 551
 552    for (GDCM_NAME_SPACE::XCoherentFileSetmap::iterator i = xcm.begin();
 553                                                   i != xcm.end();
 554                                                 ++i)
 555    {
 556       vtkGdcmReader *reader = vtkGdcmReader::New(); /// \TODO FIXME : unable to delete!
 557
 558       reader->SetFlipY(FlipY);
 559       // better user SetFileLowerLeft()
 560       /// \TODO : modify vtkGdcmReader !
 561       if (verbose)
 562                std::cout << " --- xCoherentName = [" << (*i).first << "]" << std::endl;
 563
 564       if (SortOnPosition)
 565       {
 566               if (verbose) std::cout << "SortOnPosition" << std::endl;
 567               // (will be IPPSorter, in GDCM2)
 568               s->ImagePositionPatientOrdering((*i).second);
 569               if (verbose) std::cout << "out of SortOnPosition" << std::endl;
 570       }
 571
 572       else if (SortOnOrientation)
 573       {
 574               if (verbose) std::cout << "SortOnOrientation" << std::endl;
 575             /// \TODO (?) SortOnOrientation()
 576
 577             // we still miss an algo to sort an Orientation, given by 6 cosines!
 578             //  Anything like this, in GDCM2?
 579             std::cout << "SortOnOrientation : not so easy - I(mage)O(rientation)P(atient)Sorter still missing! -" << std::endl;
 580             // have a look at SerieHelper::SplitOnOrientation() to have an idea of the mess!
 581
 582             //Better sort on the file name, right now...
 583              s->FileNameOrdering((*i).second);
 584       }
 585
 586       else if (SortOnFileName)
 587       {
 588          if (verbose) std::cout << "SortOnFileName" << std::endl;
 589          if (verbose) std::cout << "taille " << ((*i).second)->size() << std::endl;
 590
 591          s->FileNameOrdering((*i).second);
 592          if (verbose) std::cout << "Out of SortOnFileName" << std::endl;
 593       }
 594
 595       else if (SortOnTag)
 596       {
 597          if (verbose) std::cout << "SortOnTag" << std::endl;
 598          printf ("--> %04x %04x\n", SortGroup,SortElem);
 599          std::cout << "Sorry, troubles not solved yet; use SortOnUserFunction, right now!" << std::endl;
 600
 601         /*        ==> WARNING : This one has troubles; do NOT use it, right now!
 602         // a pointer to fonction cannot be casted as a pointer to member function!
 603         // Use SortOnUserFunction, instead!
 604
 605         //  if ( SortConvertToFloat )
 606          //    s->SetUserLessThanFunction( reinterpret_cast<bool (*)(gdcm13::File*, gdcm13::File*)>
 607                                                                  ( &vtkGdcm4DSplitter::CompareOnSortTagConvertToFloat));
 608         //  else
 609         //     s->SetUserLessThanFunction( reinterpret_cast<bool (*)(gdcm13::File*, gdcm13::File*)>
 610                                                                  ( &vtkGdcm4DSplitter::CompareOnSortTag));
 611
 612          // Anything like this, in GDCM2?
 613         //  s->UserOrdering((*i).second);
 614         */
 615
 616          //if (verbose) std::cout << "Out of SortOnTag" << std::endl;
 617          std::cout << "NO ordering performed  :-( " << std::endl;
 618       }
 619
 620       else if (SortOnUserFunction)
 621       {
 622           if (verbose) std::cout << "SortOnUserFunction" << std::endl;
 623           s->SetUserLessThanFunction( UserCompareFunction );
 624          // Anything like this, in GDCM2?
 625          s->UserOrdering((*i).second);
 626          if (verbose) std::cout << "Out of SortOnUserFunction" << std::endl;
 627       }
 628
 629        reader->SetCoherentFileList((*i).second);
 630        reader->Update();
 631
 632        /// \TODO : remove the following
 633        if (verbose) {
 634           std::cout << "reader->GetOutput() :" << std::endl;
 635           reader->GetOutput()->PrintSelf(std::cout, vtkIndent(2));
 636        }
 637
 638        ImageDataVector->push_back(reader->GetOutput() );
 639
 640        std::vector<vtkImageData*>::iterator it;
 641        if (verbose)
 642        for(it=ImageDataVector->begin(); it!=ImageDataVector->end(); ++it) {
 643          std::cout << "-in vtkGdcm4DSplitter--------------------------" << std::endl;
 644          (*it)->PrintSelf(std::cout, vtkIndent(2));
 645        }
 646        std::cout << std::endl;
 647    }
 648
 649    //reader->Delete();  // \TODO : fix
 650    s->Delete();
 651    f->Delete();
 652    delete l;
 653
 654    return true;
 655  }
 656
 657
 658  void vtkGdcm4DSplitter::reorgXCoherentFileSetmap(GDCM_NAME_SPACE::XCoherentFileSetmap &xcm)
 659  {
 660  /*
 661  the key of the 'XCoherentFileSetmap', is a std::string, used as if it was found in the Dicom header
 662  Normaly(?), it's suitable for almost everything ...
 663  ... but the 'Image Position Patient'.
 664  We need to order the 'XCoherentFileSetmap' (NOT the content of each XCoherentFileSet!) according to the IPP,
 665  using Jolinda Smith's algorithm.
 666  (we use a subset of the one defined in gdcm::SerieHelper)
 667
 668 */
 669
 670    ELEM e;
 671    std::vector<ELEM> vectElem;
 672
 673 /* remenber :
 674    typedef struct
 675    {
 676       std::string strIPP;
 677       double dist;
 678       GDCM_NAME_SPACE::File *file;
 679    } ELEM;
 680 */
 681
 682    bool Debug=true;
 683
 684    for (GDCM_NAME_SPACE::XCoherentFileSetmap::iterator i = xcm.begin();
 685                                                        i != xcm.end();
 686                                                      ++i)
 687    {
 688       if (verbose)
 689                std::cout << "--- xCoherentName = [" << (*i).first << "]" << std::endl;
 690
 691       e.strIPP  = (*i).first;
 692       e.file = *(((*i).second)->begin()); // all the gdcm::File of a given xcm item *have* the same IPP; first one is enough
 693       e.dist=0.0;
 694       vectElem.push_back(e);
 695    }
 696    sortVectElem(&vectElem);
 697
 698    // now, create the final std::map !
 699    // final_xcm<to_str(e.dist , xcm[e.strIPP]>
 700    // xcm = final_xcm;
 701    // check what we need to free !
 702
 703    int dist;
 704    char c_dist[100];
 705    std::string str_dist;
 706    int lgr=vectElem.size();
 707    GDCM_NAME_SPACE::XCoherentFileSetmap final_xcm;
 708    for (int i2=0; i2<lgr; i2++)
 709    {
 710       dist =  (vectElem[i2].dist*1000);
 711       sprintf(c_dist,"%010d",dist);
 712       str_dist = c_dist;
 713 /*
 714       std::cout << "dist " << vectElem[i2].dist
 715                 << " str_dist " << str_dist
 716                 << " IPP " << vectElem[i2].strIPP
 717                 << std::endl;
 718
 719 */
 720       final_xcm[str_dist] = xcm[vectElem[i2].strIPP];
 721    }
 722
 723    /// \TODO : check what needs to be cleared // JPR
 724
 725    xcm = final_xcm;
 726
 727  }
 728
 729
 730 bool vtkGdcm4DSplitter::sortVectElem(std::vector<ELEM> *fileList)
 731 {
 732 //based on Jolinda Smith's algorithm
 733 // NOTE : if you need to use Jolinda Smith's algorithm, get the one inside gdcm::SerieHelper
 734 // this one is a light version.
 735
 736 //Tags always use the same coordinate system, where "x" is left
 737 //to right, "y" is posterior to anterior, and "z" is foot to head (RAH).
 738
 739    //iop is calculated based on the file file
 740    float cosines[6];
 741    double normal[3];
 742    double ipp[3];
 743    double dist;
 744    double min = 0, max = 0;
 745    bool first = true;
 746
 747    //double ZSpacing; // useless here! // JPR
 748    bool DirectOrder = true; // remove it!
 749
 750   // ZSpacing = -1.0;  // will be updated if process doesn't fail
 751
 752    //std::multimap<double,File *> distmultimap; // JPR
 753    std::multimap<double,ELEM> distmultimap; // JPR
 754
 755    // Use a multimap to sort the distances from 0,0,0
 756    //for ( FileList::const_iterator // JPR
 757    for ( std::vector<ELEM>::iterator   // JPR
 758          it = fileList->begin();
 759          it != fileList->end(); ++it )
 760    {
 761       //gdcmDebugMacro("deal with " << (*it)->file->GetFileName() );
 762       if ( first )
 763       {
 764          (*it).file->GetImageOrientationPatient( cosines );
 765
 766    // The "Image Orientation Patient" tag gives the direction cosines
 767    // for the rows and columns for the three axes defined above.
 768    // Typical axial slices will have a value 1/0/0/0/1/0:
 769    // rows increase from left to right,
 770    // columns increase from posterior to anterior. This is your everyday
 771    // "looking up from the bottom of the head with the eyeballs up" image.
 772
 773    // The "Image Position Patient" tag gives the coordinates of the first
 774    // voxel in the image in the "RAH" coordinate system, relative to some
 775    // origin.
 776
 777    // First, calculate the slice normal from IOP :
 778
 779          // You only have to do this once for all slices in the volume. Next,
 780          // for each slice, calculate the distance along the slice normal
 781          // using the IPP ("Image Position Patient") tag.
 782          // ("dist" is initialized to zero before reading the first slice) :
 783
 784          normal[0] = cosines[1]*cosines[5] - cosines[2]*cosines[4];
 785          normal[1] = cosines[2]*cosines[3] - cosines[0]*cosines[5];
 786          normal[2] = cosines[0]*cosines[4] - cosines[1]*cosines[3];
 787
 788    // For each slice (here : the first), calculate the distance along
 789    // the slice normal using the IPP tag
 790
 791          ipp[0] = (*it).file->GetXOrigin();
 792          ipp[1] = (*it).file->GetYOrigin();
 793          ipp[2] = (*it).file->GetZOrigin();
 794
 795          dist = 0;
 796          for ( int i = 0; i < 3; ++i )
 797          {
 798             dist += normal[i]*ipp[i];
 799          }
 800
 801          //gdcmDebugMacro("dist : " << dist);
 802          distmultimap.insert(std::pair<const double,ELEM>(dist, *it));
 803
 804          max = min = dist;
 805          first = false;
 806       }
 807       else
 808       {
 809    // Next, for each slice, calculate the distance along the slice normal
 810    // using the IPP tag
 811          ipp[0] = (*it).file->GetXOrigin();
 812          ipp[1] = (*it).file->GetYOrigin();
 813          ipp[2] = (*it).file->GetZOrigin();
 814
 815          dist = 0;
 816          for ( int i = 0; i < 3; ++i )
 817          {
 818             dist += normal[i]*ipp[i];
 819          }
 820
 821          (*it).dist = dist; // JPR
 822
 823          distmultimap.insert(std::pair<const double,ELEM>(dist, *it));
 824          //gdcmDebugMacro("dist : " << dist);
 825          min = (min < dist) ? min : dist;
 826          max = (max > dist) ? max : dist;
 827       }
 828    }
 829
 830   // gdcmDebugMacro("After parsing vector, nb of elements : " << fileList->size() );
 831
 832 /*Useless here. // JPR
 833
 834    // Find out if min/max are coherent
 835    if ( min == max )
 836    {
 837      gdcmWarningMacro("Looks like all images have the exact same image position. "
 838                       << "No PositionPatientOrdering sort performed. "
 839                       << "No 'ZSpacing' calculated! ");
 840      return false;
 841    }
 842 */
 843
 844 /* Useless here, 'split' already done. // JPR
 845
 846    // Check to see if image shares a common position
 847    bool ok = true;
 848    for (std::multimap<double, File *>::iterator it2 = distmultimap.begin();
 849         it2 != distmultimap.end();
 850         ++it2)
 851    {
 852
 853       gdcmDebugMacro("Check if image shares a common position : " << ((*it2).second).file->GetFileName() );
 854
 855       if (distmultimap.count((*it2).first) != 1)
 856       {
 857          gdcmWarningMacro("File: ["
 858               << ((*it2).second->GetFileName())
 859               << "] : more than ONE file at distance: '"
 860               << (*it2).first
 861               << " (position is not unique!) "
 862               << "No PositionPatientOrdering sort performed. "
 863               << "No 'ZSpacing' calculated! ");
 864
 865          ok = false;
 866       }
 867    }
 868    if (!ok)
 869    {
 870       if (! DropDuplicatePositions)
 871          return false;
 872    }
 873
 874 */
 875
 876 // Now, we can calculate Z Spacing as the difference
 877 // between the "dist" values for the first two slices.
 878
 879 // The following (un)-commented out code is let here
 880 // to be re-used by whomsoever is interested...
 881
 882     //std::multimap<double, ELEM>::iterator it5 = distmultimap.begin();
 883     //double d1 = (*it5).first;
 884     //it5++;
 885     //double d2 = (*it5).first;
 886     //ZSpacing = d1-d2;
 887     //if (ZSpacing < 0.0)
 888     //   ZSpacing = - ZSpacing;
 889
 890    fileList->clear();  // doesn't delete list elements, only nodes
 891
 892 // Acording to user requierement, we sort direct order or reverse order.
 893    if (DirectOrder)
 894    {
 895       for (std::multimap<double, ELEM>::iterator it3 = distmultimap.begin();
 896            it3 != distmultimap.end();
 897            ++it3)
 898       {
 899          fileList->push_back( (*it3).second );
 900
 901 /*       useless here! // JPR
 902
 903          if (DropDuplicatePositions)
 904          {
 905             // ImagePositionPatientOrdering  wrong duplicates are found ???
 906             // --> fixed. See comment
 907
 908             it3 =  distmultimap.upper_bound((*it3).first); // skip all duplicates
 909            // the upper_bound function increments the iterator to the next non-duplicate entry
 910            // The for loop iteration also increments the iterator, which causes the code to skip every other image
 911            // --> decrement the iterator after the upper_bound function call
 912             it3--;
 913             if (it3 == distmultimap.end() )  // if last image, stop iterate
 914                break;
 915          }
 916 */
 917       }
 918    }
 919    else // user asked for reverse order
 920    {
 921       std::multimap<double, ELEM>::const_iterator it4;
 922       it4 = distmultimap.end();
 923       do
 924       {
 925          it4--;
 926          fileList->push_back( (*it4).second );
 927
 928 /* useless here // JPR
 929
 930          if (DropDuplicatePositions)  // skip all duplicates
 931          {
 932             // lower_bound finds the next element that is
 933             // less than or *equal to* the current value!
 934             //it4 =  distmultimap.lower_bound((*it4).first);
 935
 936            // David Feng's fix
 937            std::multimap<double, ELEM>::const_iterator itPrev = it4;
 938            while (itPrev->first == it4->first)
 939               --itPrev;
 940            it4 = itPrev;
 941
 942            if (it4 == distmultimap.begin() ) // if first image, stop iterate
 943                break;
 944          }
 945 */
 946       } while (it4 != distmultimap.begin() );
 947    }
 948
 949    distmultimap.clear();
 950
 951    return true;
 952
 953 }