]> Creatis software - gdcm.git/blobdiff - vtk/vtkGdcm4DSplitter.cxx
fix mistyping
[gdcm.git] / vtk / vtkGdcm4DSplitter.cxx
index 7e06563a749aa5f755a4ee23e62f48dd6013f802..8296e635730fb72991be2b313e6a8433cddc314c 100644 (file)
@@ -3,8 +3,8 @@
   Program:   gdcm
   Module:    $RCSfile: vtkGdcm4DSplitter.cxx,v $
   Language:  C++
-  Date:      $Date: 2011/04/04 17:01:03 $
-  Version:   $Revision: 1.6 $
+  Date:      $Date: 2011/04/15 15:23:14 $
+  Version:   $Revision: 1.12 $
                                                                                 
   Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
   l'Image). All rights reserved. See Doc/License.txt or
@@ -38,14 +38,22 @@ vtkGdcm4DSplitter
 A 'Dicom Serie' doesn't mean always the same thing :
         a given Slice along the time
         a given Volume at a given time
-Sometimes, an image within a serie is so artefacted than user decides to replace
-it by an other image.
+Sometimes, an image within a serie is so artefacted than user decides to replace it
+           by an other image.
 
 User needs to be aware, *only him* knows want he wants to do.
-vtkGdcm4DSplitter class does the job for hom
+vtkGdcm4DSplitter class does the job for him
 (despite its name, it works on 3D or 2D+T images too)
 
-User will have to specify some points
+==> To (try to) understand how the file store is organised, user is deeply encourage to use something like:
+
+PrintFile dirin=Root_Directory_holding_the_images rec > 1.txt
+
+open it with any test editor, and/or (Linux users) 
+grep a_supposed_to_be_string_of_interest 1.txt
+
+Aware user will have to specify some points :
 
 . Choose input data
 ------------------- 
@@ -56,6 +64,8 @@ User will have to specify some points
        bool setVectDirName(std::vector<std::string> &vectDirName);
 - a list of files       
        bool setVectFileName(std::vector<std::string> &vectFileName);
+- a list of gdcm::File*      
+       bool setVectGdcmFile(std::vector<GDCM_NAME_SPACE::File *> &vectGdcmFile);
 
 - Recursive directory exploration
        void setRecursive(bool recursive);
@@ -73,6 +83,20 @@ User will have to specify some points
  - UserDefined Function
         void setSortOnUserFunction (FoncComp f);
  
+ for 'true' 3D image sets :
+   - if you want to get a single 3D vtkImageData, use SplitOnOrientation -i.e. no split-
+   - if you want to get a vector of 2D vtkImageData, use SplitOnPosition  -i.e. one slice in each 'XCoherent fileset'-
+
+ for 'true' 4D multi-orientation image sets (i.e. a stack of n axial + m sagital + p coronal images, at different instants ...)
+   --> this is 5D, right?
+   (almost) nothing done, yet :
+   . use setSplitOnly()
+   . Use a first time vtkGdcm4DSplitter with setSplitOnOrientation();
+   . Get the VectGdcmFileLists (a std::vector of 'XCoherent fileset')
+   . use vtkGdcm4DSplitter, with as many setVectGdcmFile(std::vector<GDCM_NAME_SPACE::File *> &vectGdcmFile) you need 
+        one per element of std::vector<GDCM_NAME_SPACE::File *>
+        think on 'spliting' and 'sorting' it, according to your needs. 
 . Choose 'sort' criterion :
 --------------------------
 
@@ -81,16 +105,18 @@ User will have to specify some points
  - ImageOrientationPatient
        ==> Only in your dreams!
        ==> or, please, write a IOP sorter ...
- - User choosen tag
-        ==> WARNING : This one has troubles; do NOT use it, right now!
-        ==> use setSortOnUserFunction instead 
-        void setSortOnTag(unsigned short sortGroup, unsigned short sortElem);
-        void setSortConvertToFloat(bool conv)
  - UserDefined Function
         void setSortOnUserFunction (FoncComp f);
  - File name
         void setSortOnFileName()
-    
+
+- Choose 'flip' option
+----------------------
+VTK reader 'flop upside down' the images to display them correctly ?!?
+When a 'math guy' gets the image, the first thing he does is flipping it again.
+Avoid loosing CPU time using :
+      void SetFlipY(bool); // Wrong = no flip
+
 . Execute :
 -----------
         bool Go();
@@ -98,26 +124,28 @@ User will have to specify some points
 . Get the result
 ----------------
 
- -a single vtkImageData:
-        vtkImageData *GetImageData();
+// -a single vtkImageData:
+//        vtkImageData *GetImageData();
 - a vector of vtkImageData
         std::vector<vtkImageData*> *GetImageDataVector();
 
   ===================================================================== */
 
-#include "gdcmSerieHelper.h"
-
 #include "vtkGdcmReader.h"
 #include "vtkGdcm4DSplitter.h"
 #include <algorithm>
 #include "gdcmSerieHelper.h" // for ImagePositionPatientOrdering()
 #include <stdlib.h> // for atof
 
+// Constructor / Destructor
+/**
+ * \brief   Constructor from a given vtkGdcm4DSplitter
+ */
  vtkGdcm4DSplitter::vtkGdcm4DSplitter() :
                  SplitOnPosition(false), SplitOnOrientation(false), SplitOnTag(false),
-                 SplitGroup(0), SplitElem(0),
+                 SplitGroup(0), SplitElem(0), SplitConvertToFloat(false),
 
-                 SortOnPosition(false),  SortOnOrientation(false),  SortOnTag(false), 
+                 SortOnPosition(false),  SortOnOrientation(false),  SortOnTag(false), SortOnFileName(false), SortOnUserFunction(false),
                  SortGroup(0),  SortElem(0), SortConvertToFloat(false),
 
                  Recursive(false), TypeDir(0),
@@ -126,35 +154,29 @@ User will have to specify some points
  
  }
 
- std::vector<vtkImageData*> * vtkGdcm4DSplitter::GetImageDataVector() 
+/**
+ * \brief   Canonical destructor.
+ */ 
+ vtkGdcm4DSplitter::~vtkGdcm4DSplitter()
  {
- if (verbose) std::cout << "TypeDir " << TypeDir << std::endl;
-    if (TypeResult == 2)
-       return ImageDataVector;
-    else
-      if (TypeResult == 1)
-      {
-         std::vector<vtkImageData*> *t = new std::vector<vtkImageData*>; 
-         t->push_back( ImageData );
-         return t;            
-      }
-      else
-         return (std::vector<vtkImageData*>*) NULL;
+    /// \TODO : delete everything that must be! 
  }
- vtkImageData *vtkGdcm4DSplitter::GetImageData() 
- {
-    if (TypeResult == 1)
-       return ImageData;
-    else
-      if (TypeResult == 1)
-      {
-         return (*ImageDataVector)[0];      
-      }
-      else
-         return (vtkImageData*) NULL;
- }      
-       
+
+       // Locate Data to process
+       // ======================
+/**
+ * \brief sets the directories exploration mode
+ * @param recursive whether we want explore recursively the root Directory
+ */       
+void  vtkGdcm4DSplitter::setRecursive(bool recursive) 
+{ 
+   Recursive=recursive;
+}
+      
+/**
+ * \brief Sets the root Directory to get the images from
+ * @param   dirName name of the directory to deal with
+ */       
  bool vtkGdcm4DSplitter::setDirName(std::string &dirName) 
  {
     if ( ! GDCM_NAME_SPACE::DirList::IsDirectory(dirName) ) 
@@ -166,7 +188,11 @@ User will have to specify some points
     TypeDir=1;
     return true;
  }
+
+/**
+ * \brief Sets a list of Directories to get the images from
+ * @param   vectDirName vector of directory names to deal with
+ */   
  bool vtkGdcm4DSplitter::setVectDirName(std::vector<std::string> &vectDirName) 
  {
     int nbDir = vectDirName.size();
@@ -183,7 +209,11 @@ User will have to specify some points
     TypeDir=2;
     return true;
  }
+
+/**
+ * \brief Sets a list of files read
+ * @param   vectFileName vector of file names to deal with
+ */ 
  bool vtkGdcm4DSplitter::setVectFileName(std::vector<std::string> &vectFileName)
  {
     if ( vectFileName.size() == 0)
@@ -196,6 +226,160 @@ User will have to specify some points
     return true;
  }      
 
+/**
+ * \brief Sets a vector of gdcm::File *
+ * @param   vectGdcmFileName vector of gdcm::File *
+ */
+ bool vtkGdcm4DSplitter::setVectGdcmFile(GDCM_NAME_SPACE::FileList *vectGdcmFile)
+ {
+    if ( vectGdcmFile->size() == 0)
+    {
+          std::cout << "[ vectGdcmFile ] : empty list" << std::endl;
+          return false;
+    }
+    TypeDir=4;
+    VectGdcmFile = vectGdcmFile; 
+ } 
+
+
+       // Split
+       // =====
+
+/**
+ * \brief   asks for splitting  Filesets according to the Position
+ */
+       
+ void  vtkGdcm4DSplitter::setSplitOnPosition()
+ {
+    SplitOnPosition=true;
+    SplitOnOrientation=false;
+    SplitOnTag=false;
+ }
+/**
+ * \brief   asks for splitting  Filesets according to the Orientation
+ */ 
+ void  vtkGdcm4DSplitter::setSplitOnOrientation()
+ {
+    SplitOnPosition=false;
+    SplitOnOrientation=true; 
+    SplitOnTag=false;
+ }
+/**
+ * \brief   asks for splitting  Filesets according to the value of a given Tag
+ * @param   group  group number of the target Element
+ * @param   element element number of the target Element 
+ */ 
+ void  vtkGdcm4DSplitter::setSplitOnTag(unsigned short int splitGroup, unsigned short int splitElem)
+ {
+    SplitOnPosition=false;
+    SplitOnOrientation=false;
+    SplitOnTag=true;
+    SplitGroup=splitGroup;
+    SplitElem=splitElem;
+ }
+/**
+ * \brief   asks for converting to 'float' the tag values used as a splitting criteria (lexicographic order may not be suitable)
+ */   
+ void  vtkGdcm4DSplitter::setSplitConvertToFloat(bool conv) {SplitConvertToFloat=conv;}
+/**
+ * \brief   asks for splitting Filesets according to what was asked for (no sorting, no reading data)
+ */  
+ void  vtkGdcm4DSplitter::setSplitOnly(bool s)
+ {
+    SplitOnly = s;
+ }
+       // Sort
+       // ====
+
+ void  vtkGdcm4DSplitter::setSortOnPosition() 
+ {
+    SortOnPosition=true;
+    SortOnOrientation=false;
+    SortOnTag=false;
+    SortOnFileName=false;
+    SortOnUserFunction=false;
+    SortOnPosition=true;
+ }
+      // use setSortOnUserFunction, instead!
+      // void setSortOnTag(unsigned short int sortGroup, unsigned short int sortElem)
+      // {
+      //    SortOnPosition=false;
+      //    SortOnOrientation=false;
+      //    SortOnTag=true;
+      //    SortOnFileName=false;
+      //    SortOnUserFunction=false;
+      //    SortGroup=sortGroup;  SortElem=sortElem;
+      // }
+
+
+/**
+ * \brief sets a user supplied function (comparison)
+ * @param f comparison function
+ */
+ void  vtkGdcm4DSplitter::setSortOnUserFunction (FoncComp f)
+ {
+    UserCompareFunction=f;
+    SortOnPosition=false;
+    SortOnOrientation=false;
+    SortOnTag=false;
+    SortOnFileName=false;
+    SortOnUserFunction=true;
+  }
+
+
+ //  void setSortConvertToFloat(bool conv)
+ //  {
+ //     SortConvertToFloat=conv;
+ //  }
+
+/**
+ * \brief asks for sorting  the images, according to their File Name
+ */
+ void  vtkGdcm4DSplitter::setSortOnFileName()
+ {
+    SortOnPosition=false;
+    SortOnOrientation=false;
+    SortOnTag=false;
+    SortOnFileName=true;
+    SortOnUserFunction=false;
+ }
+
+
+ std::vector<vtkImageData*> * vtkGdcm4DSplitter::GetImageDataVector() 
+ { 
+    if (SplitOnly)
+       return NULL;
+
+     return ImageDataVector;
+ }
+ std::vector<GDCM_NAME_SPACE::FileList *> *vtkGdcm4DSplitter::GetVectGdcmFileLists()
+ {
+    if (SplitOnly)
+        return NULL;
+
+    GDCM_NAME_SPACE::XCoherentFileSetmap::iterator it;
+    for ( it = xcm.begin();
+          it != xcm.end();
+        ++it)
+    {
+       VectGdcmFileLists.push_back((*it).second); 
+    } 
+    return  &VectGdcmFileLists;         
+ }
+
+ vtkImageData *vtkGdcm4DSplitter::GetImageData() 
+ {
+   if (SplitOnly)
+      return NULL;
+   return (*ImageDataVector)[0]; 
+ }      
+
  bool vtkGdcm4DSplitter::CompareOnSortTagConvertToFloat(GDCM_NAME_SPACE::File *file1, GDCM_NAME_SPACE::File *file2)
  { 
   /* if (verbose) printf ("%04x %04x\n", this->SortGroup,this->SortElem);
@@ -215,7 +399,7 @@ User will have to specify some points
  
  bool vtkGdcm4DSplitter::Go()
  {
-   if (!SplitOnPosition && !SplitOnOrientation && !SplitOnTag 
+   if (!SplitOnPosition && !SplitOnOrientation && !SplitOnTag) 
    {
        ///\TODO (?) Throw an exception "Choose Splitting mode before!"
        std::cout << "Choose Splitting mode before!" << std::endl;
@@ -240,7 +424,11 @@ User will have to specify some points
 
    GDCM_NAME_SPACE::File *f;
    GDCM_NAME_SPACE::DirListType fileNames;
-   
+ //
+ // Fill  fileNames with the user supplied file names (in any)
+ // --------------------------------------
+ //
    if (TypeDir == 0 )
    {
       ///\TODO (?) Throw an exception "Set input Directory name(s) / file names  before!"
@@ -263,7 +451,7 @@ User will have to specify some points
         tmpFileNames = dirlist.GetFilenames();
         // Concat two std::vector
         //vector1.insert( vector1.end(), vector2.begin(), vector2.end() );
-       fileNames.insert( fileNames.end(), tmpFileNames.begin(), tmpFileNames.end() );
+        fileNames.insert( fileNames.end(), tmpFileNames.begin(), tmpFileNames.end() );
       }    
    }
    else if (TypeDir == 3)
@@ -271,34 +459,57 @@ User will have to specify some points
       fileNames=VectFileName;
    }  
 
+ //
+ // Fill l with the gdcm::File* corresponding to the files
+ // --------------------------------------
+ //
+
    GDCM_NAME_SPACE::FileList *l = new GDCM_NAME_SPACE::FileList; // (set of gdcm::File)
-   double floatTagvalue;  
-   // Loop on all the gdcm-readable files
-   for (GDCM_NAME_SPACE::DirListType::iterator it = fileNames.begin();
-                                    it != fileNames.end();
-                                  ++it)
+   
+   if (TypeDir == 4)
    {
-      int maxSize  = 0x7fff;         // load Elements of any length
-      f = GDCM_NAME_SPACE::File::New();
-      f->SetMaxSizeLoadEntry(maxSize);
-      f->SetFileName( *it );
-      if (f->Load())
-         l->push_back(f);
-      else 
-         std::cout << " Fail to load [" <<  *it << "]" << std::endl;          
-   }   
+   // User passed a vector of gdcm::File* 
+      l = VectGdcmFile;
+   } 
+   else
+   { 
+      double floatTagvalue;  
+      // Loop on all the gdcm-readable files
+      for (GDCM_NAME_SPACE::DirListType::iterator it = fileNames.begin();
+                                                  it != fileNames.end();
+                                                ++it)
+      {
+         int maxSize  = 0x7fff;         // load Elements of any length
+         f = GDCM_NAME_SPACE::File::New();
+         f->SetMaxSizeLoadEntry(maxSize);
+         f->SetFileName( *it );
+         if (f->Load())
+            l->push_back(f);
+         else 
+            std::cout << " Fail to load [" <<  *it << "]" << std::endl;          
+      }
+   } 
 
-   GDCM_NAME_SPACE::XCoherentFileSetmap xcm;
 
+//
+// Split the gdcm::File* set, according to user's requierements
+// ------------------------------------------------------------
+//
    if (SplitOnOrientation) 
    {
-            s->SetDropDuplicatePositions(false);
-            xcm = s->SplitOnOrientation(l);
+      s->SetDropDuplicatePositions(false);
+      xcm = s->SplitOnOrientation(l);
    }
    else if (SplitOnPosition)
    {
-            s->SetDropDuplicatePositions(true);
-            xcm = s->SplitOnPosition(l);
+      s->SetDropDuplicatePositions(true);
+      xcm = s->SplitOnPosition(l);
+
+      // the key of xcm follows lexicographical order
+      // (that may be different than the 'distance' order)
+      // we have to reorganize it!
+
+      reorgXCoherentFileSetmap(xcm);
    }
    else if (SplitOnTag) 
    {
@@ -313,38 +524,42 @@ User will have to specify some points
                 xcm = s->SplitOnTagValueConvertToFloat(l, SplitGroup, SplitElem);
             }
    }
-   
+  
    if (xcm.size() == 0)
    {
       if(verbose) std::cout << "Empty XCoherent File Set after 'split' ?!?" << std::endl;
       return false;
    }
+/*
    else if (xcm.size() == 1)
       TypeResult=1;
    else
       TypeResult=2;
+*/
 
+   if(SplitOnly)
+      return true;   
+//
+//
+// ------------------------------------------------------------
+//
    ImageDataVector = new std::vector<vtkImageData*>;
+   /// \TODO move inside the loop, or be clever using vtk!
   // vtkGdcmReader *reader = vtkGdcmReader::New(); // move inside the loop, or be clever using vtk!
    
-   for (GDCM_NAME_SPACE::XCoherentFileSetmap::iterator i = xcm.begin(); 
-                                                  i != xcm.end();
-                                                ++i)
-   {
-           if (verbose)
-               std::cout << "--- xCoherentName = [" << (*i).first << "]" << std::endl;
-   }
  // XCoherentFileSetmap map < critère de split, FileList (= std::vector de gdcm::File*) >
 
    for (GDCM_NAME_SPACE::XCoherentFileSetmap::iterator i = xcm.begin(); 
                                                   i != xcm.end();
                                                 ++i)
    {
-   
-      vtkGdcmReader *reader = vtkGdcmReader::New(); /// \FIXME : unable to delete!
-       
+      vtkGdcmReader *reader = vtkGdcmReader::New(); /// \TODO FIXME : unable to delete!
+      
+      reader->SetFlipY(FlipY);
+      // better user SetFileLowerLeft()
+      /// \TODO : modify vtkGdcmReader ! 
       if (verbose)
-               std::cout << "==========================================xCoherentName = [" << (*i).first << "]" << std::endl;
+               std::cout << " --- xCoherentName = [" << (*i).first << "]" << std::endl;
 
       if (SortOnPosition)
       {
@@ -357,7 +572,7 @@ User will have to specify some points
       else if (SortOnOrientation)
       {
               if (verbose) std::cout << "SortOnOrientation" << std::endl;
-            /// \TODO SortOnOrientation()
+            /// \TODO (?) SortOnOrientation()
       
             // we still miss an algo to sort an Orientation, given by 6 cosines!
             //  Anything like this, in GDCM2? 
@@ -387,15 +602,15 @@ User will have to specify some points
         // a pointer to fonction cannot be casted as a pointer to member function!
         // Use SortOnUserFunction, instead!
 
-         if ( SortConvertToFloat )
-            s->SetUserLessThanFunction( reinterpret_cast<bool (*)(gdcm13::File*, gdcm13::File*)> 
+        //  if ( SortConvertToFloat )
+         //    s->SetUserLessThanFunction( reinterpret_cast<bool (*)(gdcm13::File*, gdcm13::File*)> 
                                                                  ( &vtkGdcm4DSplitter::CompareOnSortTagConvertToFloat));     
-         else
-            s->SetUserLessThanFunction( reinterpret_cast<bool (*)(gdcm13::File*, gdcm13::File*)>
+        //  else
+        //     s->SetUserLessThanFunction( reinterpret_cast<bool (*)(gdcm13::File*, gdcm13::File*)>
                                                                  ( &vtkGdcm4DSplitter::CompareOnSortTag)); 
        
          // Anything like this, in GDCM2? 
-         s->UserOrdering((*i).second);
+        //  s->UserOrdering((*i).second);
         */
 
          //if (verbose) std::cout << "Out of SortOnTag" << std::endl;
@@ -415,10 +630,19 @@ User will have to specify some points
        reader->Update();
        
        /// \TODO : remove the following
-       //if (verbose) reader->GetOutput()->PrintSelf(std::cout, vtkIndent(2));
+       if (verbose) {
+          std::cout << "reader->GetOutput() :" << std::endl;
+          reader->GetOutput()->PrintSelf(std::cout, vtkIndent(2));
+       }
        
        ImageDataVector->push_back(reader->GetOutput() );
-
+       
+       std::vector<vtkImageData*>::iterator it; 
+       if (verbose)      
+       for(it=ImageDataVector->begin(); it!=ImageDataVector->end(); ++it) {
+         std::cout << "-in vtkGdcm4DSplitter--------------------------" << std::endl;
+         (*it)->PrintSelf(std::cout, vtkIndent(2));
+       }
        std::cout << std::endl;
    }
 
@@ -430,3 +654,300 @@ User will have to specify some points
    return true;
  }
 
+
+ void vtkGdcm4DSplitter::reorgXCoherentFileSetmap(GDCM_NAME_SPACE::XCoherentFileSetmap &xcm)
+ {
+ /*
+ the key of the 'XCoherentFileSetmap', is a std::string, used as if it was found in the Dicom header
+ Normaly(?), it's suitable for almost everything ...
+ ... but the 'Image Position Patient'.
+ We need to order the 'XCoherentFileSetmap' (NOT the content of each XCoherentFileSet!) according to the IPP,
+ using Jolinda Smith's algorithm.
+ (we use a subset of the one defined in gdcm::SerieHelper)
+*/
+
+   ELEM e;   
+   std::vector<ELEM> vectElem;
+
+/* remenber :
+   typedef struct 
+   {
+      std::string strIPP;
+      double dist;
+      GDCM_NAME_SPACE::File *file;
+   } ELEM;   
+*/
+
+   bool Debug=true;
+   for (GDCM_NAME_SPACE::XCoherentFileSetmap::iterator i = xcm.begin(); 
+                                                       i != xcm.end();
+                                                     ++i)
+   {
+      if (verbose)
+               std::cout << "--- xCoherentName = [" << (*i).first << "]" << std::endl;
+
+      e.strIPP  = (*i).first;
+      e.file = *(((*i).second)->begin()); // all the gdcm::File of a given xcm item *have* the same IPP; first one is enough
+      e.dist=0.0;
+      vectElem.push_back(e);   
+   }  
+   sortVectElem(&vectElem);
+   
+   // now, create the final std::map !
+   // final_xcm<to_str(e.dist , xcm[e.strIPP]>
+   // xcm = final_xcm;
+   // check what we need to free !
+
+   int dist;
+   char c_dist[100];
+   std::string str_dist;
+   int lgr=vectElem.size();
+   GDCM_NAME_SPACE::XCoherentFileSetmap final_xcm;
+   for (int i2=0; i2<lgr; i2++)
+   { 
+      dist =  (vectElem[i2].dist*1000);
+      sprintf(c_dist,"%010d",dist);
+      str_dist = c_dist;
+/*
+      std::cout << "dist " << vectElem[i2].dist 
+                << " str_dist " << str_dist
+                << " IPP " << vectElem[i2].strIPP
+                << std::endl;
+
+*/     
+      final_xcm[str_dist] = xcm[vectElem[i2].strIPP];
+   } 
+   
+   /// \TODO : check what needs to be cleared // JPR
+
+   xcm = final_xcm;
+
+ }
+
+
+bool vtkGdcm4DSplitter::sortVectElem(std::vector<ELEM> *fileList)
+{
+//based on Jolinda Smith's algorithm
+// NOTE : if you need to use Jolinda Smith's algorithm, get the one inside gdcm::SerieHelper
+// this one is a light version.
+
+//Tags always use the same coordinate system, where "x" is left
+//to right, "y" is posterior to anterior, and "z" is foot to head (RAH).
+
+   //iop is calculated based on the file file
+   float cosines[6];
+   double normal[3];
+   double ipp[3];
+   double dist;
+   double min = 0, max = 0;
+   bool first = true;
+   
+   //double ZSpacing; // useless here! // JPR
+   bool DirectOrder = true; // remove it!
+   
+  // ZSpacing = -1.0;  // will be updated if process doesn't fail
+    
+   //std::multimap<double,File *> distmultimap; // JPR
+   std::multimap<double,ELEM> distmultimap; // JPR
+   
+   // Use a multimap to sort the distances from 0,0,0
+   //for ( FileList::const_iterator // JPR
+   for ( std::vector<ELEM>::iterator   // JPR
+         it = fileList->begin();
+         it != fileList->end(); ++it )
+   {
+      //gdcmDebugMacro("deal with " << (*it)->file->GetFileName() );
+      if ( first ) 
+      {
+         (*it).file->GetImageOrientationPatient( cosines );
+
+   // The "Image Orientation Patient" tag gives the direction cosines 
+   // for the rows and columns for the three axes defined above. 
+   // Typical axial slices will have a value 1/0/0/0/1/0: 
+   // rows increase from left to right, 
+   // columns increase from posterior to anterior. This is your everyday
+   // "looking up from the bottom of the head with the eyeballs up" image. 
+   
+   // The "Image Position Patient" tag gives the coordinates of the first
+   // voxel in the image in the "RAH" coordinate system, relative to some
+   // origin.   
+
+   // First, calculate the slice normal from IOP : 
+          
+         // You only have to do this once for all slices in the volume. Next, 
+         // for each slice, calculate the distance along the slice normal 
+         // using the IPP ("Image Position Patient") tag.
+         // ("dist" is initialized to zero before reading the first slice) :
+
+         normal[0] = cosines[1]*cosines[5] - cosines[2]*cosines[4];
+         normal[1] = cosines[2]*cosines[3] - cosines[0]*cosines[5];
+         normal[2] = cosines[0]*cosines[4] - cosines[1]*cosines[3];
+
+   // For each slice (here : the first), calculate the distance along 
+   // the slice normal using the IPP tag 
+    
+         ipp[0] = (*it).file->GetXOrigin();
+         ipp[1] = (*it).file->GetYOrigin();
+         ipp[2] = (*it).file->GetZOrigin();
+
+         dist = 0;
+         for ( int i = 0; i < 3; ++i )
+         {
+            dist += normal[i]*ipp[i];
+         }
+    
+         //gdcmDebugMacro("dist : " << dist);
+         distmultimap.insert(std::pair<const double,ELEM>(dist, *it));
+
+         max = min = dist;
+         first = false;
+      }
+      else 
+      {
+   // Next, for each slice, calculate the distance along the slice normal
+   // using the IPP tag 
+         ipp[0] = (*it).file->GetXOrigin();
+         ipp[1] = (*it).file->GetYOrigin();
+         ipp[2] = (*it).file->GetZOrigin();
+
+         dist = 0;
+         for ( int i = 0; i < 3; ++i )
+         {
+            dist += normal[i]*ipp[i];
+         }
+
+         (*it).dist = dist; // JPR
+
+         distmultimap.insert(std::pair<const double,ELEM>(dist, *it));
+         //gdcmDebugMacro("dist : " << dist);
+         min = (min < dist) ? min : dist;
+         max = (max > dist) ? max : dist;
+      }
+   }
+
+  // gdcmDebugMacro("After parsing vector, nb of elements : " << fileList->size() );
+
+/*Useless here. // JPR
+
+   // Find out if min/max are coherent
+   if ( min == max )
+   {
+     gdcmWarningMacro("Looks like all images have the exact same image position. "
+                      << "No PositionPatientOrdering sort performed. "
+                      << "No 'ZSpacing' calculated! ");
+     return false;
+   }
+*/
+
+/* Useless here, 'split' already done. // JPR
+
+   // Check to see if image shares a common position
+   bool ok = true;
+   for (std::multimap<double, File *>::iterator it2 = distmultimap.begin(); 
+        it2 != distmultimap.end();
+        ++it2)
+   {
+   
+      gdcmDebugMacro("Check if image shares a common position : " << ((*it2).second).file->GetFileName() );   
+   
+      if (distmultimap.count((*it2).first) != 1)
+      {
+         gdcmWarningMacro("File: ["
+              << ((*it2).second->GetFileName())
+              << "] : more than ONE file at distance: '"
+              << (*it2).first
+              << " (position is not unique!) "
+              << "No PositionPatientOrdering sort performed. "
+              << "No 'ZSpacing' calculated! ");      
+
+         ok = false;
+      }
+   }
+   if (!ok)
+   {
+      if (! DropDuplicatePositions)
+         return false;
+   }
+
+*/
+      
+// Now, we can calculate Z Spacing as the difference
+// between the "dist" values for the first two slices.
+
+// The following (un)-commented out code is let here
+// to be re-used by whomsoever is interested...
+
+    //std::multimap<double, ELEM>::iterator it5 = distmultimap.begin();
+    //double d1 = (*it5).first;
+    //it5++;
+    //double d2 = (*it5).first;
+    //ZSpacing = d1-d2;
+    //if (ZSpacing < 0.0)
+    //   ZSpacing = - ZSpacing;
+
+   fileList->clear();  // doesn't delete list elements, only nodes
+
+// Acording to user requierement, we sort direct order or reverse order.
+   if (DirectOrder)
+   {  
+      for (std::multimap<double, ELEM>::iterator it3 = distmultimap.begin();
+           it3 != distmultimap.end();
+           ++it3)
+      {
+         fileList->push_back( (*it3).second );
+/*       useless here! // JPR
+
+         if (DropDuplicatePositions)
+         {
+            // ImagePositionPatientOrdering  wrong duplicates are found ???
+            // --> fixed. See comment
+
+            it3 =  distmultimap.upper_bound((*it3).first); // skip all duplicates
+           // the upper_bound function increments the iterator to the next non-duplicate entry
+           // The for loop iteration also increments the iterator, which causes the code to skip every other image
+           // --> decrement the iterator after the upper_bound function call
+            it3--;
+            if (it3 == distmultimap.end() )  // if last image, stop iterate
+               break;
+         }
+*/
+      }
+   }
+   else // user asked for reverse order
+   {
+      std::multimap<double, ELEM>::const_iterator it4;
+      it4 = distmultimap.end();
+      do
+      {
+         it4--;
+         fileList->push_back( (*it4).second );
+
+/* useless here // JPR
+
+         if (DropDuplicatePositions)  // skip all duplicates
+         {
+            // lower_bound finds the next element that is 
+            // less than or *equal to* the current value!
+            //it4 =  distmultimap.lower_bound((*it4).first);
+   
+           // David Feng's fix
+           std::multimap<double, ELEM>::const_iterator itPrev = it4;
+           while (itPrev->first == it4->first)
+              --itPrev;
+           it4 = itPrev;
+    
+           if (it4 == distmultimap.begin() ) // if first image, stop iterate
+               break;
+         }
+*/ 
+      } while (it4 != distmultimap.begin() );
+   }
+
+   distmultimap.clear();
+
+   return true;
+
+}