1 /*=========================================================================
4 Module: $RCSfile: vtkGdcm4DSplitter.cxx,v $
6 Date: $Date: 2011/09/20 16:09:05 $
7 Version: $Revision: 1.15 $
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.
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.
17 =========================================================================*/
19 /* Raisons ne pas utiliser itkImageSeriesReader:
21 On Wed, Feb 16, 2011 at 11:51 AM, Roger Bramon Feixas <rogerbramon@gmail.com>
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).
31 /* ====================================================================
34 3D, 2D+T, 3D+T, n*2D+T, 4D images are not always stored the same way :
35 a single 'Dicom Serie',
36 several 'Dicom series' within a single directory
37 several 'Dicom series' within several directories
38 A 'Dicom Serie' doesn't mean always the same thing :
39 a given Slice along the time
40 a given Volume at a given time
41 Sometimes, an image within a serie is so artefacted than user decides to replace it
44 User needs to be aware, *only him* knows want he wants to do.
45 vtkGdcm4DSplitter class does the job for him
46 (despite its name, it works on 3D or 2D+T images too)
48 ==> To (try to) understand how the file store is organised, user is deeply encourage to use something like:
50 PrintFile dirin=Root_Directory_holding_the_images rec > 1.txt
52 open it with any test editor, and/or (Linux users)
53 grep a_supposed_to_be_string_of_interest 1.txt
56 Aware user will have to specify some points :
62 bool setDirName(std::string &dirName);
63 - a list of directories
64 bool setVectDirName(std::vector<std::string> &vectDirName);
66 bool setVectFileName(std::vector<std::string> &vectFileName);
67 - a list of gdcm::File*
68 bool setVectGdcmFile(std::vector<GDCM_NAME_SPACE::File *> &vectGdcmFile);
70 - Recursive directory exploration
71 void setRecursive(bool recursive);
73 . Choose 'split' criterion :
74 ---------------------------
76 - ImagePositionPatient
77 void setSplitOnPosition();
78 - ImageOrientationPatient
79 void setSplitOnOrientation();
81 void setSplitOnTag(unsigned short splitGroup, unsigned short splitElem);
82 void setSplitConvertToFloat(bool conv);
83 - UserDefined Function
84 void setSortOnUserFunction (FoncComp f);
86 for 'true' 3D image sets :
87 - if you want to get a single 3D vtkImageData, use SplitOnOrientation -i.e. no split-
88 - if you want to get a vector of 2D vtkImageData, use SplitOnPosition -i.e. one slice in each 'XCoherent fileset'-
90 for 'true' 4D multi-orientation image sets (i.e. a stack of n axial + m sagital + p coronal images, at different instants ...)
91 --> this is 5D, right?
92 (almost) nothing done, yet :
94 . Use a first time vtkGdcm4DSplitter with setSplitOnOrientation();
95 . Get the VectGdcmFileLists (a std::vector of 'XCoherent fileset')
96 . use vtkGdcm4DSplitter, with as many setVectGdcmFile(std::vector<GDCM_NAME_SPACE::File *> &vectGdcmFile) you need
97 one per element of std::vector<GDCM_NAME_SPACE::File *>
98 think on 'spliting' and 'sorting' it, according to your needs.
100 . Choose 'sort' criterion :
101 --------------------------
103 - ImagePositionPatient
104 void setSortOnPosition();
105 - ImageOrientationPatient
106 ==> Only in your dreams!
107 ==> or, please, write a IOP sorter ...
108 - UserDefined Function
109 void setSortOnUserFunction (FoncComp f);
111 void setSortOnFileName()
113 - Choose 'flip' option
114 ----------------------
115 VTK reader 'flop upside down' the images to display them correctly ?!?
116 When a 'math guy' gets the image, the first thing he does is flipping it again.
117 Avoid loosing CPU time using :
118 void SetFlipY(bool); // Wrong = no flip
127 // -a single vtkImageData:
128 // vtkImageData *GetImageData();
129 - a vector of vtkImageData
130 std::vector<vtkImageData*> *GetImageDataVector();
132 ===================================================================== */
134 #include "vtkGdcmReader.h"
135 #include "vtkGdcm4DSplitter.h"
137 #include "gdcmSerieHelper.h" // for ImagePositionPatientOrdering()
138 #include <stdlib.h> // for atof
140 // Constructor / Destructor
142 * \brief Constructor from a given vtkGdcm4DSplitter
144 vtkGdcm4DSplitter::vtkGdcm4DSplitter() :
145 SplitOnPosition(false), SplitOnOrientation(false), SplitOnTag(false),
146 SplitGroup(0), SplitElem(0), SplitConvertToFloat(false), SplitOnly(false),
148 SortOnPosition(false), SortOnOrientation(false), /* SortOnTag(false),*/ SortOnFileName(false), SortOnUserFunction(false),
149 /*SortGroup(0), SortElem(0), SortConvertToFloat(false),*/
151 Recursive(false),FlipY(true), TypeDir(0),
158 * \brief Canonical destructor.
160 vtkGdcm4DSplitter::~vtkGdcm4DSplitter()
162 /// \TODO : delete everything that must be!
164 GDCM_NAME_SPACE::XCoherentFileSetmap::iterator it;
165 std::vector<GDCM_NAME_SPACE::File*>::iterator it2;
167 for ( it = xcm.begin(); // for each std::vector<gdcm::File*>
171 for ( it2 = (*it).second->begin(); // for each gdcm::File*
172 it2 != (*it).second->end();
175 (*it2)->Delete(); // delete gdcm::File
177 delete (*it).second; // delete the now empty std::vector<gdcm::File*>
180 // VectGdcmFileLists is a vector of pointers on elements of xcm : nothing to do!
184 // Locate Data to process
185 // ======================
187 * \brief sets the directories exploration mode
188 * @param recursive whether we want explore recursively the root Directory
190 void vtkGdcm4DSplitter::setRecursive(bool recursive)
196 * \brief Sets the root Directory to get the images from
197 * @param dirName name of the directory to deal with
199 bool vtkGdcm4DSplitter::setDirName(std::string &dirName)
201 if ( ! GDCM_NAME_SPACE::DirList::IsDirectory(dirName) )
203 std::cout << "[" << dirName << "] is NOT a directory" << std::endl;
212 * \brief Sets a list of Directories to get the images from
213 * @param vectDirName vector of directory names to deal with
215 bool vtkGdcm4DSplitter::setVectDirName(std::vector<std::string> &vectDirName)
217 int nbDir = vectDirName.size();
218 for (int iDir=0; iDir<nbDir; iDir++)
220 if ( ! GDCM_NAME_SPACE::DirList::IsDirectory(vectDirName[iDir]) )
222 std::cout << "[" << vectDirName[iDir] << "] is NOT a directory" << std::endl;
227 VectDirName = vectDirName;
233 * \brief Sets a list of files names to read
234 * @param vectFileName vector of file names to deal with
236 bool vtkGdcm4DSplitter::setVectFileName(std::vector<std::string> &vectFileName)
238 if ( vectFileName.size() == 0)
240 std::cout << "[ vectFileName ] : empty list" << std::endl;
243 VectFileName = vectFileName;
249 * \brief Sets an already filled std::vector of gdcm::File *
250 * @param vectGdcmFileName std::vector of gdcm::File *
253 bool vtkGdcm4DSplitter::setVectGdcmFile(GDCM_NAME_SPACE::FileList *vectGdcmFile)
255 if ( vectGdcmFile->size() == 0)
257 std::cout << "[ vectGdcmFile ] : empty list" << std::endl;
261 VectGdcmFile = vectGdcmFile;
270 * \brief asks for splitting Filesets according to the Position
272 void vtkGdcm4DSplitter::setSplitOnPosition()
274 SplitOnPosition=true;
275 SplitOnOrientation=false;
280 * \brief asks for splitting Filesets according to the Orientation
282 void vtkGdcm4DSplitter::setSplitOnOrientation()
284 SplitOnPosition=false;
285 SplitOnOrientation=true;
290 * \brief asks for splitting Filesets according to the value of a given Tag
291 * @param group group number of the target Element
292 * @param element element number of the target Element
294 void vtkGdcm4DSplitter::setSplitOnTag(unsigned short int splitGroup, unsigned short int splitElem)
296 SplitOnPosition=false;
297 SplitOnOrientation=false;
299 SplitGroup=splitGroup;
304 * \brief asks for converting to 'float' the tag values used as a splitting criteria (lexicographic order may not be suitable)
306 void vtkGdcm4DSplitter::setSplitConvertToFloat(bool conv) {SplitConvertToFloat=conv;}
309 * \brief asks for splitting Filesets according to what was asked for (no sorting, no reading data)
311 void vtkGdcm4DSplitter::setSplitOnly(bool s)
321 * \brief asks for IPP sorting each XCoherent gdcm::FILE set
323 void vtkGdcm4DSplitter::setSortOnPosition()
326 SortOnOrientation=false;
328 SortOnFileName=false;
329 SortOnUserFunction=false;
333 // use setSortOnUserFunction, instead!
334 // void setSortOnTag(unsigned short int sortGroup, unsigned short int sortElem)
336 // SortOnPosition=false;
337 // SortOnOrientation=false;
339 // SortOnFileName=false;
340 // SortOnUserFunction=false;
341 // SortGroup=sortGroup; SortElem=sortElem;
346 * \brief sets a user supplied function (comparison)
347 * @param f comparison function
349 void vtkGdcm4DSplitter::setSortOnUserFunction (FoncComp f)
351 UserCompareFunction=f;
352 SortOnPosition=false;
353 SortOnOrientation=false;
355 SortOnFileName=false;
356 SortOnUserFunction=true;
359 // void setSortConvertToFloat(bool conv)
361 // SortConvertToFloat=conv;
365 * \brief asks for sorting each XCoherent gdcm::FILE set, according to the File names
367 void vtkGdcm4DSplitter::setSortOnFileName()
369 SortOnPosition=false;
370 SortOnOrientation=false;
373 SortOnUserFunction=false;
377 * \brief returns a std::vector of gdcm::FileList* (gdcm::FileList : actually, a std::vector of gdcm::File*)
379 std::vector<GDCM_NAME_SPACE::FileList *> *vtkGdcm4DSplitter::GetVectGdcmFileLists()
382 return &VectGdcmFileLists;
384 GDCM_NAME_SPACE::XCoherentFileSetmap::iterator it;
385 for ( it = xcm.begin();
389 VectGdcmFileLists.push_back((*it).second);
391 return &VectGdcmFileLists;
395 * \brief returns a std::vector of [2D/3D, depending on what was passed] vtkImageData*
397 std::vector<vtkImageData*> * vtkGdcm4DSplitter::GetImageDataVector()
401 return ImageDataVector;
405 * \brief when user _knows_ only _one_ vtkImageData* is returned he may be interested in not getting a vector...
407 vtkImageData *vtkGdcm4DSplitter::GetImageData()
411 return (*ImageDataVector)[0];
415 // bool vtkGdcm4DSplitter::CompareOnSortTagConvertToFloat(GDCM_NAME_SPACE::File *file1, GDCM_NAME_SPACE::File *file2)
417 // return atof(file1->GetEntryString(vtkGdcm4DSplitter::SortGroup,vtkGdcm4DSplitter::SortElem).c_str()) < atof(file2->GetEntryString(vtkGdcm4DSplitter::SortGroup,vtkGdcm4DSplitter::SortElem).c_str());
418 // return atof(file1->GetEntryString(SortGroup,SortElem).c_str()) < atof(file2->GetEntryString(SortGroup,SortElem).c_str());
421 // bool vtkGdcm4DSplitter::CompareOnSortTag(GDCM_NAME_SPACE::File *file1, GDCM_NAME_SPACE::File *file2)
423 // return file1->GetEntryString(vtkGdcm4DSplitter::SortGroup,vtkGdcm4DSplitter::SortElem) < file2->GetEntryString(vtkGdcm4DSplitter::SortGroup,vtkGdcm4DSplitter::SortElem);
427 bool vtkGdcm4DSplitter::Go()
431 entree nom de directory / Vecteur de noms?
433 recuperer la liste des gdcm::File*
434 passer a SerieHelper (pas de check du Serie UID)
437 trier chaque Coherent file set
438 passer chacun a un vtkGcdmReader
439 retourner le (vecteur de) vtkImageData
442 if (!SplitOnPosition && !SplitOnOrientation && !SplitOnTag)
444 ///\TODO (?) Throw an exception "Choose Splitting mode before!"
445 std::cout << "Choose Splitting mode before!" << std::endl;
449 GDCM_NAME_SPACE::FileList *l;
451 GDCM_NAME_SPACE::SerieHelper *s;
452 s = GDCM_NAME_SPACE::SerieHelper::New();
455 // Load the gdcm::File* set, according to user's requierements
456 // ------------------------------------------------------------
458 l = getGdcmFileList();
459 std::cout << l->size() << " gdcm::File read" << std::endl;
461 // Split the gdcm::File* set, according to user's requierements
462 // ------------------------------------------------------------
464 if (SplitOnOrientation)
466 s->SetDropDuplicatePositions(false);
467 xcm = s->SplitOnOrientation(l);
469 else if (SplitOnPosition)
471 s->SetDropDuplicatePositions(true);
472 xcm = s->SplitOnPosition(l);
474 // the key of xcm follows lexicographical order
475 // (that may be different than the 'distance' order)
476 // we have to reorganize it!
478 reorgXCoherentFileSetmap(xcm);
482 s->SetDropDuplicatePositions(false);
484 // Crashes if DataElement not found
485 //std:: cout << GDCM_NAME_SPACE::Global::GetDicts()->GetDefaultPubDict()->GetEntry(groupelem[0],groupelem[1])->GetName() << std::endl;
486 if ( ! SplitConvertToFloat )
487 xcm = s->SplitOnTagValue(l, SplitGroup, SplitElem);
490 xcm = s->SplitOnTagValueConvertToFloat(l, SplitGroup, SplitElem);
497 std::cout << "Empty XCoherent File Set after 'split' ?!?" << std::endl;
502 std::cout << xcm.size() << " XCoherent entries found" << std::endl;
505 // put here, to avoid segfault when unaware user sets SplitOnly to true, and the asks for ImageDataVector
506 ImageDataVector = new std::vector<vtkImageData*>;
508 std::cout <<"SplitOnly " << SplitOnly << std::endl;
513 // ------------------------------------------------------------
515 // ImageDataVector = new std::vector<vtkImageData*>;
517 /// \TODO move inside the loop, or be clever using vtk!
518 vtkGdcmReader *reader = vtkGdcmReader::New(); // move inside the loop, or be clever using vtk!
520 // XCoherentFileSetmap map < critère de split, FileList (= std::vector de gdcm::File*) >
522 for (GDCM_NAME_SPACE::XCoherentFileSetmap::iterator i = xcm.begin();
526 //vtkGdcmReader *reader = vtkGdcmReader::New(); /// \TODO FIXME : unable to delete!
528 reader->SetFlipY(FlipY);
529 // better user SetFileLowerLeft()
530 /// \TODO : modify vtkGdcmReader !
532 std::cout << " --- xCoherentName = [" << (*i).first << "]" << std::endl;
536 if (verbose) std::cout << "SortOnPosition" << std::endl;
537 // (will be IPPSorter, in GDCM2)
538 s->ImagePositionPatientOrdering((*i).second);
539 if (verbose) std::cout << "out of SortOnPosition" << std::endl;
542 else if (SortOnOrientation)
544 if (verbose) std::cout << "SortOnOrientation" << std::endl;
545 /// \TODO (?) SortOnOrientation()
547 // we still miss an algo to sort an Orientation, given by 6 cosines!
548 // Anything like this, in GDCM2?
549 std::cout << "SortOnOrientation : not so easy - I(mage)O(rientation)P(atient)Sorter still missing! -" << std::endl;
550 // have a look at SerieHelper::SplitOnOrientation() to have an idea of the mess!
552 //Better sort on the file name, right now...
553 s->FileNameOrdering((*i).second);
556 else if (SortOnFileName)
558 if (verbose) std::cout << "SortOnFileName" << std::endl;
559 if (verbose) std::cout << "taille " << ((*i).second)->size() << std::endl;
561 s->FileNameOrdering((*i).second);
562 if (verbose) std::cout << "Out of SortOnFileName" << std::endl;
565 // else if (SortOnTag)
567 // if (verbose) std::cout << "SortOnTag" << std::endl;
568 // printf ("--> %04x %04x\n", SortGroup,SortElem);
569 // std::cout << "Sorry, troubles not solved yet; use SortOnUserFunction, right now!" << std::endl;
571 /* ==> WARNING : This one has troubles; do NOT use it, right now!
572 // a pointer to fonction cannot be casted as a pointer to member function!
573 // Use SortOnUserFunction, instead!
575 // if ( SortConvertToFloat )
576 // s->SetUserLessThanFunction( reinterpret_cast<bool (*)(gdcm13::File*, gdcm13::File*)>
577 // ( &vtkGdcm4DSplitter::CompareOnSortTagConvertToFloat));
579 // s->SetUserLessThanFunction( reinterpret_cast<bool (*)(gdcm13::File*, gdcm13::File*)>
580 ( &vtkGdcm4DSplitter::CompareOnSortTag));
582 // Anything like this, in GDCM2?
583 // s->UserOrdering((*i).second);
586 // //if (verbose) std::cout << "Out of SortOnTag" << std::endl;
587 // std::cout << "NO ordering performed :-( " << std::endl;
590 else if (SortOnUserFunction)
592 if (verbose) std::cout << "SortOnUserFunction" << std::endl;
593 s->SetUserLessThanFunction( UserCompareFunction );
594 // Anything like this, in GDCM2?
595 s->UserOrdering((*i).second);
596 if (verbose) std::cout << "Out of SortOnUserFunction" << std::endl;
599 reader->SetCoherentFileList((*i).second);
602 /// \TODO : remove the following
604 // std::cout << "reader->GetOutput() :" << std::endl;
605 // reader->GetOutput()->PrintSelf(std::cout, vtkIndent(2));
608 ImageDataVector->push_back(reader->GetOutput() );
610 std::vector<vtkImageData*>::iterator it;
611 //if (verbose) // JPR
613 for(it=ImageDataVector->begin(); it!=ImageDataVector->end(); ++it) {
614 // std::cout << "-in vtkGdcm4DSplitter --------------------------" << std::endl;
615 // (*it)->PrintSelf(std::cout, vtkIndent(2));
617 //std::cout << std::endl;
620 reader->Delete(); // \TODO : fix
630 * \brief Load the gdcm::File* set, according to user's requierements
631 * returns a std::vector of gdcm::File* (gdcm::FileList : actually, a std::vector of gdcm::File*)
635 GDCM_NAME_SPACE::FileList *vtkGdcm4DSplitter::getGdcmFileList()
638 GDCM_NAME_SPACE::File *f;
639 GDCM_NAME_SPACE::DirListType fileNames;
642 // Fill fileNames with the user supplied file names (in any)
643 // ------------------------------------------------
645 if (TypeDir == 0 ) // Nothing was set as input...
647 ///\TODO (?) Throw an exception "Set input Directory name(s) / file names before!"
648 std::cout << "Set input Directory name(s) / file names before!" << std::endl;
651 else if (TypeDir == 1) // A root directory name was set as input
653 GDCM_NAME_SPACE::DirList dirlist(DirName, Recursive); // NO recursive exploration
654 fileNames = dirlist.GetFilenames(); // all the file names
657 else if (TypeDir == 2) // a std::vector of directory names was set as input
659 int nbDir = VectDirName.size();
660 GDCM_NAME_SPACE::DirListType tmpFileNames;
661 for (int iDir=0; iDir<nbDir; iDir++)
663 GDCM_NAME_SPACE::DirList dirlist(VectDirName[iDir], Recursive);
664 tmpFileNames = dirlist.GetFilenames();
665 // Concat two std::vector
666 //vector1.insert( vector1.end(), vector2.begin(), vector2.end() );
667 fileNames.insert( fileNames.end(), tmpFileNames.begin(), tmpFileNames.end() );
670 else if (TypeDir == 3) // a list of files names was set as input
672 fileNames=VectFileName;
676 // Fill l with the gdcm::File* corresponding to the files
677 // --------------------------------------
680 GDCM_NAME_SPACE::FileList *l = new GDCM_NAME_SPACE::FileList; // (set of gdcm::File*)
682 if (TypeDir == 4) // an already existing std::vector of gdcm::File* was set as input
688 double floatTagvalue;
689 // Loop on all the gdcm-readable files
690 for (GDCM_NAME_SPACE::DirListType::iterator it = fileNames.begin();
691 it != fileNames.end();
694 int maxSize = 0x7fff; // load Elements of any length
695 /// \TODO ? gdcm::File are never free'd
696 f = GDCM_NAME_SPACE::File::New();
697 f->SetMaxSizeLoadEntry(maxSize);
698 f->SetFileName( *it );
702 std::cout << " Fail to load [" << *it << "]" << std::endl;
710 void vtkGdcm4DSplitter::reorgXCoherentFileSetmap(GDCM_NAME_SPACE::XCoherentFileSetmap &xcm)
713 the key of the 'XCoherentFileSetmap', is a std::string, used as if it was found in the Dicom header
714 Normaly(?), it's suitable for almost everything ...
715 ... but the 'Image Position Patient'.
716 We need to order the 'XCoherentFileSetmap' (NOT the content of each XCoherentFileSet!) according to the IPP,
717 using Jolinda Smith's algorithm.
718 (we use a subset of the one defined in gdcm::SerieHelper)
723 std::vector<ELEM> vectElem;
725 /* just to remember :
730 GDCM_NAME_SPACE::File *file;
736 for (GDCM_NAME_SPACE::XCoherentFileSetmap::iterator i = xcm.begin();
741 std::cout << "--- xCoherentName = [" << (*i).first << "]" << std::endl;
743 e.strIPP = (*i).first;
744 e.file = *(((*i).second)->begin()); // all the gdcm::File of a given xcm item *have* the same IPP; first one is enough
746 vectElem.push_back(e);
748 sortVectElem(&vectElem);
750 // now, create the final std::map !
751 // final_xcm<to_str(e.dist , xcm[e.strIPP]>
753 // check what we need to free !
757 std::string str_dist;
758 int lgr=vectElem.size();
759 GDCM_NAME_SPACE::XCoherentFileSetmap final_xcm;
760 for (int i2=0; i2<lgr; i2++)
762 dist = (vectElem[i2].dist*1000);
763 sprintf(c_dist,"%010d",dist);
766 std::cout << "dist " << vectElem[i2].dist
767 << " str_dist " << str_dist
768 << " IPP " << vectElem[i2].strIPP
772 final_xcm[str_dist] = xcm[vectElem[i2].strIPP];
775 /// \TODO : check what needs to be cleared // JPR
782 bool vtkGdcm4DSplitter::sortVectElem(std::vector<ELEM> *fileList)
784 //based on Jolinda Smith's algorithm
786 // NOTE : if you need to use Jolinda Smith's algorithm, get the one inside gdcm::SerieHelper
787 // this one is a light version.
789 //Tags always use the same coordinate system, where "x" is left
790 //to right, "y" is posterior to anterior, and "z" is foot to head (RAH).
792 //iop is calculated based on the file file
797 double min = 0, max = 0;
800 //double ZSpacing; // useless here! // JPR
801 bool DirectOrder = true; // remove it!
803 // ZSpacing = -1.0; // will be updated if process doesn't fail
805 //std::multimap<double,File *> distmultimap; // JPR
806 std::multimap<double,ELEM> distmultimap; // JPR
808 // Use a multimap to sort the distances from 0,0,0
809 //for ( FileList::const_iterator // JPR
810 for ( std::vector<ELEM>::iterator // JPR
811 it = fileList->begin();
812 it != fileList->end(); ++it )
814 //gdcmDebugMacro("deal with " << (*it)->file->GetFileName() );
817 (*it).file->GetImageOrientationPatient( cosines );
819 // The "Image Orientation Patient" tag gives the direction cosines
820 // for the rows and columns for the three axes defined above.
821 // Typical axial slices will have a value 1/0/0/0/1/0:
822 // rows increase from left to right,
823 // columns increase from posterior to anterior. This is your everyday
824 // "looking up from the bottom of the head with the eyeballs up" image.
826 // The "Image Position Patient" tag gives the coordinates of the first
827 // voxel in the image in the "RAH" coordinate system, relative to some
830 // First, calculate the slice normal from IOP :
832 // You only have to do this once for all slices in the volume. Next,
833 // for each slice, calculate the distance along the slice normal
834 // using the IPP ("Image Position Patient") tag.
835 // ("dist" is initialized to zero before reading the first slice) :
837 normal[0] = cosines[1]*cosines[5] - cosines[2]*cosines[4];
838 normal[1] = cosines[2]*cosines[3] - cosines[0]*cosines[5];
839 normal[2] = cosines[0]*cosines[4] - cosines[1]*cosines[3];
841 // For each slice (here : the first), calculate the distance along
842 // the slice normal using the IPP tag
844 ipp[0] = (*it).file->GetXOrigin();
845 ipp[1] = (*it).file->GetYOrigin();
846 ipp[2] = (*it).file->GetZOrigin();
849 for ( int i = 0; i < 3; ++i )
851 dist += normal[i]*ipp[i];
854 //gdcmDebugMacro("dist : " << dist);
855 distmultimap.insert(std::pair<const double,ELEM>(dist, *it));
862 // Next, for each slice, calculate the distance along the slice normal
864 ipp[0] = (*it).file->GetXOrigin();
865 ipp[1] = (*it).file->GetYOrigin();
866 ipp[2] = (*it).file->GetZOrigin();
869 for ( int i = 0; i < 3; ++i )
871 dist += normal[i]*ipp[i];
874 (*it).dist = dist; // JPR
876 distmultimap.insert(std::pair<const double,ELEM>(dist, *it));
877 //gdcmDebugMacro("dist : " << dist);
878 min = (min < dist) ? min : dist;
879 max = (max > dist) ? max : dist;
883 // gdcmDebugMacro("After parsing vector, nb of elements : " << fileList->size() );
885 /*Useless here. // JPR
887 // Find out if min/max are coherent
890 gdcmWarningMacro("Looks like all images have the exact same image position. "
891 << "No PositionPatientOrdering sort performed. "
892 << "No 'ZSpacing' calculated! ");
897 /* Useless here, 'split' already done. // JPR
899 // Check to see if image shares a common position
901 for (std::multimap<double, File *>::iterator it2 = distmultimap.begin();
902 it2 != distmultimap.end();
906 gdcmDebugMacro("Check if image shares a common position : " << ((*it2).second).file->GetFileName() );
908 if (distmultimap.count((*it2).first) != 1)
910 gdcmWarningMacro("File: ["
911 << ((*it2).second->GetFileName())
912 << "] : more than ONE file at distance: '"
914 << " (position is not unique!) "
915 << "No PositionPatientOrdering sort performed. "
916 << "No 'ZSpacing' calculated! ");
923 if (! DropDuplicatePositions)
929 // Now, we can calculate Z Spacing as the difference
930 // between the "dist" values for the first two slices.
932 // The following (un)-commented out code is let here
933 // to be re-used by whomsoever is interested...
935 //std::multimap<double, ELEM>::iterator it5 = distmultimap.begin();
936 //double d1 = (*it5).first;
938 //double d2 = (*it5).first;
940 //if (ZSpacing < 0.0)
941 // ZSpacing = - ZSpacing;
943 fileList->clear(); // doesn't delete list elements, only nodes
945 // Acording to user requierement, we sort direct order or reverse order.
948 for (std::multimap<double, ELEM>::iterator it3 = distmultimap.begin();
949 it3 != distmultimap.end();
952 fileList->push_back( (*it3).second );
954 /* useless here! // JPR
956 if (DropDuplicatePositions)
958 // ImagePositionPatientOrdering wrong duplicates are found ???
959 // --> fixed. See comment
961 it3 = distmultimap.upper_bound((*it3).first); // skip all duplicates
962 // the upper_bound function increments the iterator to the next non-duplicate entry
963 // The for loop iteration also increments the iterator, which causes the code to skip every other image
964 // --> decrement the iterator after the upper_bound function call
966 if (it3 == distmultimap.end() ) // if last image, stop iterate
972 else // user asked for reverse order
974 std::multimap<double, ELEM>::const_iterator it4;
975 it4 = distmultimap.end();
979 fileList->push_back( (*it4).second );
981 /* useless here // JPR
983 if (DropDuplicatePositions) // skip all duplicates
985 // lower_bound finds the next element that is
986 // less than or *equal to* the current value!
987 //it4 = distmultimap.lower_bound((*it4).first);
990 std::multimap<double, ELEM>::const_iterator itPrev = it4;
991 while (itPrev->first == it4->first)
995 if (it4 == distmultimap.begin() ) // if first image, stop iterate
999 } while (it4 != distmultimap.begin() );
1002 distmultimap.clear();