src/gdcmSerieHeader.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: gdcmSerieHeader.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2005/02/01 10:29:56 $
   7   Version:   $Revision: 1.18 $
   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 #include "gdcmSerieHeader.h"
  20 #include "gdcmDirList.h"
  21 #include "gdcmFile.h"
  22 #include "gdcmDebug.h"
  23
  24 #include <math.h>
  25 #include <algorithm>
  26 #include <vector>
  27
  28 namespace gdcm
  29 {
  30 typedef std::list<File* > GdcmFileList;
  31 typedef std::vector<File* > GdcmFileVector;
  32
  33 //-----------------------------------------------------------------------------
  34 // Constructor / Destructor
  35 /**
  36  * \brief   Constructor from a given SerieHeader
  37  */
  38 SerieHeader::SerieHeader()
  39 {
  40    //TODO
  41    //CoherentGdcmFileList.clear();
  42 }
  43
  44 /**
  45  * \brief   Canonical destructor.
  46  */
  47 SerieHeader::~SerieHeader()
  48 {
  49    // For all the Coherent File lists of the gdcm::Serie
  50    GdcmFileList *l = GetFirstCoherentFileList();
  51    while (l)
  52    {
  53       // For all the files of a Coherent File list
  54       for (GdcmFileList::iterator it = l->begin();
  55                                   it != l->end();
  56                                 ++it)
  57       {
  58          delete *it;
  59       }
  60       l->clear();
  61       l = GetNextCoherentFileList();
  62    }
  63 }
  64
  65 //-----------------------------------------------------------------------------
  66 // Public
  67 /**
  68  * \brief add a gdcm::File to the list corresponding to its Serie UID
  69  * @param   filename Name of the file to deal with
  70  */
  71 void SerieHeader::AddFileName(std::string const &filename)
  72 {
  73    //directly use string and not const char*:
  74    File *header = new File( filename );
  75    if( header->IsReadable() )
  76    {
  77       // 0020 000e UI REL Series Instance UID
  78       std::string uid =  header->GetEntryValue (0x0020, 0x000e);
  79       // if uid == GDCM_UNFOUND then consistently we should find GDCM_UNFOUND
  80       // no need here to do anything special
  81
  82       if ( CoherentGdcmFileListHT.count(uid) == 0 )
  83       {
  84          gdcmVerboseMacro(" New Serie UID :[" << uid << "]");
  85          // create a std::list in 'uid' position
  86          CoherentGdcmFileListHT[uid] = new GdcmFileList;
  87       }
  88       // Current Serie UID and DICOM header seems to match add the file:
  89       CoherentGdcmFileListHT[uid]->push_back( header );
  90    }
  91    else
  92    {
  93       gdcmVerboseMacro("Could not read file: " << filename );
  94       delete header;
  95    }
  96 }
  97
  98 /**
  99  * \brief Sets the root Directory
 100  * @param   dir Name of the directory to deal with
 101  * @param recursive whether we want explore recursively the Directory
 102  */
 103 void SerieHeader::SetDirectory(std::string const &dir, bool recursive)
 104 {
 105    DirList dirList(dir, recursive); // OS specific
 106
 107    DirListType filenames_list = dirList.GetFilenames();
 108    for( DirListType::const_iterator it = filenames_list.begin();
 109         it != filenames_list.end(); ++it)
 110    {
 111       AddFileName( *it );
 112    }
 113 }
 114
 115 /**
 116  * \brief Sorts the given File List
 117  * \warning This could be implemented in a 'Strategy Pattern' approach
 118  *          But as I don't know how to do it, I leave it this way
 119  *          BTW, this is also a Strategy, I don't know this is the best approach :)
 120  */
 121 void SerieHeader::OrderGdcmFileList(GdcmFileList *CoherentGdcmFileList)
 122 {
 123    if( ImagePositionPatientOrdering( CoherentGdcmFileList ) )
 124    {
 125       return ;
 126    }
 127    else if( ImageNumberOrdering(CoherentGdcmFileList ) )
 128    {
 129       return ;
 130    }
 131    else
 132    {
 133       FileNameOrdering(CoherentGdcmFileList );
 134    }
 135 }
 136
 137 /**
 138  * \brief   Get the first List while visiting the CoherentFileListHT
 139  * @return  The first GdcmFileList if found, otherwhise NULL
 140  */
 141  std::list<File* > *SerieHeader::GetFirstCoherentFileList()
 142 // Why doesn't it compile ?!?
 143 //GdcmFileList *SerieHeader::GetFirstCoherentFileList()
 144 {
 145    ItListHt = CoherentGdcmFileListHT.begin();
 146    if( ItListHt != CoherentGdcmFileListHT.end() )
 147       return ItListHt->second;
 148    return NULL;
 149 }
 150
 151 /**
 152  * \brief   Get the next List while visiting the CoherentFileListHT
 153  * \note : meaningfull only if GetFirstCoherentFileList already called
 154  * @return  The next GdcmFileList if found, otherwhise NULL
 155  */
 156 std::list<File* > *SerieHeader::GetNextCoherentFileList()
 157 {
 158    gdcmAssertMacro (ItListHt != CoherentGdcmFileListHT.end());
 159
 160    ++ItListHt;
 161    if ( ItListHt != CoherentGdcmFileListHT.end() )
 162       return ItListHt->second;
 163    return NULL;
 164 }
 165
 166 /**
 167  * \brief   Get the Coherent Files list according to its Serie UID
 168  * @param SerieUID SerieUID
 169  * \return  pointer to the Coherent Filseslist if found, otherwhise NULL
 170  */
 171 GdcmFileList *SerieHeader::GetCoherentFileList(std::string SerieUID)
 172 {
 173    if ( CoherentGdcmFileListHT.count(SerieUID) == 0 )
 174       return 0;
 175    return CoherentGdcmFileListHT[SerieUID];
 176 }
 177
 178
 179 //-----------------------------------------------------------------------------
 180 // Protected
 181
 182 //-----------------------------------------------------------------------------
 183 // Private
 184 /**
 185  * \brief sorts the images, according to their Patient Position
 186  *  We may order, considering :
 187  *   -# Image Position Patient
 188  *   -# Image Number
 189  *   -# More to come :-)
 190  * @param CoherentGdcmFileList Coherent File list (same Serie UID) to sort
 191  * @return false only if the header is bugged !
 192  */
 193 bool SerieHeader::ImagePositionPatientOrdering(
 194                                        GdcmFileList *CoherentGdcmFileList )
 195 //based on Jolinda's algorithm
 196 {
 197    //iop is calculated based on the file file
 198    float cosines[6];
 199    float normal[3];
 200    float ipp[3];
 201    float dist;
 202    float min = 0, max = 0;
 203    bool first = true;
 204    int n=0;
 205    std::vector<float> distlist;
 206
 207    //!\todo rewrite this for loop.
 208    for ( GdcmFileList::const_iterator
 209          it = CoherentGdcmFileList->begin();
 210          it != CoherentGdcmFileList->end(); ++it )
 211    {
 212       if( first )
 213       {
 214          (*it)->GetImageOrientationPatient( cosines );
 215
 216          // You only have to do this once for all slices in the volume. Next,
 217          // for each slice, calculate the distance along the slice normal
 218          // using the IPP tag ("dist" is initialized to zero before reading
 219          // the first slice) :
 220          normal[0] = cosines[1]*cosines[5] - cosines[2]*cosines[4];
 221          normal[1] = cosines[2]*cosines[3] - cosines[0]*cosines[5];
 222          normal[2] = cosines[0]*cosines[4] - cosines[1]*cosines[3];
 223
 224          ipp[0] = (*it)->GetXOrigin();
 225          ipp[1] = (*it)->GetYOrigin();
 226          ipp[2] = (*it)->GetZOrigin();
 227
 228          dist = 0;
 229          for ( int i = 0; i < 3; ++i )
 230          {
 231             dist += normal[i]*ipp[i];
 232          }
 233
 234          if( dist == 0 )
 235          {
 236             return false;
 237          }
 238
 239          distlist.push_back( dist );
 240
 241          max = min = dist;
 242          first = false;
 243       }
 244       else
 245       {
 246          ipp[0] = (*it)->GetXOrigin();
 247          ipp[1] = (*it)->GetYOrigin();
 248          ipp[2] = (*it)->GetZOrigin();
 249
 250          dist = 0;
 251          for ( int i = 0; i < 3; ++i )
 252          {
 253             dist += normal[i]*ipp[i];
 254          }
 255
 256          if( dist == 0 )
 257          {
 258             return false;
 259          }
 260
 261          distlist.push_back( dist );
 262
 263          min = (min < dist) ? min : dist;
 264          max = (max > dist) ? max : dist;
 265       }
 266       ++n;
 267    }
 268
 269    // Then I order the slices according to the value "dist". Finally, once
 270    // I've read in all the slices, I calculate the z-spacing as the difference
 271    // between the "dist" values for the first two slices.
 272    GdcmFileVector CoherentGdcmFileVector(n);
 273    // CoherentGdcmFileVector.reserve( n );
 274    CoherentGdcmFileVector.resize( n );
 275    // gdcmAssertMacro( CoherentGdcmFileVector.capacity() >= n );
 276
 277    float step = (max - min)/(n - 1);
 278    int pos;
 279    n = 0;
 280
 281    //VC++ don't understand what scope is !! it -> it2
 282    for (GdcmFileList::const_iterator it2  = CoherentGdcmFileList->begin();
 283         it2 != CoherentGdcmFileList->end(); ++it2, ++n)
 284    {
 285       //2*n sort algo !!
 286       //Assumption: all files are present (no one missing)
 287       pos = (int)( fabs( (distlist[n]-min)/step) + .5 );
 288
 289       // a Dicom 'Serie' may contain scout views
 290       // and images may have differents directions
 291       // -> More than one may have the same 'pos'
 292       // Sorting has then NO meaning !
 293       if (CoherentGdcmFileVector[pos]==NULL)
 294          CoherentGdcmFileVector[pos] = *it2;
 295       else
 296          return false;
 297    }
 298
 299    CoherentGdcmFileList->clear();  // doesn't delete list elements, only node
 300
 301    //VC++ don't understand what scope is !! it -> it3
 302    for (GdcmFileVector::const_iterator it3  = CoherentGdcmFileVector.begin();
 303         it3 != CoherentGdcmFileVector.end(); ++it3)
 304    {
 305       CoherentGdcmFileList->push_back( *it3 );
 306    }
 307
 308    distlist.clear();
 309    CoherentGdcmFileVector.clear();
 310
 311    return true;
 312 }
 313
 314 /**
 315  * \brief sorts the images, according to their Image Number
 316  * @param CoherentGdcmFileList Coherent File list (same Serie UID) to sort
 317  * @return false only if the header is bugged !
 318  */
 319 bool SerieHeader::ImageNumberOrdering(GdcmFileList *CoherentGdcmFileList)
 320 {
 321    int min, max, pos;
 322    int n = 0;//CoherentGdcmFileList.size() is a O(N) operation
 323    unsigned char *partition;
 324
 325    GdcmFileList::const_iterator it = CoherentGdcmFileList->begin();
 326    min = max = (*it)->GetImageNumber();
 327
 328    for (; it != CoherentGdcmFileList->end(); ++it, ++n)
 329    {
 330       pos = (*it)->GetImageNumber();
 331
 332       //else
 333       min = (min < pos) ? min : pos;
 334       max = (max > pos) ? max : pos;
 335    }
 336
 337    // Find out if sorting worked:
 338    if( min == max || max == 0 || max > (n+min)) return false;
 339
 340    //bzeros(partition, n); //This function is deprecated, better use memset.
 341    partition = new unsigned char[n];
 342    memset(partition, 0, n);
 343
 344    GdcmFileVector CoherentGdcmFileVector(n);
 345
 346    //VC++ don't understand what scope is !! it -> it2
 347    for (GdcmFileList::const_iterator it2 = CoherentGdcmFileList->begin();
 348         it2 != CoherentGdcmFileList->end(); ++it2)
 349    {
 350       pos = (*it2)->GetImageNumber();
 351       CoherentGdcmFileVector[pos - min] = *it2;
 352       partition[pos - min]++;
 353    }
 354
 355    unsigned char mult = 1;
 356    for( int i=0; i<n ; i++ )
 357    {
 358       mult *= partition[i];
 359    }
 360
 361    //VC++ don't understand what scope is !! it -> it3
 362    CoherentGdcmFileList->clear();  // doesn't delete list elements, only node
 363    for ( GdcmFileVector::const_iterator it3 = CoherentGdcmFileVector.begin();
 364          it3 != CoherentGdcmFileVector.end(); ++it3 )
 365    {
 366       CoherentGdcmFileList->push_back( *it3 );
 367    }
 368    CoherentGdcmFileVector.clear();
 369
 370    delete[] partition;
 371
 372    return mult != 0;
 373 }
 374
 375 /**
 376  * \brief sorts the images, according to their File Name
 377  * @param CoherentGdcmFileList Coherent File list (same Serie UID) to sort
 378  * @return false only if the header is bugged !
 379  */
 380 bool SerieHeader::FileNameOrdering(GdcmFileList *)
 381 {
 382    //using the sort
 383    //sort(CoherentGdcmFileList.begin(), CoherentGdcmFileList.end());
 384    return true;
 385 }
 386
 387 //-----------------------------------------------------------------------------
 388 // Print
 389 /**
 390  * \brief   Canonical printer.
 391  */
 392 void SerieHeader::Print()
 393 {
 394    // For all the Coherent File lists of the gdcm::Serie
 395    CoherentFileListmap::iterator itl = CoherentGdcmFileListHT.begin();
 396    if ( itl == CoherentGdcmFileListHT.end() )
 397    {
 398       gdcmVerboseMacro( "No Coherent File list found" );
 399       return;
 400    }
 401    while (itl != CoherentGdcmFileListHT.end())
 402    {
 403       std::cout << "Serie UID :[" << itl->first << "]" << std::endl;
 404
 405       // For all the files of a Coherent File list
 406       for (GdcmFileList::iterator it =  (itl->second)->begin();
 407                                   it != (itl->second)->end();
 408                                 ++it)
 409       {
 410          std::cout << " --- " << (*it)->GetFileName() << std::endl;
 411       }
 412       ++itl;
 413    }
 414 }
 415
 416 //-----------------------------------------------------------------------------
 417 } // end namespace gdcm