src/gdcmHeaderHelper.cxx

   1 // gdcmHeaderHelper.cxx
   2 //-----------------------------------------------------------------------------
   3 #include "gdcmHeaderHelper.h"
   4 #include "gdcmDirList.h"
   5
   6 #include "gdcmDebug.h"
   7 #include <math.h>
   8 #include <algorithm>
   9 #include <vector>
  10
  11 //-----------------------------------------------------------------------------
  12 /**
  13  * \ingroup gdcmHeaderHelper
  14  * \brief   constructor
  15  */
  16 gdcmHeaderHelper::gdcmHeaderHelper() : gdcmHeader( ) {
  17
  18 }
  19
  20 /**
  21  * \ingroup gdcmHeaderHelper
  22  * \brief   constructor
  23  * @param   InFilename Name of the file to deal with
  24  * @param   exception_on_error
  25  * @param   enable_sequences = true to allow the header
  26  *          to be parsed *inside* the SeQuences,
  27  *          when they have an actual length
  28  * @param   ignore_shadow = true if user wants to skip shadow groups
  29  *           during parsing, to save memory space
  30  */
  31 gdcmHeaderHelper::gdcmHeaderHelper(const char *InFilename,
  32                                    bool exception_on_error,
  33                                    bool enable_sequences,
  34                                    bool ignore_shadow)
  35                  : gdcmHeader( InFilename,
  36                                exception_on_error,
  37                                enable_sequences,
  38                                ignore_shadow)
  39 {
  40 }
  41
  42 //-----------------------------------------------------------------------------
  43 // Print
  44
  45 //-----------------------------------------------------------------------------
  46 // Public
  47 /**
  48  * \ingroup gdcmHeaderHelper
  49  * \brief   Returns the size (in bytes) of a single pixel of data.
  50  * @return  The size in bytes of a single pixel of data.
  51  *
  52  */
  53 int gdcmHeaderHelper::GetPixelSize() {
  54
  55      // 0028 0100 US IMG Bits Allocated
  56      // (in order no to be messed up by old RGB images)
  57    if (gdcmHeader::GetEntryByNumber(0x0028,0x0100) == "24")
  58       return 3;
  59
  60    std::string PixelType = GetPixelType();
  61    if (PixelType == "8U"  || PixelType == "8S")
  62       return 1;
  63    if (PixelType == "16U" || PixelType == "16S")
  64       return 2;
  65    if (PixelType == "32U" || PixelType == "32S")
  66       return 4;
  67    if (PixelType == "FD") // to help unfortunate users to manage DOUBLE
  68       return 8;
  69    dbg.Verbose(0, "gdcmHeader::GetPixelSize: Unknown pixel type");
  70    return 0;
  71 }
  72
  73 /**
  74  * \ingroup gdcmHeaderHelper
  75  * \brief   Build the Pixel Type of the image.
  76  *          Possible values are:
  77  *          - 8U  unsigned  8 bit,
  78  *          - 8S    signed  8 bit,
  79  *          - 16U unsigned 16 bit,
  80  *          - 16S   signed 16 bit,
  81  *          - 32U unsigned 32 bit,
  82  *          - 32S   signed 32 bit,
  83  *          - FD    Double,
  84  * \warning 12 bit images appear as 16 bit.
  85  *          24 bit images appear as 8 bit
  86  *          64 bit means 'DOUBLE' images
  87  *                 (no DOUBLE images in kosher DICOM,
  88  *                  but so usefull for people that miss them ;-)
  89  * @return
  90  */
  91 std::string gdcmHeaderHelper::GetPixelType() {
  92    std::string BitsAlloc;
  93    BitsAlloc = GetEntryByNumber(0x0028, 0x0100);
  94    if (BitsAlloc == GDCM_UNFOUND) { // Bits Allocated
  95       dbg.Verbose(0, "gdcmHeader::GetPixelType: unfound Bits Allocated");
  96       BitsAlloc = std::string("16");
  97    }
  98    if (BitsAlloc == "12")           // It will be unpacked
  99       BitsAlloc = std::string("16");
 100    else if (BitsAlloc == "24")      // (in order no to be messed up
 101       BitsAlloc = std::string("8"); // by old RGB images)
 102
 103    std::string Signed;
 104    Signed = GetEntryByNumber(0x0028, 0x0103);
 105    if (Signed == GDCM_UNFOUND) { // "Pixel Representation"
 106       dbg.Verbose(0, "gdcmHeader::GetPixelType: unfound Pixel Representation");
 107       BitsAlloc = std::string("0");
 108    }
 109    if (BitsAlloc == "64") // to help users that want to deal with DOUBLE
 110       return("FD");
 111
 112    if (Signed == "0")
 113       Signed = std::string("U");
 114    else
 115       Signed = std::string("S");
 116
 117    return( BitsAlloc + Signed);
 118 }
 119
 120 /**
 121   * \ingroup gdcmHeaderHelper
 122   * \brief gets the info from 0028,0030 : Pixel Spacing
 123   *             else 1.0
 124   * @return X dimension of a pixel
 125   */
 126 float gdcmHeaderHelper::GetXSpacing() {
 127     float xspacing, yspacing;
 128     std::string StrSpacing = GetEntryByNumber(0x0028,0x0030);
 129
 130    if (StrSpacing == GDCM_UNFOUND) {
 131       dbg.Verbose(0, "gdcmHeader::GetXSpacing: unfound Pixel Spacing (0028,0030)");
 132       return 1.;
 133     }
 134   int nbValues;
 135   if( (nbValues = sscanf( StrSpacing.c_str(), "%f\\%f", &yspacing, &xspacing)) != 2) {
 136     if (nbValues==1)  // if single value is found, xspacing is defaulted to yspacing
 137        return yspacing;
 138   }
 139   if (xspacing == 0.) {
 140     dbg.Verbose(0, "gdcmHeader::GetYSpacing: gdcmData/CT-MONO2-8-abdo.dcm problem");
 141     // seems to be a bug in the header ...
 142     sscanf( StrSpacing.c_str(), "%f\\0\\%f", &yspacing, &xspacing);
 143   }
 144   return xspacing;
 145 }
 146
 147 /**
 148   * \ingroup gdcmHeaderHelper
 149   * \brief gets the info from 0028,0030 : Pixel Spacing
 150   *             else 1.0
 151   * @return Y dimension of a pixel
 152   */
 153 float gdcmHeaderHelper::GetYSpacing() {
 154    float yspacing;
 155    std::string StrSpacing = GetEntryByNumber(0x0028,0x0030);
 156
 157    if (StrSpacing == GDCM_UNFOUND) {
 158       dbg.Verbose(0, "gdcmHeader::GetYSpacing: unfound Pixel Spacing (0028,0030)");
 159       return 1.;
 160     }
 161   sscanf( StrSpacing.c_str(), "%f", &yspacing);
 162   return yspacing;
 163 }
 164
 165 /**
 166   *\ingroup gdcmHeaderHelper
 167   *\brief gets the info from 0018,0088 : Space Between Slices
 168   *                else from 0018,0050 : Slice Thickness
 169    *                else 1.0
 170   * @return Z dimension of a voxel-to be
 171   */
 172 float gdcmHeaderHelper::GetZSpacing() {
 173    // Spacing Between Slices : distance entre le milieu de chaque coupe
 174    // Les coupes peuvent etre :
 175    //   jointives     (Spacing between Slices = Slice Thickness)
 176    //   chevauchantes (Spacing between Slices < Slice Thickness)
 177    //   disjointes    (Spacing between Slices > Slice Thickness)
 178    // Slice Thickness : epaisseur de tissus sur laquelle est acquis le signal
 179    //   ca interesse le physicien de l'IRM, pas le visualisateur de volumes ...
 180    //   Si le Spacing Between Slices est absent,
 181    //   on suppose que les coupes sont jointives
 182
 183    std::string StrSpacingBSlices = GetEntryByNumber(0x0018,0x0088);
 184
 185    if (StrSpacingBSlices == GDCM_UNFOUND) {
 186       dbg.Verbose(0, "gdcmHeader::GetZSpacing: unfound StrSpacingBSlices");
 187       std::string StrSliceThickness = GetEntryByNumber(0x0018,0x0050);
 188       if (StrSliceThickness == GDCM_UNFOUND)
 189          return 1.;
 190       else
 191          // if no 'Spacing Between Slices' is found,
 192          // we assume slices join together
 193          // (no overlapping, no interslice gap)
 194          // if they don't, we're fucked up
 195          return atof(StrSliceThickness.c_str());
 196    } else {
 197       return atof(StrSpacingBSlices.c_str());
 198    }
 199 }
 200
 201 /**
 202   *\ingroup gdcmHeaderHelper
 203   *\brief gets the info from 0028,1052 : Rescale Intercept
 204   * @return Rescale Intercept
 205  */
 206 float gdcmHeaderHelper::GetRescaleIntercept() {
 207   float resInter = 0.;
 208   std::string StrRescInter = GetEntryByNumber(0x0028,0x1052); //0028 1052 DS IMG Rescale Intercept
 209   if (StrRescInter != GDCM_UNFOUND) {
 210       if( sscanf( StrRescInter.c_str(), "%f", &resInter) != 1) {
 211          dbg.Verbose(0, "gdcmHeader::GetRescaleIntercept: Rescale Slope is empty");
 212            // bug in the element 0x0028,0x1052
 213       }
 214    }
 215   return resInter;
 216 }
 217
 218 /**
 219   *\ingroup gdcmHeaderHelper
 220   *\brief gets the info from 0028,1053 : Rescale Slope
 221   * @return Rescale Slope
 222  */
 223  float gdcmHeaderHelper::GetRescaleSlope() {
 224   float resSlope = 1.;
 225   std::string StrRescSlope = GetEntryByNumber(0x0028,0x1053); //0028 1053 DS IMG Rescale Slope
 226   if (StrRescSlope != GDCM_UNFOUND) {
 227       if( sscanf( StrRescSlope.c_str(), "%f", &resSlope) != 1) {
 228          dbg.Verbose(0, "gdcmHeader::GetRescaleSlope: Rescale Slope is empty");
 229            // bug in the element 0x0028,0x1053
 230       }
 231    }
 232    return resSlope;
 233 }
 234
 235 /**
 236   * \ingroup gdcmHeaderHelper
 237   * \brief This function is intended to user who doesn't want
 238   *   to have to manage a LUT and expects to get an RBG Pixel image
 239   *   (or a monochrome one ...)
 240   * \warning to be used with GetImagePixels()
 241   * @return 1 if Gray level, 3 if Color (RGB, YBR or PALETTE COLOR)
 242   */
 243 int gdcmHeaderHelper::GetNumberOfScalarComponents() {
 244    if (GetSamplesPerPixel() ==3)
 245       return 3;
 246
 247      // 0028 0100 US IMG Bits Allocated
 248      // (in order no to be messed up by old RGB images)
 249    if (gdcmHeader::GetEntryByNumber(0x0028,0x0100) == "24")
 250       return 3;
 251
 252    std::string PhotometricInterpretation =
 253                   gdcmHeader::GetEntryByNumber(0x0028,0x0004);
 254
 255    if ( ( PhotometricInterpretation == "PALETTE COLOR ") ) {
 256       if (HasLUT())   // PALETTE COLOR is NOT enough
 257          return 3;
 258       else
 259          return 1;
 260    }
 261
 262       //beware of trailing space at end of string
 263    if (PhotometricInterpretation.find(GDCM_UNFOUND) <
 264                            PhotometricInterpretation.length() ||
 265        PhotometricInterpretation.find("MONOCHROME1") <
 266                            PhotometricInterpretation.length() ||
 267        PhotometricInterpretation.find("MONOCHROME2") <
 268                            PhotometricInterpretation.length() )
 269        return 1;
 270     else
 271     // we assume that *all* kinds of YBR are dealt with
 272       return 3;
 273 }
 274
 275 /**
 276   * \ingroup gdcmHeaderHelper
 277   * \brief This function is intended to user that DOESN'T want
 278   *  to get RGB pixels image when it's stored as a PALETTE COLOR image
 279   *   - the (vtk) user is supposed to know how deal with LUTs -
 280   * \warning to be used with GetImagePixelsRaw()
 281   * @return 1 if Gray level, 3 if Color (RGB or YBR - NOT 'PALETTE COLOR' -)
 282   */
 283 int gdcmHeaderHelper::GetNumberOfScalarComponentsRaw() {
 284
 285      // 0028 0100 US IMG Bits Allocated
 286      // (in order no to be messed up by old RGB images)
 287    if (gdcmHeader::GetEntryByNumber(0x0028,0x0100) == "24")
 288       return 3;
 289
 290     // we assume that *all* kinds of YBR are dealt with
 291       return GetSamplesPerPixel();
 292 }
 293
 294 /**
 295   *\ingroup gdcmHeaderHelper
 296   *\brief gets the info from 0020,000d : Study Instance UID
 297   *\todo ? : return the ACR-NEMA element value if DICOM one is not found
 298   * @return Study Instance UID
 299  */
 300  std::string gdcmHeaderHelper::GetStudyUID(){
 301   return GetEntryByNumber(0x0020,0x000d); //0020 000d UI REL Study Instance UID
 302 }
 303
 304 /**
 305   *\ingroup gdcmHeaderHelper
 306   *\brief gets the info from 0020,000e : Series Instance UID
 307   *\todo ? : return the ACR-NEMA element value if DICOM one is not found
 308   * @return Series Instance UID
 309  */
 310  std::string gdcmHeaderHelper::GetSeriesUID(){
 311   return GetEntryByNumber(0x0020,0x000e); //0020 000e UI REL Series Instance UID
 312 }
 313
 314 /**
 315   *\ingroup gdcmHeaderHelper
 316   *\brief gets the info from 0008,0016 : SOP Class UID
 317   *\todo ? : return the ACR-NEMA element value if DICOM one is not found
 318   * @return SOP Class UID
 319  */
 320  std::string gdcmHeaderHelper::GetClassUID(){
 321   return GetEntryByNumber(0x0008,0x0016); //0008 0016 UI ID SOP Class UID
 322 }
 323
 324 /**
 325   *\ingroup gdcmHeaderHelper
 326   *\brief gets the info from 0008,0018 : SOP Instance UID
 327   *\todo ? : return the ACR-NEMA element value if DICOM one is not found
 328   * @return SOP Instance UID
 329  */
 330  std::string gdcmHeaderHelper::GetInstanceUID(){
 331   return GetEntryByNumber(0x0008,0x0018); //0008 0018 UI ID SOP Instance UID
 332 }
 333 //
 334 // --------------  Remember ! ----------------------------------
 335 //
 336 // Image Position Patient                              (0020,0032):
 337 // If not found (ACR_NEMA) we try Image Position       (0020,0030)
 338 // If not found (ACR-NEMA), we consider Slice Location (0020,1041)
 339 //                                   or Location       (0020,0050)
 340 // as the Z coordinate,
 341 // 0. for all the coordinates if nothing is found
 342
 343 // TODO : find a way to inform the caller nothing was found
 344 // TODO : How to tell the caller a wrong number of values was found?
 345 //
 346 // ---------------------------------------------------------------
 347 //
 348
 349 /**
 350   * \ingroup gdcmHeaderHelper
 351   * \brief gets the info from 0020,0032 : Image Position Patient
 352   *                 else from 0020,0030 : Image Position (RET)
 353   *                 else 0.
 354   * @return up-left image corner X position
 355   */
 356
 357 float gdcmHeaderHelper::GetXOrigin() {
 358     float xImPos, yImPos, zImPos;
 359     std::string StrImPos = GetEntryByNumber(0x0020,0x0032);
 360
 361     if (StrImPos == GDCM_UNFOUND) {
 362        dbg.Verbose(0, "gdcmHeader::GetXImagePosition: unfound Image Position Patient (0020,0032)");
 363        StrImPos = GetEntryByNumber(0x0020,0x0030); // For ACR-NEMA images
 364        if (StrImPos == GDCM_UNFOUND) {
 365           dbg.Verbose(0, "gdcmHeader::GetXImagePosition: unfound Image Position (RET) (0020,0030)");
 366           // How to tell the caller nothing was found ?
 367          return 0.;
 368        }
 369      }
 370    if( sscanf( StrImPos.c_str(), "%f\\%f\\%f", &xImPos, &yImPos, &zImPos) != 3)
 371      return 0.;
 372    return xImPos;
 373 }
 374
 375 /**
 376   * \ingroup gdcmHeaderHelper
 377   * \brief gets the info from 0020,0032 : Image Position Patient
 378   *                 else from 0020,0030 : Image Position (RET)
 379   *                 else 0.
 380   * @return up-left image corner Y position
 381   */
 382 float gdcmHeaderHelper::GetYOrigin() {
 383     float xImPos, yImPos, zImPos;
 384     std::string StrImPos = GetEntryByNumber(0x0020,0x0032);
 385
 386     if (StrImPos == GDCM_UNFOUND) {
 387        dbg.Verbose(0, "gdcmHeader::GetYImagePosition: unfound Image Position Patient (0020,0032)");
 388        StrImPos = GetEntryByNumber(0x0020,0x0030); // For ACR-NEMA images
 389        if (StrImPos == GDCM_UNFOUND) {
 390           dbg.Verbose(0, "gdcmHeader::GetYImagePosition: unfound Image Position (RET) (0020,0030)");
 391           // How to tell the caller nothing was found ?
 392            return 0.;
 393        }
 394      }
 395    if( sscanf( StrImPos.c_str(), "%f\\%f\\%f", &xImPos, &yImPos, &zImPos) != 3)
 396      return 0.;
 397    return yImPos;
 398 }
 399
 400 /**
 401   * \ingroup gdcmHeaderHelper
 402   * \brief gets the info from 0020,0032 : Image Position Patient
 403   * \               else from 0020,0030 : Image Position (RET)
 404   * \               else from 0020,1041 : Slice Location
 405   * \               else from 0020,0050 : Location
 406   * \               else 0.
 407   * @return up-left image corner Z position
 408   */
 409 float gdcmHeaderHelper::GetZOrigin() {
 410    float xImPos, yImPos, zImPos;
 411    std::string StrImPos = GetEntryByNumber(0x0020,0x0032);
 412    if (StrImPos != GDCM_UNFOUND) {
 413       if( sscanf( StrImPos.c_str(), "%f\\%f\\%f", &xImPos, &yImPos, &zImPos) != 3) {
 414          dbg.Verbose(0, "gdcmHeader::GetZImagePosition: wrong Image Position Patient (0020,0032)");
 415          return 0.;  // bug in the element 0x0020,0x0032
 416       } else {
 417          return zImPos;
 418       }
 419    }
 420    StrImPos = GetEntryByNumber(0x0020,0x0030); // For ACR-NEMA images
 421    if (StrImPos != GDCM_UNFOUND) {
 422       if( sscanf( StrImPos.c_str(), "%f\\%f\\%f", &xImPos, &yImPos, &zImPos) != 3) {
 423          dbg.Verbose(0, "gdcmHeader::GetZImagePosition: wrong Image Position (RET) (0020,0030)");
 424          return 0.;  // bug in the element 0x0020,0x0032
 425       } else {
 426          return zImPos;
 427       }
 428    }
 429    std::string StrSliceLocation = GetEntryByNumber(0x0020,0x1041);// for *very* old ACR-NEMA images
 430    if (StrSliceLocation != GDCM_UNFOUND) {
 431       if( sscanf( StrSliceLocation.c_str(), "%f", &zImPos) !=1) {
 432          dbg.Verbose(0, "gdcmHeader::GetZImagePosition: wrong Slice Location (0020,1041)");
 433          return 0.;  // bug in the element 0x0020,0x1041
 434       } else {
 435          return zImPos;
 436       }
 437    }
 438    dbg.Verbose(0, "gdcmHeader::GetZImagePosition: unfound Slice Location (0020,1041)");
 439    std::string StrLocation = GetEntryByNumber(0x0020,0x0050);
 440    if (StrLocation != GDCM_UNFOUND) {
 441       if( sscanf( StrLocation.c_str(), "%f", &zImPos) !=1) {
 442          dbg.Verbose(0, "gdcmHeader::GetZImagePosition: wrong Location (0020,0050)");
 443          return 0.;  // bug in the element 0x0020,0x0050
 444       } else {
 445          return zImPos;
 446       }
 447    }
 448    dbg.Verbose(0, "gdcmHeader::GetYImagePosition: unfound Location (0020,0050)");
 449    return 0.; // Hopeless
 450 }
 451
 452 /**
 453   * \ingroup gdcmHeaderHelper
 454   * \brief gets the info from 0020,0013 : Image Number
 455   * \               else 0.
 456   * @return image number
 457   */
 458 int gdcmHeaderHelper::GetImageNumber() {
 459   //The function i atoi() takes the address of an area of memory as parameter and converts
 460   //the string stored at that location to an integer using the external decimal to internal
 461   //binary conversion rules. This may be preferable to sscanf() since atoi() is a much smaller,
 462   // simpler and faster function. sscanf() can do all possible conversions whereas atoi() can
 463   //only do single decimal integer conversions.
 464   std::string StrImNumber = GetEntryByNumber(0x0020,0x0013); //0020 0013 IS REL Image Number
 465   if (StrImNumber != GDCM_UNFOUND) {
 466     return atoi( StrImNumber.c_str() );
 467   }
 468   return 0;   //Hopeless
 469 }
 470
 471 /**
 472   * \ingroup gdcmHeaderHelper
 473   * \brief gets the info from 0008,0060 : Modality
 474   * @return Modality Type
 475   */
 476 ModalityType gdcmHeaderHelper::GetModality(void) {
 477   std::string StrModality = GetEntryByNumber(0x0008,0x0060); //0008 0060 CS ID Modality
 478   if (StrModality != GDCM_UNFOUND) {
 479          if ( StrModality.find("AU") < StrModality.length()) return AU;
 480     else if ( StrModality.find("AS") < StrModality.length()) return AS;
 481     else if ( StrModality.find("BI") < StrModality.length()) return BI;
 482     else if ( StrModality.find("CF") < StrModality.length()) return CF;
 483     else if ( StrModality.find("CP") < StrModality.length()) return CP;
 484     else if ( StrModality.find("CR") < StrModality.length()) return CR;
 485     else if ( StrModality.find("CT") < StrModality.length()) return CT;
 486     else if ( StrModality.find("CS") < StrModality.length()) return CS;
 487     else if ( StrModality.find("DD") < StrModality.length()) return DD;
 488     else if ( StrModality.find("DF") < StrModality.length()) return DF;
 489     else if ( StrModality.find("DG") < StrModality.length()) return DG;
 490     else if ( StrModality.find("DM") < StrModality.length()) return DM;
 491     else if ( StrModality.find("DS") < StrModality.length()) return DS;
 492     else if ( StrModality.find("DX") < StrModality.length()) return DX;
 493     else if ( StrModality.find("ECG") < StrModality.length()) return ECG;
 494     else if ( StrModality.find("EPS") < StrModality.length()) return EPS;
 495     else if ( StrModality.find("FA") < StrModality.length()) return FA;
 496     else if ( StrModality.find("FS") < StrModality.length()) return FS;
 497     else if ( StrModality.find("HC") < StrModality.length()) return HC;
 498     else if ( StrModality.find("HD") < StrModality.length()) return HD;
 499     else if ( StrModality.find("LP") < StrModality.length()) return LP;
 500     else if ( StrModality.find("LS") < StrModality.length()) return LS;
 501     else if ( StrModality.find("MA") < StrModality.length()) return MA;
 502     else if ( StrModality.find("MR") < StrModality.length()) return MR;
 503     else if ( StrModality.find("NM") < StrModality.length()) return NM;
 504     else if ( StrModality.find("OT") < StrModality.length()) return OT;
 505     else if ( StrModality.find("PT") < StrModality.length()) return PT;
 506     else if ( StrModality.find("RF") < StrModality.length()) return RF;
 507     else if ( StrModality.find("RG") < StrModality.length()) return RG;
 508     else if ( StrModality.find("RTDOSE")  < StrModality.length()) return RTDOSE;
 509     else if ( StrModality.find("RTIMAGE") < StrModality.length()) return RTIMAGE;
 510     else if ( StrModality.find("RTPLAN")  < StrModality.length()) return RTPLAN;
 511     else if ( StrModality.find("RTSTRUCT")< StrModality.length()) return RTSTRUCT;
 512     else if ( StrModality.find("SM") < StrModality.length()) return SM;
 513     else if ( StrModality.find("ST") < StrModality.length()) return ST;
 514     else if ( StrModality.find("TG") < StrModality.length()) return TG;
 515     else if ( StrModality.find("US") < StrModality.length()) return US;
 516     else if ( StrModality.find("VF") < StrModality.length()) return VF;
 517     else if ( StrModality.find("XA") < StrModality.length()) return XA;
 518     else if ( StrModality.find("XC") < StrModality.length()) return XC;
 519
 520     else
 521     {
 522       //throw error return value ???
 523       // specified <> unknow in our database
 524       return Unknow;
 525     }
 526   }
 527   return Unknow;
 528 }
 529
 530 /**
 531   * \ingroup gdcmHeaderHelper
 532   * \brief gets the info from 0020,0037 : Image Orientation Patient
 533   * @param iop adress of the (6)float aray to receive values
 534   * @return cosines of image orientation patient
 535   */
 536 void gdcmHeaderHelper::GetImageOrientationPatient( float* iop ) {
 537
 538   //iop is supposed to be float[6]
 539   iop[0] = iop[1] = iop[2] = iop[3] = iop[4] = iop[5] = 0;
 540
 541   std::string StrImOriPat = GetEntryByNumber(0x0020,0x0037); // 0020 0037 DS REL Image Orientation (Patient)
 542   if (StrImOriPat != GDCM_UNFOUND) {
 543     if( sscanf( StrImOriPat.c_str(), "%f\\%f\\%f\\%f\\%f\\%f",
 544             &iop[0], &iop[1], &iop[2], &iop[3], &iop[4], &iop[5]) != 6) {
 545          dbg.Verbose(0, "gdcmHeader::GetImageOrientationPatient: wrong Image Orientation Patient (0020,0037)");
 546          return ;  // bug in the element 0x0020,0x0037
 547     }
 548     else
 549       return ;
 550   }
 551
 552   //For ACR-NEMA
 553   StrImOriPat = GetEntryByNumber(0x0020,0x0035); //0020 0035 DS REL Image Orientation (RET)
 554   if (StrImOriPat != GDCM_UNFOUND) {
 555     if( sscanf( StrImOriPat.c_str(), "%f\\%f\\%f\\%f\\%f\\%f",
 556             &iop[0], &iop[1], &iop[2], &iop[3], &iop[4], &iop[5]) != 6) {
 557          dbg.Verbose(0, "gdcmHeader::GetImageOrientationPatient: wrong Image Orientation Patient (0020,0035)");
 558          return ;  // bug in the element 0x0020,0x0035
 559     }
 560     else
 561       return ;
 562   }
 563 }
 564
 565 //-----------------------------------------------------------------------------
 566 // Protected
 567
 568 //-----------------------------------------------------------------------------
 569 // Private
 570
 571 //-----------------------------------------------------------------------------
 572
 573
 574
 575 //-----------------------------------------------------------------------------
 576 // gdcmSerieHeaderHelper
 577 //-----------------------------------------------------------------------------
 578 // Constructor / Destructor
 579 gdcmSerieHeaderHelper::~gdcmSerieHeaderHelper(){
 580   //! \todo
 581   for (std::list<gdcmHeaderHelper*>::iterator it  = CoherentGdcmFileList.begin();
 582         it != CoherentGdcmFileList.end(); it++)
 583   {
 584     delete *it;
 585   }
 586   CoherentGdcmFileList.clear();
 587 }
 588
 589 //-----------------------------------------------------------------------------
 590 // Print
 591
 592 //-----------------------------------------------------------------------------
 593 // Public
 594 /**
 595  * \ingroup gdcmHeaderHelper
 596  * \brief add a gdcmFile to the list based on file name
 597  * @param   filename Name of the file to deal with
 598  */
 599 void gdcmSerieHeaderHelper::AddFileName(std::string filename) {
 600   gdcmHeaderHelper *GdcmFile = new gdcmHeaderHelper( filename.c_str() );
 601   this->CoherentGdcmFileList.push_back( GdcmFile );
 602 }
 603
 604 /**
 605  * \ingroup gdcmHeaderHelper
 606  * \brief add a gdcmFile to the list
 607  * @param   file gdcmHeaderHelper to add
 608  */
 609 void gdcmSerieHeaderHelper::AddGdcmFile(gdcmHeaderHelper *file){
 610   this->CoherentGdcmFileList.push_back( file );
 611 }
 612
 613 /**
 614  * \ingroup gdcmHeaderHelper
 615  * \brief Sets the Directory
 616  * @param   dir Name of the directory to deal with
 617  */
 618 void gdcmSerieHeaderHelper::SetDirectory(std::string dir){
 619   gdcmDirList filenames_list(dir);  //OS specific
 620
 621   for(gdcmDirList::iterator it = filenames_list.begin();
 622       it !=filenames_list.end(); it++)
 623   {
 624     gdcmHeaderHelper *file = new gdcmHeaderHelper( it->c_str() );
 625     this->CoherentGdcmFileList.push_back( file );
 626   }
 627 }
 628
 629 /**
 630  * \ingroup gdcmHeaderHelper
 631  * \brief Sorts the File List
 632  * \warning This could be implemented in a 'Strategy Pattern' approach
 633  *          But as I don't know how to do it, I leave it this way
 634  *          BTW, this is also a Strategy, I don't know this is the best approach :)
 635 */
 636 void gdcmSerieHeaderHelper::OrderGdcmFileList(){
 637   if( ImagePositionPatientOrdering() ) {
 638     return ;
 639   }
 640   else if( ImageNumberOrdering() ) {
 641     return ;
 642   } else  {
 643     FileNameOrdering();
 644   }
 645 }
 646
 647 /**
 648  * \ingroup gdcmHeaderHelper
 649  * \brief Gets the *coherent* File List
 650  * @return the *coherent* File List
 651 */
 652 std::list<gdcmHeaderHelper*> &gdcmSerieHeaderHelper::GetGdcmFileList() {
 653   return CoherentGdcmFileList;
 654 }
 655
 656 //-----------------------------------------------------------------------------
 657 // Protected
 658
 659 //-----------------------------------------------------------------------------
 660 // Private
 661 /**
 662  * \ingroup gdcmHeaderHelper
 663  * \brief sorts the images, according to their Patient Position
 664  *  We may order, considering :
 665  *   -# Image Number
 666  *   -# Image Position Patient
 667  *   -# More to come :)
 668  * @return false only if the header is bugged !
 669  */
 670 bool gdcmSerieHeaderHelper::ImagePositionPatientOrdering()
 671 //based on Jolinda's algorithm
 672 {
 673   //iop is calculated based on the file file
 674   float *cosines = new float[6];
 675   float normal[3];
 676   float ipp[3];
 677   float dist;
 678   float min, max;
 679   bool first = true;
 680   int n=0;
 681   std::vector<float> distlist;
 682
 683   //!\todo rewrite this for loop.
 684   for (std::list<gdcmHeaderHelper*>::iterator it  = CoherentGdcmFileList.begin();
 685         it != CoherentGdcmFileList.end(); it++)
 686   {
 687     if(first) {
 688       (*it)->GetImageOrientationPatient(cosines);
 689
 690       //You only have to do this once for all slices in the volume. Next, for
 691       //each slice, calculate the distance along the slice normal using the IPP
 692       //tag ("dist" is initialized to zero before reading the first slice) :
 693       normal[0] = cosines[1]*cosines[5] - cosines[2]*cosines[4];
 694       normal[1] = cosines[2]*cosines[3] - cosines[0]*cosines[5];
 695       normal[2] = cosines[0]*cosines[4] - cosines[1]*cosines[3];
 696
 697       ipp[0] = (*it)->GetXOrigin();
 698       ipp[1] = (*it)->GetYOrigin();
 699       ipp[2] = (*it)->GetZOrigin();
 700
 701       dist = 0;
 702       for (int i = 0; i < 3; ++i)
 703           dist += normal[i]*ipp[i];
 704
 705       if( dist == 0 )
 706       {
 707         delete[] cosines;
 708         return false;
 709       }
 710
 711       distlist.push_back( dist );
 712
 713       max = min = dist;
 714       first = false;
 715     }
 716     else {
 717       ipp[0] = (*it)->GetXOrigin();
 718       ipp[1] = (*it)->GetYOrigin();
 719       ipp[2] = (*it)->GetZOrigin();
 720
 721       dist = 0;
 722       for (int i = 0; i < 3; ++i)
 723           dist += normal[i]*ipp[i];
 724
 725       if( dist == 0 )
 726       {
 727         delete[] cosines;
 728         return false;
 729       }
 730
 731       distlist.push_back( dist );
 732
 733       min = (min < dist) ? min : dist;
 734       max = (max > dist) ? max : dist;
 735     }
 736     n++;
 737   }
 738
 739     //Then I order the slices according to the value "dist". Finally, once
 740     //I've read in all the slices, I calculate the z-spacing as the difference
 741     //between the "dist" values for the first two slices.
 742     std::vector<gdcmHeaderHelper*> CoherentGdcmFileVector(n);
 743     //CoherentGdcmFileVector.reserve( n );
 744     CoherentGdcmFileVector.resize( n );
 745     //assert( CoherentGdcmFileVector.capacity() >= n );
 746
 747     float step = (max - min)/(n - 1);
 748     int pos;
 749     n = 0;
 750
 751     //VC++ don't understand what scope is !! it -> it2
 752     for (std::list<gdcmHeaderHelper*>::iterator it2  = CoherentGdcmFileList.begin();
 753         it2 != CoherentGdcmFileList.end(); it2++, n++)
 754     {
 755       //2*n sort algo !!
 756       //Assumption: all files are present (no one missing)
 757       pos = (int)( fabs( (distlist[n]-min)/step) + .5 );
 758
 759       CoherentGdcmFileVector[pos] = *it2;
 760     }
 761
 762   CoherentGdcmFileList.clear();  //this doesn't delete list's element, node only
 763
 764   //VC++ don't understand what scope is !! it -> it3
 765   for (std::vector<gdcmHeaderHelper*>::iterator it3  = CoherentGdcmFileVector.begin();
 766         it3 != CoherentGdcmFileVector.end(); it3++)
 767   {
 768     CoherentGdcmFileList.push_back( *it3 );
 769   }
 770
 771   distlist.clear();
 772   CoherentGdcmFileVector.clear();
 773   delete[] cosines;
 774
 775   return true;
 776 }
 777
 778 /**
 779  * \ingroup gdcmHeaderHelper
 780  * \brief sorts the images, according to their Image Number
 781  * @return false only if the header is bugged !
 782  */
 783
 784 bool gdcmSerieHeaderHelper::ImageNumberOrdering() {
 785   int min, max, pos;
 786   int n = 0;//CoherentGdcmFileList.size() is a O(N) operation !!
 787   unsigned char *partition;
 788
 789   std::list<gdcmHeaderHelper*>::iterator it  = CoherentGdcmFileList.begin();
 790   min = max = (*it)->GetImageNumber();
 791
 792   for (; it != CoherentGdcmFileList.end(); it++, n++)
 793   {
 794     pos = (*it)->GetImageNumber();
 795
 796     //else
 797     min = (min < pos) ? min : pos;
 798   }
 799
 800   //bzeros(partition, n); //Cette fonction est déconseillée, utilisez plutôt memset.
 801   partition = new unsigned char[n];
 802   memset(partition, 0, n);
 803
 804   std::vector<gdcmHeaderHelper*> CoherentGdcmFileVector(n);
 805
 806   //VC++ don't understand what scope is !! it -> it2
 807   for (std::list<gdcmHeaderHelper*>::iterator it2  = CoherentGdcmFileList.begin();
 808         it2 != CoherentGdcmFileList.end(); it2++)
 809   {
 810     pos = (*it2)->GetImageNumber();
 811     CoherentGdcmFileVector[pos - min] = *it2;
 812     partition[pos - min]++;
 813   }
 814
 815   unsigned char mult = 1;
 816   for(int i=0; i<n ; i++)
 817   {
 818     mult *= partition[i];
 819   }
 820
 821   //VC++ don't understand what scope is !! it -> it3
 822   CoherentGdcmFileList.clear();  //this doesn't delete list's element, node only
 823   for (std::vector<gdcmHeaderHelper*>::iterator it3  = CoherentGdcmFileVector.begin();
 824         it3 != CoherentGdcmFileVector.end(); it3++)
 825   {
 826     CoherentGdcmFileList.push_back( *it3 );
 827   }
 828   CoherentGdcmFileVector.clear();
 829
 830   delete[] partition;
 831   return (mult!=0);
 832 }
 833
 834
 835 /**
 836  * \ingroup gdcmHeaderHelper
 837  * \brief sorts the images, according to their File Name
 838  * @return false only if the header is bugged !
 839  */
 840  bool gdcmSerieHeaderHelper::FileNameOrdering() {
 841   //using the sort
 842   //sort(CoherentGdcmFileList.begin(), CoherentGdcmFileList.end());
 843   return true;
 844 }
 845
 846 //-----------------------------------------------------------------------------