vtk/vtkGdcm4DSplitter.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: vtkGdcm4DSplitter.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2011/04/21 09:14:31 $
   7   Version:   $Revision: 1.14 $
   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       void 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), SplitOnly(false),
 147
 148                  SortOnPosition(false),  SortOnOrientation(false), /* SortOnTag(false),*/ SortOnFileName(false), SortOnUserFunction(false),
 149                  /*SortGroup(0),  SortElem(0), SortConvertToFloat(false),*/
 150
 151                  Recursive(false),FlipY(true), 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 names to 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 an already filled std::vector of gdcm::File *
 231  * @param   vectGdcmFileName std::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 criterion
 247        // ===============
 248        //
 249
 250 /**
 251  * \brief   asks for splitting  Filesets according to the Position
 252  */
 253  void  vtkGdcm4DSplitter::setSplitOnPosition()
 254  {
 255     SplitOnPosition=true;
 256     SplitOnOrientation=false;
 257     SplitOnTag=false;
 258  }
 259
 260 /**
 261  * \brief   asks for splitting  Filesets according to the Orientation
 262  */
 263  void  vtkGdcm4DSplitter::setSplitOnOrientation()
 264  {
 265     SplitOnPosition=false;
 266     SplitOnOrientation=true;
 267     SplitOnTag=false;
 268  }
 269
 270 /**
 271  * \brief   asks for splitting  Filesets according to the value of a given Tag
 272  * @param   group  group number of the target Element
 273  * @param   element element number of the target Element
 274  */
 275  void  vtkGdcm4DSplitter::setSplitOnTag(unsigned short int splitGroup, unsigned short int splitElem)
 276  {
 277     SplitOnPosition=false;
 278     SplitOnOrientation=false;
 279     SplitOnTag=true;
 280     SplitGroup=splitGroup;
 281     SplitElem=splitElem;
 282  }
 283
 284 /**
 285  * \brief   asks for converting to 'float' the tag values used as a splitting criteria (lexicographic order may not be suitable)
 286  */
 287  void  vtkGdcm4DSplitter::setSplitConvertToFloat(bool conv) {SplitConvertToFloat=conv;}
 288
 289 /**
 290  * \brief   asks for splitting Filesets according to what was asked for (no sorting, no reading data)
 291  */
 292  void  vtkGdcm4DSplitter::setSplitOnly(bool s)
 293  {
 294     SplitOnly = s;
 295  }
 296
 297        // Sort criterion
 298        // ==============
 299        //
 300
 301 /**
 302  * \brief   asks for IPP sorting each XCoherent gdcm::FILE set
 303  */
 304  void  vtkGdcm4DSplitter::setSortOnPosition()
 305  {
 306     SortOnPosition=true;
 307     SortOnOrientation=false;
 308     //SortOnTag=false;
 309     SortOnFileName=false;
 310     SortOnUserFunction=false;
 311     SortOnPosition=true;
 312  }
 313
 314  // use setSortOnUserFunction, instead!
 315  // void setSortOnTag(unsigned short int sortGroup, unsigned short int sortElem)
 316  // {
 317  //    SortOnPosition=false;
 318  //    SortOnOrientation=false;
 319  //    SortOnTag=true;
 320  //    SortOnFileName=false;
 321  //    SortOnUserFunction=false;
 322  //    SortGroup=sortGroup;  SortElem=sortElem;
 323  // }
 324
 325
 326 /**
 327  * \brief sets a user supplied function (comparison)
 328  * @param f comparison function
 329  */
 330  void  vtkGdcm4DSplitter::setSortOnUserFunction (FoncComp f)
 331  {
 332     UserCompareFunction=f;
 333     SortOnPosition=false;
 334     SortOnOrientation=false;
 335     //SortOnTag=false;
 336     SortOnFileName=false;
 337     SortOnUserFunction=true;
 338   }
 339
 340  //  void setSortConvertToFloat(bool conv)
 341  //  {
 342  //     SortConvertToFloat=conv;
 343  //  }
 344
 345 /**
 346  * \brief asks for sorting each XCoherent gdcm::FILE set, according to the File names
 347  */
 348  void  vtkGdcm4DSplitter::setSortOnFileName()
 349  {
 350     SortOnPosition=false;
 351     SortOnOrientation=false;
 352     //SortOnTag=false;
 353     SortOnFileName=true;
 354     SortOnUserFunction=false;
 355  }
 356
 357 /**
 358  * \brief returns a std::vector of gdcm::FileList* (gdcm::FileList : actually, a std::vector of gdcm::File*)
 359  */
 360  std::vector<GDCM_NAME_SPACE::FileList *> *vtkGdcm4DSplitter::GetVectGdcmFileLists()
 361  {
 362     if (SplitOnly)
 363         return &VectGdcmFileLists;
 364
 365     GDCM_NAME_SPACE::XCoherentFileSetmap::iterator it;
 366     for ( it = xcm.begin();
 367           it != xcm.end();
 368         ++it)
 369     {
 370        VectGdcmFileLists.push_back((*it).second);
 371     }
 372     return  &VectGdcmFileLists;
 373  }
 374
 375 /**
 376  * \brief returns a std::vector of [2D/3D, depending on what was passed] vtkImageData*
 377  */
 378  std::vector<vtkImageData*> * vtkGdcm4DSplitter::GetImageDataVector()
 379  {
 380     if (SplitOnly)
 381        return NULL;
 382      return ImageDataVector;
 383  }
 384
 385 /**
 386  * \brief when user _knows_ only _one_ vtkImageData* is returned he may be interested in not getting a vector...
 387  */
 388  vtkImageData *vtkGdcm4DSplitter::GetImageData()
 389  {
 390    if (SplitOnly)
 391       return NULL;
 392    return (*ImageDataVector)[0];
 393  }
 394
 395
 396 // bool vtkGdcm4DSplitter::CompareOnSortTagConvertToFloat(GDCM_NAME_SPACE::File *file1, GDCM_NAME_SPACE::File *file2)
 397 // {
 398 //   return atof(file1->GetEntryString(vtkGdcm4DSplitter::SortGroup,vtkGdcm4DSplitter::SortElem).c_str()) < atof(file2->GetEntryString(vtkGdcm4DSplitter::SortGroup,vtkGdcm4DSplitter::SortElem).c_str());
 399 //   return atof(file1->GetEntryString(SortGroup,SortElem).c_str()) < atof(file2->GetEntryString(SortGroup,SortElem).c_str());
 400 // }
 401
 402 // bool vtkGdcm4DSplitter::CompareOnSortTag(GDCM_NAME_SPACE::File *file1, GDCM_NAME_SPACE::File *file2)
 403 // {
 404 //   return file1->GetEntryString(vtkGdcm4DSplitter::SortGroup,vtkGdcm4DSplitter::SortElem) < file2->GetEntryString(vtkGdcm4DSplitter::SortGroup,vtkGdcm4DSplitter::SortElem);
 405 // }
 406
 407
 408  bool vtkGdcm4DSplitter::Go()
 409  {
 410    /// How To :
 411    /*
 412    entree nom de directory / Vecteur de noms?
 413    recursif oui/non
 414    recuperer la liste des gdcm::File*
 415    passer a SerieHelper (pas de check du Serie UID)
 416    set critere de split
 417
 418    trier chaque Coherent file set
 419    passer chacun a un vtkGcdmReader
 420    retourner le (vecteur de) vtkImageData
 421    */
 422
 423    if (!SplitOnPosition && !SplitOnOrientation && !SplitOnTag)
 424    {
 425        ///\TODO (?) Throw an exception "Choose Splitting mode before!"
 426        std::cout << "Choose Splitting mode before!" << std::endl;
 427        return false;
 428    }
 429
 430    GDCM_NAME_SPACE::FileList *l;
 431
 432    GDCM_NAME_SPACE::SerieHelper *s;
 433    s = GDCM_NAME_SPACE::SerieHelper::New();
 434
 435 //
 436 // Load the gdcm::File* set, according to user's requierements
 437 // ------------------------------------------------------------
 438 //
 439    l = getGdcmFileList();
 440 std::cout << l->size() << " gdcm::File read" << std::endl;
 441 //
 442 // Split the gdcm::File* set, according to user's requierements
 443 // ------------------------------------------------------------
 444 //
 445    if (SplitOnOrientation)
 446    {
 447       s->SetDropDuplicatePositions(false);
 448       xcm = s->SplitOnOrientation(l);
 449    }
 450    else if (SplitOnPosition)
 451    {
 452       s->SetDropDuplicatePositions(true);
 453       xcm = s->SplitOnPosition(l);
 454
 455       // the key of xcm follows lexicographical order
 456       // (that may be different than the 'distance' order)
 457       // we have to reorganize it!
 458
 459       reorgXCoherentFileSetmap(xcm);
 460    }
 461    else if (SplitOnTag)
 462    {
 463          s->SetDropDuplicatePositions(false);
 464
 465          // Crashes if DataElement not found
 466          //std:: cout << GDCM_NAME_SPACE::Global::GetDicts()->GetDefaultPubDict()->GetEntry(groupelem[0],groupelem[1])->GetName() << std::endl;
 467             if ( ! SplitConvertToFloat )
 468                xcm = s->SplitOnTagValue(l, SplitGroup, SplitElem);
 469             else
 470             {
 471                xcm = s->SplitOnTagValueConvertToFloat(l, SplitGroup, SplitElem);
 472             }
 473   }
 474
 475    if (xcm.size() == 0)
 476    {
 477       if(verbose)
 478          std::cout << "Empty XCoherent File Set after 'split' ?!?" << std::endl;
 479       return false;
 480    }
 481    else
 482    {
 483       std::cout << xcm.size() << " XCoherent entries found" << std::endl;
 484    }
 485
 486 // put here, to avoid segfault when unaware user sets SplitOnly to true, and the asks for ImageDataVector
 487 ImageDataVector = new std::vector<vtkImageData*>;
 488
 489 std::cout <<"SplitOnly " << SplitOnly << std::endl;
 490    if(SplitOnly)
 491       return true;
 492 //
 493 //
 494 // ------------------------------------------------------------
 495 //
 496 //   ImageDataVector = new std::vector<vtkImageData*>;
 497
 498    /// \TODO move inside the loop, or be clever using vtk!
 499   // vtkGdcmReader *reader = vtkGdcmReader::New(); // move inside the loop, or be clever using vtk!
 500
 501  // XCoherentFileSetmap map < critère de split, FileList (= std::vector de gdcm::File*) >
 502
 503    for (GDCM_NAME_SPACE::XCoherentFileSetmap::iterator i = xcm.begin();
 504                                                   i != xcm.end();
 505                                                 ++i)
 506    {
 507       vtkGdcmReader *reader = vtkGdcmReader::New(); /// \TODO FIXME : unable to delete!
 508
 509       reader->SetFlipY(FlipY);
 510       // better user SetFileLowerLeft()
 511       /// \TODO : modify vtkGdcmReader !
 512       if (verbose)
 513                std::cout << " --- xCoherentName = [" << (*i).first << "]" << std::endl;
 514
 515       if (SortOnPosition)
 516       {
 517               if (verbose) std::cout << "SortOnPosition" << std::endl;
 518               // (will be IPPSorter, in GDCM2)
 519               s->ImagePositionPatientOrdering((*i).second);
 520               if (verbose) std::cout << "out of SortOnPosition" << std::endl;
 521       }
 522
 523       else if (SortOnOrientation)
 524       {
 525               if (verbose) std::cout << "SortOnOrientation" << std::endl;
 526             /// \TODO (?) SortOnOrientation()
 527
 528             // we still miss an algo to sort an Orientation, given by 6 cosines!
 529             //  Anything like this, in GDCM2?
 530             std::cout << "SortOnOrientation : not so easy - I(mage)O(rientation)P(atient)Sorter still missing! -" << std::endl;
 531             // have a look at SerieHelper::SplitOnOrientation() to have an idea of the mess!
 532
 533             //Better sort on the file name, right now...
 534              s->FileNameOrdering((*i).second);
 535       }
 536
 537       else if (SortOnFileName)
 538       {
 539          if (verbose) std::cout << "SortOnFileName" << std::endl;
 540          if (verbose) std::cout << "taille " << ((*i).second)->size() << std::endl;
 541
 542          s->FileNameOrdering((*i).second);
 543          if (verbose) std::cout << "Out of SortOnFileName" << std::endl;
 544       }
 545
 546 //      else if (SortOnTag)
 547 //      {
 548 //         if (verbose) std::cout << "SortOnTag" << std::endl;
 549 //         printf ("--> %04x %04x\n", SortGroup,SortElem);
 550 //         std::cout << "Sorry, troubles not solved yet; use SortOnUserFunction, right now!" << std::endl;
 551
 552         /*        ==> WARNING : This one has troubles; do NOT use it, right now!
 553         // a pointer to fonction cannot be casted as a pointer to member function!
 554         // Use SortOnUserFunction, instead!
 555
 556         //  if ( SortConvertToFloat )
 557         //    s->SetUserLessThanFunction( reinterpret_cast<bool (*)(gdcm13::File*, gdcm13::File*)>
 558         //                                                         ( &vtkGdcm4DSplitter::CompareOnSortTagConvertToFloat));
 559         //  else
 560         //     s->SetUserLessThanFunction( reinterpret_cast<bool (*)(gdcm13::File*, gdcm13::File*)>
 561                                                                     ( &vtkGdcm4DSplitter::CompareOnSortTag));
 562
 563          // Anything like this, in GDCM2?
 564         //  s->UserOrdering((*i).second);
 565         */
 566
 567 //         //if (verbose) std::cout << "Out of SortOnTag" << std::endl;
 568 //         std::cout << "NO ordering performed  :-( " << std::endl;
 569 //      }
 570
 571       else if (SortOnUserFunction)
 572       {
 573           if (verbose) std::cout << "SortOnUserFunction" << std::endl;
 574           s->SetUserLessThanFunction( UserCompareFunction );
 575          // Anything like this, in GDCM2?
 576          s->UserOrdering((*i).second);
 577          if (verbose) std::cout << "Out of SortOnUserFunction" << std::endl;
 578       }
 579
 580        reader->SetCoherentFileList((*i).second);
 581        reader->Update();
 582
 583        /// \TODO : remove the following
 584        if (verbose) {
 585          // std::cout << "reader->GetOutput() :" << std::endl;
 586          // reader->GetOutput()->PrintSelf(std::cout, vtkIndent(2));
 587        }
 588
 589        ImageDataVector->push_back(reader->GetOutput() );
 590
 591        std::vector<vtkImageData*>::iterator it;
 592        //if (verbose) // JPR
 593
 594        for(it=ImageDataVector->begin(); it!=ImageDataVector->end(); ++it) {
 595         // std::cout << "-in vtkGdcm4DSplitter --------------------------" << std::endl;
 596         // (*it)->PrintSelf(std::cout, vtkIndent(2));
 597        }
 598        //std::cout << std::endl;
 599    }
 600
 601    //reader->Delete();  // \TODO : fix
 602    s->Delete();
 603   // f->Delete();
 604    delete l;
 605
 606    return true;
 607  }
 608
 609
 610 /**
 611  * \brief Load the gdcm::File* set, according to user's requierements
 612  * returns a std::vector of gdcm::File* (gdcm::FileList : actually, a std::vector of gdcm::File*)
 613  */
 614
 615
 616 GDCM_NAME_SPACE::FileList *vtkGdcm4DSplitter::getGdcmFileList()
 617 {
 618
 619   GDCM_NAME_SPACE::File *f;
 620   GDCM_NAME_SPACE::DirListType fileNames;
 621
 622  //
 623  // Fill fileNames with the user supplied file names (in any)
 624  // ------------------------------------------------
 625  //
 626    if (TypeDir == 0 )  // Nothing was set as input...
 627    {
 628       ///\TODO (?) Throw an exception "Set input Directory name(s) / file names  before!"
 629       std::cout << "Set input Directory name(s) / file names  before!" << std::endl;
 630       return false;
 631    }
 632    else if (TypeDir == 1) // A root directory name was set as input
 633    {
 634       GDCM_NAME_SPACE::DirList dirlist(DirName, Recursive); // NO recursive exploration
 635       fileNames = dirlist.GetFilenames(); // all the file names
 636    }
 637
 638    else if (TypeDir == 2) // a std::vector of directory names was set as input
 639    {
 640       int nbDir = VectDirName.size();
 641       GDCM_NAME_SPACE::DirListType tmpFileNames;
 642       for (int iDir=0; iDir<nbDir; iDir++)
 643       {
 644         GDCM_NAME_SPACE::DirList dirlist(VectDirName[iDir], Recursive);
 645         tmpFileNames = dirlist.GetFilenames();
 646         // Concat two std::vector
 647         //vector1.insert( vector1.end(), vector2.begin(), vector2.end() );
 648         fileNames.insert( fileNames.end(), tmpFileNames.begin(), tmpFileNames.end() );
 649       }
 650    }
 651    else if (TypeDir == 3) // a list of files names was set as input
 652    {
 653       fileNames=VectFileName;
 654    }
 655
 656  //
 657  // Fill l with the gdcm::File* corresponding to the files
 658  // --------------------------------------
 659  //
 660
 661    GDCM_NAME_SPACE::FileList *l = new GDCM_NAME_SPACE::FileList; // (set of gdcm::File*)
 662
 663    if (TypeDir == 4) // an already existing std::vector of gdcm::File* was set as input
 664    {
 665       l = VectGdcmFile;
 666    }
 667    else
 668    {
 669       double floatTagvalue;
 670       // Loop on all the gdcm-readable files
 671       for (GDCM_NAME_SPACE::DirListType::iterator it = fileNames.begin();
 672                                                   it != fileNames.end();
 673                                                 ++it)
 674       {
 675          int maxSize  = 0x7fff;         // load Elements of any length
 676          f = GDCM_NAME_SPACE::File::New();
 677          f->SetMaxSizeLoadEntry(maxSize);
 678          f->SetFileName( *it );
 679          if (f->Load())
 680             l->push_back(f);
 681          else
 682             std::cout << " Fail to load [" <<  *it << "]" << std::endl;
 683       }
 684    }
 685    return l;
 686 }
 687
 688
 689
 690  void vtkGdcm4DSplitter::reorgXCoherentFileSetmap(GDCM_NAME_SPACE::XCoherentFileSetmap &xcm)
 691  {
 692  /*
 693  the key of the 'XCoherentFileSetmap', is a std::string, used as if it was found in the Dicom header
 694  Normaly(?), it's suitable for almost everything ...
 695  ... but the 'Image Position Patient'.
 696  We need to order the 'XCoherentFileSetmap' (NOT the content of each XCoherentFileSet!) according to the IPP,
 697  using Jolinda Smith's algorithm.
 698  (we use a subset of the one defined in gdcm::SerieHelper)
 699
 700 */
 701
 702    ELEM e;
 703    std::vector<ELEM> vectElem;
 704
 705 /* just to remember :
 706    typedef struct
 707    {
 708       std::string strIPP;
 709       double dist;
 710       GDCM_NAME_SPACE::File *file;
 711    } ELEM;
 712 */
 713
 714    bool Debug=true;
 715
 716    for (GDCM_NAME_SPACE::XCoherentFileSetmap::iterator i = xcm.begin();
 717                                                        i != xcm.end();
 718                                                      ++i)
 719    {
 720       if (verbose)
 721                std::cout << "--- xCoherentName = [" << (*i).first << "]" << std::endl;
 722
 723       e.strIPP  = (*i).first;
 724       e.file = *(((*i).second)->begin()); // all the gdcm::File of a given xcm item *have* the same IPP; first one is enough
 725       e.dist=0.0;
 726       vectElem.push_back(e);
 727    }
 728    sortVectElem(&vectElem);
 729
 730    // now, create the final std::map !
 731    // final_xcm<to_str(e.dist , xcm[e.strIPP]>
 732    // xcm = final_xcm;
 733    // check what we need to free !
 734
 735    int dist;
 736    char c_dist[100];
 737    std::string str_dist;
 738    int lgr=vectElem.size();
 739    GDCM_NAME_SPACE::XCoherentFileSetmap final_xcm;
 740    for (int i2=0; i2<lgr; i2++)
 741    {
 742       dist =  (vectElem[i2].dist*1000);
 743       sprintf(c_dist,"%010d",dist);
 744       str_dist = c_dist;
 745 /*
 746       std::cout << "dist " << vectElem[i2].dist
 747                 << " str_dist " << str_dist
 748                 << " IPP " << vectElem[i2].strIPP
 749                 << std::endl;
 750
 751 */
 752       final_xcm[str_dist] = xcm[vectElem[i2].strIPP];
 753    }
 754
 755    /// \TODO : check what needs to be cleared // JPR
 756
 757    xcm = final_xcm;
 758
 759  }
 760
 761
 762 bool vtkGdcm4DSplitter::sortVectElem(std::vector<ELEM> *fileList)
 763 {
 764 //based on Jolinda Smith's algorithm
 765 //
 766 // NOTE : if you need to use Jolinda Smith's algorithm, get the one inside gdcm::SerieHelper
 767 // this one is a light version.
 768
 769 //Tags always use the same coordinate system, where "x" is left
 770 //to right, "y" is posterior to anterior, and "z" is foot to head (RAH).
 771
 772    //iop is calculated based on the file file
 773    float cosines[6];
 774    double normal[3];
 775    double ipp[3];
 776    double dist;
 777    double min = 0, max = 0;
 778    bool first = true;
 779
 780    //double ZSpacing; // useless here! // JPR
 781    bool DirectOrder = true; // remove it!
 782
 783   // ZSpacing = -1.0;  // will be updated if process doesn't fail
 784
 785    //std::multimap<double,File *> distmultimap; // JPR
 786    std::multimap<double,ELEM> distmultimap; // JPR
 787
 788    // Use a multimap to sort the distances from 0,0,0
 789    //for ( FileList::const_iterator // JPR
 790    for ( std::vector<ELEM>::iterator   // JPR
 791          it = fileList->begin();
 792          it != fileList->end(); ++it )
 793    {
 794       //gdcmDebugMacro("deal with " << (*it)->file->GetFileName() );
 795       if ( first )
 796       {
 797          (*it).file->GetImageOrientationPatient( cosines );
 798
 799    // The "Image Orientation Patient" tag gives the direction cosines
 800    // for the rows and columns for the three axes defined above.
 801    // Typical axial slices will have a value 1/0/0/0/1/0:
 802    // rows increase from left to right,
 803    // columns increase from posterior to anterior. This is your everyday
 804    // "looking up from the bottom of the head with the eyeballs up" image.
 805
 806    // The "Image Position Patient" tag gives the coordinates of the first
 807    // voxel in the image in the "RAH" coordinate system, relative to some
 808    // origin.
 809
 810    // First, calculate the slice normal from IOP :
 811
 812          // You only have to do this once for all slices in the volume. Next,
 813          // for each slice, calculate the distance along the slice normal
 814          // using the IPP ("Image Position Patient") tag.
 815          // ("dist" is initialized to zero before reading the first slice) :
 816
 817          normal[0] = cosines[1]*cosines[5] - cosines[2]*cosines[4];
 818          normal[1] = cosines[2]*cosines[3] - cosines[0]*cosines[5];
 819          normal[2] = cosines[0]*cosines[4] - cosines[1]*cosines[3];
 820
 821    // For each slice (here : the first), calculate the distance along
 822    // the slice normal using the IPP tag
 823
 824          ipp[0] = (*it).file->GetXOrigin();
 825          ipp[1] = (*it).file->GetYOrigin();
 826          ipp[2] = (*it).file->GetZOrigin();
 827
 828          dist = 0;
 829          for ( int i = 0; i < 3; ++i )
 830          {
 831             dist += normal[i]*ipp[i];
 832          }
 833
 834          //gdcmDebugMacro("dist : " << dist);
 835          distmultimap.insert(std::pair<const double,ELEM>(dist, *it));
 836
 837          max = min = dist;
 838          first = false;
 839       }
 840       else
 841       {
 842    // Next, for each slice, calculate the distance along the slice normal
 843    // using the IPP tag
 844          ipp[0] = (*it).file->GetXOrigin();
 845          ipp[1] = (*it).file->GetYOrigin();
 846          ipp[2] = (*it).file->GetZOrigin();
 847
 848          dist = 0;
 849          for ( int i = 0; i < 3; ++i )
 850          {
 851             dist += normal[i]*ipp[i];
 852          }
 853
 854          (*it).dist = dist; // JPR
 855
 856          distmultimap.insert(std::pair<const double,ELEM>(dist, *it));
 857          //gdcmDebugMacro("dist : " << dist);
 858          min = (min < dist) ? min : dist;
 859          max = (max > dist) ? max : dist;
 860       }
 861    }
 862
 863   // gdcmDebugMacro("After parsing vector, nb of elements : " << fileList->size() );
 864
 865 /*Useless here. // JPR
 866
 867    // Find out if min/max are coherent
 868    if ( min == max )
 869    {
 870      gdcmWarningMacro("Looks like all images have the exact same image position. "
 871                       << "No PositionPatientOrdering sort performed. "
 872                       << "No 'ZSpacing' calculated! ");
 873      return false;
 874    }
 875 */
 876
 877 /* Useless here, 'split' already done. // JPR
 878
 879    // Check to see if image shares a common position
 880    bool ok = true;
 881    for (std::multimap<double, File *>::iterator it2 = distmultimap.begin();
 882         it2 != distmultimap.end();
 883         ++it2)
 884    {
 885
 886       gdcmDebugMacro("Check if image shares a common position : " << ((*it2).second).file->GetFileName() );
 887
 888       if (distmultimap.count((*it2).first) != 1)
 889       {
 890          gdcmWarningMacro("File: ["
 891               << ((*it2).second->GetFileName())
 892               << "] : more than ONE file at distance: '"
 893               << (*it2).first
 894               << " (position is not unique!) "
 895               << "No PositionPatientOrdering sort performed. "
 896               << "No 'ZSpacing' calculated! ");
 897
 898          ok = false;
 899       }
 900    }
 901    if (!ok)
 902    {
 903       if (! DropDuplicatePositions)
 904          return false;
 905    }
 906
 907 */
 908
 909 // Now, we can calculate Z Spacing as the difference
 910 // between the "dist" values for the first two slices.
 911
 912 // The following (un)-commented out code is let here
 913 // to be re-used by whomsoever is interested...
 914
 915     //std::multimap<double, ELEM>::iterator it5 = distmultimap.begin();
 916     //double d1 = (*it5).first;
 917     //it5++;
 918     //double d2 = (*it5).first;
 919     //ZSpacing = d1-d2;
 920     //if (ZSpacing < 0.0)
 921     //   ZSpacing = - ZSpacing;
 922
 923    fileList->clear();  // doesn't delete list elements, only nodes
 924
 925 // Acording to user requierement, we sort direct order or reverse order.
 926    if (DirectOrder)
 927    {
 928       for (std::multimap<double, ELEM>::iterator it3 = distmultimap.begin();
 929            it3 != distmultimap.end();
 930            ++it3)
 931       {
 932          fileList->push_back( (*it3).second );
 933
 934 /*       useless here! // JPR
 935
 936          if (DropDuplicatePositions)
 937          {
 938             // ImagePositionPatientOrdering  wrong duplicates are found ???
 939             // --> fixed. See comment
 940
 941             it3 =  distmultimap.upper_bound((*it3).first); // skip all duplicates
 942            // the upper_bound function increments the iterator to the next non-duplicate entry
 943            // The for loop iteration also increments the iterator, which causes the code to skip every other image
 944            // --> decrement the iterator after the upper_bound function call
 945             it3--;
 946             if (it3 == distmultimap.end() )  // if last image, stop iterate
 947                break;
 948          }
 949 */
 950       }
 951    }
 952    else // user asked for reverse order
 953    {
 954       std::multimap<double, ELEM>::const_iterator it4;
 955       it4 = distmultimap.end();
 956       do
 957       {
 958          it4--;
 959          fileList->push_back( (*it4).second );
 960
 961 /* useless here // JPR
 962
 963          if (DropDuplicatePositions)  // skip all duplicates
 964          {
 965             // lower_bound finds the next element that is
 966             // less than or *equal to* the current value!
 967             //it4 =  distmultimap.lower_bound((*it4).first);
 968
 969            // David Feng's fix
 970            std::multimap<double, ELEM>::const_iterator itPrev = it4;
 971            while (itPrev->first == it4->first)
 972               --itPrev;
 973            it4 = itPrev;
 974
 975            if (it4 == distmultimap.begin() ) // if first image, stop iterate
 976                break;
 977          }
 978 */
 979       } while (it4 != distmultimap.begin() );
 980    }
 981
 982    distmultimap.clear();
 983
 984    return true;
 985
 986 }