vtk/vtkGdcm4DSplitter.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: vtkGdcm4DSplitter.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2011/04/20 15:03:54 $
   7   Version:   $Revision: 1.13 $
   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 NULL;
 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 std::cout <<"SplitOnly " << SplitOnly << std::endl;
 487    if(SplitOnly)
 488       return true;
 489 //
 490 //
 491 // ------------------------------------------------------------
 492 //
 493    ImageDataVector = new std::vector<vtkImageData*>;
 494
 495    /// \TODO move inside the loop, or be clever using vtk!
 496   // vtkGdcmReader *reader = vtkGdcmReader::New(); // move inside the loop, or be clever using vtk!
 497
 498  // XCoherentFileSetmap map < critère de split, FileList (= std::vector de gdcm::File*) >
 499
 500    for (GDCM_NAME_SPACE::XCoherentFileSetmap::iterator i = xcm.begin();
 501                                                   i != xcm.end();
 502                                                 ++i)
 503    {
 504       vtkGdcmReader *reader = vtkGdcmReader::New(); /// \TODO FIXME : unable to delete!
 505
 506       reader->SetFlipY(FlipY);
 507       // better user SetFileLowerLeft()
 508       /// \TODO : modify vtkGdcmReader !
 509       if (verbose)
 510                std::cout << " --- xCoherentName = [" << (*i).first << "]" << std::endl;
 511
 512       if (SortOnPosition)
 513       {
 514               if (verbose) std::cout << "SortOnPosition" << std::endl;
 515               // (will be IPPSorter, in GDCM2)
 516               s->ImagePositionPatientOrdering((*i).second);
 517               if (verbose) std::cout << "out of SortOnPosition" << std::endl;
 518       }
 519
 520       else if (SortOnOrientation)
 521       {
 522               if (verbose) std::cout << "SortOnOrientation" << std::endl;
 523             /// \TODO (?) SortOnOrientation()
 524
 525             // we still miss an algo to sort an Orientation, given by 6 cosines!
 526             //  Anything like this, in GDCM2?
 527             std::cout << "SortOnOrientation : not so easy - I(mage)O(rientation)P(atient)Sorter still missing! -" << std::endl;
 528             // have a look at SerieHelper::SplitOnOrientation() to have an idea of the mess!
 529
 530             //Better sort on the file name, right now...
 531              s->FileNameOrdering((*i).second);
 532       }
 533
 534       else if (SortOnFileName)
 535       {
 536          if (verbose) std::cout << "SortOnFileName" << std::endl;
 537          if (verbose) std::cout << "taille " << ((*i).second)->size() << std::endl;
 538
 539          s->FileNameOrdering((*i).second);
 540          if (verbose) std::cout << "Out of SortOnFileName" << std::endl;
 541       }
 542
 543 //      else if (SortOnTag)
 544 //      {
 545 //         if (verbose) std::cout << "SortOnTag" << std::endl;
 546 //         printf ("--> %04x %04x\n", SortGroup,SortElem);
 547 //         std::cout << "Sorry, troubles not solved yet; use SortOnUserFunction, right now!" << std::endl;
 548
 549         /*        ==> WARNING : This one has troubles; do NOT use it, right now!
 550         // a pointer to fonction cannot be casted as a pointer to member function!
 551         // Use SortOnUserFunction, instead!
 552
 553         //  if ( SortConvertToFloat )
 554         //    s->SetUserLessThanFunction( reinterpret_cast<bool (*)(gdcm13::File*, gdcm13::File*)>
 555         //                                                         ( &vtkGdcm4DSplitter::CompareOnSortTagConvertToFloat));
 556         //  else
 557         //     s->SetUserLessThanFunction( reinterpret_cast<bool (*)(gdcm13::File*, gdcm13::File*)>
 558                                                                     ( &vtkGdcm4DSplitter::CompareOnSortTag));
 559
 560          // Anything like this, in GDCM2?
 561         //  s->UserOrdering((*i).second);
 562         */
 563
 564 //         //if (verbose) std::cout << "Out of SortOnTag" << std::endl;
 565 //         std::cout << "NO ordering performed  :-( " << std::endl;
 566 //      }
 567
 568       else if (SortOnUserFunction)
 569       {
 570           if (verbose) std::cout << "SortOnUserFunction" << std::endl;
 571           s->SetUserLessThanFunction( UserCompareFunction );
 572          // Anything like this, in GDCM2?
 573          s->UserOrdering((*i).second);
 574          if (verbose) std::cout << "Out of SortOnUserFunction" << std::endl;
 575       }
 576
 577        reader->SetCoherentFileList((*i).second);
 578        reader->Update();
 579
 580        /// \TODO : remove the following
 581        if (verbose) {
 582          // std::cout << "reader->GetOutput() :" << std::endl;
 583          // reader->GetOutput()->PrintSelf(std::cout, vtkIndent(2));
 584        }
 585
 586        ImageDataVector->push_back(reader->GetOutput() );
 587
 588        std::vector<vtkImageData*>::iterator it;
 589        //if (verbose) // JPR
 590
 591        for(it=ImageDataVector->begin(); it!=ImageDataVector->end(); ++it) {
 592         // std::cout << "-in vtkGdcm4DSplitter --------------------------" << std::endl;
 593         // (*it)->PrintSelf(std::cout, vtkIndent(2));
 594        }
 595        //std::cout << std::endl;
 596    }
 597
 598    //reader->Delete();  // \TODO : fix
 599    s->Delete();
 600   // f->Delete();
 601    delete l;
 602
 603    return true;
 604  }
 605
 606
 607 /**
 608  * \brief Load the gdcm::File* set, according to user's requierements
 609  * returns a std::vector of gdcm::File* (gdcm::FileList : actually, a std::vector of gdcm::File*)
 610  */
 611
 612
 613 GDCM_NAME_SPACE::FileList *vtkGdcm4DSplitter::getGdcmFileList()
 614 {
 615
 616   GDCM_NAME_SPACE::File *f;
 617   GDCM_NAME_SPACE::DirListType fileNames;
 618
 619  //
 620  // Fill fileNames with the user supplied file names (in any)
 621  // ------------------------------------------------
 622  //
 623    if (TypeDir == 0 )  // Nothing was set as input...
 624    {
 625       ///\TODO (?) Throw an exception "Set input Directory name(s) / file names  before!"
 626       std::cout << "Set input Directory name(s) / file names  before!" << std::endl;
 627       return false;
 628    }
 629    else if (TypeDir == 1) // A root directory name was set as input
 630    {
 631       GDCM_NAME_SPACE::DirList dirlist(DirName, Recursive); // NO recursive exploration
 632       fileNames = dirlist.GetFilenames(); // all the file names
 633    }
 634
 635    else if (TypeDir == 2) // a std::vector of directory names was set as input
 636    {
 637       int nbDir = VectDirName.size();
 638       GDCM_NAME_SPACE::DirListType tmpFileNames;
 639       for (int iDir=0; iDir<nbDir; iDir++)
 640       {
 641         GDCM_NAME_SPACE::DirList dirlist(VectDirName[iDir], Recursive);
 642         tmpFileNames = dirlist.GetFilenames();
 643         // Concat two std::vector
 644         //vector1.insert( vector1.end(), vector2.begin(), vector2.end() );
 645         fileNames.insert( fileNames.end(), tmpFileNames.begin(), tmpFileNames.end() );
 646       }
 647    }
 648    else if (TypeDir == 3) // a list of files names was set as input
 649    {
 650       fileNames=VectFileName;
 651    }
 652
 653  //
 654  // Fill l with the gdcm::File* corresponding to the files
 655  // --------------------------------------
 656  //
 657
 658    GDCM_NAME_SPACE::FileList *l = new GDCM_NAME_SPACE::FileList; // (set of gdcm::File*)
 659
 660    if (TypeDir == 4) // an already existing std::vector of gdcm::File* was set as input
 661    {
 662       l = VectGdcmFile;
 663    }
 664    else
 665    {
 666       double floatTagvalue;
 667       // Loop on all the gdcm-readable files
 668       for (GDCM_NAME_SPACE::DirListType::iterator it = fileNames.begin();
 669                                                   it != fileNames.end();
 670                                                 ++it)
 671       {
 672          int maxSize  = 0x7fff;         // load Elements of any length
 673          f = GDCM_NAME_SPACE::File::New();
 674          f->SetMaxSizeLoadEntry(maxSize);
 675          f->SetFileName( *it );
 676          if (f->Load())
 677             l->push_back(f);
 678          else
 679             std::cout << " Fail to load [" <<  *it << "]" << std::endl;
 680       }
 681    }
 682    return l;
 683 }
 684
 685
 686
 687  void vtkGdcm4DSplitter::reorgXCoherentFileSetmap(GDCM_NAME_SPACE::XCoherentFileSetmap &xcm)
 688  {
 689  /*
 690  the key of the 'XCoherentFileSetmap', is a std::string, used as if it was found in the Dicom header
 691  Normaly(?), it's suitable for almost everything ...
 692  ... but the 'Image Position Patient'.
 693  We need to order the 'XCoherentFileSetmap' (NOT the content of each XCoherentFileSet!) according to the IPP,
 694  using Jolinda Smith's algorithm.
 695  (we use a subset of the one defined in gdcm::SerieHelper)
 696
 697 */
 698
 699    ELEM e;
 700    std::vector<ELEM> vectElem;
 701
 702 /* just to remember :
 703    typedef struct
 704    {
 705       std::string strIPP;
 706       double dist;
 707       GDCM_NAME_SPACE::File *file;
 708    } ELEM;
 709 */
 710
 711    bool Debug=true;
 712
 713    for (GDCM_NAME_SPACE::XCoherentFileSetmap::iterator i = xcm.begin();
 714                                                        i != xcm.end();
 715                                                      ++i)
 716    {
 717       if (verbose)
 718                std::cout << "--- xCoherentName = [" << (*i).first << "]" << std::endl;
 719
 720       e.strIPP  = (*i).first;
 721       e.file = *(((*i).second)->begin()); // all the gdcm::File of a given xcm item *have* the same IPP; first one is enough
 722       e.dist=0.0;
 723       vectElem.push_back(e);
 724    }
 725    sortVectElem(&vectElem);
 726
 727    // now, create the final std::map !
 728    // final_xcm<to_str(e.dist , xcm[e.strIPP]>
 729    // xcm = final_xcm;
 730    // check what we need to free !
 731
 732    int dist;
 733    char c_dist[100];
 734    std::string str_dist;
 735    int lgr=vectElem.size();
 736    GDCM_NAME_SPACE::XCoherentFileSetmap final_xcm;
 737    for (int i2=0; i2<lgr; i2++)
 738    {
 739       dist =  (vectElem[i2].dist*1000);
 740       sprintf(c_dist,"%010d",dist);
 741       str_dist = c_dist;
 742 /*
 743       std::cout << "dist " << vectElem[i2].dist
 744                 << " str_dist " << str_dist
 745                 << " IPP " << vectElem[i2].strIPP
 746                 << std::endl;
 747
 748 */
 749       final_xcm[str_dist] = xcm[vectElem[i2].strIPP];
 750    }
 751
 752    /// \TODO : check what needs to be cleared // JPR
 753
 754    xcm = final_xcm;
 755
 756  }
 757
 758
 759 bool vtkGdcm4DSplitter::sortVectElem(std::vector<ELEM> *fileList)
 760 {
 761 //based on Jolinda Smith's algorithm
 762 //
 763 // NOTE : if you need to use Jolinda Smith's algorithm, get the one inside gdcm::SerieHelper
 764 // this one is a light version.
 765
 766 //Tags always use the same coordinate system, where "x" is left
 767 //to right, "y" is posterior to anterior, and "z" is foot to head (RAH).
 768
 769    //iop is calculated based on the file file
 770    float cosines[6];
 771    double normal[3];
 772    double ipp[3];
 773    double dist;
 774    double min = 0, max = 0;
 775    bool first = true;
 776
 777    //double ZSpacing; // useless here! // JPR
 778    bool DirectOrder = true; // remove it!
 779
 780   // ZSpacing = -1.0;  // will be updated if process doesn't fail
 781
 782    //std::multimap<double,File *> distmultimap; // JPR
 783    std::multimap<double,ELEM> distmultimap; // JPR
 784
 785    // Use a multimap to sort the distances from 0,0,0
 786    //for ( FileList::const_iterator // JPR
 787    for ( std::vector<ELEM>::iterator   // JPR
 788          it = fileList->begin();
 789          it != fileList->end(); ++it )
 790    {
 791       //gdcmDebugMacro("deal with " << (*it)->file->GetFileName() );
 792       if ( first )
 793       {
 794          (*it).file->GetImageOrientationPatient( cosines );
 795
 796    // The "Image Orientation Patient" tag gives the direction cosines
 797    // for the rows and columns for the three axes defined above.
 798    // Typical axial slices will have a value 1/0/0/0/1/0:
 799    // rows increase from left to right,
 800    // columns increase from posterior to anterior. This is your everyday
 801    // "looking up from the bottom of the head with the eyeballs up" image.
 802
 803    // The "Image Position Patient" tag gives the coordinates of the first
 804    // voxel in the image in the "RAH" coordinate system, relative to some
 805    // origin.
 806
 807    // First, calculate the slice normal from IOP :
 808
 809          // You only have to do this once for all slices in the volume. Next,
 810          // for each slice, calculate the distance along the slice normal
 811          // using the IPP ("Image Position Patient") tag.
 812          // ("dist" is initialized to zero before reading the first slice) :
 813
 814          normal[0] = cosines[1]*cosines[5] - cosines[2]*cosines[4];
 815          normal[1] = cosines[2]*cosines[3] - cosines[0]*cosines[5];
 816          normal[2] = cosines[0]*cosines[4] - cosines[1]*cosines[3];
 817
 818    // For each slice (here : the first), calculate the distance along
 819    // the slice normal using the IPP tag
 820
 821          ipp[0] = (*it).file->GetXOrigin();
 822          ipp[1] = (*it).file->GetYOrigin();
 823          ipp[2] = (*it).file->GetZOrigin();
 824
 825          dist = 0;
 826          for ( int i = 0; i < 3; ++i )
 827          {
 828             dist += normal[i]*ipp[i];
 829          }
 830
 831          //gdcmDebugMacro("dist : " << dist);
 832          distmultimap.insert(std::pair<const double,ELEM>(dist, *it));
 833
 834          max = min = dist;
 835          first = false;
 836       }
 837       else
 838       {
 839    // Next, for each slice, calculate the distance along the slice normal
 840    // using the IPP tag
 841          ipp[0] = (*it).file->GetXOrigin();
 842          ipp[1] = (*it).file->GetYOrigin();
 843          ipp[2] = (*it).file->GetZOrigin();
 844
 845          dist = 0;
 846          for ( int i = 0; i < 3; ++i )
 847          {
 848             dist += normal[i]*ipp[i];
 849          }
 850
 851          (*it).dist = dist; // JPR
 852
 853          distmultimap.insert(std::pair<const double,ELEM>(dist, *it));
 854          //gdcmDebugMacro("dist : " << dist);
 855          min = (min < dist) ? min : dist;
 856          max = (max > dist) ? max : dist;
 857       }
 858    }
 859
 860   // gdcmDebugMacro("After parsing vector, nb of elements : " << fileList->size() );
 861
 862 /*Useless here. // JPR
 863
 864    // Find out if min/max are coherent
 865    if ( min == max )
 866    {
 867      gdcmWarningMacro("Looks like all images have the exact same image position. "
 868                       << "No PositionPatientOrdering sort performed. "
 869                       << "No 'ZSpacing' calculated! ");
 870      return false;
 871    }
 872 */
 873
 874 /* Useless here, 'split' already done. // JPR
 875
 876    // Check to see if image shares a common position
 877    bool ok = true;
 878    for (std::multimap<double, File *>::iterator it2 = distmultimap.begin();
 879         it2 != distmultimap.end();
 880         ++it2)
 881    {
 882
 883       gdcmDebugMacro("Check if image shares a common position : " << ((*it2).second).file->GetFileName() );
 884
 885       if (distmultimap.count((*it2).first) != 1)
 886       {
 887          gdcmWarningMacro("File: ["
 888               << ((*it2).second->GetFileName())
 889               << "] : more than ONE file at distance: '"
 890               << (*it2).first
 891               << " (position is not unique!) "
 892               << "No PositionPatientOrdering sort performed. "
 893               << "No 'ZSpacing' calculated! ");
 894
 895          ok = false;
 896       }
 897    }
 898    if (!ok)
 899    {
 900       if (! DropDuplicatePositions)
 901          return false;
 902    }
 903
 904 */
 905
 906 // Now, we can calculate Z Spacing as the difference
 907 // between the "dist" values for the first two slices.
 908
 909 // The following (un)-commented out code is let here
 910 // to be re-used by whomsoever is interested...
 911
 912     //std::multimap<double, ELEM>::iterator it5 = distmultimap.begin();
 913     //double d1 = (*it5).first;
 914     //it5++;
 915     //double d2 = (*it5).first;
 916     //ZSpacing = d1-d2;
 917     //if (ZSpacing < 0.0)
 918     //   ZSpacing = - ZSpacing;
 919
 920    fileList->clear();  // doesn't delete list elements, only nodes
 921
 922 // Acording to user requierement, we sort direct order or reverse order.
 923    if (DirectOrder)
 924    {
 925       for (std::multimap<double, ELEM>::iterator it3 = distmultimap.begin();
 926            it3 != distmultimap.end();
 927            ++it3)
 928       {
 929          fileList->push_back( (*it3).second );
 930
 931 /*       useless here! // JPR
 932
 933          if (DropDuplicatePositions)
 934          {
 935             // ImagePositionPatientOrdering  wrong duplicates are found ???
 936             // --> fixed. See comment
 937
 938             it3 =  distmultimap.upper_bound((*it3).first); // skip all duplicates
 939            // the upper_bound function increments the iterator to the next non-duplicate entry
 940            // The for loop iteration also increments the iterator, which causes the code to skip every other image
 941            // --> decrement the iterator after the upper_bound function call
 942             it3--;
 943             if (it3 == distmultimap.end() )  // if last image, stop iterate
 944                break;
 945          }
 946 */
 947       }
 948    }
 949    else // user asked for reverse order
 950    {
 951       std::multimap<double, ELEM>::const_iterator it4;
 952       it4 = distmultimap.end();
 953       do
 954       {
 955          it4--;
 956          fileList->push_back( (*it4).second );
 957
 958 /* useless here // JPR
 959
 960          if (DropDuplicatePositions)  // skip all duplicates
 961          {
 962             // lower_bound finds the next element that is
 963             // less than or *equal to* the current value!
 964             //it4 =  distmultimap.lower_bound((*it4).first);
 965
 966            // David Feng's fix
 967            std::multimap<double, ELEM>::const_iterator itPrev = it4;
 968            while (itPrev->first == it4->first)
 969               --itPrev;
 970            it4 = itPrev;
 971
 972            if (it4 == distmultimap.begin() ) // if first image, stop iterate
 973                break;
 974          }
 975 */
 976       } while (it4 != distmultimap.begin() );
 977    }
 978
 979    distmultimap.clear();
 980
 981    return true;
 982
 983 }