vtk/vtkGdcm4DSplitter.cxx

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