vtk/vtkGdcm4DSplitter.cxx

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