src/gdcmSerieHeader.cxx

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