]> Creatis software - gdcm.git/blob - vtk/vtkGdcm4DSplitter.cxx
6413a24768db186fb14411b26a5d00cb6d9b2365
[gdcm.git] / 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 }