1 /*=========================================================================
4 Module: $RCSfile: vtkGdcm4DSplitter.cxx,v $
6 Date: $Date: 2011/04/11 11:28:31 $
7 Version: $Revision: 1.9 $
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
44 User needs to be aware, *only him* knows want he wants to do.
45 vtkGdcm4DSplitter class does the job for hom
46 (despite its name, it works on 3D or 2D+T images too)
48 User will have to specify some points
54 bool setDirName(std::string &dirName);
55 - a list of directories
56 bool setVectDirName(std::vector<std::string> &vectDirName);
58 bool setVectFileName(std::vector<std::string> &vectFileName);
60 - Recursive directory exploration
61 void setRecursive(bool recursive);
63 . Choose 'split' criterion :
64 ---------------------------
66 - ImagePositionPatient
67 void setSplitOnPosition();
68 - ImageOrientationPatient
69 void setSplitOnOrientation();
71 void setSplitOnTag(unsigned short splitGroup, unsigned short splitElem);
72 void setSplitConvertToFloat(bool conv);
73 - UserDefined Function
74 void setSortOnUserFunction (FoncComp f);
76 for 'true' 3D image sets :
77 - if you want to get a single 3D vtkImageData, use SplitOnOrientation -i.e. no split-
78 - if you want to get a vector of 2D vtkImageData, use SplitOnPosition -i.e. one slice in each 'XCoherent filesite'-
80 for 'true' 4D multi-orientation image sets (i.e. a stack of axial + sagital + coronal images, at different instants ...)
81 --> this is 5D, right?
85 . Choose 'sort' criterion :
86 --------------------------
88 - ImagePositionPatient
89 void setSortOnPosition();
90 - ImageOrientationPatient
91 ==> Only in your dreams!
92 ==> or, please, write a IOP sorter ...
93 - UserDefined Function
94 void setSortOnUserFunction (FoncComp f);
96 void setSortOnFileName()
105 // -a single vtkImageData:
106 // vtkImageData *GetImageData();
107 - a vector of vtkImageData
108 std::vector<vtkImageData*> *GetImageDataVector();
110 ===================================================================== */
112 #include "gdcmSerieHelper.h"
114 #include "vtkGdcmReader.h"
115 #include "vtkGdcm4DSplitter.h"
117 #include "gdcmSerieHelper.h" // for ImagePositionPatientOrdering()
118 #include <stdlib.h> // for atof
120 vtkGdcm4DSplitter::vtkGdcm4DSplitter() :
121 SplitOnPosition(false), SplitOnOrientation(false), SplitOnTag(false),
122 SplitGroup(0), SplitElem(0),
124 SortOnPosition(false), SortOnOrientation(false), SortOnTag(false),
125 SortGroup(0), SortElem(0), SortConvertToFloat(false),
127 Recursive(false), TypeDir(0),
133 std::vector<vtkImageData*> * vtkGdcm4DSplitter::GetImageDataVector()
136 if (verbose) std::cout << "GetImageDataVector : TypeResult " << TypeResult << std::endl;
138 return ImageDataVector;
142 std::vector<vtkImageData*> *t = new std::vector<vtkImageData*>;
143 t->push_back( ImageData );
147 return (std::vector<vtkImageData*>*) NULL;
149 return ImageDataVector;
152 vtkImageData *vtkGdcm4DSplitter::GetImageData()
155 if (verbose) std::cout << "GetImageData : TypeResult " << TypeResult << std::endl;
161 return (*ImageDataVector)[0];
164 return (vtkImageData*) NULL;
166 return (*ImageDataVector)[0];
169 bool vtkGdcm4DSplitter::setDirName(std::string &dirName)
171 if ( ! GDCM_NAME_SPACE::DirList::IsDirectory(dirName) )
173 std::cout << "[" << dirName << "] is NOT a directory" << std::endl;
181 bool vtkGdcm4DSplitter::setVectDirName(std::vector<std::string> &vectDirName)
183 int nbDir = vectDirName.size();
184 for (int iDir=0; iDir<nbDir; iDir++)
186 if ( ! GDCM_NAME_SPACE::DirList::IsDirectory(vectDirName[iDir]) )
188 std::cout << "[" << vectDirName[iDir] << "] is NOT a directory" << std::endl;
193 VectDirName = vectDirName;
198 bool vtkGdcm4DSplitter::setVectFileName(std::vector<std::string> &vectFileName)
200 if ( vectFileName.size() == 0)
202 std::cout << "[ vectFileName ] : empty list" << std::endl;
205 VectFileName = vectFileName;
210 bool vtkGdcm4DSplitter::CompareOnSortTagConvertToFloat(GDCM_NAME_SPACE::File *file1, GDCM_NAME_SPACE::File *file2)
212 /* if (verbose) printf ("%04x %04x\n", this->SortGroup,this->SortElem);
213 if (verbose) std :: cout << file1->GetEntryString(SortGroup,SortElem).c_str() << " : "
214 << atof(file1->GetEntryString(SortGroup,SortElem).c_str())
217 // return atof(file1->GetEntryString(vtkGdcm4DSplitter::SortGroup,vtkGdcm4DSplitter::SortElem).c_str()) < atof(file2->GetEntryString(vtkGdcm4DSplitter::SortGroup,vtkGdcm4DSplitter::SortElem).c_str());
218 return atof(file1->GetEntryString(SortGroup,SortElem).c_str()) < atof(file2->GetEntryString(SortGroup,SortElem).c_str());
221 bool vtkGdcm4DSplitter::CompareOnSortTag(GDCM_NAME_SPACE::File *file1, GDCM_NAME_SPACE::File *file2)
223 return file1->GetEntryString(vtkGdcm4DSplitter::SortGroup,vtkGdcm4DSplitter::SortElem) < file2->GetEntryString(vtkGdcm4DSplitter::SortGroup,vtkGdcm4DSplitter::SortElem);
227 bool vtkGdcm4DSplitter::Go()
229 if (!SplitOnPosition && !SplitOnOrientation && !SplitOnTag )
231 ///\TODO (?) Throw an exception "Choose Splitting mode before!"
232 std::cout << "Choose Splitting mode before!" << std::endl;
238 entree nom de directory / Vecteur de noms?
240 recuperer la liste des gdcm::File*
241 passer a SerieHelper (pas de check du Serie UID)
244 trier chaque Coherent file set
245 passer chacun a un vtkGcdmReader
246 retourner le (vecteur de) vtkImageData
249 GDCM_NAME_SPACE::SerieHelper *s;
250 s = GDCM_NAME_SPACE::SerieHelper::New();
252 GDCM_NAME_SPACE::File *f;
253 GDCM_NAME_SPACE::DirListType fileNames;
257 ///\TODO (?) Throw an exception "Set input Directory name(s) / file names before!"
258 std::cout << "Set input Directory name(s) / file names before!" << std::endl;
261 else if (TypeDir == 1)
263 GDCM_NAME_SPACE::DirList dirlist(DirName, Recursive); // NO recursive exploration
264 fileNames = dirlist.GetFilenames(); // all the file names
267 else if (TypeDir == 2)
269 int nbDir = VectDirName.size();
270 GDCM_NAME_SPACE::DirListType tmpFileNames;
271 for (int iDir=0; iDir<nbDir; iDir++)
273 GDCM_NAME_SPACE::DirList dirlist(VectDirName[iDir], Recursive);
274 tmpFileNames = dirlist.GetFilenames();
275 // Concat two std::vector
276 //vector1.insert( vector1.end(), vector2.begin(), vector2.end() );
277 fileNames.insert( fileNames.end(), tmpFileNames.begin(), tmpFileNames.end() );
280 else if (TypeDir == 3)
282 fileNames=VectFileName;
285 GDCM_NAME_SPACE::FileList *l = new GDCM_NAME_SPACE::FileList; // (set of gdcm::File)
286 double floatTagvalue;
287 // Loop on all the gdcm-readable files
288 for (GDCM_NAME_SPACE::DirListType::iterator it = fileNames.begin();
289 it != fileNames.end();
292 int maxSize = 0x7fff; // load Elements of any length
293 f = GDCM_NAME_SPACE::File::New();
294 f->SetMaxSizeLoadEntry(maxSize);
295 f->SetFileName( *it );
299 std::cout << " Fail to load [" << *it << "]" << std::endl;
302 GDCM_NAME_SPACE::XCoherentFileSetmap xcm;
304 if (SplitOnOrientation)
306 s->SetDropDuplicatePositions(false);
307 xcm = s->SplitOnOrientation(l);
309 else if (SplitOnPosition)
311 s->SetDropDuplicatePositions(true);
312 xcm = s->SplitOnPosition(l);
313 // reorg the std::map xcm according to position // JPR
314 // the key of xcm follows lexicographical order
315 // (that may be different than the 'distance' order)
316 // we have to reorganize it!
317 reorgXCoherentFileSetmap(xcm);
321 s->SetDropDuplicatePositions(false);
323 // Crashes if DataElement not found
324 //std:: cout << GDCM_NAME_SPACE::Global::GetDicts()->GetDefaultPubDict()->GetEntry(groupelem[0],groupelem[1])->GetName() << std::endl;
325 if ( ! SplitConvertToFloat )
326 xcm = s->SplitOnTagValue(l, SplitGroup, SplitElem);
329 xcm = s->SplitOnTagValueConvertToFloat(l, SplitGroup, SplitElem);
335 if(verbose) std::cout << "Empty XCoherent File Set after 'split' ?!?" << std::endl;
338 else if (xcm.size() == 1)
343 ImageDataVector = new std::vector<vtkImageData*>;
344 // vtkGdcmReader *reader = vtkGdcmReader::New(); // move inside the loop, or be clever using vtk!
346 for (GDCM_NAME_SPACE::XCoherentFileSetmap::iterator i = xcm.begin();
351 std::cout << "--- xCoherentName = [" << (*i).first << "]" << std::endl;
354 // XCoherentFileSetmap map < critère de split, FileList (= std::vector de gdcm::File*) >
356 for (GDCM_NAME_SPACE::XCoherentFileSetmap::iterator i = xcm.begin();
361 vtkGdcmReader *reader = vtkGdcmReader::New(); /// \FIXME : unable to delete!
364 std::cout << "==========================================xCoherentName = [" << (*i).first << "]" << std::endl;
368 if (verbose) std::cout << "SortOnPosition" << std::endl;
369 // (will be IPPSorter, in GDCM2)
370 s->ImagePositionPatientOrdering((*i).second);
371 if (verbose) std::cout << "out of SortOnPosition" << std::endl;
374 else if (SortOnOrientation)
376 if (verbose) std::cout << "SortOnOrientation" << std::endl;
377 /// \TODO SortOnOrientation()
379 // we still miss an algo to sort an Orientation, given by 6 cosines!
380 // Anything like this, in GDCM2?
381 std::cout << "SortOnOrientation : not so easy - I(mage)O(rientation)P(atient)Sorter still missing! -" << std::endl;
382 // have a look at SerieHelper::SplitOnOrientation() to have an idea of the mess!
384 //Better sort on the file name, right now...
385 s->FileNameOrdering((*i).second);
388 else if (SortOnFileName)
390 if (verbose) std::cout << "SortOnFileName" << std::endl;
391 if (verbose) std::cout << "taille " << ((*i).second)->size() << std::endl;
393 s->FileNameOrdering((*i).second);
394 if (verbose) std::cout << "Out of SortOnFileName" << std::endl;
399 if (verbose) std::cout << "SortOnTag" << std::endl;
400 printf ("--> %04x %04x\n", SortGroup,SortElem);
401 std::cout << "Sorry, troubles not solved yet; use SortOnUserFunction, right now!" << std::endl;
403 /* ==> WARNING : This one has troubles; do NOT use it, right now!
404 // a pointer to fonction cannot be casted as a pointer to member function!
405 // Use SortOnUserFunction, instead!
407 if ( SortConvertToFloat )
408 s->SetUserLessThanFunction( reinterpret_cast<bool (*)(gdcm13::File*, gdcm13::File*)>
409 ( &vtkGdcm4DSplitter::CompareOnSortTagConvertToFloat));
411 s->SetUserLessThanFunction( reinterpret_cast<bool (*)(gdcm13::File*, gdcm13::File*)>
412 ( &vtkGdcm4DSplitter::CompareOnSortTag));
414 // Anything like this, in GDCM2?
415 s->UserOrdering((*i).second);
418 //if (verbose) std::cout << "Out of SortOnTag" << std::endl;
419 std::cout << "NO ordering performed :-( " << std::endl;
422 else if (SortOnUserFunction)
424 if (verbose) std::cout << "SortOnUserFunction" << std::endl;
425 s->SetUserLessThanFunction( UserCompareFunction );
426 // Anything like this, in GDCM2?
427 s->UserOrdering((*i).second);
428 if (verbose) std::cout << "Out of SortOnUserFunction" << std::endl;
431 reader->SetCoherentFileList((*i).second);
434 /// \TODO : remove the following
436 std::cout << "reader->GetOutput() :" << std::endl;
437 reader->GetOutput()->PrintSelf(std::cout, vtkIndent(2));
440 ImageDataVector->push_back(reader->GetOutput() );
442 std::vector<vtkImageData*>::iterator it;
444 for(it=ImageDataVector->begin(); it!=ImageDataVector->end(); ++it) {
445 std::cout << "-in vtkGdcm4DSplitter--------------------------" << std::endl;
446 (*it)->PrintSelf(std::cout, vtkIndent(2));
448 std::cout << std::endl;
451 //reader->Delete(); // \TODO : fix
460 void vtkGdcm4DSplitter::reorgXCoherentFileSetmap(GDCM_NAME_SPACE::XCoherentFileSetmap &xcm)
463 std::vector<ELEM> vectElem;
469 GDCM_NAME_SPACE::File *file;
475 for (GDCM_NAME_SPACE::XCoherentFileSetmap::iterator i = xcm.begin();
480 std::cout << "--- xCoherentName = [" << (*i).first << "]" << std::endl;
482 e.strIPP = (*i).first;
483 e.file = *(((*i).second)->begin()); // all the gdcm::File of a given xcm item *have* the same IPP; first one is enough
485 vectElem.push_back(e);
487 sortVectElem(&vectElem);
489 // now, create the final std::map !
490 // final_xcm<to_str(e.dist , xcm[e.strIPP]>
492 // check what we need to free !
496 std::string str_dist;
497 int lgr=vectElem.size();
498 GDCM_NAME_SPACE::XCoherentFileSetmap final_xcm;
499 for (int i2=0; i2<lgr; i2++)
501 dist = (vectElem[i2].dist*1000);
502 sprintf(c_dist,"%010d",dist);
505 std::cout << "dist " << vectElem[i2].dist
506 << " str_dist " << str_dist
507 << " IPP " << vectElem[i2].strIPP
511 final_xcm[str_dist] = xcm[vectElem[i2].strIPP];
514 /// \TODO : check what needs to be cleared // JPR
521 bool vtkGdcm4DSplitter::sortVectElem(std::vector<ELEM> *fileList)
523 //based on Jolinda Smith's algorithm
525 //Tags always use the same coordinate system, where "x" is left
526 //to right, "y" is posterior to anterior, and "z" is foot to head (RAH).
528 //iop is calculated based on the file file
533 double min = 0, max = 0;
536 double ZSpacing; // useless here! // JPR
537 bool DirectOrder = true; // remove it!
539 ZSpacing = -1.0; // will be updated if process doesn't fail
541 //std::multimap<double,File *> distmultimap; // JPR
542 std::multimap<double,ELEM> distmultimap; // JPR
544 // Use a multimap to sort the distances from 0,0,0
545 //for ( FileList::const_iterator // JPR
546 for ( std::vector<ELEM>::iterator // JPR
547 it = fileList->begin();
548 it != fileList->end(); ++it )
550 //gdcmDebugMacro("deal with " << (*it)->file->GetFileName() );
553 (*it).file->GetImageOrientationPatient( cosines );
555 // The "Image Orientation Patient" tag gives the direction cosines
556 // for the rows and columns for the three axes defined above.
557 // Typical axial slices will have a value 1/0/0/0/1/0:
558 // rows increase from left to right,
559 // columns increase from posterior to anterior. This is your everyday
560 // "looking up from the bottom of the head with the eyeballs up" image.
562 // The "Image Position Patient" tag gives the coordinates of the first
563 // voxel in the image in the "RAH" coordinate system, relative to some
566 // First, calculate the slice normal from IOP :
568 // You only have to do this once for all slices in the volume. Next,
569 // for each slice, calculate the distance along the slice normal
570 // using the IPP ("Image Position Patient") tag.
571 // ("dist" is initialized to zero before reading the first slice) :
573 normal[0] = cosines[1]*cosines[5] - cosines[2]*cosines[4];
574 normal[1] = cosines[2]*cosines[3] - cosines[0]*cosines[5];
575 normal[2] = cosines[0]*cosines[4] - cosines[1]*cosines[3];
577 // For each slice (here : the first), calculate the distance along
578 // the slice normal using the IPP tag
580 ipp[0] = (*it).file->GetXOrigin();
581 ipp[1] = (*it).file->GetYOrigin();
582 ipp[2] = (*it).file->GetZOrigin();
585 for ( int i = 0; i < 3; ++i )
587 dist += normal[i]*ipp[i];
590 //gdcmDebugMacro("dist : " << dist);
591 distmultimap.insert(std::pair<const double,ELEM>(dist, *it));
598 // Next, for each slice, calculate the distance along the slice normal
600 ipp[0] = (*it).file->GetXOrigin();
601 ipp[1] = (*it).file->GetYOrigin();
602 ipp[2] = (*it).file->GetZOrigin();
605 for ( int i = 0; i < 3; ++i )
607 dist += normal[i]*ipp[i];
610 (*it).dist = dist; // JPR
612 distmultimap.insert(std::pair<const double,ELEM>(dist, *it));
613 //gdcmDebugMacro("dist : " << dist);
614 min = (min < dist) ? min : dist;
615 max = (max > dist) ? max : dist;
619 // gdcmDebugMacro("After parsing vector, nb of elements : " << fileList->size() );
621 /*Useless here. // JPR
623 // Find out if min/max are coherent
626 gdcmWarningMacro("Looks like all images have the exact same image position. "
627 << "No PositionPatientOrdering sort performed. "
628 << "No 'ZSpacing' calculated! ");
633 /* Useless here, 'split' already done. // JPR
635 // Check to see if image shares a common position
637 for (std::multimap<double, File *>::iterator it2 = distmultimap.begin();
638 it2 != distmultimap.end();
642 gdcmDebugMacro("Check if image shares a common position : " << ((*it2).second).file->GetFileName() );
644 if (distmultimap.count((*it2).first) != 1)
646 gdcmWarningMacro("File: ["
647 << ((*it2).second->GetFileName())
648 << "] : more than ONE file at distance: '"
650 << " (position is not unique!) "
651 << "No PositionPatientOrdering sort performed. "
652 << "No 'ZSpacing' calculated! ");
659 if (! DropDuplicatePositions)
665 // Now, we can calculate Z Spacing as the difference
666 // between the "dist" values for the first two slices.
668 // The following (un)-commented out code is let here
669 // to be re-used by whomsoever is interested...
671 std::multimap<double, ELEM>::iterator it5 = distmultimap.begin();
672 double d1 = (*it5).first;
674 double d2 = (*it5).first;
677 ZSpacing = - ZSpacing;
679 fileList->clear(); // doesn't delete list elements, only nodes
681 // Acording to user requierement, we sort direct order or reverse order.
684 for (std::multimap<double, ELEM>::iterator it3 = distmultimap.begin();
685 it3 != distmultimap.end();
688 fileList->push_back( (*it3).second );
690 /* useless here! // JPR
692 if (DropDuplicatePositions)
694 // ImagePositionPatientOrdering wrong duplicates are found ???
695 // --> fixed. See comment
697 it3 = distmultimap.upper_bound((*it3).first); // skip all duplicates
698 // the upper_bound function increments the iterator to the next non-duplicate entry
699 // The for loop iteration also increments the iterator, which causes the code to skip every other image
700 // --> decrement the iterator after the upper_bound function call
702 if (it3 == distmultimap.end() ) // if last image, stop iterate
708 else // user asked for reverse order
710 std::multimap<double, ELEM>::const_iterator it4;
711 it4 = distmultimap.end();
715 fileList->push_back( (*it4).second );
717 /* useless here // JPR
719 if (DropDuplicatePositions) // skip all duplicates
721 // lower_bound finds the next element that is
722 // less than or *equal to* the current value!
723 //it4 = distmultimap.lower_bound((*it4).first);
726 std::multimap<double, ELEM>::const_iterator itPrev = it4;
727 while (itPrev->first == it4->first)
731 if (it4 == distmultimap.begin() ) // if first image, stop iterate
735 } while (it4 != distmultimap.begin() );
738 distmultimap.clear();