vtk/vtkGdcm4DSplitter.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: vtkGdcm4DSplitter.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2011/04/08 00:11:36 $
   7   Version:   $Revision: 1.8 $
   8
   9   Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
  10   l'Image). All rights reserved. See Doc/License.txt or
  11   http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
  12
  13      This software is distributed WITHOUT ANY WARRANTY; without even
  14      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  15      PURPOSE.  See the above copyright notices for more information.
  16
  17 =========================================================================*/
  18
  19 /* Raisons ne pas utiliser itkImageSeriesReader:
  20
  21 On Wed, Feb 16, 2011 at 11:51 AM, Roger Bramon Feixas <rogerbramon@gmail.com>
  22     Hi,
  23     I'm developing with ITK 3.20 + GDCM 2.0.17 + VTK 5.6 and I've noticed
  24     itkImageSeriesReader is ~2x slower than vtkGDCMImageReader (from GDCM2).
  25     I compared both codes and I think the difference is the extra copy which
  26     itkImageSeriesReader makes from ImageFileReader's output to its own output
  27     (ImageSeriesReader::GenerateData() line 393).
  28 */
  29
  30
  31 /* ====================================================================
  32 vtkGdcm4DSplitter
  33
  34 3D, 2D+T, 3D+T, n*2D+T, 4D images are not always stored the same way :
  35         a single 'Dicom Serie',
  36         several 'Dicom series' within a single directory
  37         several 'Dicom series' within several directories
  38 A 'Dicom Serie' doesn't mean always the same thing :
  39         a given Slice along the time
  40         a given Volume at a given time
  41 Sometimes, an image within a serie is so artefacted than user decides to replace
  42 it by an other image.
  43
  44 User needs to be aware, *only him* knows want he wants to do.
  45 vtkGdcm4DSplitter class does the job for hom
  46 (despite its name, it works on 3D or 2D+T images too)
  47
  48 User will have to specify some points
  49
  50 . Choose input data
  51 -------------------
  52
  53 - a single directory
  54        bool setDirName(std::string &dirName);
  55 - a list of directories
  56        bool setVectDirName(std::vector<std::string> &vectDirName);
  57 - a list of files
  58        bool setVectFileName(std::vector<std::string> &vectFileName);
  59
  60 - Recursive directory exploration
  61        void setRecursive(bool recursive);
  62
  63 . Choose 'split' criterion :
  64 ---------------------------
  65
  66  - ImagePositionPatient
  67         void setSplitOnPosition();
  68  - ImageOrientationPatient
  69         void setSplitOnOrientation();
  70  - User choosen tag
  71         void setSplitOnTag(unsigned short splitGroup, unsigned short splitElem);
  72         void setSplitConvertToFloat(bool conv);
  73  - UserDefined Function
  74         void setSortOnUserFunction (FoncComp f);
  75
  76  for 'true' 3D image sets :
  77    - if you want to get a single 3D vtkImageData, use SplitOnOrientation -i.e. no split-
  78    - if you want to get a vector of 2D vtkImageData, use SplitOnPosition  -i.e. one slice in each 'XCoherent filesite'-
  79
  80  for 'true' 4D multi-orientation image sets (i.e. a stack of axial + sagital + coronal images, at different instants ...)
  81    --> this is 5D, right?
  82    Nothing done, yet.
  83
  84
  85 . Choose 'sort' criterion :
  86 --------------------------
  87
  88  - ImagePositionPatient
  89         void setSortOnPosition();
  90  - ImageOrientationPatient
  91        ==> Only in your dreams!
  92        ==> or, please, write a IOP sorter ...
  93  - UserDefined Function
  94         void setSortOnUserFunction (FoncComp f);
  95  - File name
  96         void setSortOnFileName()
  97
  98 . Execute :
  99 -----------
 100         bool Go();
 101
 102 . Get the result
 103 ----------------
 104
 105 // -a single vtkImageData:
 106 //        vtkImageData *GetImageData();
 107 - a vector of vtkImageData
 108         std::vector<vtkImageData*> *GetImageDataVector();
 109
 110   ===================================================================== */
 111
 112 #include "gdcmSerieHelper.h"
 113
 114 #include "vtkGdcmReader.h"
 115 #include "vtkGdcm4DSplitter.h"
 116 #include <algorithm>
 117 #include "gdcmSerieHelper.h" // for ImagePositionPatientOrdering()
 118 #include <stdlib.h> // for atof
 119
 120  vtkGdcm4DSplitter::vtkGdcm4DSplitter() :
 121                  SplitOnPosition(false), SplitOnOrientation(false), SplitOnTag(false),
 122                  SplitGroup(0), SplitElem(0),
 123
 124                  SortOnPosition(false),  SortOnOrientation(false),  SortOnTag(false),
 125                  SortGroup(0),  SortElem(0), SortConvertToFloat(false),
 126
 127                  Recursive(false), TypeDir(0),
 128                  verbose(false)
 129  {
 130
 131  }
 132
 133  std::vector<vtkImageData*> * vtkGdcm4DSplitter::GetImageDataVector()
 134  {
 135 /*
 136  if (verbose) std::cout << "GetImageDataVector : TypeResult " << TypeResult << std::endl;
 137     if (TypeResult == 2)
 138        return ImageDataVector;
 139     else
 140       if (TypeResult == 1)
 141       {
 142          std::vector<vtkImageData*> *t = new std::vector<vtkImageData*>;
 143          t->push_back( ImageData );
 144          return t;
 145       }
 146       else
 147          return (std::vector<vtkImageData*>*) NULL;
 148 */
 149      return ImageDataVector;
 150  }
 151
 152  vtkImageData *vtkGdcm4DSplitter::GetImageData()
 153  {
 154  /*
 155   if (verbose) std::cout << "GetImageData : TypeResult " << TypeResult << std::endl;
 156     if (TypeResult == 1)
 157        return ImageData;
 158     else
 159       if (TypeResult == 2)
 160       {
 161          return (*ImageDataVector)[0];
 162       }
 163       else
 164          return (vtkImageData*) NULL;
 165 */
 166    return (*ImageDataVector)[0];
 167  }
 168
 169  bool vtkGdcm4DSplitter::setDirName(std::string &dirName)
 170  {
 171     if ( ! GDCM_NAME_SPACE::DirList::IsDirectory(dirName) )
 172     {
 173        std::cout << "[" << dirName << "] is NOT a directory" << std::endl;
 174        return false;
 175     }
 176     DirName = dirName;
 177     TypeDir=1;
 178     return true;
 179  }
 180
 181  bool vtkGdcm4DSplitter::setVectDirName(std::vector<std::string> &vectDirName)
 182  {
 183     int nbDir = vectDirName.size();
 184     for (int iDir=0; iDir<nbDir; iDir++)
 185     {
 186        if ( ! GDCM_NAME_SPACE::DirList::IsDirectory(vectDirName[iDir]) )
 187        {
 188           std::cout << "[" << vectDirName[iDir] << "] is NOT a directory" << std::endl;
 189           return false;
 190        }
 191     }
 192
 193     VectDirName = vectDirName;
 194     TypeDir=2;
 195     return true;
 196  }
 197
 198  bool vtkGdcm4DSplitter::setVectFileName(std::vector<std::string> &vectFileName)
 199  {
 200     if ( vectFileName.size() == 0)
 201     {
 202           std::cout << "[ vectFileName ] : empty list" << std::endl;
 203           return false;
 204     }
 205     VectFileName = vectFileName;
 206     TypeDir=3;
 207     return true;
 208  }
 209
 210  bool vtkGdcm4DSplitter::CompareOnSortTagConvertToFloat(GDCM_NAME_SPACE::File *file1, GDCM_NAME_SPACE::File *file2)
 211  {
 212   /* if (verbose) printf ("%04x %04x\n", this->SortGroup,this->SortElem);
 213      if (verbose) std :: cout << file1->GetEntryString(SortGroup,SortElem).c_str() << " : "
 214                             << atof(file1->GetEntryString(SortGroup,SortElem).c_str())
 215                             << std::endl;
 216 */
 217 //   return atof(file1->GetEntryString(vtkGdcm4DSplitter::SortGroup,vtkGdcm4DSplitter::SortElem).c_str()) < atof(file2->GetEntryString(vtkGdcm4DSplitter::SortGroup,vtkGdcm4DSplitter::SortElem).c_str());
 218    return atof(file1->GetEntryString(SortGroup,SortElem).c_str()) < atof(file2->GetEntryString(SortGroup,SortElem).c_str());
 219  }
 220
 221  bool vtkGdcm4DSplitter::CompareOnSortTag(GDCM_NAME_SPACE::File *file1, GDCM_NAME_SPACE::File *file2)
 222  {
 223    return file1->GetEntryString(vtkGdcm4DSplitter::SortGroup,vtkGdcm4DSplitter::SortElem) < file2->GetEntryString(vtkGdcm4DSplitter::SortGroup,vtkGdcm4DSplitter::SortElem);
 224  }
 225
 226
 227  bool vtkGdcm4DSplitter::Go()
 228  {
 229    if (!SplitOnPosition && !SplitOnOrientation && !SplitOnTag )
 230    {
 231        ///\TODO (?) Throw an exception "Choose Splitting mode before!"
 232        std::cout << "Choose Splitting mode before!" << std::endl;
 233        return false;
 234    }
 235
 236    /// How To :
 237    /*
 238    entree nom de directory / Vecteur de noms?
 239    recursif oui/non
 240    recuperer la liste des gdcm::File*
 241    passer a SerieHelper (pas de check du Serie UID)
 242    set critere de split
 243
 244    trier chaque Coherent file set
 245    passer chacun a un vtkGcdmReader
 246    retourner le (vecteur de) vtkImageData
 247    */
 248
 249    GDCM_NAME_SPACE::SerieHelper *s;
 250    s = GDCM_NAME_SPACE::SerieHelper::New();
 251
 252    GDCM_NAME_SPACE::File *f;
 253    GDCM_NAME_SPACE::DirListType fileNames;
 254
 255    if (TypeDir == 0 )
 256    {
 257       ///\TODO (?) Throw an exception "Set input Directory name(s) / file names  before!"
 258       std::cout << "Set input Directory name(s) / file names  before!" << std::endl;
 259       return false;
 260    }
 261    else if (TypeDir == 1)
 262    {
 263       GDCM_NAME_SPACE::DirList dirlist(DirName, Recursive); // NO recursive exploration
 264       fileNames = dirlist.GetFilenames(); // all the file names
 265    }
 266
 267    else if (TypeDir == 2)
 268    {
 269       int nbDir = VectDirName.size();
 270       GDCM_NAME_SPACE::DirListType tmpFileNames;
 271       for (int iDir=0; iDir<nbDir; iDir++)
 272       {
 273         GDCM_NAME_SPACE::DirList dirlist(VectDirName[iDir], Recursive);
 274         tmpFileNames = dirlist.GetFilenames();
 275         // Concat two std::vector
 276         //vector1.insert( vector1.end(), vector2.begin(), vector2.end() );
 277        fileNames.insert( fileNames.end(), tmpFileNames.begin(), tmpFileNames.end() );
 278       }
 279    }
 280    else if (TypeDir == 3)
 281    {
 282       fileNames=VectFileName;
 283    }
 284
 285    GDCM_NAME_SPACE::FileList *l = new GDCM_NAME_SPACE::FileList; // (set of gdcm::File)
 286    double floatTagvalue;
 287    // Loop on all the gdcm-readable files
 288    for (GDCM_NAME_SPACE::DirListType::iterator it = fileNames.begin();
 289                                     it != fileNames.end();
 290                                   ++it)
 291    {
 292       int maxSize  = 0x7fff;         // load Elements of any length
 293       f = GDCM_NAME_SPACE::File::New();
 294       f->SetMaxSizeLoadEntry(maxSize);
 295       f->SetFileName( *it );
 296       if (f->Load())
 297          l->push_back(f);
 298       else
 299          std::cout << " Fail to load [" <<  *it << "]" << std::endl;
 300    }
 301
 302    GDCM_NAME_SPACE::XCoherentFileSetmap xcm;
 303
 304    if (SplitOnOrientation)
 305    {
 306             s->SetDropDuplicatePositions(false);
 307             xcm = s->SplitOnOrientation(l);
 308    }
 309    else if (SplitOnPosition)
 310    {
 311             s->SetDropDuplicatePositions(true);
 312             xcm = s->SplitOnPosition(l);
 313    }
 314    else if (SplitOnTag)
 315    {
 316          s->SetDropDuplicatePositions(false);
 317
 318          // Crashes if DataElement not found
 319          //std:: cout << GDCM_NAME_SPACE::Global::GetDicts()->GetDefaultPubDict()->GetEntry(groupelem[0],groupelem[1])->GetName() << std::endl;
 320             if ( ! SplitConvertToFloat )
 321                xcm = s->SplitOnTagValue(l, SplitGroup, SplitElem);
 322             else
 323             {
 324                 xcm = s->SplitOnTagValueConvertToFloat(l, SplitGroup, SplitElem);
 325             }
 326    }
 327
 328    if (xcm.size() == 0)
 329    {
 330       if(verbose) std::cout << "Empty XCoherent File Set after 'split' ?!?" << std::endl;
 331       return false;
 332    }
 333    else if (xcm.size() == 1)
 334       TypeResult=1;
 335    else
 336       TypeResult=2;
 337
 338    ImageDataVector = new std::vector<vtkImageData*>;
 339   // vtkGdcmReader *reader = vtkGdcmReader::New(); // move inside the loop, or be clever using vtk!
 340
 341    for (GDCM_NAME_SPACE::XCoherentFileSetmap::iterator i = xcm.begin();
 342                                                   i != xcm.end();
 343                                                 ++i)
 344    {
 345            if (verbose)
 346                std::cout << "--- xCoherentName = [" << (*i).first << "]" << std::endl;
 347    }
 348  // XCoherentFileSetmap map < critère de split, FileList (= std::vector de gdcm::File*) >
 349
 350    for (GDCM_NAME_SPACE::XCoherentFileSetmap::iterator i = xcm.begin();
 351                                                   i != xcm.end();
 352                                                 ++i)
 353    {
 354
 355       vtkGdcmReader *reader = vtkGdcmReader::New(); /// \FIXME : unable to delete!
 356
 357       if (verbose)
 358                std::cout << "==========================================xCoherentName = [" << (*i).first << "]" << std::endl;
 359
 360       if (SortOnPosition)
 361       {
 362               if (verbose) std::cout << "SortOnPosition" << std::endl;
 363               // (will be IPPSorter, in GDCM2)
 364               s->ImagePositionPatientOrdering((*i).second);
 365               if (verbose) std::cout << "out of SortOnPosition" << std::endl;
 366       }
 367
 368       else if (SortOnOrientation)
 369       {
 370               if (verbose) std::cout << "SortOnOrientation" << std::endl;
 371             /// \TODO SortOnOrientation()
 372
 373             // we still miss an algo to sort an Orientation, given by 6 cosines!
 374             //  Anything like this, in GDCM2?
 375             std::cout << "SortOnOrientation : not so easy - I(mage)O(rientation)P(atient)Sorter still missing! -" << std::endl;
 376             // have a look at SerieHelper::SplitOnOrientation() to have an idea of the mess!
 377
 378             //Better sort on the file name, right now...
 379              s->FileNameOrdering((*i).second);
 380       }
 381
 382       else if (SortOnFileName)
 383       {
 384          if (verbose) std::cout << "SortOnFileName" << std::endl;
 385          if (verbose) std::cout << "taille " << ((*i).second)->size() << std::endl;
 386
 387          s->FileNameOrdering((*i).second);
 388          if (verbose) std::cout << "Out of SortOnFileName" << std::endl;
 389       }
 390
 391       else if (SortOnTag)
 392       {
 393          if (verbose) std::cout << "SortOnTag" << std::endl;
 394          printf ("--> %04x %04x\n", SortGroup,SortElem);
 395          std::cout << "Sorry, troubles not solved yet; use SortOnUserFunction, right now!" << std::endl;
 396
 397         /*        ==> WARNING : This one has troubles; do NOT use it, right now!
 398         // a pointer to fonction cannot be casted as a pointer to member function!
 399         // Use SortOnUserFunction, instead!
 400
 401          if ( SortConvertToFloat )
 402             s->SetUserLessThanFunction( reinterpret_cast<bool (*)(gdcm13::File*, gdcm13::File*)>
 403                                                                  ( &vtkGdcm4DSplitter::CompareOnSortTagConvertToFloat));
 404          else
 405             s->SetUserLessThanFunction( reinterpret_cast<bool (*)(gdcm13::File*, gdcm13::File*)>
 406                                                                  ( &vtkGdcm4DSplitter::CompareOnSortTag));
 407
 408          // Anything like this, in GDCM2?
 409          s->UserOrdering((*i).second);
 410         */
 411
 412          //if (verbose) std::cout << "Out of SortOnTag" << std::endl;
 413          std::cout << "NO ordering performed  :-( " << std::endl;
 414       }
 415
 416       else if (SortOnUserFunction)
 417       {
 418           if (verbose) std::cout << "SortOnUserFunction" << std::endl;
 419           s->SetUserLessThanFunction( UserCompareFunction );
 420          // Anything like this, in GDCM2?
 421          s->UserOrdering((*i).second);
 422          if (verbose) std::cout << "Out of SortOnUserFunction" << std::endl;
 423       }
 424
 425        reader->SetCoherentFileList((*i).second);
 426        reader->Update();
 427
 428        /// \TODO : remove the following
 429        if (verbose) {
 430           std::cout << "reader->GetOutput() :" << std::endl;
 431           reader->GetOutput()->PrintSelf(std::cout, vtkIndent(2));
 432        }
 433
 434        ImageDataVector->push_back(reader->GetOutput() );
 435
 436        std::vector<vtkImageData*>::iterator it;
 437        if (verbose)
 438        for(it=ImageDataVector->begin(); it!=ImageDataVector->end(); ++it) {
 439          std::cout << "-in vtkGdcm4DSplitter--------------------------" << std::endl;
 440          (*it)->PrintSelf(std::cout, vtkIndent(2));
 441        }
 442        std::cout << std::endl;
 443    }
 444
 445    //reader->Delete();  // \TODO : fix
 446    s->Delete();
 447    f->Delete();
 448    delete l;
 449
 450    return true;
 451  }
 452