]> Creatis software - gdcm.git/blob - src/gdcmSerieHelper.cxx
re indent
[gdcm.git] / src / gdcmSerieHelper.cxx
1 /*=========================================================================
2                                                                                 
3   Program:   gdcm
4   Module:    $RCSfile: gdcmSerieHelper.cxx,v $
5   Language:  C++
6   Date:      $Date: 2011/04/04 17:33:58 $
7   Version:   $Revision: 1.73 $
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 #include "gdcmSerieHelper.h"
20 #include "gdcmDirList.h"
21 #include "gdcmFile.h"
22 //#include "gdcmDictEntry.h" // for TranslateToKey : no more !
23 #include "gdcmDebug.h"
24 #include "gdcmUtil.h"
25
26 #include <math.h>
27 #include <vector>
28 #include <algorithm>
29 #include <map>
30 #include <stdio.h>  //for sscanf
31
32 namespace GDCM_NAME_SPACE
33 {
34
35 //-----------------------------------------------------------------------------
36 // Constructor / Destructor
37 /**
38  * \brief   Constructor from a given SerieHelper
39  */
40 SerieHelper::SerieHelper()
41 {
42    m_UseSeriesDetails = false;
43    ClearAll();
44    UserLessThanFunction = 0;
45    DirectOrder = true;
46    DropDuplicatePositions = false;   
47 }
48
49 /**
50  * \brief   Canonical destructor.
51  */
52 SerieHelper::~SerieHelper()
53 {
54    ClearAll();
55 }
56
57 /**
58  * \brief  Preventively, clear everything at constructor time.
59  *         ( use it at destructor time.)
60  */
61 void SerieHelper::ClearAll()
62 {
63    // For all the 'Single SerieUID' Filesets that may already exist 
64    FileList *l = GetFirstSingleSerieUIDFileSet();
65    while (l)
66    { 
67       // For all the GDCM_NAME_SPACE::File of a File set
68       for (GDCM_NAME_SPACE::FileList::iterator it  = l->begin();
69                                     it != l->end(); 
70                                   ++it)
71       {
72          (*it)->Delete(); // remove each entry
73       }
74       l->clear();
75       delete l;     // remove the container
76       l = GetNextSingleSerieUIDFileSet();
77    }
78    // Need to clear that too:
79    SingleSerieUIDFileSetHT.clear();
80 }
81
82 //-----------------------------------------------------------------------------
83
84 //-----------------------------------------------------------------------------
85
86 // Public
87 /**
88  * \brief add a GDCM_NAME_SPACE::File to the Fileset corresponding to its Serie UID
89  * @param   filename Name of the file to deal with
90  */
91 void SerieHelper::AddFileName(std::string const &filename)
92 {
93    // Create a DICOM file
94    File *header = File::New();
95    header->SetLoadMode(LoadMode);
96    header->SetFileName( filename ); 
97    header->Load();
98
99    if ( header->IsReadable() )
100    {
101       if ( !AddFile( header ) )
102       {
103          // at least one rule was unmatched we need to deallocate the file:
104          header->Delete();
105       }
106    }
107    else
108    {
109       gdcmWarningMacro("Could not read file: " << filename );
110       header->Delete();
111    }
112 }
113
114 /**
115  * \brief add a GDCM_NAME_SPACE::File to the first (and supposed to be unique) file set
116  *        of the GDCM_NAME_SPACE::SerieHelper.
117  * \warning : this method should be used by aware users only!
118  *           Passing a GDCM_NAME_SPACE::File* has the same effect than passing a file name!
119  * \todo : decide which one is wrong (the method, or the commentary)!
120  *           the following comment doesn't match the method :-(
121  *            User is supposed to know the files he want to deal with
122  *           and consider them they belong to the same Serie
123  *           (even if their Serie UID is different)
124  *           user will probabely OrderFileList() this list (actually, ordering
125  *           user choosen gdm::File is the sole interest of this method)
126  *           Moreover, using vtkGdcmReader::SetCoherentFileList() will avoid
127  *           vtkGdcmReader parsing twice the same files. 
128  *           *no* coherence check is performed, but those specified
129  *           by SerieHelper::AddRestriction()
130  * @param   header GDCM_NAME_SPACE::File* of the file to deal with
131  * @return  true if file was added, false if file was rejected
132  */
133 bool SerieHelper::AddFile(File *header)
134 {
135    int allrules = 1;
136    // First step the user has defined a set of rules for the DICOM 
137    // he is looking for.
138    // make sure the file correspond to his set of rules:
139
140    std::string s;
141    for(SerieExRestrictions::iterator it2 = ExRestrictions.begin();
142      it2 != ExRestrictions.end();
143      ++it2)
144    {
145       const ExRule &r = *it2;
146       s = header->GetEntryString( r.group, r.elem );
147       if ( !Util::CompareDicomString(s, r.value.c_str(), r.op) )
148       {
149          // Argh ! This rule is unmatched; let's just quit
150          allrules = 0;
151          break;
152       }
153    }
154
155    if ( allrules ) // all rules are respected:
156    {
157       // Allright! we have a found a DICOM that matches the user expectation. 
158       // Let's add it to the specific 'id' which by default is uid (Serie UID)
159       // but can be `refined` by user with more paramater (see AddRestriction(g,e))
160  
161       std::string id = CreateUniqueSeriesIdentifier( header );
162       // if id == GDCM_UNFOUND then consistently we should find GDCM_UNFOUND
163       // no need here to do anything special
164
165       if ( SingleSerieUIDFileSetHT.count(id) == 0 )
166       {
167          gdcmDebugMacro(" New/gdcmSerieHelper.cxx Serie UID :[" << id << "]");
168          // create a std::list in 'id' position
169          SingleSerieUIDFileSetHT[id] = new FileList;
170       }
171       // Current Serie UID and DICOM header seems to match add the file:
172       SingleSerieUIDFileSetHT[id]->push_back( header );
173    }
174    else
175    {
176       // one rule not matched, tell user:
177       return false;
178    }
179    return true;
180 }
181
182 /**
183  * \brief add a rule for restricting a DICOM file to be in the serie we are
184  * trying to find. For example you can select only the DICOM files from a
185  * directory which would have a particular EchoTime==4.0.
186  * This method is a user level, value is not required to be formatted as a DICOM
187  * string
188  * \todo find a trick to allow user to say if he wants the Rectrictions 
189  *       to be *ored* (and not only *anded*)
190  * @param   key  Target tag we want restrict on a given value
191  * @param value value to be checked to exclude File
192  * @param op  operator we want to use to check
193  */
194 void SerieHelper::AddRestriction(TagKey const &key, 
195                                  std::string const &value, int op)
196 {
197    ExRule r;
198    r.group = key[0];
199    r.elem  = key[1];
200    r.value = value;
201    r.op    = op;
202    ExRestrictions.push_back( r ); 
203 }
204
205 void SerieHelper::AddRestriction(TagKey const &key)
206 {
207   ExRule r;
208   r.group = key[0];
209   r.elem  = key[1];
210   ExRefine.push_back( r );
211 }
212
213 //#ifndef GDCM_LEGACY_REMOVE
214 /**
215  * \brief add a rule for restricting a DICOM file to be in the serie we are
216  * trying to find. For example you can select only the DICOM files from a
217  * directory which would have a particular EchoTime==4.0.
218  * This method is a user level, value is not required to be formatted as a DICOM
219  * string
220  * @param   group tag group number we want restrict on a given value
221  * @param   elem  tag element number we want restrict on a given value 
222  * @param value value to be checked to exclude File
223  * @param op  operator we want to use to check
224  * @deprecated use : AddRestriction(TagKey const &key, 
225  *                                 std::string const &value, int op);
226  */
227
228 void SerieHelper::AddRestriction(uint16_t group, uint16_t elem, 
229                                  std::string const &value, int op)
230 {
231   TagKey t(group, elem);
232   AddRestriction(t, value, op);
233 }
234
235 //#endif
236
237 /**
238  * \brief add an extra  'SerieDetail' for building a 'Serie Identifier'
239  *        that ensures (hope so) File consistency (Series Instance UID doesn't)
240  * @param   group tag group number we want restrict on a given value
241  * @param   elem  tag element number we want restrict on a given value
242  * @param  convert wether we want 'convertion', to allow further ordering
243  *         e.g : 100 would be *before* 20; 000020.00 vs 00100.00 : OK 
244  */
245 void SerieHelper::AddSeriesDetail(uint16_t group, uint16_t elem, bool convert)
246 {   
247    ExDetail d;
248    d.group   = group;
249    d.elem    = elem;
250    d.convert = convert;
251    ExDetails.push_back( d ); 
252 }
253 /**
254  * \brief Sets the root Directory
255  * @param   dir Name of the directory to deal with
256  * @param recursive whether we want explore recursively the root Directory
257  */
258 void SerieHelper::SetDirectory(std::string const &dir, bool recursive)
259 {
260    DirList dirList(dir, recursive); // OS specific
261
262    DirListType filenames_list = dirList.GetFilenames();
263    for( DirListType::const_iterator it = filenames_list.begin(); 
264         it != filenames_list.end(); ++it)
265    {
266      // std::cout << "-----------------------------filename [" << *it << "]"
267      //           << std::endl;
268       gdcmDebugMacro("filename [" << *it << "]" );
269       AddFileName( *it );
270    }
271 }
272
273 /**
274  * \brief Sets the DicomDirSerie
275  * @param   se DicomDirSerie to deal with
276  */
277 void SerieHelper::SetDicomDirSerie(DicomDirSerie *se)
278 {
279    DirList dirList(se);
280   
281    DirListType filenames_list = dirList.GetFilenames();
282    for( DirListType::const_iterator it = filenames_list.begin(); 
283         it != filenames_list.end(); ++it)
284    {
285       AddFileName( *it );
286    }
287 }
288
289 /**
290  * \brief Sorts the given Fileset
291  * \warning This could be implemented in a 'Strategy Pattern' approach
292  *          But as I don't know how to do it, I leave it this way
293  *          BTW, this is also a Strategy, I don't know this is 
294  *          the best approach :)
295  */
296 void SerieHelper::OrderFileList(FileList *fileSet)
297 {
298    // Only computed during ImagePositionPatientOrdering
299    // (need to sort the FileList using IPP and IOP !)
300    ZSpacing = -1.0;
301    
302    if ( SerieHelper::UserLessThanFunction )
303    {
304       gdcmDebugMacro("Use UserLessThanFunction");
305       UserOrdering( fileSet );
306       return; 
307    }
308    else if ( ImagePositionPatientOrdering( fileSet ) )
309    { 
310       gdcmDebugMacro("ImagePositionPatientOrdering succeeded");
311       return ;
312    }
313    else if ( ImageNumberOrdering(fileSet ) )
314    {
315       gdcmDebugMacro("ImageNumberOrdering succeeded");   
316       return ;
317    }
318    else  
319    {
320       gdcmDebugMacro("Use FileNameOrdering");    
321       FileNameOrdering(fileSet );
322    }
323 }
324
325 /**
326  * \brief Elementary coherence checking of the files with the same Serie UID
327  * Only sizes and pixel type are checked right now ...
328  */ 
329 bool SerieHelper::IsCoherent(FileList *fileSet)
330 {
331    if(fileSet->size() == 1)
332    return true;
333
334    FileList::const_iterator it = fileSet->begin();
335
336    int nX =               (*it)->GetXSize();
337    int nY =               (*it)->GetYSize();
338    int pixelSize =        (*it)->GetPixelSize();
339    bool signedPixelData = (*it)->IsSignedPixelData();
340    it ++;
341    for ( ;
342          it != fileSet->end();
343        ++it)
344    {
345       if ( (*it)->GetXSize() != nX )
346          return false;
347       if ( (*it)->GetYSize() != nY )
348          return false;
349       if ( (*it)->GetPixelSize() != pixelSize )
350          return false;
351       if ( (*it)->IsSignedPixelData() != signedPixelData )
352          return false;
353       // probabely more is to be checked (?)
354    }
355    return true;
356 }
357
358 //#ifndef GDCM_LEGACY_REMOVE
359 /**
360  * \brief   accessor (DEPRECATED :  use GetFirstSingleSerieUIDFileSet )
361  *          Warning : 'coherent' means here they have the same Serie UID
362  * @return  The first FileList if found, otherwhise NULL
363  */
364  /*
365 FileList *SerieHelper::GetFirstCoherentFileList()
366 {
367    ItFileSetHt = SingleSerieUIDFileSetHT.begin();
368    if ( ItFileSetHt != SingleSerieUIDFileSetHT.end() )
369       return ItFileSetHt->second;
370    return NULL;
371 }
372 */
373 /**
374  * \brief   accessor (DEPRECATED :  use GetNextSingleSerieUIDFileSet )
375  *          Warning : 'coherent' means here they have the same Serie UID
376  * \note : meaningfull only if GetFirstCoherentFileList() already called 
377  * @return  The next FileList if found, otherwhise NULL
378  */
379  /*
380 FileList *SerieHelper::GetNextCoherentFileList()
381 {
382    gdcmAssertMacro (ItFileSetHt != SingleSerieUIDFileSetHT.end());
383   
384    ++ItFileSetHt;
385    if ( ItFileSetHt != SingleSerieUIDFileSetHT.end() )
386       return ItFileSetHt->second;
387    return NULL;
388 }
389 */
390
391 /**
392  * \brief   accessor (DEPRECATED :  use GetSingleSerieUIDFileSet )
393   *          Warning : 'coherent' means here they have the same Serie UID
394  * @param SerieUID SerieUID
395  * \return  pointer to the FileList if found, otherwhise NULL
396  */
397  /*
398 FileList *SerieHelper::GetCoherentFileList(std::string SerieUID)
399 {
400    if ( SingleSerieUIDFileSetHT.count(SerieUID) == 0 )
401       return 0;     
402    return SingleSerieUIDFileSetHT[SerieUID];
403 }
404 */
405 //#endif
406
407
408 /**
409  * \brief   Get the first Fileset while visiting the SingleSerieUIDFileSetmap
410  * @return  The first FileList (SingleSerieUIDFileSet) if found, otherwhise 0
411  */
412 FileList *SerieHelper::GetFirstSingleSerieUIDFileSet()
413 {
414    ItFileSetHt = SingleSerieUIDFileSetHT.begin();
415    if ( ItFileSetHt != SingleSerieUIDFileSetHT.end() )
416       return ItFileSetHt->second;
417    return NULL;
418 }
419
420 /**
421  * \brief   Get the next Fileset while visiting the SingleSerieUIDFileSetmap
422  * \note : meaningfull only if GetNextSingleSerieUIDFileSet() already called 
423  * @return  The next FileList (SingleSerieUIDFileSet) if found, otherwhise 0
424  */
425 FileList *SerieHelper::GetNextSingleSerieUIDFileSet()
426 {
427    gdcmAssertMacro (ItFileSetHt != SingleSerieUIDFileSetHT.end());
428   
429    ++ItFileSetHt;
430    if ( ItFileSetHt != SingleSerieUIDFileSetHT.end() )
431       return ItFileSetHt->second;
432    return NULL;
433 }
434
435 /**
436  * \brief   Get the SingleSerieUIDFileSet according to its Serie UID
437  * @param SerieUID SerieUID to retrieve
438  * \return pointer to the FileList (SingleSerieUIDFileSet) if found, otherwhise 0
439  */
440 FileList *SerieHelper::GetSingleSerieUIDFileSet(std::string SerieUID)
441 {
442    if ( SingleSerieUIDFileSetHT.count(SerieUID) == 0 )
443       return 0;     
444    return SingleSerieUIDFileSetHT[SerieUID];
445 }
446
447 /**
448  * \brief   Splits a Single SerieUID Fileset according to the Orientations
449  * @param fileSet File Set to be splitted
450  * \return  std::map of 'Xcoherent' File sets
451  */
452
453 XCoherentFileSetmap SerieHelper::SplitOnOrientation(FileList *fileSet)
454 {
455    XCoherentFileSetmap CoherentFileSet;
456
457    int nb = fileSet->size();
458    if (nb == 0 ) {
459       gdcmWarningMacro("Empty FileList passed to SplitOnOrientation");
460       return CoherentFileSet;
461    }
462
463    float iop[6];
464    std::string strOrient;
465    std::ostringstream ossOrient;
466
467    FileList::const_iterator it = fileSet->begin();
468    for ( ;
469          it != fileSet->end();
470        ++it)
471    {     
472       // Information is in :      
473       // 0020 0037 : Image Orientation (Patient) or
474       // 0020 0035 : Image Orientation (RET)
475
476       // Let's build again the 'cosines' string, to be sure of its format      
477       (*it)->GetImageOrientationPatient(iop);
478
479       ossOrient << iop[0];      
480       for (int i = 1; i < 6; i++)
481       {
482         ossOrient << "\\";
483         ossOrient << iop[i]; 
484       }      
485       strOrient = ossOrient.str();
486       ossOrient.str("");
487       if ( CoherentFileSet.count(strOrient) == 0 )
488       {
489          gdcmDebugMacro(" New Orientation :[" << strOrient << "]");
490          // create a File set in 'orientation' position
491          CoherentFileSet[strOrient] = new FileList;
492          gdcmDebugMacro(" CoherentFileSet[strOrient]" << strOrient << "created");
493       }
494       // Current Orientation and DICOM header match; add the file:
495       CoherentFileSet[strOrient]->push_back( (*it) );
496       gdcmDebugMacro(" CoherentFileSet[strOrient]" << "pushed back")    
497    }
498    return CoherentFileSet;
499 }
500
501 /**
502  * \brief   Splits a 'Single SerieUID' Fileset according to the Positions
503  * @param fileSet File Set to be splitted
504  * \return  std::map of 'Xcoherent' File sets
505  */
506
507 XCoherentFileSetmap SerieHelper::SplitOnPosition(FileList *fileSet)
508 {
509    XCoherentFileSetmap CoherentFileSet;
510
511    int nb = fileSet->size();
512    if (nb == 0 ) {
513       gdcmWarningMacro("Empty FileList passed to SplitOnPosition");
514       return CoherentFileSet;
515    }
516    float pos[3];
517    std::string strImPos;  // read on disc
518    std::ostringstream ossPosition;
519    std::string strPosition; // re computed
520    FileList::const_iterator it = fileSet->begin();
521    for ( ;
522          it != fileSet->end();
523        ++it)
524    {     
525       // Information is in :      
526       // 0020,0032 : Image Position Patient
527       // 0020,0030 : Image Position (RET)
528
529       strImPos = (*it)->GetEntryString(0x0020,0x0032);
530       if ( strImPos == GDCM_UNFOUND)
531       {
532          gdcmWarningMacro( "Unfound Image Position Patient (0020,0032)");
533          strImPos = (*it)->GetEntryString(0x0020,0x0030); // For ACR-NEMA images
534          if ( strImPos == GDCM_UNFOUND )
535          {
536             gdcmWarningMacro( "Unfound Image Position (RET) (0020,0030)");
537             // User wants to split on the 'Position'
538             // No 'Position' info found.
539             // We return an empty Htable !
540             return CoherentFileSet;
541          }  
542       }
543
544       if ( sscanf( strImPos.c_str(), "%f \\%f \\%f ", 
545                                               &pos[0], &pos[1], &pos[2]) != 3 )
546       {
547             gdcmWarningMacro( "Wrong number for Position : ["
548                        << strImPos << "]" );
549             return CoherentFileSet;
550       }
551
552       // Let's build again the 'position' string, to be sure of it's format      
553
554       ossPosition << pos[0];      
555       for (int i = 1; i < 3; i++)
556       {
557         ossPosition << "\\";
558         ossPosition << pos[i]; 
559       }      
560       strPosition = ossPosition.str();
561       ossPosition.str("");
562
563       if ( CoherentFileSet.count(strPosition) == 0 )
564       {
565          gdcmDebugMacro(" New Position :[" << strPosition << "]");
566          // create a File set in 'position' position
567          CoherentFileSet[strPosition] = new FileList;
568       }
569       // Current Position and DICOM header match; add the file:
570       CoherentFileSet[strPosition]->push_back( (*it) );
571    }   
572    return CoherentFileSet;
573 }
574
575 /**
576  * \brief   Splits a 'Single SerieUID' File set Coherent according to the
577  *          value of a given Tag
578  * @param fileSet File Set to be splitted
579  * @param   group  group number of the target Element
580  * @param   element element number of the target Element
581  * \return  std::map of 'Xcoherent' File sets
582  */
583
584 XCoherentFileSetmap SerieHelper::SplitOnTagValue(FileList *fileSet,
585                                                uint16_t group, uint16_t element)
586 {
587    XCoherentFileSetmap CoherentFileSet;
588
589    int nb = fileSet->size();
590    if (nb == 0 ) {
591       gdcmWarningMacro("Empty FileList passed to SplitOnTagValue");
592       return CoherentFileSet;
593    }
594
595    std::string strTagValue;  // read on disc
596
597    FileList::const_iterator it = fileSet->begin();
598    //it ++;
599    for ( ;
600          it != fileSet->end();
601        ++it)
602    {     
603       // Information is in :      
604       // 0020,0032 : Image Position Patient
605       // 0020,0030 : Image Position (RET)
606
607       strTagValue = (*it)->GetEntryString(group,element);
608
609       if ( CoherentFileSet.count(strTagValue) == 0 )
610       {
611          gdcmDebugMacro("  :[" << strTagValue << "]");
612          // create a File set in 'position' position
613          CoherentFileSet[strTagValue] = new FileList;
614       }
615       // Current Tag value and DICOM header match; add the file:
616       CoherentFileSet[strTagValue]->push_back( (*it) );
617    }
618    return CoherentFileSet;
619 }
620
621
622 /**
623  * \brief   Splits a 'Single SerieUID' File set Coherent according to the
624  *          value of a given Tag
625  * @param fileSet File Set to be splitted
626  * @param   group  group number of the target Element
627  * @param   element element number of the target Element
628  * \return  std::map of 'Xcoherent' File sets
629  */
630
631 XCoherentFileSetmap SerieHelper::SplitOnTagValueConvertToFloat(FileList *fileSet,
632                                                uint16_t group, uint16_t element)
633 {
634    XCoherentFileSetmap CoherentFileSet;
635
636    int nb = fileSet->size();
637    if (nb == 0 ) {
638       gdcmWarningMacro("Empty FileList passed to SplitOnTagValue");
639       return CoherentFileSet;
640    }
641
642    std::string strTagValue;  // read on disc
643
644    double      dTagValue;
645    int         iTagValue;
646    char        cTagValue[11];
647    std::string sTagValue;
648
649    FileList::const_iterator it = fileSet->begin();
650    //it ++;
651    for ( ;
652          it != fileSet->end();
653        ++it)
654    {
655       /// \TODO : find a trick to create a string whose value follows lexicographical order 
656
657       strTagValue = (*it)->GetEntryString(group,element);
658       dTagValue = atof( strTagValue.c_str());
659       iTagValue = dTagValue*10000.;
660       sprintf(cTagValue, "%010d", iTagValue);
661       strTagValue = cTagValue;
662       
663       if ( CoherentFileSet.count(strTagValue) == 0 )
664       {
665          gdcmDebugMacro("  :[" << strTagValue << "]");
666          // create a File set in 'position' position
667          CoherentFileSet[strTagValue] = new FileList;
668       }
669       // Current Tag value and DICOM header match; add the file:
670       CoherentFileSet[strTagValue]->push_back( (*it) );
671    }
672    return CoherentFileSet;
673 }
674
675 //-----------------------------------------------------------------------------
676 // Protected
677
678 //-----------------------------------------------------------------------------
679 // Private
680 /**
681  * \brief sorts the images, according to their Patient Position.
682  *  As a side effect, it computes the ZSpacing, according to Jolinda Smith's
683  *  algorithm. (get it with double GetZSpacing() !)
684  *  We may order, considering :
685  *   -# Image Position Patient
686  *   -# Image Number
687  *   -# file name
688  *   -# More to come :-)
689  * \note : FileList = std::vector<File* >
690  * @param fileList Coherent File list (same Serie UID) to sort
691  * @return false only if the header is bugged !
692  */
693 bool SerieHelper::ImagePositionPatientOrdering( FileList *fileList )
694 //based on Jolinda Smith's algorithm
695 {
696 //Tags always use the same coordinate system, where "x" is left
697 //to right, "y" is posterior to anterior, and "z" is foot to head (RAH).
698
699    //iop is calculated based on the file file
700    float cosines[6];
701    double normal[3];
702    double ipp[3];
703    double dist;
704    double min = 0, max = 0;
705    bool first = true;
706    ZSpacing = -1.0;  // will be updated if process doesn't fail
707
708    gdcmDebugMacro("============================================DropDuplicatePositions : " << DropDuplicatePositions );
709     
710    std::multimap<double,File *> distmultimap;
711    // Use a multimap to sort the distances from 0,0,0
712    for ( FileList::const_iterator 
713          it = fileList->begin();
714          it != fileList->end(); ++it )
715    {
716       gdcmDebugMacro("deal with " << (*it)->GetFileName() );
717       if ( first ) 
718       {
719          (*it)->GetImageOrientationPatient( cosines );
720
721    // The "Image Orientation Patient" tag gives the direction cosines 
722    // for the rows and columns for the three axes defined above. 
723    // Typical axial slices will have a value 1/0/0/0/1/0: 
724    // rows increase from left to right, 
725    // columns increase from posterior to anterior. This is your everyday
726    // "looking up from the bottom of the head with the eyeballs up" image. 
727    
728    // The "Image Position Patient" tag gives the coordinates of the first
729    // voxel in the image in the "RAH" coordinate system, relative to some
730    // origin.   
731
732    // First, calculate the slice normal from IOP : 
733           
734          // You only have to do this once for all slices in the volume. Next, 
735          // for each slice, calculate the distance along the slice normal 
736          // using the IPP ("Image Position Patient") tag.
737          // ("dist" is initialized to zero before reading the first slice) :
738          normal[0] = cosines[1]*cosines[5] - cosines[2]*cosines[4];
739          normal[1] = cosines[2]*cosines[3] - cosines[0]*cosines[5];
740          normal[2] = cosines[0]*cosines[4] - cosines[1]*cosines[3];
741
742    // For each slice (here : the first), calculate the distance along 
743    // the slice normal using the IPP tag 
744     
745          ipp[0] = (*it)->GetXOrigin();
746          ipp[1] = (*it)->GetYOrigin();
747          ipp[2] = (*it)->GetZOrigin();
748
749          dist = 0;
750          for ( int i = 0; i < 3; ++i )
751          {
752             dist += normal[i]*ipp[i];
753          }
754     
755          gdcmDebugMacro("dist : " << dist);
756          distmultimap.insert(std::pair<const double,File *>(dist, *it));
757
758          max = min = dist;
759          first = false;
760       }
761       else 
762       {
763    // Next, for each slice, calculate the distance along the slice normal
764    // using the IPP tag 
765          ipp[0] = (*it)->GetXOrigin();
766          ipp[1] = (*it)->GetYOrigin();
767          ipp[2] = (*it)->GetZOrigin();
768
769          dist = 0;
770          for ( int i = 0; i < 3; ++i )
771          {
772             dist += normal[i]*ipp[i];
773          }
774
775          distmultimap.insert(std::pair<const double,File *>(dist, *it));
776          gdcmDebugMacro("dist : " << dist);
777          min = (min < dist) ? min : dist;
778          max = (max > dist) ? max : dist;
779       }
780    }
781
782    gdcmDebugMacro("After parsing vector, nb of elements : " << fileList->size() );
783
784    // Find out if min/max are coherent
785    if ( min == max )
786    {
787      gdcmWarningMacro("Looks like all images have the exact same image position. "
788                       << "No PositionPatientOrdering sort performed. "
789                       << "No 'ZSpacing' calculated! ");
790      return false;
791    }
792
793    // Check to see if image shares a common position
794    bool ok = true;
795    for (std::multimap<double, File *>::iterator it2 = distmultimap.begin();
796         it2 != distmultimap.end();
797         ++it2)
798    {
799    
800       gdcmDebugMacro("Check if image shares a common position : " << (*it2).second->GetFileName() );   
801    
802       if (distmultimap.count((*it2).first) != 1)
803       {
804          gdcmWarningMacro("File: ["
805               << ((*it2).second->GetFileName())
806               << "] : more than ONE file at distance: '"
807               << (*it2).first
808               << " (position is not unique!) "
809               << "No PositionPatientOrdering sort performed. "
810               << "No 'ZSpacing' calculated! ");      
811
812          ok = false;
813       }
814    }
815    if (!ok)
816    {
817       if (! DropDuplicatePositions)
818          return false;
819    }
820       
821 // Now, we can calculate Z Spacing as the difference
822 // between the "dist" values for the first two slices.
823
824 // The following (un)-commented out code is let here
825 // to be re-used by whomsoever is interested...
826
827     std::multimap<double, File *>::iterator it5 = distmultimap.begin();
828     double d1 = (*it5).first;
829     it5++;
830     double d2 = (*it5).first;
831     ZSpacing = d1-d2;
832     if (ZSpacing < 0.0)
833        ZSpacing = - ZSpacing;
834
835    fileList->clear();  // doesn't delete list elements, only nodes
836
837 // Acording to user requierement, we sort direct order or reverse order.
838    if (DirectOrder)
839    {  
840       for (std::multimap<double, File *>::iterator it3 = distmultimap.begin();
841            it3 != distmultimap.end();
842            ++it3)
843       {
844          fileList->push_back( (*it3).second );
845          if (DropDuplicatePositions)
846          {
847             // ImagePositionPatientOrdering  wrong duplicates are found ???
848             // --> fixed. See comment
849
850             it3 =  distmultimap.upper_bound((*it3).first); // skip all duplicates
851            // the upper_bound function increments the iterator to the next non-duplicate entry
852            // The for loop iteration also increments the iterator, which causes the code to skip every other image
853            // --> decrement the iterator after the upper_bound function call
854             it3--;
855             if (it3 == distmultimap.end() )  // if last image, stop iterate
856                break;
857          }
858       }
859    }
860    else // user asked for reverse order
861    {
862       std::multimap<double, File *>::const_iterator it4;
863       it4 = distmultimap.end();
864       do
865       {
866          it4--;
867          fileList->push_back( (*it4).second );
868          if (DropDuplicatePositions)  // skip all duplicates
869          {
870             // lower_bound finds the next element that is 
871             // less than or *equal to* the current value!
872             //it4 =  distmultimap.lower_bound((*it4).first);
873    
874            // David Feng's fix
875            std::multimap<double, File *>::const_iterator itPrev = it4;
876            while (itPrev->first == it4->first)
877               --itPrev;
878            it4 = itPrev;
879     
880            if (it4 == distmultimap.begin() ) // if first image, stop iterate
881                break;
882          } 
883       } while (it4 != distmultimap.begin() );
884    }
885
886    distmultimap.clear();
887
888    return true;
889 }
890
891 bool SerieHelper::ImageNumberLessThan(File *file1, File *file2)
892 {
893   return file1->GetImageNumber() < file2->GetImageNumber();
894 }
895
896 bool SerieHelper::ImageNumberGreaterThan(File *file1, File *file2)
897 {
898   return file1->GetImageNumber() > file2->GetImageNumber();
899 }
900
901 /**
902  * \brief sorts the images, according to their Image Number
903  * \note Works only on bona fide files  (i.e image number is a character string
904  *                                      corresponding to an integer)
905  *             within a bona fide serie (i.e image numbers are consecutive)
906  * @param fileList File set (same Serie UID) to sort 
907  * @return false if non bona fide stuff encountered
908  */
909 bool SerieHelper::ImageNumberOrdering(FileList *fileList) 
910 {
911    int min, max, pos;
912    int n = fileList->size();
913
914    FileList::const_iterator it = fileList->begin();
915    min = max = (*it)->GetImageNumber();
916
917    for (; it != fileList->end(); ++it, ++n)
918    {
919       pos = (*it)->GetImageNumber();
920       min = (min < pos) ? min : pos;
921       max = (max > pos) ? max : pos;
922    }
923
924    // Find out if image numbers are coherent (consecutive)
925    if ( min == max || max == 0 || max >= (n+min) )
926    {
927       gdcmWarningMacro( " 'Image numbers' not coherent. "
928                         << " No ImageNumberOrdering sort performed.");
929       return false;
930    }
931    if (DirectOrder)
932       Sort(fileList,SerieHelper::ImageNumberLessThan);
933    else
934       Sort(fileList,SerieHelper::ImageNumberGreaterThan);
935
936    return true;
937 }
938
939 bool SerieHelper::FileNameLessThan(File *file1, File *file2)
940 {
941    return file1->GetFileName() < file2->GetFileName();
942 }
943
944 bool SerieHelper::FileNameGreaterThan(File *file1, File *file2)
945 {
946    return file1->GetFileName() > file2->GetFileName();
947 }
948 /**
949  * \brief sorts the images, according to their File Name
950  * @param fileList Coherent File list (same Serie UID) to sort
951  * @return false only if the header is bugged !
952  */
953 bool SerieHelper::FileNameOrdering(FileList *fileList)
954 {
955    if (DirectOrder) 
956       Sort(fileList,SerieHelper::FileNameLessThan);
957    else
958       Sort(fileList,SerieHelper::FileNameGreaterThan);   
959
960    return true;
961 }
962
963 /**
964  * \brief sorts the images, according to user supplied function
965  * @param fileList Coherent File list (same Serie UID) to sort
966  * @return false only if the header is bugged !
967  */
968 bool SerieHelper::UserOrdering(FileList *fileList)
969 {
970    Sort(fileList,SerieHelper::UserLessThanFunction);   
971    if (!DirectOrder) 
972    {
973       std::reverse(fileList->begin(), fileList->end());
974    }
975    return true;
976 }
977
978 //-----------------------------------------------------------------------------
979 // Print
980 /**
981  * \brief   Canonical printer.
982  */
983 void SerieHelper::Print(std::ostream &os, std::string const &indent)
984 {
985    // For all the Coherent File lists of the GDCM_NAME_SPACE::Serie
986    SingleSerieUIDFileSetmap::iterator itl = SingleSerieUIDFileSetHT.begin();
987    if ( itl == SingleSerieUIDFileSetHT.end() )
988    {
989       gdcmWarningMacro( "No SingleSerieUID File set found" );
990       return;
991    }
992    while (itl != SingleSerieUIDFileSetHT.end())
993    { 
994       os << "Serie UID :[" << itl->first << "]" << std::endl;
995
996       // For all the files of a SingleSerieUID File set
997       for (FileList::iterator it =  (itl->second)->begin();
998                               it != (itl->second)->end(); 
999                             ++it)
1000       {
1001          os << indent << " --- " << (*it)->GetFileName() << std::endl;
1002       }
1003       ++itl;
1004    }
1005 }
1006
1007 void SerieHelper::CreateDefaultUniqueSeriesIdentifier()
1008 {
1009    // If the user requests, additional information can be appended
1010    // to the SeriesUID to further differentiate volumes in the DICOM
1011    // objects being processed.
1012  
1013    // 0020 0011 Series Number
1014    // A scout scan prior to a CT volume scan can share the same
1015    //   SeriesUID, but they will sometimes have a different Series Number
1016    AddRestriction( TagKey(0x0020, 0x0011) );
1017    
1018    // 0018 0024 Sequence Name
1019    // For T1-map and phase-contrast MRA, the different flip angles and
1020    //   directions are only distinguished by the Sequence Name
1021    AddRestriction( TagKey(0x0018, 0x0024) );
1022    
1023    // 0018 0050 Slice Thickness
1024    // On some CT systems, scout scans and subsequence volume scans will
1025    //   have the same SeriesUID and Series Number - YET the slice 
1026    //   thickness will differ from the scout slice and the volume slices.
1027    AddRestriction( TagKey(0x0018, 0x0050));
1028    
1029    // 0028 0010 Rows
1030    // If the 2D images in a sequence don't have the same number of rows,
1031    //   then it is difficult to reconstruct them into a 3D volume.
1032    AddRestriction( TagKey(0x0028, 0x0010));
1033    
1034    // 0028 0011 Columns
1035    // If the 2D images in a sequence don't have the same number of columns,
1036    //   then it is difficult to reconstruct them into a 3D volume.
1037    AddRestriction( TagKey(0x0028, 0x0011));
1038 }
1039
1040 /**
1041  * \brief Heuristics to *try* to build a Serie Identifier that would ensure
1042  *        all the images are coherent.
1043  *
1044  * By default, uses the SeriesUID.  If UseSeriesDetails(true) has been called,
1045  *         then additional identifying information is used.
1046  *  We allow user to add his own critierions, using AddSeriesDetail
1047  *        (he knows more than we do about his images!)
1048  *        ex : in tagging series, the only pertinent tag is
1049  *        0018|1312 [In-plane Phase Encoding Direction] value : ROW/COLUMN
1050  * @param inFile GDCM_NAME_SPACE::File we want to build a Serie Identifier for.
1051  * @return the SeriesIdentifier
1052  */
1053 std::string SerieHelper::CreateUniqueSeriesIdentifier( File *inFile )
1054 {
1055    if( inFile->IsReadable() )
1056    {
1057     // 0020 000e UI REL Series Instance UID
1058     std::string uid = inFile->GetEntryString (0x0020, 0x000e);
1059     std::string id = uid.c_str();
1060     if(m_UseSeriesDetails)
1061       {
1062       for(SerieExRestrictions::iterator it2 = ExRefine.begin();
1063         it2 != ExRefine.end();
1064         ++it2)
1065         {
1066         const ExRule &r = *it2;
1067         std::string s = inFile->GetEntryString( r.group, r.elem );
1068         if( s == GDCM_UNFOUND )
1069           {
1070           s = "";
1071           }
1072         if( id == uid && !s.empty() )
1073           {
1074           id += "."; // add separator
1075           }
1076         id += s;
1077         }
1078       }
1079     // Eliminate non-alnum characters, including whitespace...
1080     //   that may have been introduced by concats.
1081     unsigned int s_size = id.size();
1082     for(unsigned int i=0; i<s_size; i++)
1083     {
1084      while(i<s_size
1085        && !( id[i] == '.' || id[i] == '%' || id[i] == '_'
1086          || (id[i] >= '+' && id[i] <= '-')       
1087          || (id[i] >= 'a' && id[i] <= 'z')
1088          || (id[i] >= '0' && id[i] <= '9')
1089          || (id[i] >= 'A' && id[i] <= 'Z')))
1090       {
1091          id.replace(i, 1, "_");  // ImagePositionPatient related stuff will be more human readable
1092       }
1093    }
1094    // deal with Dicom strings trailing '\0' 
1095     if(s_size && id[s_size-1] == '_')
1096       id.erase(s_size-1, 1);
1097     return id;
1098     }
1099   else // Could not open inFile
1100     {
1101     gdcmWarningMacro("Could not parse series info.");
1102     std::string id = GDCM_UNFOUND;
1103     return id;
1104     }
1105 }
1106
1107 /**
1108  * \brief Allow user to build is own File Identifier (to be able to sort
1109  *        temporal series just as he wants)
1110  *        Criterions will be set with AddSeriesDetail.
1111  *        (Maybe the method should be moved elsewhere 
1112  *       -File class? FileHelper class?-
1113  * @return FileIdentifier (Tokenizable on '%%%'. Hope it's enough !)
1114  */
1115 std::string SerieHelper::CreateUserDefinedFileIdentifier( File *inFile )
1116 {
1117   //     Deal with all user supplied tags.
1118   //      (user knows more than we do about his images!)
1119   
1120    double converted;
1121    std::string id;
1122    std::string s; 
1123    char charConverted[17]; 
1124    
1125    for(SeriesExDetails::iterator it2 = ExDetails.begin();
1126       it2 != ExDetails.end();
1127       ++it2)
1128    {
1129       const ExDetail &r = *it2;
1130       s = inFile->GetEntryString( r.group, r.elem );
1131       if (s == "") // avoid troubles when empty string is found
1132          s = "-";
1133
1134       // User is allowed to ask for 'convertion', to allow further ordering
1135       // e.g : 100 would be *before* 20; 000020.00 vs 00100.00 : OK
1136       if (it2->convert)
1137       {
1138          if ( s != GDCM_UNFOUND) // Don't convert unfound fields !
1139          {
1140             converted = atof(s.c_str());
1141             // probabely something much more complicated is possible, 
1142             // using C++ features
1143             /// \todo check the behaviour when there are >0 and <0 numbers
1144             sprintf(charConverted, "%016.6f",converted);
1145             s = charConverted;
1146          }
1147       }
1148       // Eliminate non-alphanum characters, including whitespace.
1149       unsigned int s_size = s.size();
1150       if(s_size == 0)
1151       { // to avoid further troubles when wild anonymization was performed
1152          s = "a";
1153       }
1154       else
1155       {
1156          for(unsigned int i=0; i<s_size; i++)
1157          {
1158             while(i<s_size
1159                && !( s[i] == '.' || s[i] == '%' || s[i] == '_'
1160                  || (s[i] >= '+' && s[i] <= '-')       
1161                  || (s[i] >= 'a' && s[i] <= 'z')
1162                  || (s[i] >= '0' && s[i] <= '9')
1163                  || (s[i] >= 'A' && s[i] <= 'Z')))
1164             {
1165                s.replace(i, 1, "_");  // ImagePositionPatient related stuff will be more human readable
1166             }
1167          }
1168          // deal with Dicom strings trailing '\0'
1169          if(s[s_size-1] == '_')
1170             s.erase(s_size-1, 1);
1171       }
1172       id += s.c_str();
1173       id += "%%%"; // make the FileIdentifier Tokenizable
1174    }
1175    id += inFile->GetFileName();
1176    id += "%%%"; 
1177    return id;             
1178 }
1179
1180 //-----------------------------------------------------------------------------
1181 // Sort
1182 /**
1183  * \brief   Sort FileList.
1184  */
1185 void SerieHelper::Sort(FileList *fileList, bool (*pt2Func)( File *file1, File *file2) )
1186 {
1187  std::sort(fileList->begin(), fileList->end(), pt2Func );
1188 }
1189
1190 /*
1191 #ifndef GDCM_LEGACY_REMOVE
1192 bool SerieHelper::AddGdcmFile(File* header)
1193 {
1194   return AddFile(header);
1195 }
1196 #endif
1197 */
1198
1199 //-----------------------------------------------------------------------------
1200 } // end namespace gdcm