vtk/vtkGdcm4DSplitter.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: vtkGdcm4DSplitter.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2011/04/11 11:28:31 $
   7   Version:   $Revision: 1.9 $
   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
  42 it 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 hom
  46 (despite its name, it works on 3D or 2D+T images too)
  47
  48 User will have to specify some points
  49
  50 . Choose input data
  51 -------------------
  52
  53 - a single directory
  54        bool setDirName(std::string &dirName);
  55 - a list of directories
  56        bool setVectDirName(std::vector<std::string> &vectDirName);
  57 - a list of files
  58        bool setVectFileName(std::vector<std::string> &vectFileName);
  59
  60 - Recursive directory exploration
  61        void setRecursive(bool recursive);
  62
  63 . Choose 'split' criterion :
  64 ---------------------------
  65
  66  - ImagePositionPatient
  67         void setSplitOnPosition();
  68  - ImageOrientationPatient
  69         void setSplitOnOrientation();
  70  - User choosen tag
  71         void setSplitOnTag(unsigned short splitGroup, unsigned short splitElem);
  72         void setSplitConvertToFloat(bool conv);
  73  - UserDefined Function
  74         void setSortOnUserFunction (FoncComp f);
  75
  76  for 'true' 3D image sets :
  77    - if you want to get a single 3D vtkImageData, use SplitOnOrientation -i.e. no split-
  78    - if you want to get a vector of 2D vtkImageData, use SplitOnPosition  -i.e. one slice in each 'XCoherent filesite'-
  79
  80  for 'true' 4D multi-orientation image sets (i.e. a stack of axial + sagital + coronal images, at different instants ...)
  81    --> this is 5D, right?
  82    Nothing done, yet.
  83
  84
  85 . Choose 'sort' criterion :
  86 --------------------------
  87
  88  - ImagePositionPatient
  89         void setSortOnPosition();
  90  - ImageOrientationPatient
  91        ==> Only in your dreams!
  92        ==> or, please, write a IOP sorter ...
  93  - UserDefined Function
  94         void setSortOnUserFunction (FoncComp f);
  95  - File name
  96         void setSortOnFileName()
  97
  98 . Execute :
  99 -----------
 100         bool Go();
 101
 102 . Get the result
 103 ----------------
 104
 105 // -a single vtkImageData:
 106 //        vtkImageData *GetImageData();
 107 - a vector of vtkImageData
 108         std::vector<vtkImageData*> *GetImageDataVector();
 109
 110   ===================================================================== */
 111
 112 #include "gdcmSerieHelper.h"
 113
 114 #include "vtkGdcmReader.h"
 115 #include "vtkGdcm4DSplitter.h"
 116 #include <algorithm>
 117 #include "gdcmSerieHelper.h" // for ImagePositionPatientOrdering()
 118 #include <stdlib.h> // for atof
 119
 120  vtkGdcm4DSplitter::vtkGdcm4DSplitter() :
 121                  SplitOnPosition(false), SplitOnOrientation(false), SplitOnTag(false),
 122                  SplitGroup(0), SplitElem(0),
 123
 124                  SortOnPosition(false),  SortOnOrientation(false),  SortOnTag(false),
 125                  SortGroup(0),  SortElem(0), SortConvertToFloat(false),
 126
 127                  Recursive(false), TypeDir(0),
 128                  verbose(false)
 129  {
 130
 131  }
 132
 133  std::vector<vtkImageData*> * vtkGdcm4DSplitter::GetImageDataVector()
 134  {
 135 /*
 136  if (verbose) std::cout << "GetImageDataVector : TypeResult " << TypeResult << std::endl;
 137     if (TypeResult == 2)
 138        return ImageDataVector;
 139     else
 140       if (TypeResult == 1)
 141       {
 142          std::vector<vtkImageData*> *t = new std::vector<vtkImageData*>;
 143          t->push_back( ImageData );
 144          return t;
 145       }
 146       else
 147          return (std::vector<vtkImageData*>*) NULL;
 148 */
 149      return ImageDataVector;
 150  }
 151
 152  vtkImageData *vtkGdcm4DSplitter::GetImageData()
 153  {
 154  /*
 155   if (verbose) std::cout << "GetImageData : TypeResult " << TypeResult << std::endl;
 156     if (TypeResult == 1)
 157        return ImageData;
 158     else
 159       if (TypeResult == 2)
 160       {
 161          return (*ImageDataVector)[0];
 162       }
 163       else
 164          return (vtkImageData*) NULL;
 165 */
 166    return (*ImageDataVector)[0];
 167  }
 168
 169  bool vtkGdcm4DSplitter::setDirName(std::string &dirName)
 170  {
 171     if ( ! GDCM_NAME_SPACE::DirList::IsDirectory(dirName) )
 172     {
 173        std::cout << "[" << dirName << "] is NOT a directory" << std::endl;
 174        return false;
 175     }
 176     DirName = dirName;
 177     TypeDir=1;
 178     return true;
 179  }
 180
 181  bool vtkGdcm4DSplitter::setVectDirName(std::vector<std::string> &vectDirName)
 182  {
 183     int nbDir = vectDirName.size();
 184     for (int iDir=0; iDir<nbDir; iDir++)
 185     {
 186        if ( ! GDCM_NAME_SPACE::DirList::IsDirectory(vectDirName[iDir]) )
 187        {
 188           std::cout << "[" << vectDirName[iDir] << "] is NOT a directory" << std::endl;
 189           return false;
 190        }
 191     }
 192
 193     VectDirName = vectDirName;
 194     TypeDir=2;
 195     return true;
 196  }
 197
 198  bool vtkGdcm4DSplitter::setVectFileName(std::vector<std::string> &vectFileName)
 199  {
 200     if ( vectFileName.size() == 0)
 201     {
 202           std::cout << "[ vectFileName ] : empty list" << std::endl;
 203           return false;
 204     }
 205     VectFileName = vectFileName;
 206     TypeDir=3;
 207     return true;
 208  }
 209
 210  bool vtkGdcm4DSplitter::CompareOnSortTagConvertToFloat(GDCM_NAME_SPACE::File *file1, GDCM_NAME_SPACE::File *file2)
 211  {
 212   /* if (verbose) printf ("%04x %04x\n", this->SortGroup,this->SortElem);
 213      if (verbose) std :: cout << file1->GetEntryString(SortGroup,SortElem).c_str() << " : "
 214                             << atof(file1->GetEntryString(SortGroup,SortElem).c_str())
 215                             << std::endl;
 216 */
 217 //   return atof(file1->GetEntryString(vtkGdcm4DSplitter::SortGroup,vtkGdcm4DSplitter::SortElem).c_str()) < atof(file2->GetEntryString(vtkGdcm4DSplitter::SortGroup,vtkGdcm4DSplitter::SortElem).c_str());
 218    return atof(file1->GetEntryString(SortGroup,SortElem).c_str()) < atof(file2->GetEntryString(SortGroup,SortElem).c_str());
 219  }
 220
 221  bool vtkGdcm4DSplitter::CompareOnSortTag(GDCM_NAME_SPACE::File *file1, GDCM_NAME_SPACE::File *file2)
 222  {
 223    return file1->GetEntryString(vtkGdcm4DSplitter::SortGroup,vtkGdcm4DSplitter::SortElem) < file2->GetEntryString(vtkGdcm4DSplitter::SortGroup,vtkGdcm4DSplitter::SortElem);
 224  }
 225
 226
 227  bool vtkGdcm4DSplitter::Go()
 228  {
 229    if (!SplitOnPosition && !SplitOnOrientation && !SplitOnTag )
 230    {
 231        ///\TODO (?) Throw an exception "Choose Splitting mode before!"
 232        std::cout << "Choose Splitting mode before!" << std::endl;
 233        return false;
 234    }
 235
 236    /// How To :
 237    /*
 238    entree nom de directory / Vecteur de noms?
 239    recursif oui/non
 240    recuperer la liste des gdcm::File*
 241    passer a SerieHelper (pas de check du Serie UID)
 242    set critere de split
 243
 244    trier chaque Coherent file set
 245    passer chacun a un vtkGcdmReader
 246    retourner le (vecteur de) vtkImageData
 247    */
 248
 249    GDCM_NAME_SPACE::SerieHelper *s;
 250    s = GDCM_NAME_SPACE::SerieHelper::New();
 251
 252    GDCM_NAME_SPACE::File *f;
 253    GDCM_NAME_SPACE::DirListType fileNames;
 254
 255    if (TypeDir == 0 )
 256    {
 257       ///\TODO (?) Throw an exception "Set input Directory name(s) / file names  before!"
 258       std::cout << "Set input Directory name(s) / file names  before!" << std::endl;
 259       return false;
 260    }
 261    else if (TypeDir == 1)
 262    {
 263       GDCM_NAME_SPACE::DirList dirlist(DirName, Recursive); // NO recursive exploration
 264       fileNames = dirlist.GetFilenames(); // all the file names
 265    }
 266
 267    else if (TypeDir == 2)
 268    {
 269       int nbDir = VectDirName.size();
 270       GDCM_NAME_SPACE::DirListType tmpFileNames;
 271       for (int iDir=0; iDir<nbDir; iDir++)
 272       {
 273         GDCM_NAME_SPACE::DirList dirlist(VectDirName[iDir], Recursive);
 274         tmpFileNames = dirlist.GetFilenames();
 275         // Concat two std::vector
 276         //vector1.insert( vector1.end(), vector2.begin(), vector2.end() );
 277        fileNames.insert( fileNames.end(), tmpFileNames.begin(), tmpFileNames.end() );
 278       }
 279    }
 280    else if (TypeDir == 3)
 281    {
 282       fileNames=VectFileName;
 283    }
 284
 285    GDCM_NAME_SPACE::FileList *l = new GDCM_NAME_SPACE::FileList; // (set of gdcm::File)
 286    double floatTagvalue;
 287    // Loop on all the gdcm-readable files
 288    for (GDCM_NAME_SPACE::DirListType::iterator it = fileNames.begin();
 289                                     it != fileNames.end();
 290                                   ++it)
 291    {
 292       int maxSize  = 0x7fff;         // load Elements of any length
 293       f = GDCM_NAME_SPACE::File::New();
 294       f->SetMaxSizeLoadEntry(maxSize);
 295       f->SetFileName( *it );
 296       if (f->Load())
 297          l->push_back(f);
 298       else
 299          std::cout << " Fail to load [" <<  *it << "]" << std::endl;
 300    }
 301
 302    GDCM_NAME_SPACE::XCoherentFileSetmap xcm;
 303
 304    if (SplitOnOrientation)
 305    {
 306       s->SetDropDuplicatePositions(false);
 307       xcm = s->SplitOnOrientation(l);
 308    }
 309    else if (SplitOnPosition)
 310    {
 311       s->SetDropDuplicatePositions(true);
 312       xcm = s->SplitOnPosition(l);
 313       // reorg the std::map xcm according to position // JPR
 314       // the key of xcm follows lexicographical order
 315       // (that may be different than the 'distance' order)
 316       // we have to reorganize it!
 317       reorgXCoherentFileSetmap(xcm);
 318    }
 319    else if (SplitOnTag)
 320    {
 321          s->SetDropDuplicatePositions(false);
 322
 323          // Crashes if DataElement not found
 324          //std:: cout << GDCM_NAME_SPACE::Global::GetDicts()->GetDefaultPubDict()->GetEntry(groupelem[0],groupelem[1])->GetName() << std::endl;
 325             if ( ! SplitConvertToFloat )
 326                xcm = s->SplitOnTagValue(l, SplitGroup, SplitElem);
 327             else
 328             {
 329                 xcm = s->SplitOnTagValueConvertToFloat(l, SplitGroup, SplitElem);
 330             }
 331    }
 332
 333    if (xcm.size() == 0)
 334    {
 335       if(verbose) std::cout << "Empty XCoherent File Set after 'split' ?!?" << std::endl;
 336       return false;
 337    }
 338    else if (xcm.size() == 1)
 339       TypeResult=1;
 340    else
 341       TypeResult=2;
 342
 343    ImageDataVector = new std::vector<vtkImageData*>;
 344   // vtkGdcmReader *reader = vtkGdcmReader::New(); // move inside the loop, or be clever using vtk!
 345
 346    for (GDCM_NAME_SPACE::XCoherentFileSetmap::iterator i = xcm.begin();
 347                                                   i != xcm.end();
 348                                                 ++i)
 349    {
 350            if (verbose)
 351                std::cout << "--- xCoherentName = [" << (*i).first << "]" << std::endl;
 352    }
 353
 354  // XCoherentFileSetmap map < critère de split, FileList (= std::vector de gdcm::File*) >
 355
 356    for (GDCM_NAME_SPACE::XCoherentFileSetmap::iterator i = xcm.begin();
 357                                                   i != xcm.end();
 358                                                 ++i)
 359    {
 360
 361       vtkGdcmReader *reader = vtkGdcmReader::New(); /// \FIXME : unable to delete!
 362
 363       if (verbose)
 364                std::cout << "==========================================xCoherentName = [" << (*i).first << "]" << std::endl;
 365
 366       if (SortOnPosition)
 367       {
 368               if (verbose) std::cout << "SortOnPosition" << std::endl;
 369               // (will be IPPSorter, in GDCM2)
 370               s->ImagePositionPatientOrdering((*i).second);
 371               if (verbose) std::cout << "out of SortOnPosition" << std::endl;
 372       }
 373
 374       else if (SortOnOrientation)
 375       {
 376               if (verbose) std::cout << "SortOnOrientation" << std::endl;
 377             /// \TODO SortOnOrientation()
 378
 379             // we still miss an algo to sort an Orientation, given by 6 cosines!
 380             //  Anything like this, in GDCM2?
 381             std::cout << "SortOnOrientation : not so easy - I(mage)O(rientation)P(atient)Sorter still missing! -" << std::endl;
 382             // have a look at SerieHelper::SplitOnOrientation() to have an idea of the mess!
 383
 384             //Better sort on the file name, right now...
 385              s->FileNameOrdering((*i).second);
 386       }
 387
 388       else if (SortOnFileName)
 389       {
 390          if (verbose) std::cout << "SortOnFileName" << std::endl;
 391          if (verbose) std::cout << "taille " << ((*i).second)->size() << std::endl;
 392
 393          s->FileNameOrdering((*i).second);
 394          if (verbose) std::cout << "Out of SortOnFileName" << std::endl;
 395       }
 396
 397       else if (SortOnTag)
 398       {
 399          if (verbose) std::cout << "SortOnTag" << std::endl;
 400          printf ("--> %04x %04x\n", SortGroup,SortElem);
 401          std::cout << "Sorry, troubles not solved yet; use SortOnUserFunction, right now!" << std::endl;
 402
 403         /*        ==> WARNING : This one has troubles; do NOT use it, right now!
 404         // a pointer to fonction cannot be casted as a pointer to member function!
 405         // Use SortOnUserFunction, instead!
 406
 407          if ( SortConvertToFloat )
 408             s->SetUserLessThanFunction( reinterpret_cast<bool (*)(gdcm13::File*, gdcm13::File*)>
 409                                                                  ( &vtkGdcm4DSplitter::CompareOnSortTagConvertToFloat));
 410          else
 411             s->SetUserLessThanFunction( reinterpret_cast<bool (*)(gdcm13::File*, gdcm13::File*)>
 412                                                                  ( &vtkGdcm4DSplitter::CompareOnSortTag));
 413
 414          // Anything like this, in GDCM2?
 415          s->UserOrdering((*i).second);
 416         */
 417
 418          //if (verbose) std::cout << "Out of SortOnTag" << std::endl;
 419          std::cout << "NO ordering performed  :-( " << std::endl;
 420       }
 421
 422       else if (SortOnUserFunction)
 423       {
 424           if (verbose) std::cout << "SortOnUserFunction" << std::endl;
 425           s->SetUserLessThanFunction( UserCompareFunction );
 426          // Anything like this, in GDCM2?
 427          s->UserOrdering((*i).second);
 428          if (verbose) std::cout << "Out of SortOnUserFunction" << std::endl;
 429       }
 430
 431        reader->SetCoherentFileList((*i).second);
 432        reader->Update();
 433
 434        /// \TODO : remove the following
 435        if (verbose) {
 436           std::cout << "reader->GetOutput() :" << std::endl;
 437           reader->GetOutput()->PrintSelf(std::cout, vtkIndent(2));
 438        }
 439
 440        ImageDataVector->push_back(reader->GetOutput() );
 441
 442        std::vector<vtkImageData*>::iterator it;
 443        if (verbose)
 444        for(it=ImageDataVector->begin(); it!=ImageDataVector->end(); ++it) {
 445          std::cout << "-in vtkGdcm4DSplitter--------------------------" << std::endl;
 446          (*it)->PrintSelf(std::cout, vtkIndent(2));
 447        }
 448        std::cout << std::endl;
 449    }
 450
 451    //reader->Delete();  // \TODO : fix
 452    s->Delete();
 453    f->Delete();
 454    delete l;
 455
 456    return true;
 457  }
 458
 459
 460  void vtkGdcm4DSplitter::reorgXCoherentFileSetmap(GDCM_NAME_SPACE::XCoherentFileSetmap &xcm)
 461  {
 462    ELEM e;
 463    std::vector<ELEM> vectElem;
 464 /*
 465    typedef struct
 466    {
 467       std::string strIPP;
 468       double dist;
 469       GDCM_NAME_SPACE::File *file;
 470    } ELEM;
 471 */
 472
 473    bool Debug=true;
 474
 475    for (GDCM_NAME_SPACE::XCoherentFileSetmap::iterator i = xcm.begin();
 476                                                   i != xcm.end();
 477                                                 ++i)
 478    {
 479       if (verbose)
 480                std::cout << "--- xCoherentName = [" << (*i).first << "]" << std::endl;
 481
 482       e.strIPP  = (*i).first;
 483       e.file = *(((*i).second)->begin()); // all the gdcm::File of a given xcm item *have* the same IPP; first one is enough
 484       e.dist=0.0;
 485       vectElem.push_back(e);
 486    }
 487    sortVectElem(&vectElem);
 488
 489    // now, create the final std::map !
 490    // final_xcm<to_str(e.dist , xcm[e.strIPP]>
 491    // xcm = final_xcm;
 492    // check what we need to free !
 493
 494    int dist;
 495    char c_dist[100];
 496    std::string str_dist;
 497    int lgr=vectElem.size();
 498    GDCM_NAME_SPACE::XCoherentFileSetmap final_xcm;
 499    for (int i2=0; i2<lgr; i2++)
 500    {
 501       dist =  (vectElem[i2].dist*1000);
 502       sprintf(c_dist,"%010d",dist);
 503       str_dist = c_dist;
 504 /*
 505       std::cout << "dist " << vectElem[i2].dist
 506                 << " str_dist " << str_dist
 507                 << " IPP " << vectElem[i2].strIPP
 508                 << std::endl;
 509
 510 */
 511       final_xcm[str_dist] = xcm[vectElem[i2].strIPP];
 512    }
 513
 514    /// \TODO : check what needs to be cleared // JPR
 515
 516    xcm = final_xcm;
 517
 518  }
 519
 520
 521 bool vtkGdcm4DSplitter::sortVectElem(std::vector<ELEM> *fileList)
 522 {
 523 //based on Jolinda Smith's algorithm
 524
 525 //Tags always use the same coordinate system, where "x" is left
 526 //to right, "y" is posterior to anterior, and "z" is foot to head (RAH).
 527
 528    //iop is calculated based on the file file
 529    float cosines[6];
 530    double normal[3];
 531    double ipp[3];
 532    double dist;
 533    double min = 0, max = 0;
 534    bool first = true;
 535
 536    double ZSpacing; // useless here! // JPR
 537    bool DirectOrder = true; // remove it!
 538
 539    ZSpacing = -1.0;  // will be updated if process doesn't fail
 540
 541    //std::multimap<double,File *> distmultimap; // JPR
 542    std::multimap<double,ELEM> distmultimap; // JPR
 543
 544    // Use a multimap to sort the distances from 0,0,0
 545    //for ( FileList::const_iterator // JPR
 546    for ( std::vector<ELEM>::iterator   // JPR
 547          it = fileList->begin();
 548          it != fileList->end(); ++it )
 549    {
 550       //gdcmDebugMacro("deal with " << (*it)->file->GetFileName() );
 551       if ( first )
 552       {
 553          (*it).file->GetImageOrientationPatient( cosines );
 554
 555    // The "Image Orientation Patient" tag gives the direction cosines
 556    // for the rows and columns for the three axes defined above.
 557    // Typical axial slices will have a value 1/0/0/0/1/0:
 558    // rows increase from left to right,
 559    // columns increase from posterior to anterior. This is your everyday
 560    // "looking up from the bottom of the head with the eyeballs up" image.
 561
 562    // The "Image Position Patient" tag gives the coordinates of the first
 563    // voxel in the image in the "RAH" coordinate system, relative to some
 564    // origin.
 565
 566    // First, calculate the slice normal from IOP :
 567
 568          // You only have to do this once for all slices in the volume. Next,
 569          // for each slice, calculate the distance along the slice normal
 570          // using the IPP ("Image Position Patient") tag.
 571          // ("dist" is initialized to zero before reading the first slice) :
 572
 573          normal[0] = cosines[1]*cosines[5] - cosines[2]*cosines[4];
 574          normal[1] = cosines[2]*cosines[3] - cosines[0]*cosines[5];
 575          normal[2] = cosines[0]*cosines[4] - cosines[1]*cosines[3];
 576
 577    // For each slice (here : the first), calculate the distance along
 578    // the slice normal using the IPP tag
 579
 580          ipp[0] = (*it).file->GetXOrigin();
 581          ipp[1] = (*it).file->GetYOrigin();
 582          ipp[2] = (*it).file->GetZOrigin();
 583
 584          dist = 0;
 585          for ( int i = 0; i < 3; ++i )
 586          {
 587             dist += normal[i]*ipp[i];
 588          }
 589
 590          //gdcmDebugMacro("dist : " << dist);
 591          distmultimap.insert(std::pair<const double,ELEM>(dist, *it));
 592
 593          max = min = dist;
 594          first = false;
 595       }
 596       else
 597       {
 598    // Next, for each slice, calculate the distance along the slice normal
 599    // using the IPP tag
 600          ipp[0] = (*it).file->GetXOrigin();
 601          ipp[1] = (*it).file->GetYOrigin();
 602          ipp[2] = (*it).file->GetZOrigin();
 603
 604          dist = 0;
 605          for ( int i = 0; i < 3; ++i )
 606          {
 607             dist += normal[i]*ipp[i];
 608          }
 609
 610          (*it).dist = dist; // JPR
 611
 612          distmultimap.insert(std::pair<const double,ELEM>(dist, *it));
 613          //gdcmDebugMacro("dist : " << dist);
 614          min = (min < dist) ? min : dist;
 615          max = (max > dist) ? max : dist;
 616       }
 617    }
 618
 619   // gdcmDebugMacro("After parsing vector, nb of elements : " << fileList->size() );
 620
 621 /*Useless here. // JPR
 622
 623    // Find out if min/max are coherent
 624    if ( min == max )
 625    {
 626      gdcmWarningMacro("Looks like all images have the exact same image position. "
 627                       << "No PositionPatientOrdering sort performed. "
 628                       << "No 'ZSpacing' calculated! ");
 629      return false;
 630    }
 631 */
 632
 633 /* Useless here, 'split' already done. // JPR
 634
 635    // Check to see if image shares a common position
 636    bool ok = true;
 637    for (std::multimap<double, File *>::iterator it2 = distmultimap.begin();
 638         it2 != distmultimap.end();
 639         ++it2)
 640    {
 641
 642       gdcmDebugMacro("Check if image shares a common position : " << ((*it2).second).file->GetFileName() );
 643
 644       if (distmultimap.count((*it2).first) != 1)
 645       {
 646          gdcmWarningMacro("File: ["
 647               << ((*it2).second->GetFileName())
 648               << "] : more than ONE file at distance: '"
 649               << (*it2).first
 650               << " (position is not unique!) "
 651               << "No PositionPatientOrdering sort performed. "
 652               << "No 'ZSpacing' calculated! ");
 653
 654          ok = false;
 655       }
 656    }
 657    if (!ok)
 658    {
 659       if (! DropDuplicatePositions)
 660          return false;
 661    }
 662
 663 */
 664
 665 // Now, we can calculate Z Spacing as the difference
 666 // between the "dist" values for the first two slices.
 667
 668 // The following (un)-commented out code is let here
 669 // to be re-used by whomsoever is interested...
 670
 671     std::multimap<double, ELEM>::iterator it5 = distmultimap.begin();
 672     double d1 = (*it5).first;
 673     it5++;
 674     double d2 = (*it5).first;
 675     ZSpacing = d1-d2;
 676     if (ZSpacing < 0.0)
 677        ZSpacing = - ZSpacing;
 678
 679    fileList->clear();  // doesn't delete list elements, only nodes
 680
 681 // Acording to user requierement, we sort direct order or reverse order.
 682    if (DirectOrder)
 683    {
 684       for (std::multimap<double, ELEM>::iterator it3 = distmultimap.begin();
 685            it3 != distmultimap.end();
 686            ++it3)
 687       {
 688          fileList->push_back( (*it3).second );
 689
 690 /*       useless here! // JPR
 691
 692          if (DropDuplicatePositions)
 693          {
 694             // ImagePositionPatientOrdering  wrong duplicates are found ???
 695             // --> fixed. See comment
 696
 697             it3 =  distmultimap.upper_bound((*it3).first); // skip all duplicates
 698            // the upper_bound function increments the iterator to the next non-duplicate entry
 699            // The for loop iteration also increments the iterator, which causes the code to skip every other image
 700            // --> decrement the iterator after the upper_bound function call
 701             it3--;
 702             if (it3 == distmultimap.end() )  // if last image, stop iterate
 703                break;
 704          }
 705 */
 706       }
 707    }
 708    else // user asked for reverse order
 709    {
 710       std::multimap<double, ELEM>::const_iterator it4;
 711       it4 = distmultimap.end();
 712       do
 713       {
 714          it4--;
 715          fileList->push_back( (*it4).second );
 716
 717 /* useless here // JPR
 718
 719          if (DropDuplicatePositions)  // skip all duplicates
 720          {
 721             // lower_bound finds the next element that is
 722             // less than or *equal to* the current value!
 723             //it4 =  distmultimap.lower_bound((*it4).first);
 724
 725            // David Feng's fix
 726            std::multimap<double, ELEM>::const_iterator itPrev = it4;
 727            while (itPrev->first == it4->first)
 728               --itPrev;
 729            it4 = itPrev;
 730
 731            if (it4 == distmultimap.begin() ) // if first image, stop iterate
 732                break;
 733          }
 734 */
 735       } while (it4 != distmultimap.begin() );
 736    }
 737
 738    distmultimap.clear();
 739
 740    return true;
 741
 742 }