vtk/vtkGdcm4DSplitter.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: vtkGdcm4DSplitter.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2011/04/04 17:01:03 $
   7   Version:   $Revision: 1.6 $
   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 . Choose 'sort' criterion :
  77 --------------------------
  78
  79  - ImagePositionPatient
  80         void setSortOnPosition();
  81  - ImageOrientationPatient
  82        ==> Only in your dreams!
  83        ==> or, please, write a IOP sorter ...
  84  - User choosen tag
  85         ==> WARNING : This one has troubles; do NOT use it, right now!
  86         ==> use setSortOnUserFunction instead
  87         void setSortOnTag(unsigned short sortGroup, unsigned short sortElem);
  88         void setSortConvertToFloat(bool conv)
  89  - UserDefined Function
  90         void setSortOnUserFunction (FoncComp f);
  91  - File name
  92         void setSortOnFileName()
  93
  94 . Execute :
  95 -----------
  96         bool Go();
  97
  98 . Get the result
  99 ----------------
 100
 101  -a single vtkImageData:
 102         vtkImageData *GetImageData();
 103 - a vector of vtkImageData
 104         std::vector<vtkImageData*> *GetImageDataVector();
 105
 106   ===================================================================== */
 107
 108 #include "gdcmSerieHelper.h"
 109
 110 #include "vtkGdcmReader.h"
 111 #include "vtkGdcm4DSplitter.h"
 112 #include <algorithm>
 113 #include "gdcmSerieHelper.h" // for ImagePositionPatientOrdering()
 114 #include <stdlib.h> // for atof
 115
 116  vtkGdcm4DSplitter::vtkGdcm4DSplitter() :
 117                  SplitOnPosition(false), SplitOnOrientation(false), SplitOnTag(false),
 118                  SplitGroup(0), SplitElem(0),
 119
 120                  SortOnPosition(false),  SortOnOrientation(false),  SortOnTag(false),
 121                  SortGroup(0),  SortElem(0), SortConvertToFloat(false),
 122
 123                  Recursive(false), TypeDir(0),
 124                  verbose(false)
 125  {
 126
 127  }
 128
 129  std::vector<vtkImageData*> * vtkGdcm4DSplitter::GetImageDataVector()
 130  {
 131  if (verbose) std::cout << "TypeDir " << TypeDir << std::endl;
 132     if (TypeResult == 2)
 133        return ImageDataVector;
 134     else
 135       if (TypeResult == 1)
 136       {
 137          std::vector<vtkImageData*> *t = new std::vector<vtkImageData*>;
 138          t->push_back( ImageData );
 139          return t;
 140       }
 141       else
 142          return (std::vector<vtkImageData*>*) NULL;
 143  }
 144
 145  vtkImageData *vtkGdcm4DSplitter::GetImageData()
 146  {
 147     if (TypeResult == 1)
 148        return ImageData;
 149     else
 150       if (TypeResult == 1)
 151       {
 152          return (*ImageDataVector)[0];
 153       }
 154       else
 155          return (vtkImageData*) NULL;
 156  }
 157
 158  bool vtkGdcm4DSplitter::setDirName(std::string &dirName)
 159  {
 160     if ( ! GDCM_NAME_SPACE::DirList::IsDirectory(dirName) )
 161     {
 162        std::cout << "[" << dirName << "] is NOT a directory" << std::endl;
 163        return false;
 164     }
 165     DirName = dirName;
 166     TypeDir=1;
 167     return true;
 168  }
 169
 170  bool vtkGdcm4DSplitter::setVectDirName(std::vector<std::string> &vectDirName)
 171  {
 172     int nbDir = vectDirName.size();
 173     for (int iDir=0; iDir<nbDir; iDir++)
 174     {
 175        if ( ! GDCM_NAME_SPACE::DirList::IsDirectory(vectDirName[iDir]) )
 176        {
 177           std::cout << "[" << vectDirName[iDir] << "] is NOT a directory" << std::endl;
 178           return false;
 179        }
 180     }
 181
 182     VectDirName = vectDirName;
 183     TypeDir=2;
 184     return true;
 185  }
 186
 187  bool vtkGdcm4DSplitter::setVectFileName(std::vector<std::string> &vectFileName)
 188  {
 189     if ( vectFileName.size() == 0)
 190     {
 191           std::cout << "[ vectFileName ] : empty list" << std::endl;
 192           return false;
 193     }
 194     VectFileName = vectFileName;
 195     TypeDir=3;
 196     return true;
 197  }
 198
 199  bool vtkGdcm4DSplitter::CompareOnSortTagConvertToFloat(GDCM_NAME_SPACE::File *file1, GDCM_NAME_SPACE::File *file2)
 200  {
 201   /* if (verbose) printf ("%04x %04x\n", this->SortGroup,this->SortElem);
 202      if (verbose) std :: cout << file1->GetEntryString(SortGroup,SortElem).c_str() << " : "
 203                             << atof(file1->GetEntryString(SortGroup,SortElem).c_str())
 204                             << std::endl;
 205 */
 206 //   return atof(file1->GetEntryString(vtkGdcm4DSplitter::SortGroup,vtkGdcm4DSplitter::SortElem).c_str()) < atof(file2->GetEntryString(vtkGdcm4DSplitter::SortGroup,vtkGdcm4DSplitter::SortElem).c_str());
 207    return atof(file1->GetEntryString(SortGroup,SortElem).c_str()) < atof(file2->GetEntryString(SortGroup,SortElem).c_str());
 208  }
 209
 210  bool vtkGdcm4DSplitter::CompareOnSortTag(GDCM_NAME_SPACE::File *file1, GDCM_NAME_SPACE::File *file2)
 211  {
 212    return file1->GetEntryString(vtkGdcm4DSplitter::SortGroup,vtkGdcm4DSplitter::SortElem) < file2->GetEntryString(vtkGdcm4DSplitter::SortGroup,vtkGdcm4DSplitter::SortElem);
 213  }
 214
 215
 216  bool vtkGdcm4DSplitter::Go()
 217  {
 218    if (!SplitOnPosition && !SplitOnOrientation && !SplitOnTag )
 219    {
 220        ///\TODO (?) Throw an exception "Choose Splitting mode before!"
 221        std::cout << "Choose Splitting mode before!" << std::endl;
 222        return false;
 223    }
 224
 225    /// How To :
 226    /*
 227    entree nom de directory / Vecteur de noms?
 228    recursif oui/non
 229    recuperer la liste des gdcm::File*
 230    passer a SerieHelper (pas de check du Serie UID)
 231    set critere de split
 232
 233    trier chaque Coherent file set
 234    passer chacun a un vtkGcdmReader
 235    retourner le (vecteur de) vtkImageData
 236    */
 237
 238    GDCM_NAME_SPACE::SerieHelper *s;
 239    s = GDCM_NAME_SPACE::SerieHelper::New();
 240
 241    GDCM_NAME_SPACE::File *f;
 242    GDCM_NAME_SPACE::DirListType fileNames;
 243
 244    if (TypeDir == 0 )
 245    {
 246       ///\TODO (?) Throw an exception "Set input Directory name(s) / file names  before!"
 247       std::cout << "Set input Directory name(s) / file names  before!" << std::endl;
 248       return false;
 249    }
 250    else if (TypeDir == 1)
 251    {
 252       GDCM_NAME_SPACE::DirList dirlist(DirName, Recursive); // NO recursive exploration
 253       fileNames = dirlist.GetFilenames(); // all the file names
 254    }
 255
 256    else if (TypeDir == 2)
 257    {
 258       int nbDir = VectDirName.size();
 259       GDCM_NAME_SPACE::DirListType tmpFileNames;
 260       for (int iDir=0; iDir<nbDir; iDir++)
 261       {
 262         GDCM_NAME_SPACE::DirList dirlist(VectDirName[iDir], Recursive);
 263         tmpFileNames = dirlist.GetFilenames();
 264         // Concat two std::vector
 265         //vector1.insert( vector1.end(), vector2.begin(), vector2.end() );
 266        fileNames.insert( fileNames.end(), tmpFileNames.begin(), tmpFileNames.end() );
 267       }
 268    }
 269    else if (TypeDir == 3)
 270    {
 271       fileNames=VectFileName;
 272    }
 273
 274    GDCM_NAME_SPACE::FileList *l = new GDCM_NAME_SPACE::FileList; // (set of gdcm::File)
 275    double floatTagvalue;
 276    // Loop on all the gdcm-readable files
 277    for (GDCM_NAME_SPACE::DirListType::iterator it = fileNames.begin();
 278                                     it != fileNames.end();
 279                                   ++it)
 280    {
 281       int maxSize  = 0x7fff;         // load Elements of any length
 282       f = GDCM_NAME_SPACE::File::New();
 283       f->SetMaxSizeLoadEntry(maxSize);
 284       f->SetFileName( *it );
 285       if (f->Load())
 286          l->push_back(f);
 287       else
 288          std::cout << " Fail to load [" <<  *it << "]" << std::endl;
 289    }
 290
 291    GDCM_NAME_SPACE::XCoherentFileSetmap xcm;
 292
 293    if (SplitOnOrientation)
 294    {
 295             s->SetDropDuplicatePositions(false);
 296             xcm = s->SplitOnOrientation(l);
 297    }
 298    else if (SplitOnPosition)
 299    {
 300             s->SetDropDuplicatePositions(true);
 301             xcm = s->SplitOnPosition(l);
 302    }
 303    else if (SplitOnTag)
 304    {
 305          s->SetDropDuplicatePositions(false);
 306
 307          // Crashes if DataElement not found
 308          //std:: cout << GDCM_NAME_SPACE::Global::GetDicts()->GetDefaultPubDict()->GetEntry(groupelem[0],groupelem[1])->GetName() << std::endl;
 309             if ( ! SplitConvertToFloat )
 310                xcm = s->SplitOnTagValue(l, SplitGroup, SplitElem);
 311             else
 312             {
 313                 xcm = s->SplitOnTagValueConvertToFloat(l, SplitGroup, SplitElem);
 314             }
 315    }
 316
 317    if (xcm.size() == 0)
 318    {
 319       if(verbose) std::cout << "Empty XCoherent File Set after 'split' ?!?" << std::endl;
 320       return false;
 321    }
 322    else if (xcm.size() == 1)
 323       TypeResult=1;
 324    else
 325       TypeResult=2;
 326
 327    ImageDataVector = new std::vector<vtkImageData*>;
 328   // vtkGdcmReader *reader = vtkGdcmReader::New(); // move inside the loop, or be clever using vtk!
 329
 330    for (GDCM_NAME_SPACE::XCoherentFileSetmap::iterator i = xcm.begin();
 331                                                   i != xcm.end();
 332                                                 ++i)
 333    {
 334            if (verbose)
 335                std::cout << "--- xCoherentName = [" << (*i).first << "]" << std::endl;
 336    }
 337  // XCoherentFileSetmap map < critère de split, FileList (= std::vector de gdcm::File*) >
 338
 339    for (GDCM_NAME_SPACE::XCoherentFileSetmap::iterator i = xcm.begin();
 340                                                   i != xcm.end();
 341                                                 ++i)
 342    {
 343
 344       vtkGdcmReader *reader = vtkGdcmReader::New(); /// \FIXME : unable to delete!
 345
 346       if (verbose)
 347                std::cout << "==========================================xCoherentName = [" << (*i).first << "]" << std::endl;
 348
 349       if (SortOnPosition)
 350       {
 351               if (verbose) std::cout << "SortOnPosition" << std::endl;
 352               // (will be IPPSorter, in GDCM2)
 353               s->ImagePositionPatientOrdering((*i).second);
 354               if (verbose) std::cout << "out of SortOnPosition" << std::endl;
 355       }
 356
 357       else if (SortOnOrientation)
 358       {
 359               if (verbose) std::cout << "SortOnOrientation" << std::endl;
 360             /// \TODO SortOnOrientation()
 361
 362             // we still miss an algo to sort an Orientation, given by 6 cosines!
 363             //  Anything like this, in GDCM2?
 364             std::cout << "SortOnOrientation : not so easy - I(mage)O(rientation)P(atient)Sorter still missing! -" << std::endl;
 365             // have a look at SerieHelper::SplitOnOrientation() to have an idea of the mess!
 366
 367             //Better sort on the file name, right now...
 368              s->FileNameOrdering((*i).second);
 369       }
 370
 371       else if (SortOnFileName)
 372       {
 373          if (verbose) std::cout << "SortOnFileName" << std::endl;
 374          if (verbose) std::cout << "taille " << ((*i).second)->size() << std::endl;
 375
 376          s->FileNameOrdering((*i).second);
 377          if (verbose) std::cout << "Out of SortOnFileName" << std::endl;
 378       }
 379
 380       else if (SortOnTag)
 381       {
 382          if (verbose) std::cout << "SortOnTag" << std::endl;
 383          printf ("--> %04x %04x\n", SortGroup,SortElem);
 384          std::cout << "Sorry, troubles not solved yet; use SortOnUserFunction, right now!" << std::endl;
 385
 386         /*        ==> WARNING : This one has troubles; do NOT use it, right now!
 387         // a pointer to fonction cannot be casted as a pointer to member function!
 388         // Use SortOnUserFunction, instead!
 389
 390          if ( SortConvertToFloat )
 391             s->SetUserLessThanFunction( reinterpret_cast<bool (*)(gdcm13::File*, gdcm13::File*)>
 392                                                                  ( &vtkGdcm4DSplitter::CompareOnSortTagConvertToFloat));
 393          else
 394             s->SetUserLessThanFunction( reinterpret_cast<bool (*)(gdcm13::File*, gdcm13::File*)>
 395                                                                  ( &vtkGdcm4DSplitter::CompareOnSortTag));
 396
 397          // Anything like this, in GDCM2?
 398          s->UserOrdering((*i).second);
 399         */
 400
 401          //if (verbose) std::cout << "Out of SortOnTag" << std::endl;
 402          std::cout << "NO ordering performed  :-( " << std::endl;
 403       }
 404
 405       else if (SortOnUserFunction)
 406       {
 407           if (verbose) std::cout << "SortOnUserFunction" << std::endl;
 408           s->SetUserLessThanFunction( UserCompareFunction );
 409          // Anything like this, in GDCM2?
 410          s->UserOrdering((*i).second);
 411          if (verbose) std::cout << "Out of SortOnUserFunction" << std::endl;
 412       }
 413
 414        reader->SetCoherentFileList((*i).second);
 415        reader->Update();
 416
 417        /// \TODO : remove the following
 418        //if (verbose) reader->GetOutput()->PrintSelf(std::cout, vtkIndent(2));
 419
 420        ImageDataVector->push_back(reader->GetOutput() );
 421
 422        std::cout << std::endl;
 423    }
 424
 425    //reader->Delete();  // \TODO : fix
 426    s->Delete();
 427    f->Delete();
 428    delete l;
 429
 430    return true;
 431  }
 432