src/gdcmSerieHeader.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: gdcmSerieHeader.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2005/01/24 16:10:53 $
   7   Version:   $Revision: 1.13 $
   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
  31 typedef std::vector<File* > GdcmFileVector;
  32 //-----------------------------------------------------------------------------
  33 // Constructor / Destructor
  34 SerieHeader::SerieHeader()
  35 {
  36    CoherentGdcmFileList.clear();
  37    // Later will contain: 0020 000e UI REL Series Instance UID
  38    CurrentSerieUID = "";
  39 }
  40
  41 SerieHeader::~SerieHeader()
  42 {
  43    /// \todo
  44    for ( GdcmFileList::const_iterator it = CoherentGdcmFileList.begin();
  45          it != CoherentGdcmFileList.end(); ++it)
  46    {
  47       delete *it;
  48    }
  49    CoherentGdcmFileList.clear();
  50 }
  51
  52 //-----------------------------------------------------------------------------
  53 // Print
  54
  55 //-----------------------------------------------------------------------------
  56 // Public
  57 /**
  58  * \brief add a File to the list based on file name
  59  * @param   filename Name of the file to deal with
  60  */
  61 void SerieHeader::AddFileName(std::string const &filename)
  62 {
  63    //directly use string and not const char*:
  64    File *header = new File( filename );
  65    if( header->IsReadable() )
  66    {
  67       // 0020 000e UI REL Series Instance UID
  68       std::string uid =  header->GetEntryValue (0x0020, 0x000e);
  69       // if uid == GDCM_UNFOUND then consistenly we should find GDCM_UNFOUND
  70       // no need here to do anything special
  71       if( CurrentSerieUID == "" )
  72       {
  73          // Set the current one
  74          CurrentSerieUID = uid;
  75       }
  76       if( CurrentSerieUID == uid )
  77       {
  78          // Current Serie UID and DICOM header seems to match add the file:
  79          CoherentGdcmFileList.push_back( header );
  80       }
  81       else
  82       {
  83          gdcmVerboseMacro("Wrong Serie Instance UID should be:" << CurrentSerieUID );
  84       }
  85    }
  86    else
  87    {
  88       gdcmVerboseMacro("Could not read file: " << filename );
  89       delete header;
  90    }
  91 }
  92
  93 /**
  94  * \brief Sets the Directory
  95  * @param   dir Name of the directory to deal with
  96  */
  97 void SerieHeader::SetDirectory(std::string const &dir)
  98 {
  99    CurrentSerieUID = ""; //Reset previous Serie Instance UID
 100    DirList dirList(dir);  //OS specific
 101
 102    DirListType filenames_list = dirList.GetFilenames();
 103    for( DirListType::const_iterator it = filenames_list.begin();
 104         it != filenames_list.end(); ++it)
 105    {
 106       AddFileName( *it );
 107    }
 108 }
 109
 110 /**
 111  * \brief Sorts the File List
 112  * \warning This could be implemented in a 'Strategy Pattern' approach
 113  *          But as I don't know how to do it, I leave it this way
 114  *          BTW, this is also a Strategy, I don't know this is the best approach :)
 115  */
 116 void SerieHeader::OrderGdcmFileList()
 117 {
 118    if( ImagePositionPatientOrdering() )
 119    {
 120       return ;
 121    }
 122    else if( ImageNumberOrdering() )
 123    {
 124       return ;
 125    }
 126    else
 127    {
 128       FileNameOrdering();
 129    }
 130 }
 131
 132 //-----------------------------------------------------------------------------
 133 // Protected
 134
 135 //-----------------------------------------------------------------------------
 136 // Private
 137 /**
 138  * \brief sorts the images, according to their Patient Position
 139  *  We may order, considering :
 140  *   -# Image Number
 141  *   -# Image Position Patient
 142  *   -# More to come :)
 143  * @return false only if the header is bugged !
 144  */
 145 bool SerieHeader::ImagePositionPatientOrdering()
 146 //based on Jolinda's algorithm
 147 {
 148    //iop is calculated based on the file file
 149    float cosines[6];
 150    float normal[3];
 151    float ipp[3];
 152    float dist;
 153    float min = 0, max = 0;
 154    bool first = true;
 155    int n=0;
 156    std::vector<float> distlist;
 157
 158    //!\todo rewrite this for loop.
 159    for ( GdcmFileList::const_iterator
 160          it = CoherentGdcmFileList.begin();
 161          it != CoherentGdcmFileList.end(); ++it )
 162    {
 163       if( first )
 164       {
 165          (*it)->GetImageOrientationPatient( cosines );
 166
 167          //You only have to do this once for all slices in the volume. Next,
 168          // for each slice, calculate the distance along the slice normal
 169          // using the IPP tag ("dist" is initialized to zero before reading
 170          // the first slice) :
 171          normal[0] = cosines[1]*cosines[5] - cosines[2]*cosines[4];
 172          normal[1] = cosines[2]*cosines[3] - cosines[0]*cosines[5];
 173          normal[2] = cosines[0]*cosines[4] - cosines[1]*cosines[3];
 174
 175          ipp[0] = (*it)->GetXOrigin();
 176          ipp[1] = (*it)->GetYOrigin();
 177          ipp[2] = (*it)->GetZOrigin();
 178
 179          dist = 0;
 180          for ( int i = 0; i < 3; ++i )
 181          {
 182             dist += normal[i]*ipp[i];
 183          }
 184
 185          if( dist == 0 )
 186          {
 187             return false;
 188          }
 189
 190          distlist.push_back( dist );
 191
 192          max = min = dist;
 193          first = false;
 194       }
 195       else
 196       {
 197          ipp[0] = (*it)->GetXOrigin();
 198          ipp[1] = (*it)->GetYOrigin();
 199          ipp[2] = (*it)->GetZOrigin();
 200
 201          dist = 0;
 202          for ( int i = 0; i < 3; ++i )
 203          {
 204             dist += normal[i]*ipp[i];
 205          }
 206
 207          if( dist == 0 )
 208          {
 209             return false;
 210          }
 211
 212          distlist.push_back( dist );
 213
 214          min = (min < dist) ? min : dist;
 215          max = (max > dist) ? max : dist;
 216       }
 217       ++n;
 218    }
 219
 220    // Then I order the slices according to the value "dist". Finally, once
 221    // I've read in all the slices, I calculate the z-spacing as the difference
 222    // between the "dist" values for the first two slices.
 223    GdcmFileVector CoherentGdcmFileVector(n);
 224    // CoherentGdcmFileVector.reserve( n );
 225    CoherentGdcmFileVector.resize( n );
 226    // gdcmAssertMacro( CoherentGdcmFileVector.capacity() >= n );
 227
 228    float step = (max - min)/(n - 1);
 229    int pos;
 230    n = 0;
 231
 232    //VC++ don't understand what scope is !! it -> it2
 233    for (GdcmFileList::const_iterator it2  = CoherentGdcmFileList.begin();
 234         it2 != CoherentGdcmFileList.end(); ++it2, ++n)
 235    {
 236       //2*n sort algo !!
 237       //Assumption: all files are present (no one missing)
 238       pos = (int)( fabs( (distlist[n]-min)/step) + .5 );
 239
 240       // a Dicom 'Serie' may contain scout views
 241       // and images may have differents directions
 242       // -> More than one may have the same 'pos'
 243       // Sorting has then NO meaning !
 244       if (CoherentGdcmFileVector[pos]==NULL)
 245          CoherentGdcmFileVector[pos] = *it2;
 246       else
 247          return false;
 248    }
 249
 250    CoherentGdcmFileList.clear();  //this doesn't delete list's element, node only
 251
 252    //VC++ don't understand what scope is !! it -> it3
 253    for (GdcmFileVector::const_iterator it3  = CoherentGdcmFileVector.begin();
 254         it3 != CoherentGdcmFileVector.end(); ++it3)
 255    {
 256       CoherentGdcmFileList.push_back( *it3 );
 257    }
 258
 259    distlist.clear();
 260    CoherentGdcmFileVector.clear();
 261
 262    return true;
 263 }
 264
 265 /**
 266  * \brief sorts the images, according to their Image Number
 267  * @return false only if the header is bugged !
 268  */
 269
 270 bool SerieHeader::ImageNumberOrdering()
 271 {
 272    int min, max, pos;
 273    int n = 0;//CoherentGdcmFileList.size() is a O(N) operation
 274    unsigned char *partition;
 275
 276    GdcmFileList::const_iterator it = CoherentGdcmFileList.begin();
 277    min = max = (*it)->GetImageNumber();
 278
 279    for (; it != CoherentGdcmFileList.end(); ++it, ++n)
 280    {
 281       pos = (*it)->GetImageNumber();
 282
 283       //else
 284       min = (min < pos) ? min : pos;
 285       max = (max > pos) ? max : pos;
 286    }
 287
 288    // Find out if sorting worked:
 289    if( min == max || max == 0 || max > (n+min)) return false;
 290
 291    //bzeros(partition, n); //This function is deprecated, better use memset.
 292    partition = new unsigned char[n];
 293    memset(partition, 0, n);
 294
 295    GdcmFileVector CoherentGdcmFileVector(n);
 296
 297    //VC++ don't understand what scope is !! it -> it2
 298    for (GdcmFileList::const_iterator it2 = CoherentGdcmFileList.begin();
 299         it2 != CoherentGdcmFileList.end(); ++it2)
 300    {
 301       pos = (*it2)->GetImageNumber();
 302       CoherentGdcmFileVector[pos - min] = *it2;
 303       partition[pos - min]++;
 304    }
 305
 306    unsigned char mult = 1;
 307    for( int i=0; i<n ; i++ )
 308    {
 309       mult *= partition[i];
 310    }
 311
 312    //VC++ don't understand what scope is !! it -> it3
 313    CoherentGdcmFileList.clear();  //this doesn't delete list's element, node only
 314    for ( GdcmFileVector::const_iterator it3 = CoherentGdcmFileVector.begin();
 315          it3 != CoherentGdcmFileVector.end(); ++it3 )
 316    {
 317       CoherentGdcmFileList.push_back( *it3 );
 318    }
 319    CoherentGdcmFileVector.clear();
 320
 321    delete[] partition;
 322
 323    return mult != 0;
 324 }
 325
 326
 327 /**
 328  * \brief sorts the images, according to their File Name
 329  * @return false only if the header is bugged !
 330  */
 331 bool SerieHeader::FileNameOrdering()
 332 {
 333    //using the sort
 334    //sort(CoherentGdcmFileList.begin(), CoherentGdcmFileList.end());
 335    return true;
 336 }
 337
 338 } // end namespace gdcm
 339 //-----------------------------------------------------------------------------