vtk/vtkGdcm4DSplitter.cxx

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