vtk/vtkGdcm4DSplitter.cxx

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