vtk/vtkGdcm4DSplitter.cxx

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