src/gdcmFile.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: gdcmFile.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2005/07/21 14:01:49 $
   7   Version:   $Revision: 1.255 $
   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 //
  20 // --------------  Remember ! ----------------------------------
  21 //
  22 // Image Position Patient                              (0020,0032):
  23 // If not found (ACR_NEMA) we try Image Position       (0020,0030)
  24 // If not found (ACR-NEMA), we consider Slice Location (0020,1041)
  25 //                                   or Location       (0020,0050)
  26 //                                   as the Z coordinate,
  27 // 0. for all the coordinates if nothing is found
  28 //
  29 // Image Position (Patient) (0020,0032) VM=3 What is it used for?
  30 // -->
  31 //  The attribute Patient Orientation (0020,0020) from the General Image Module
  32 // is of type 2C and has the condition Required if image does not require
  33 // Image Orientation (0020,0037) and Image Position (0020,0032).
  34 // However, if the image does require the attributes
  35 // - Image Orientation (Patient) (0020,0037), VM=6
  36 // - Image Position Patient (0020,0032), VM=3
  37 // then attribute Patient Orientation (0020,0020) should not be present
  38 //  in the images.
  39 //
  40 // Remember also :
  41 // Patient Position (0018,5100) values : HFP   = Head First-Prone
  42 //                                       HFS   = Head First-Supine
  43 //                                       HFDR  = Head First-Decubitus Right
  44 //                                       HFDL = Head First-Decubitus Left
  45 //                                       FFDR = Feet First-Decubitus Right
  46 //                                       FFDL = Feet First-Decubitus Left
  47 //                                       FFP  = Feet First-Prone
  48 //                                       FFS  = Feet First-Supine
  49 //                    can also find      SEMIERECT
  50 //                                       SUPINE
  51 // CS 2 Patient Orientation (0020 0020)
  52 //               When the coordinates of the image
  53 //               are always present, this field is almost never used.
  54 //               Better we don't tust it too much ...
  55 //               Found Values are :      L\P
  56 //                                       L\FP
  57 //                                       P\F
  58 //                                       L\F
  59 //                                       P\FR
  60 //                                       R\F
  61 //
  62 // (0020|0037) [Image Orientation (Patient)] [1\0\0\0\1\0 ]
  63
  64
  65 // ---------------------------------------------------------------
  66 //
  67 #include "gdcmFile.h"
  68 #include "gdcmGlobal.h"
  69 #include "gdcmUtil.h"
  70 #include "gdcmDebug.h"
  71 #include "gdcmTS.h"
  72 #include "gdcmValEntry.h"
  73 #include "gdcmBinEntry.h"
  74 #include "gdcmSeqEntry.h"
  75 #include "gdcmRLEFramesInfo.h"
  76 #include "gdcmJPEGFragmentsInfo.h"
  77
  78 #include <vector>
  79 #include <stdio.h> //sscanf
  80 #include <stdlib.h> // for atoi
  81
  82 namespace gdcm
  83 {
  84 //-----------------------------------------------------------------------------
  85 // Constructor / Destructor
  86
  87 /**
  88  * \brief Constructor used when we want to generate dicom files from scratch
  89  */
  90 File::File():
  91    Document()
  92 {
  93    RLEInfo  = new RLEFramesInfo;
  94    JPEGInfo = new JPEGFragmentsInfo;
  95    GrPixel  = 0x7fe0;  // to avoid further troubles
  96    NumPixel = 0x0010;
  97 }
  98
  99
 100 /**
 101  * \brief   Canonical destructor.
 102  */
 103 File::~File ()
 104 {
 105    if ( RLEInfo )
 106       delete RLEInfo;
 107    if ( JPEGInfo )
 108       delete JPEGInfo;
 109 }
 110
 111 //-----------------------------------------------------------------------------
 112 // Public
 113 /**
 114  * \brief   Loader
 115  * @return false if file cannot be open or no swap info was found,
 116  *         or no tag was found.
 117  */
 118 bool File::Load( )
 119 {
 120    if ( ! this->Document::Load( ) )
 121       return false;
 122
 123     return DoTheLoadingJob( );
 124 }
 125
 126 /**
 127  * \brief   Does the Loading Job (internal use only)
 128  * @return false if file cannot be open or no swap info was found,
 129  *         or no tag was found.
 130  */
 131 bool File::DoTheLoadingJob( )
 132 {
 133
 134    // for some ACR-NEMA images GrPixel, NumPixel is *not* 7fe0,0010
 135    // We may encounter the 'RETired' (0x0028, 0x0200) tag
 136    // (Image Location") . This entry contains the number of
 137    // the group that contains the pixel data (hence the "Pixel Data"
 138    // is found by indirection through the "Image Location").
 139    // Inside the group pointed by "Image Location" the searched element
 140    // is conventionally the element 0x0010 (when the norm is respected).
 141    // When the "Image Location" is missing we default to group 0x7fe0.
 142    // Note: this IS the right place for the code
 143
 144    // Image Location
 145    const std::string &imgLocation = GetEntryValue(0x0028, 0x0200);
 146    if ( imgLocation == GDCM_UNFOUND )
 147    {
 148       // default value
 149       GrPixel = 0x7fe0;
 150    }
 151    else
 152    {
 153       GrPixel = (uint16_t) atoi( imgLocation.c_str() );
 154    }
 155
 156    // sometimes Image Location value doesn't follow
 157    // the supposed processor endianness.
 158    // see gdcmData/cr172241.dcm
 159    if ( GrPixel == 0xe07f )
 160    {
 161       GrPixel = 0x7fe0;
 162    }
 163
 164    if ( GrPixel != 0x7fe0 )
 165    {
 166       // This is a kludge for old dirty Philips imager.
 167       NumPixel = 0x1010;
 168    }
 169    else
 170    {
 171       NumPixel = 0x0010;
 172    }
 173
 174    // Now, we know GrPixel and NumPixel.
 175    // Let's create a VirtualDictEntry to allow a further VR modification
 176    // and force VR to match with BitsAllocated.
 177    DocEntry *entry = GetDocEntry(GrPixel, NumPixel);
 178    if ( entry != 0 )
 179    {
 180       // Compute the RLE or JPEG info
 181       OpenFile();
 182       const std::string &ts = GetTransferSyntax();
 183       Fp->seekg( entry->GetOffset(), std::ios::beg );
 184       if ( Global::GetTS()->IsRLELossless(ts) )
 185          ComputeRLEInfo();
 186       else if ( Global::GetTS()->IsJPEG(ts) )
 187          ComputeJPEGFragmentInfo();
 188       CloseFile();
 189
 190       // Create a new BinEntry to change the DictEntry
 191       // The changed DictEntry will have
 192       // - a correct PixelVR OB or OW)
 193       // - the name to "Pixel Data"
 194       BinEntry *oldEntry = dynamic_cast<BinEntry *>(entry);
 195       if (oldEntry)
 196       {
 197          std::string PixelVR;
 198          // 8 bits allocated is a 'O Bytes' , as well as 24 (old ACR-NEMA RGB)
 199          // more than 8 (i.e 12, 16) is a 'O Words'
 200          if ( GetBitsAllocated() == 8 || GetBitsAllocated() == 24 )
 201             PixelVR = "OB";
 202          else
 203             PixelVR = "OW";
 204
 205          // Change only made if usefull
 206          if ( PixelVR != oldEntry->GetVR() )
 207          {
 208             DictEntry* newDict = NewVirtualDictEntry(GrPixel,NumPixel,
 209                                                      PixelVR,"1","Pixel Data");
 210
 211             BinEntry *newEntry = new BinEntry(newDict);
 212             newEntry->Copy(entry);
 213             newEntry->SetBinArea(oldEntry->GetBinArea(),oldEntry->IsSelfArea());
 214             oldEntry->SetSelfArea(false);
 215
 216             RemoveEntry(oldEntry);
 217             AddEntry(newEntry);
 218          }
 219       }
 220    }
 221    return true;
 222 }
 223 /**
 224  * \brief  This predicate, based on hopefully reasonable heuristics,
 225  *         decides whether or not the current File was properly parsed
 226  *         and contains the mandatory information for being considered as
 227  *         a well formed and usable Dicom/Acr File.
 228  * @return true when File is the one of a reasonable Dicom/Acr file,
 229  *         false otherwise.
 230  */
 231 bool File::IsReadable()
 232 {
 233    if ( !Document::IsReadable() )
 234    {
 235       return false;
 236    }
 237
 238    const std::string &res = GetEntryValue(0x0028, 0x0005);
 239    if ( res != GDCM_UNFOUND && atoi(res.c_str()) > 4 )
 240    {
 241       gdcmWarningMacro("Wrong Image Dimensions" << res);
 242       return false; // Image Dimensions
 243    }
 244    bool b0028_0100 = true;
 245    if ( !GetDocEntry(0x0028, 0x0100) )
 246    {
 247       gdcmWarningMacro("Bits Allocated (0028|0100) not found");
 248       //return false; // "Bits Allocated"
 249       b0028_0100 = false;
 250    }
 251    bool b0028_0101 = true;
 252    if ( !GetDocEntry(0x0028, 0x0101) )
 253    {
 254       gdcmWarningMacro("Bits Stored (0028|0101) not found");
 255       //return false; // "Bits Stored"
 256       b0028_0101 = false;
 257    }
 258    bool b0028_0102 = true;
 259    if ( !GetDocEntry(0x0028, 0x0102) )
 260    {
 261       gdcmWarningMacro("Hight Bit (0028|0102) not found");
 262       //return false; // "High Bit"
 263       b0028_0102 = false;
 264    }
 265    bool b0028_0103 = true;
 266    if ( !GetDocEntry(0x0028, 0x0103) )
 267    {
 268       gdcmWarningMacro("Pixel Representation (0028|0103) not found");
 269       //return false; // "Pixel Representation" i.e. 'Sign' ( 0 : unsigned, 1 : signed)
 270       b0028_0103 = false;
 271    }
 272
 273    if ( !b0028_0100 && !b0028_0101 && !b0028_0102 && !b0028_0103)
 274    {
 275       gdcmWarningMacro("Too much mandatory Tags missing !");
 276       return false;
 277    }
 278
 279    if ( !GetDocEntry(GrPixel, NumPixel) )
 280    {
 281       gdcmWarningMacro("Pixel Dicom Element " << std::hex <<
 282                         GrPixel << "|" << NumPixel << "not found");
 283       return false; // Pixel Dicom Element not found :-(
 284    }
 285    return true;
 286 }
 287
 288 /**
 289  * \brief gets the info from 0020,0013 : Image Number else 0.
 290  * @return image number
 291  */
 292 int File::GetImageNumber()
 293 {
 294    //0020 0013 : Image Number
 295    std::string strImNumber = GetEntryValue(0x0020,0x0013);
 296    if ( strImNumber != GDCM_UNFOUND )
 297    {
 298       return atoi( strImNumber.c_str() );
 299    }
 300    return 0;   //Hopeless
 301 }
 302
 303 /**
 304  * \brief gets the info from 0008,0060 : Modality
 305  * @return Modality Type
 306  */
 307 ModalityType File::GetModality()
 308 {
 309    // 0008 0060 : Modality
 310    std::string strModality = GetEntryValue(0x0008,0x0060);
 311    if ( strModality != GDCM_UNFOUND )
 312    {
 313            if ( strModality.find("AU")  < strModality.length()) return AU;
 314       else if ( strModality.find("AS")  < strModality.length()) return AS;
 315       else if ( strModality.find("BI")  < strModality.length()) return BI;
 316       else if ( strModality.find("CF")  < strModality.length()) return CF;
 317       else if ( strModality.find("CP")  < strModality.length()) return CP;
 318       else if ( strModality.find("CR")  < strModality.length()) return CR;
 319       else if ( strModality.find("CT")  < strModality.length()) return CT;
 320       else if ( strModality.find("CS")  < strModality.length()) return CS;
 321       else if ( strModality.find("DD")  < strModality.length()) return DD;
 322       else if ( strModality.find("DF")  < strModality.length()) return DF;
 323       else if ( strModality.find("DG")  < strModality.length()) return DG;
 324       else if ( strModality.find("DM")  < strModality.length()) return DM;
 325       else if ( strModality.find("DS")  < strModality.length()) return DS;
 326       else if ( strModality.find("DX")  < strModality.length()) return DX;
 327       else if ( strModality.find("ECG") < strModality.length()) return ECG;
 328       else if ( strModality.find("EPS") < strModality.length()) return EPS;
 329       else if ( strModality.find("FA")  < strModality.length()) return FA;
 330       else if ( strModality.find("FS")  < strModality.length()) return FS;
 331       else if ( strModality.find("HC")  < strModality.length()) return HC;
 332       else if ( strModality.find("HD")  < strModality.length()) return HD;
 333       else if ( strModality.find("LP")  < strModality.length()) return LP;
 334       else if ( strModality.find("LS")  < strModality.length()) return LS;
 335       else if ( strModality.find("MA")  < strModality.length()) return MA;
 336       else if ( strModality.find("MR")  < strModality.length()) return MR;
 337       else if ( strModality.find("NM")  < strModality.length()) return NM;
 338       else if ( strModality.find("OT")  < strModality.length()) return OT;
 339       else if ( strModality.find("PT")  < strModality.length()) return PT;
 340       else if ( strModality.find("RF")  < strModality.length()) return RF;
 341       else if ( strModality.find("RG")  < strModality.length()) return RG;
 342       else if ( strModality.find("RTDOSE")
 343                                         < strModality.length()) return RTDOSE;
 344       else if ( strModality.find("RTIMAGE")
 345                                         < strModality.length()) return RTIMAGE;
 346       else if ( strModality.find("RTPLAN")
 347                                         < strModality.length()) return RTPLAN;
 348       else if ( strModality.find("RTSTRUCT")
 349                                         < strModality.length()) return RTSTRUCT;
 350       else if ( strModality.find("SM")  < strModality.length()) return SM;
 351       else if ( strModality.find("ST")  < strModality.length()) return ST;
 352       else if ( strModality.find("TG")  < strModality.length()) return TG;
 353       else if ( strModality.find("US")  < strModality.length()) return US;
 354       else if ( strModality.find("VF")  < strModality.length()) return VF;
 355       else if ( strModality.find("XA")  < strModality.length()) return XA;
 356       else if ( strModality.find("XC")  < strModality.length()) return XC;
 357
 358       else
 359       {
 360          /// \todo throw error return value ???
 361          /// specified <> unknown in our database
 362          return Unknow;
 363       }
 364    }
 365    return Unknow;
 366 }
 367
 368 /**
 369  * \brief   Retrieve the number of columns of image.
 370  * @return  The encountered size when found, 0 by default.
 371  *          0 means the file is NOT USABLE. The caller will have to check
 372  */
 373 int File::GetXSize()
 374 {
 375    const std::string &strSize = GetEntryValue(0x0028,0x0011);
 376    if ( strSize == GDCM_UNFOUND )
 377    {
 378       return 0;
 379    }
 380    return atoi( strSize.c_str() );
 381 }
 382
 383 /**
 384  * \brief   Retrieve the number of lines of image.
 385  * \warning The defaulted value is 1 as opposed to File::GetXSize()
 386  * @return  The encountered size when found, 1 by default
 387  *          (The ACR-NEMA file contains a Signal, not an Image).
 388  */
 389 int File::GetYSize()
 390 {
 391    const std::string &strSize = GetEntryValue(0x0028,0x0010);
 392    if ( strSize != GDCM_UNFOUND )
 393    {
 394       return atoi( strSize.c_str() );
 395    }
 396    if ( IsDicomV3() )
 397    {
 398       return 0;
 399    }
 400
 401    // The Rows (0028,0010) entry was optional for ACR/NEMA.
 402    // (at least some images didn't have it.)
 403    // It might hence be a signal (1D image). So we default to 1:
 404    return 1;
 405 }
 406
 407 /**
 408  * \brief   Retrieve the number of planes of volume or the number
 409  *          of frames of a multiframe.
 410  * \warning When present we consider the "Number of Frames" as the third
 411  *          dimension. When missing we consider the third dimension as
 412  *          being the ACR-NEMA "Planes" tag content.
 413  * @return  The encountered size when found, 1 by default (single image).
 414  */
 415 int File::GetZSize()
 416 {
 417    // Both  DicomV3 and ACR/Nema consider the "Number of Frames"
 418    // as the third dimension.
 419    const std::string &strSize = GetEntryValue(0x0028,0x0008);
 420    if ( strSize != GDCM_UNFOUND )
 421    {
 422       return atoi( strSize.c_str() );
 423    }
 424
 425    // We then consider the "Planes" entry as the third dimension
 426    const std::string &strSize2 = GetEntryValue(0x0028,0x0012);
 427    if ( strSize2 != GDCM_UNFOUND )
 428    {
 429       return atoi( strSize2.c_str() );
 430    }
 431    return 1;
 432 }
 433
 434 /**
 435   * \brief gets the info from 0018,1164 : ImagerPixelSpacing
 436   *                      then 0028,0030 : Pixel Spacing
 437   *             else 1.0
 438   * @return X dimension of a pixel
 439   */
 440 float File::GetXSpacing()
 441 {
 442    float xspacing = 1.0;
 443    float yspacing = 1.0;
 444    int nbValues;
 445
 446    // To follow David Clunie's advice, we first check ImagerPixelSpacing
 447
 448    const std::string &strImagerPixelSpacing = GetEntryValue(0x0018,0x1164);
 449    if ( strImagerPixelSpacing != GDCM_UNFOUND )
 450    {
 451       if ( ( nbValues = sscanf( strImagerPixelSpacing.c_str(),
 452             "%f\\%f", &yspacing, &xspacing)) != 2 )
 453       {
 454          // if no values, xspacing is set to 1.0
 455          if ( nbValues == 0 )
 456             xspacing = 1.0;
 457          // if single value is found, xspacing is defaulted to yspacing
 458          if ( nbValues == 1 )
 459             xspacing = yspacing;
 460
 461          if ( xspacing == 0.0 )
 462             xspacing = 1.0;
 463
 464          return xspacing;
 465       }
 466    }
 467
 468    const std::string &strSpacing = GetEntryValue(0x0028,0x0030);
 469
 470    if ( strSpacing == GDCM_UNFOUND )
 471    {
 472       gdcmWarningMacro( "Unfound Pixel Spacing (0028,0030)" );
 473       return 1.;
 474    }
 475
 476    if ( ( nbValues = sscanf( strSpacing.c_str(),
 477          "%f \\%f ", &yspacing, &xspacing)) != 2 )
 478    {
 479       // if no values, xspacing is set to 1.0
 480       if ( nbValues == 0 )
 481          xspacing = 1.0;
 482       // if single value is found, xspacing is defaulted to yspacing
 483       if ( nbValues == 1 )
 484          xspacing = yspacing;
 485
 486       if ( xspacing == 0.0 )
 487          xspacing = 1.0;
 488
 489       return xspacing;
 490    }
 491
 492    // to avoid troubles with David Clunie's-like images (at least one)
 493    if ( xspacing == 0. && yspacing == 0.)
 494       return 1.;
 495
 496    if ( xspacing == 0.)
 497    {
 498       gdcmWarningMacro("gdcmData/CT-MONO2-8-abdo.dcm-like problem");
 499       // seems to be a bug in the header ...
 500       nbValues = sscanf( strSpacing.c_str(), "%f \\0\\%f ", &yspacing, &xspacing);
 501       gdcmAssertMacro( nbValues == 2 );
 502    }
 503
 504    return xspacing;
 505 }
 506
 507 /**
 508   * \brief gets the info from 0018,1164 : ImagerPixelSpacing
 509   *               then from   0028,0030 : Pixel Spacing
 510   *             else 1.0
 511   * @return Y dimension of a pixel
 512   */
 513 float File::GetYSpacing()
 514 {
 515    float yspacing = 1.;
 516    int nbValues;
 517    // To follow David Clunie's advice, we first check ImagerPixelSpacing
 518
 519    const std::string &strImagerPixelSpacing = GetEntryValue(0x0018,0x1164);
 520    if ( strImagerPixelSpacing != GDCM_UNFOUND )
 521    {
 522       nbValues = sscanf( strImagerPixelSpacing.c_str(), "%f", &yspacing);
 523
 524    // if sscanf cannot read any float value, it won't affect yspacing
 525       if ( nbValues == 0 )
 526             yspacing = 1.0;
 527
 528       if ( yspacing == 0.0 )
 529             yspacing = 1.0;
 530
 531       return yspacing;
 532    }
 533
 534    std::string strSpacing = GetEntryValue(0x0028,0x0030);
 535    if ( strSpacing == GDCM_UNFOUND )
 536    {
 537       gdcmWarningMacro("Unfound Pixel Spacing (0028,0030)");
 538       return 1.;
 539     }
 540
 541    // if sscanf cannot read any float value, it won't affect yspacing
 542    nbValues = sscanf( strSpacing.c_str(), "%f", &yspacing);
 543
 544    // if no values, yspacing is set to 1.0
 545    if ( nbValues == 0 )
 546       yspacing = 1.0;
 547
 548    if ( yspacing == 0.0 )
 549       yspacing = 1.0;
 550
 551    return yspacing;
 552 }
 553
 554 /**
 555  * \brief gets the info from 0018,0088 : Space Between Slices
 556  *                 else from 0018,0050 : Slice Thickness
 557  *                 else 1.0
 558  * @return Z dimension of a voxel-to be
 559  */
 560 float File::GetZSpacing()
 561 {
 562    // Spacing Between Slices : distance between the middle of 2 slices
 563    // Slices may be :
 564    //   jointives     (Spacing between Slices = Slice Thickness)
 565    //   overlapping   (Spacing between Slices < Slice Thickness)
 566    //   disjointes    (Spacing between Slices > Slice Thickness)
 567    // Slice Thickness : epaisseur de tissus sur laquelle est acquis le signal
 568    //   It only concerns the MRI guys, not people wanting to visualize volumes
 569    //   If Spacing Between Slices is missing,
 570    //   we suppose slices joint together
 571
 572    const std::string &strSpacingBSlices = GetEntryValue(0x0018,0x0088);
 573
 574    if ( strSpacingBSlices == GDCM_UNFOUND )
 575    {
 576       gdcmWarningMacro("Unfound Spacing Between Slices (0018,0088)");
 577       const std::string &strSliceThickness = GetEntryValue(0x0018,0x0050);
 578       if ( strSliceThickness == GDCM_UNFOUND )
 579       {
 580          gdcmWarningMacro("Unfound Slice Thickness (0018,0050)");
 581          return 1.;
 582       }
 583       else
 584       {
 585          // if no 'Spacing Between Slices' is found,
 586          // we assume slices join together
 587          // (no overlapping, no interslice gap)
 588          // if they don't, we're fucked up
 589          return (float)atof( strSliceThickness.c_str() );
 590       }
 591    }
 592    //else
 593    return (float)atof( strSpacingBSlices.c_str() );
 594 }
 595
 596 /**
 597  * \brief gets the info from 0020,0032 : Image Position Patient
 598  *                 else from 0020,0030 : Image Position (RET)
 599  *                 else 0.
 600  * @return up-left image corner X position
 601  */
 602 float File::GetXOrigin()
 603 {
 604    float xImPos, yImPos, zImPos;
 605    std::string strImPos = GetEntryValue(0x0020,0x0032);
 606
 607    if ( strImPos == GDCM_UNFOUND )
 608    {
 609       gdcmWarningMacro( "Unfound Image Position Patient (0020,0032)");
 610       strImPos = GetEntryValue(0x0020,0x0030); // For ACR-NEMA images
 611       if ( strImPos == GDCM_UNFOUND )
 612       {
 613          gdcmWarningMacro( "Unfound Image Position (RET) (0020,0030)");
 614          return 0.;
 615       }
 616    }
 617
 618    if ( sscanf( strImPos.c_str(), "%f \\%f \\%f ", &xImPos, &yImPos, &zImPos) != 3 )
 619    {
 620       return 0.;
 621    }
 622
 623    return xImPos;
 624 }
 625
 626 /**
 627  * \brief gets the info from 0020,0032 : Image Position Patient
 628  *                 else from 0020,0030 : Image Position (RET)
 629  *                 else 0.
 630  * @return up-left image corner Y position
 631  */
 632 float File::GetYOrigin()
 633 {
 634    float xImPos, yImPos, zImPos;
 635    std::string strImPos = GetEntryValue(0x0020,0x0032);
 636
 637    if ( strImPos == GDCM_UNFOUND)
 638    {
 639       gdcmWarningMacro( "Unfound Image Position Patient (0020,0032)");
 640       strImPos = GetEntryValue(0x0020,0x0030); // For ACR-NEMA images
 641       if ( strImPos == GDCM_UNFOUND )
 642       {
 643          gdcmWarningMacro( "Unfound Image Position (RET) (0020,0030)");
 644          return 0.;
 645       }
 646    }
 647
 648    if ( sscanf( strImPos.c_str(), "%f \\%f \\%f ", &xImPos, &yImPos, &zImPos) != 3 )
 649    {
 650       return 0.;
 651    }
 652
 653    return yImPos;
 654 }
 655
 656 /**
 657  * \brief gets the info from 0020,0032 : Image Position Patient
 658  *                 else from 0020,0030 : Image Position (RET)
 659  *                 else from 0020,1041 : Slice Location
 660  *                 else from 0020,0050 : Location
 661  *                 else 0.
 662  * @return up-left image corner Z position
 663  */
 664 float File::GetZOrigin()
 665 {
 666    float xImPos, yImPos, zImPos;
 667    std::string strImPos = GetEntryValue(0x0020,0x0032);
 668
 669    if ( strImPos != GDCM_UNFOUND )
 670    {
 671       if ( sscanf( strImPos.c_str(), "%f \\%f \\%f ", &xImPos, &yImPos, &zImPos) != 3)
 672       {
 673          gdcmWarningMacro( "Wrong Image Position Patient (0020,0032)");
 674          return 0.;  // bug in the element 0x0020,0x0032
 675       }
 676       else
 677       {
 678          return zImPos;
 679       }
 680    }
 681
 682    strImPos = GetEntryValue(0x0020,0x0030); // For ACR-NEMA images
 683    if ( strImPos != GDCM_UNFOUND )
 684    {
 685       if ( sscanf( strImPos.c_str(),
 686           "%f \\%f \\%f ", &xImPos, &yImPos, &zImPos ) != 3 )
 687       {
 688          gdcmWarningMacro( "Wrong Image Position (RET) (0020,0030)");
 689          return 0.;  // bug in the element 0x0020,0x0032
 690       }
 691       else
 692       {
 693          return zImPos;
 694       }
 695    }
 696
 697    // for *very* old ACR-NEMA images
 698    std::string strSliceLocation = GetEntryValue(0x0020,0x1041);
 699    if ( strSliceLocation != GDCM_UNFOUND )
 700    {
 701       if ( sscanf( strSliceLocation.c_str(), "%f ", &zImPos) != 1)
 702       {
 703          gdcmWarningMacro( "Wrong Slice Location (0020,1041)");
 704          return 0.;  // bug in the element 0x0020,0x1041
 705       }
 706       else
 707       {
 708          return zImPos;
 709       }
 710    }
 711    gdcmWarningMacro( "Unfound Slice Location (0020,1041)");
 712
 713    std::string strLocation = GetEntryValue(0x0020,0x0050);
 714    if ( strLocation != GDCM_UNFOUND )
 715    {
 716       if ( sscanf( strLocation.c_str(), "%f ", &zImPos) != 1 )
 717       {
 718          gdcmWarningMacro( "Wrong Location (0020,0050)");
 719          return 0.;  // bug in the element 0x0020,0x0050
 720       }
 721       else
 722       {
 723          return zImPos;
 724       }
 725    }
 726    gdcmWarningMacro( "Unfound Location (0020,0050)");
 727
 728    return 0.; // Hopeless
 729 }
 730
 731 /**
 732   * \brief gets the info from 0020,0037 : Image Orientation Patient
 733   * (needed to organize DICOM files based on their x,y,z position)
 734   * @param iop adress of the (6)float array to receive values
 735   * @return cosines of image orientation patient
 736   */
 737 bool File::GetImageOrientationPatient( float iop[6] )
 738 {
 739    std::string strImOriPat;
 740    //iop is supposed to be float[6]
 741    iop[0] = iop[1] = iop[2] = iop[3] = iop[4] = iop[5] = 0.;
 742
 743    // 0020 0037 DS REL Image Orientation (Patient)
 744    if ( (strImOriPat = GetEntryValue(0x0020,0x0037)) != GDCM_UNFOUND )
 745    {
 746       if ( sscanf( strImOriPat.c_str(), "%f \\ %f \\%f \\%f \\%f \\%f ",
 747           &iop[0], &iop[1], &iop[2], &iop[3], &iop[4], &iop[5]) != 6 )
 748       {
 749          gdcmWarningMacro( "Wrong Image Orientation Patient (0020,0037). Less than 6 values were found." );
 750          return false;
 751       }
 752    }
 753    //For ACR-NEMA
 754    // 0020 0035 DS REL Image Orientation (RET)
 755    else if ( (strImOriPat = GetEntryValue(0x0020,0x0035)) != GDCM_UNFOUND )
 756    {
 757       if ( sscanf( strImOriPat.c_str(), "%f \\ %f \\%f \\%f \\%f \\%f ",
 758           &iop[0], &iop[1], &iop[2], &iop[3], &iop[4], &iop[5]) != 6 )
 759       {
 760          gdcmWarningMacro( "wrong Image Orientation Patient (0020,0035). Less than 6 values were found." );
 761          return false;
 762       }
 763    }
 764    return true;
 765 }
 766
 767 /**
 768  * \brief   Retrieve the number of Bits Stored (actually used)
 769  *          (as opposed to number of Bits Allocated)
 770  * @return  The encountered number of Bits Stored, 0 by default.
 771  *          0 means the file is NOT USABLE. The caller has to check it !
 772  */
 773 int File::GetBitsStored()
 774 {
 775    std::string strSize = GetEntryValue( 0x0028, 0x0101 );
 776    if ( strSize == GDCM_UNFOUND )
 777    {
 778       gdcmWarningMacro("(0028,0101) is supposed to be mandatory");
 779       return 0;  // It's supposed to be mandatory
 780                  // the caller will have to check
 781    }
 782    return atoi( strSize.c_str() );
 783 }
 784
 785 /**
 786  * \brief   Retrieve the number of Bits Allocated
 787  *          (8, 12 -compacted ACR-NEMA files-, 16, ...)
 788  * @return  The encountered number of Bits Allocated, 0 by default.
 789  *          0 means the file is NOT USABLE. The caller has to check it !
 790  */
 791 int File::GetBitsAllocated()
 792 {
 793    std::string strSize = GetEntryValue(0x0028,0x0100);
 794    if ( strSize == GDCM_UNFOUND  )
 795    {
 796       gdcmWarningMacro( "(0028,0100) is supposed to be mandatory");
 797       return 0; // It's supposed to be mandatory
 798                 // the caller will have to check
 799    }
 800    return atoi( strSize.c_str() );
 801 }
 802
 803 /**
 804  * \brief   Retrieve the high bit position.
 805  * \warning The method defaults to 0 when information is missing.
 806  *          The responsability of checking this value is left to the caller.
 807  * @return  The high bit position when present. 0 when missing.
 808  */
 809 int File::GetHighBitPosition()
 810 {
 811    std::string strSize = GetEntryValue( 0x0028, 0x0102 );
 812    if ( strSize == GDCM_UNFOUND )
 813    {
 814       gdcmWarningMacro( "(0028,0102) is supposed to be mandatory");
 815       return 0;
 816    }
 817    return atoi( strSize.c_str() );
 818 }
 819
 820 /**
 821  * \brief   Retrieve the number of Samples Per Pixel
 822  *          (1 : gray level, 3 : RGB/YBR -1 or 3 Planes-)
 823  * @return  The encountered number of Samples Per Pixel, 1 by default.
 824  *          (we assume Gray level Pixels)
 825  */
 826 int File::GetSamplesPerPixel()
 827 {
 828    const std::string &strSize = GetEntryValue(0x0028,0x0002);
 829    if ( strSize == GDCM_UNFOUND )
 830    {
 831       gdcmWarningMacro( "(0028,0002) is supposed to be mandatory");
 832       return 1; // Well, it's supposed to be mandatory ...
 833                 // but sometimes it's missing : *we* assume Gray pixels
 834    }
 835    return atoi( strSize.c_str() );
 836 }
 837
 838 /**
 839  * \brief   Retrieve the Planar Configuration for RGB images
 840  *          (0 : RGB Pixels , 1 : R Plane + G Plane + B Plane)
 841  * @return  The encountered Planar Configuration, 0 by default.
 842  */
 843 int File::GetPlanarConfiguration()
 844 {
 845    std::string strSize = GetEntryValue(0x0028,0x0006);
 846    if ( strSize == GDCM_UNFOUND )
 847    {
 848       gdcmWarningMacro( "Not found : Planar Configuration (0028,0006)");
 849       return 0;
 850    }
 851    return atoi( strSize.c_str() );
 852 }
 853
 854 /**
 855  * \brief   Return the size (in bytes) of a single pixel of data.
 856  * @return  The size in bytes of a single pixel of data; 0 by default
 857  *          0 means the file is NOT USABLE; the caller will have to check
 858  */
 859 int File::GetPixelSize()
 860 {
 861    // 0028 0100 US IMG Bits Allocated
 862    // (in order no to be messed up by old ACR-NEMA RGB images)
 863    //   if (File::GetEntryValue(0x0028,0x0100) == "24")
 864    //      return 3;
 865
 866    std::string pixelType = GetPixelType();
 867    if ( pixelType ==  "8U" || pixelType == "8S" )
 868    {
 869       return 1;
 870    }
 871    if ( pixelType == "16U" || pixelType == "16S")
 872    {
 873       return 2;
 874    }
 875    if ( pixelType == "32U" || pixelType == "32S")
 876    {
 877       return 4;
 878    }
 879    if ( pixelType == "FD" )
 880    {
 881       return 8;
 882    }
 883    gdcmWarningMacro( "Unknown pixel type");
 884    return 0;
 885 }
 886
 887 /**
 888  * \brief   Build the Pixel Type of the image.
 889  *          Possible values are:
 890  *          - 8U  unsigned  8 bit,
 891  *          - 8S    signed  8 bit,
 892  *          - 16U unsigned 16 bit,
 893  *          - 16S   signed 16 bit,
 894  *          - 32U unsigned 32 bit,
 895  *          - 32S   signed 32 bit,
 896  *          - FD floating double 64 bits (Not kosher DICOM, but so usefull!)
 897  * \warning 12 bit images appear as 16 bit.
 898  *          24 bit images appear as 8 bit + photochromatic interp ="RGB "
 899  *                                        + Planar Configuration = 0
 900  * @return  0S if nothing found. NOT USABLE file. The caller has to check
 901  */
 902 std::string File::GetPixelType()
 903 {
 904    std::string bitsAlloc = GetEntryValue(0x0028, 0x0100); // Bits Allocated
 905    if ( bitsAlloc == GDCM_UNFOUND )
 906    {
 907       gdcmWarningMacro( "Missing  Bits Allocated (0028,0100)");
 908       bitsAlloc = "16"; // default and arbitrary value, not to polute the output
 909    }
 910
 911    if ( bitsAlloc == "64" )
 912    {
 913       return "FD";
 914    }
 915    else if ( bitsAlloc == "12" )
 916    {
 917       // It will be unpacked
 918       bitsAlloc = "16";
 919    }
 920    else if ( bitsAlloc == "24" )
 921    {
 922       // (in order no to be messed up
 923       bitsAlloc = "8";  // by old RGB images)
 924    }
 925
 926    std::string sign = GetEntryValue(0x0028, 0x0103);//"Pixel Representation"
 927
 928    if (sign == GDCM_UNFOUND )
 929    {
 930       gdcmWarningMacro( "Missing Pixel Representation (0028,0103)");
 931       sign = "U"; // default and arbitrary value, not to polute the output
 932    }
 933    else if ( sign == "0" )
 934    {
 935       sign = "U";
 936    }
 937    else
 938    {
 939       sign = "S";
 940    }
 941    return bitsAlloc + sign;
 942 }
 943
 944 /**
 945  * \brief   Check whether the pixels are signed (1) or UNsigned (0) data.
 946  * \warning The method defaults to false (UNsigned) when tag 0028|0103
 947  *          is missing.
 948  *          The responsability of checking this value is left to the caller
 949  *          (NO transformation is performed on the pixels to make then >0)
 950  * @return  True when signed, false when UNsigned
 951  */
 952 bool File::IsSignedPixelData()
 953 {
 954    std::string strSign = GetEntryValue( 0x0028, 0x0103 );
 955    if ( strSign == GDCM_UNFOUND )
 956    {
 957       gdcmWarningMacro( "(0028,0103) is supposed to be mandatory");
 958       return false;
 959    }
 960    int sign = atoi( strSign.c_str() );
 961    if ( sign == 0 )
 962    {
 963       return false;
 964    }
 965    return true;
 966 }
 967
 968 /**
 969  * \brief   Check whether this a monochrome picture (gray levels) or not,
 970  *          using "Photometric Interpretation" tag (0x0028,0x0004).
 971  * @return  true when "MONOCHROME1" or "MONOCHROME2". False otherwise.
 972  */
 973 bool File::IsMonochrome()
 974 {
 975    const std::string &PhotometricInterp = GetEntryValue( 0x0028, 0x0004 );
 976    if (  Util::DicomStringEqual(PhotometricInterp, "MONOCHROME1")
 977       || Util::DicomStringEqual(PhotometricInterp, "MONOCHROME2") )
 978    {
 979       return true;
 980    }
 981    if ( PhotometricInterp == GDCM_UNFOUND )
 982    {
 983       gdcmWarningMacro( "Not found : Photometric Interpretation (0028,0004)");
 984    }
 985    return false;
 986 }
 987
 988 /**
 989  * \brief   Check whether this a MONOCHROME1 picture (high values = dark)
 990  *            or not using "Photometric Interpretation" tag (0x0028,0x0004).
 991  * @return  true when "MONOCHROME1" . False otherwise.
 992  */
 993 bool File::IsMonochrome1()
 994 {
 995    const std::string &PhotometricInterp = GetEntryValue( 0x0028, 0x0004 );
 996    if (  Util::DicomStringEqual(PhotometricInterp, "MONOCHROME1") )
 997    {
 998       return true;
 999    }
1000    if ( PhotometricInterp == GDCM_UNFOUND )
1001    {
1002       gdcmWarningMacro( "Not found : Photometric Interpretation (0028,0004)");
1003    }
1004    return false;
1005 }
1006
1007 /**
1008  * \brief   Check whether this a "PALETTE COLOR" picture or not by accessing
1009  *          the "Photometric Interpretation" tag ( 0x0028, 0x0004 ).
1010  * @return  true when "PALETTE COLOR". False otherwise.
1011  */
1012 bool File::IsPaletteColor()
1013 {
1014    std::string PhotometricInterp = GetEntryValue( 0x0028, 0x0004 );
1015    if (   PhotometricInterp == "PALETTE COLOR " )
1016    {
1017       return true;
1018    }
1019    if ( PhotometricInterp == GDCM_UNFOUND )
1020    {
1021       gdcmWarningMacro( "Not found : Palette color (0028,0004)");
1022    }
1023    return false;
1024 }
1025
1026 /**
1027  * \brief   Check whether this a "YBR_FULL" color picture or not by accessing
1028  *          the "Photometric Interpretation" tag ( 0x0028, 0x0004 ).
1029  * @return  true when "YBR_FULL". False otherwise.
1030  */
1031 bool File::IsYBRFull()
1032 {
1033    std::string PhotometricInterp = GetEntryValue( 0x0028, 0x0004 );
1034    if (   PhotometricInterp == "YBR_FULL" )
1035    {
1036       return true;
1037    }
1038    if ( PhotometricInterp == GDCM_UNFOUND )
1039    {
1040       gdcmWarningMacro( "Not found : YBR Full (0028,0004)");
1041    }
1042    return false;
1043 }
1044
1045 /**
1046   * \brief tells us if LUT are used
1047   * \warning Right now, 'Segmented xxx Palette Color Lookup Table Data'
1048   *          are NOT considered as LUT, since nobody knows
1049   *          how to deal with them
1050   *          Please warn me if you know sbdy that *does* know ... jprx
1051   * @return true if LUT Descriptors and LUT Tables were found
1052   */
1053 bool File::HasLUT()
1054 {
1055    // Check the presence of the LUT Descriptors, and LUT Tables
1056    // LutDescriptorRed
1057    if ( !GetDocEntry(0x0028,0x1101) )
1058    {
1059       return false;
1060    }
1061    // LutDescriptorGreen
1062    if ( !GetDocEntry(0x0028,0x1102) )
1063    {
1064       return false;
1065    }
1066    // LutDescriptorBlue
1067    if ( !GetDocEntry(0x0028,0x1103) )
1068    {
1069       return false;
1070    }
1071    // Red Palette Color Lookup Table Data
1072    if ( !GetDocEntry(0x0028,0x1201) )
1073    {
1074       return false;
1075    }
1076    // Green Palette Color Lookup Table Data
1077    if ( !GetDocEntry(0x0028,0x1202) )
1078    {
1079       return false;
1080    }
1081    // Blue Palette Color Lookup Table Data
1082    if ( !GetDocEntry(0x0028,0x1203) )
1083    {
1084       return false;
1085    }
1086
1087    // FIXME : (0x0028,0x3006) : LUT Data (CTX dependent)
1088    //         NOT taken into account, but we don't know how to use it ...
1089    return true;
1090 }
1091
1092 /**
1093   * \brief gets the info from 0028,1101 : Lookup Table Desc-Red
1094   *             else 0
1095   * @return Lookup Table number of Bits , 0 by default
1096   *          when (0028,0004),Photometric Interpretation = [PALETTE COLOR ]
1097   * @ return bit number of each LUT item
1098   */
1099 int File::GetLUTNbits()
1100 {
1101    std::vector<std::string> tokens;
1102    int lutNbits;
1103
1104    //Just hope Lookup Table Desc-Red = Lookup Table Desc-Red
1105    //                                = Lookup Table Desc-Blue
1106    // Consistency already checked in GetLUTLength
1107    std::string lutDescription = GetEntryValue(0x0028,0x1101);
1108    if ( lutDescription == GDCM_UNFOUND )
1109    {
1110       return 0;
1111    }
1112
1113    tokens.clear(); // clean any previous value
1114    Util::Tokenize ( lutDescription, tokens, "\\" );
1115    //LutLength=atoi(tokens[0].c_str());
1116    //LutDepth=atoi(tokens[1].c_str());
1117
1118    lutNbits = atoi( tokens[2].c_str() );
1119    tokens.clear();
1120
1121    return lutNbits;
1122 }
1123
1124 /**
1125  *\brief gets the info from 0028,1052 : Rescale Intercept
1126  * @return Rescale Intercept
1127  */
1128 float File::GetRescaleIntercept()
1129 {
1130    float resInter = 0.;
1131    /// 0028 1052 DS IMG Rescale Intercept
1132    const std::string &strRescInter = GetEntryValue(0x0028,0x1052);
1133    if ( strRescInter != GDCM_UNFOUND )
1134    {
1135       if ( sscanf( strRescInter.c_str(), "%f ", &resInter) != 1 )
1136       {
1137          // bug in the element 0x0028,0x1052
1138          gdcmWarningMacro( "Rescale Intercept (0028,1052) is empty." );
1139       }
1140    }
1141
1142    return resInter;
1143 }
1144
1145 /**
1146  *\brief   gets the info from 0028,1053 : Rescale Slope
1147  * @return Rescale Slope
1148  */
1149 float File::GetRescaleSlope()
1150 {
1151    float resSlope = 1.;
1152    //0028 1053 DS IMG Rescale Slope
1153    std::string strRescSlope = GetEntryValue(0x0028,0x1053);
1154    if ( strRescSlope != GDCM_UNFOUND )
1155    {
1156       if ( sscanf( strRescSlope.c_str(), "%f ", &resSlope) != 1 )
1157       {
1158          // bug in the element 0x0028,0x1053
1159          gdcmWarningMacro( "Rescale Slope (0028,1053) is empty.");
1160       }
1161    }
1162
1163    return resSlope;
1164 }
1165
1166 /**
1167  * \brief This function is intended to user who doesn't want
1168  *   to have to manage a LUT and expects to get an RBG Pixel image
1169  *   (or a monochrome one ...)
1170  * \warning to be used with GetImagePixels()
1171  * @return 1 if Gray level, 3 if Color (RGB, YBR, *or PALETTE COLOR*)
1172  */
1173 int File::GetNumberOfScalarComponents()
1174 {
1175    if ( GetSamplesPerPixel() == 3 )
1176    {
1177       return 3;
1178    }
1179
1180    // 0028 0100 US IMG Bits Allocated
1181    // (in order no to be messed up by old RGB images)
1182    if ( GetEntryValue(0x0028,0x0100) == "24" )
1183    {
1184       return 3;
1185    }
1186
1187    std::string strPhotometricInterpretation = GetEntryValue(0x0028,0x0004);
1188
1189    if ( ( strPhotometricInterpretation == "PALETTE COLOR ") )
1190    {
1191       if ( HasLUT() )// PALETTE COLOR is NOT enough
1192       {
1193          return 3;
1194       }
1195       else
1196       {
1197          return 1;
1198       }
1199    }
1200
1201    // beware of trailing space at end of string
1202    // DICOM tags are never of odd length
1203    if ( strPhotometricInterpretation == GDCM_UNFOUND   ||
1204         Util::DicomStringEqual(strPhotometricInterpretation, "MONOCHROME1") ||
1205         Util::DicomStringEqual(strPhotometricInterpretation, "MONOCHROME2") )
1206    {
1207       return 1;
1208    }
1209    else
1210    {
1211       // we assume that *all* kinds of YBR are dealt with
1212       return 3;
1213    }
1214 }
1215
1216 /**
1217  * \brief This function is intended to user that DOESN'T want
1218  *  to get RGB pixels image when it's stored as a PALETTE COLOR image
1219  *   - the (vtk) user is supposed to know how deal with LUTs -
1220  * \warning to be used with GetImagePixelsRaw()
1221  * @return 1 if Gray level, 3 if Color (RGB or YBR - NOT 'PALETTE COLOR' -)
1222  */
1223 int File::GetNumberOfScalarComponentsRaw()
1224 {
1225    // 0028 0100 US IMG Bits Allocated
1226    // (in order no to be messed up by old RGB images)
1227    if ( File::GetEntryValue(0x0028,0x0100) == "24" )
1228    {
1229       return 3;
1230    }
1231
1232    // we assume that *all* kinds of YBR are dealt with
1233    return GetSamplesPerPixel();
1234 }
1235
1236 /**
1237  * \brief   Recover the offset (from the beginning of the file)
1238  *          of *image* pixels (not *icone image* pixels, if any !)
1239  * @return Pixel Offset
1240  */
1241 size_t File::GetPixelOffset()
1242 {
1243    DocEntry *pxlElement = GetDocEntry(GrPixel, NumPixel);
1244    if ( pxlElement )
1245    {
1246       return pxlElement->GetOffset();
1247    }
1248    else
1249    {
1250       gdcmDebugMacro( "Big trouble : Pixel Element ("
1251                       << std::hex << GrPixel<<","<< NumPixel<< ") NOT found" );
1252       return 0;
1253    }
1254 }
1255
1256 /**
1257  * \brief   Recover the pixel area length (in Bytes)
1258  * @return Pixel Element Length, as stored in the header
1259  *         (NOT the memory space necessary to hold the Pixels
1260  *          -in case of embeded compressed image-)
1261  *         0 : NOT USABLE file. The caller has to check.
1262  */
1263 size_t File::GetPixelAreaLength()
1264 {
1265    DocEntry *pxlElement = GetDocEntry(GrPixel, NumPixel);
1266    if ( pxlElement )
1267    {
1268       return pxlElement->GetLength();
1269    }
1270    else
1271    {
1272       gdcmDebugMacro( "Big trouble : Pixel Element ("
1273                       << std::hex << GrPixel<<","<< NumPixel<< ") NOT found" );
1274       return 0;
1275    }
1276 }
1277
1278 /**
1279  * \brief Adds the characteristics of a new element we want to anonymize
1280  * @param   group  Group number of the target tag.
1281  * @param   elem Element number of the target tag.
1282  * @param   value new value (string) to substitute with
1283  */
1284 void File::AddAnonymizeElement (uint16_t group, uint16_t elem,
1285                                 std::string const &value)
1286 {
1287    Element el;
1288    el.Group = group;
1289    el.Elem  = elem;
1290    el.Value = value;
1291    AnonymizeList.push_back(el);
1292 }
1293
1294 /**
1295  * \brief Overwrites in the file the values of the DicomElements
1296  *       held in the list
1297  */
1298 void File::AnonymizeNoLoad()
1299 {
1300    std::fstream *fp = new std::fstream(Filename.c_str(),
1301                               std::ios::in | std::ios::out | std::ios::binary);
1302    gdcm::DocEntry *d;
1303    uint32_t offset;
1304    uint32_t lgth;
1305    uint32_t valLgth = 0;
1306    std::string *spaces;
1307    for (ListElements::iterator it = AnonymizeList.begin();
1308                                it != AnonymizeList.end();
1309                              ++it)
1310    {
1311       d = GetDocEntry( (*it).Group, (*it).Elem);
1312
1313       if ( d == NULL)
1314          continue;
1315
1316          if ( dynamic_cast<SeqEntry *>(d) )
1317          {
1318             gdcmWarningMacro( "You cannot 'Anonymize a SeqEntry ");
1319             continue;
1320          }
1321
1322       offset = d->GetOffset();
1323       lgth =   d->GetLength();
1324       if (valLgth < lgth)
1325       {
1326          spaces = new std::string( lgth-valLgth, ' ');
1327          (*it).Value = (*it).Value + *spaces;
1328          delete spaces;
1329       }
1330       fp->seekp( offset, std::ios::beg );
1331       fp->write( (*it).Value.c_str(), lgth );
1332
1333    }
1334    fp->close();
1335    delete fp;
1336 }
1337
1338 /**
1339  * \brief anonymize a File (remove Patient's personal info passed with
1340  *        AddAnonymizeElement()
1341  * \note You cannot Anonymize a BinEntry (to be fixed)
1342  */
1343 bool File::AnonymizeFile()
1344 {
1345    // If Anonymisation list is empty, let's perform some basic anonymization
1346    if ( AnonymizeList.begin() == AnonymizeList.end() )
1347    {
1348       // If exist, replace by spaces
1349       SetValEntry ("  ",0x0010, 0x2154); // Telephone
1350       SetValEntry ("  ",0x0010, 0x1040); // Adress
1351       SetValEntry ("  ",0x0010, 0x0020); // Patient ID
1352
1353       DocEntry* patientNameHE = GetDocEntry (0x0010, 0x0010);
1354
1355       if ( patientNameHE ) // we replace it by Study Instance UID (why not ?)
1356       {
1357          std::string studyInstanceUID =  GetEntryValue (0x0020, 0x000d);
1358          if ( studyInstanceUID != GDCM_UNFOUND )
1359          {
1360             SetValEntry(studyInstanceUID, 0x0010, 0x0010);
1361          }
1362          else
1363          {
1364             SetValEntry("anonymised", 0x0010, 0x0010);
1365          }
1366       }
1367    }
1368    else
1369    {
1370       gdcm::DocEntry *d;
1371       for (ListElements::iterator it = AnonymizeList.begin();
1372                                   it != AnonymizeList.end();
1373                                 ++it)
1374       {
1375          d = GetDocEntry( (*it).Group, (*it).Elem);
1376
1377          if ( d == NULL)
1378             continue;
1379
1380          if ( dynamic_cast<SeqEntry *>(d) )
1381          {
1382             gdcmWarningMacro( "You cannot 'Anonymize' a SeqEntry ");
1383             continue;
1384          }
1385
1386          if ( dynamic_cast<BinEntry *>(d) )
1387          {
1388             gdcmWarningMacro( "To 'Anonymize' a BinEntry, better use AnonymizeNoLoad (FIXME) ");
1389             continue;
1390          }
1391          else
1392             SetValEntry ((*it).Value, (*it).Group, (*it).Elem);
1393       }
1394 }
1395
1396   // In order to make definitively impossible any further identification
1397   // remove or replace all the stuff that contains a Date
1398
1399 //0008 0012 DA ID Instance Creation Date
1400 //0008 0020 DA ID Study Date
1401 //0008 0021 DA ID Series Date
1402 //0008 0022 DA ID Acquisition Date
1403 //0008 0023 DA ID Content Date
1404 //0008 0024 DA ID Overlay Date
1405 //0008 0025 DA ID Curve Date
1406 //0008 002a DT ID Acquisition Datetime
1407 //0018 9074 DT ACQ Frame Acquisition Datetime
1408 //0018 9151 DT ACQ Frame Reference Datetime
1409 //0018 a002 DT ACQ Contribution Date Time
1410 //0020 3403 SH REL Modified Image Date (RET)
1411 //0032 0032 DA SDY Study Verified Date
1412 //0032 0034 DA SDY Study Read Date
1413 //0032 1000 DA SDY Scheduled Study Start Date
1414 //0032 1010 DA SDY Scheduled Study Stop Date
1415 //0032 1040 DA SDY Study Arrival Date
1416 //0032 1050 DA SDY Study Completion Date
1417 //0038 001a DA VIS Scheduled Admission Date
1418 //0038 001c DA VIS Scheduled Discharge Date
1419 //0038 0020 DA VIS Admitting Date
1420 //0038 0030 DA VIS Discharge Date
1421 //0040 0002 DA PRC Scheduled Procedure Step Start Date
1422 //0040 0004 DA PRC Scheduled Procedure Step End Date
1423 //0040 0244 DA PRC Performed Procedure Step Start Date
1424 //0040 0250 DA PRC Performed Procedure Step End Date
1425 //0040 2004 DA PRC Issue Date of Imaging Service Request
1426 //0040 4005 DT PRC Scheduled Procedure Step Start Date and Time
1427 //0040 4011 DT PRC Expected Completion Date and Time
1428 //0040 a030 DT PRC Verification Date Time
1429 //0040 a032 DT PRC Observation Date Time
1430 //0040 a120 DT PRC DateTime
1431 //0040 a121 DA PRC Date
1432 //0040 a13a DT PRC Referenced Datetime
1433 //0070 0082 DA ??? Presentation Creation Date
1434 //0100 0420 DT ??? SOP Autorization Date and Time
1435 //0400 0105 DT ??? Digital Signature DateTime
1436 //2100 0040 DA PJ Creation Date
1437 //3006 0008 DA SSET Structure Set Date
1438 //3008 0024 DA ??? Treatment Control Point Date
1439 //3008 0054 DA ??? First Treatment Date
1440 //3008 0056 DA ??? Most Recent Treatment Date
1441 //3008 0162 DA ??? Safe Position Exit Date
1442 //3008 0166 DA ??? Safe Position Return Date
1443 //3008 0250 DA ??? Treatment Date
1444 //300a 0006 DA RT RT Plan Date
1445 //300a 022c DA RT Air Kerma Rate Reference Date
1446 //300e 0004 DA RT Review Date
1447
1448    return true;
1449 }
1450
1451 /**
1452  * \brief Performs some consistency checking on various 'File related'
1453  *       (as opposed to 'DicomDir related') entries
1454  *       then writes in a file all the (Dicom Elements) included the Pixels
1455  * @param fileName file name to write to
1456  * @param writetype type of the file to be written
1457  *          (ACR, ExplicitVR, ImplicitVR)
1458  */
1459 bool File::Write(std::string fileName, FileType writetype)
1460 {
1461    std::ofstream *fp = new std::ofstream(fileName.c_str(),
1462                                          std::ios::out | std::ios::binary);
1463    if (*fp == NULL)
1464    {
1465       gdcmWarningMacro("Failed to open (write) File: " << fileName.c_str());
1466       return false;
1467    }
1468
1469    // Entry : 0002|0000 = group length -> recalculated
1470    ValEntry*e0000 = GetValEntry(0x0002,0x0000);
1471    if ( e0000 )
1472    {
1473       std::ostringstream sLen;
1474       sLen << ComputeGroup0002Length(writetype);
1475       e0000->SetValue(sLen.str());
1476    }
1477
1478    int i_lgPix = GetEntryLength(GrPixel, NumPixel);
1479    if (i_lgPix != -2)
1480    {
1481       // no (GrPixel, NumPixel) element
1482       std::string s_lgPix = Util::Format("%d", i_lgPix+12);
1483       s_lgPix = Util::DicomString( s_lgPix.c_str() );
1484       InsertValEntry(s_lgPix,GrPixel, 0x0000);
1485    }
1486
1487    Document::WriteContent(fp, writetype);
1488
1489    fp->close();
1490    delete fp;
1491
1492    return true;
1493 }
1494
1495 //-----------------------------------------------------------------------------
1496 // Protected
1497
1498
1499 //-----------------------------------------------------------------------------
1500 // Private
1501 /**
1502  * \brief Parse pixel data from disk of [multi-]fragment RLE encoding.
1503  *        Compute the RLE extra information and store it in \ref RLEInfo
1504  *        for later pixel retrieval usage.
1505  */
1506 void File::ComputeRLEInfo()
1507 {
1508    std::string ts = GetTransferSyntax();
1509    if ( !Global::GetTS()->IsRLELossless(ts) )
1510    {
1511       return;
1512    }
1513
1514    // Encoded pixel data: for the time being we are only concerned with
1515    // Jpeg or RLE Pixel data encodings.
1516    // As stated in PS 3.5-2003, section 8.2 p44:
1517    // "If sent in Encapsulated Format (i.e. other than the Native Format) the
1518    //  value representation OB is used".
1519    // Hence we expect an OB value representation. Concerning OB VR,
1520    // the section PS 3.5-2003, section A.4.c p 58-59, states:
1521    // "For the Value Representations OB and OW, the encoding shall meet the
1522    //   following specifications depending on the Data element tag:"
1523    //   [...snip...]
1524    //    - the first item in the sequence of items before the encoded pixel
1525    //      data stream shall be basic offset table item. The basic offset table
1526    //      item value, however, is not required to be present"
1527    ReadAndSkipEncapsulatedBasicOffsetTable();
1528
1529    // Encapsulated RLE Compressed Images (see PS 3.5-2003, Annex G)
1530    // Loop on the individual frame[s] and store the information
1531    // on the RLE fragments in a RLEFramesInfo.
1532    // Note: - when only a single frame is present, this is a
1533    //         classical image.
1534    //       - when more than one frame are present, then we are in
1535    //         the case of a multi-frame image.
1536    long frameLength;
1537    while ( (frameLength = ReadTagLength(0xfffe, 0xe000)) != 0 )
1538    {
1539       // Parse the RLE Header and store the corresponding RLE Segment
1540       // Offset Table information on fragments of this current Frame.
1541       // Note that the fragment pixels themselves are not loaded
1542       // (but just skipped).
1543       long frameOffset = Fp->tellg();
1544
1545       uint32_t nbRleSegments = ReadInt32();
1546       if ( nbRleSegments > 16 )
1547       {
1548          // There should be at most 15 segments (refer to RLEFrame class)
1549          gdcmWarningMacro( "Too many segments.");
1550       }
1551
1552       uint32_t rleSegmentOffsetTable[16];
1553       for( int k = 1; k <= 15; k++ )
1554       {
1555          rleSegmentOffsetTable[k] = ReadInt32();
1556       }
1557
1558       // Deduce from both RLE Header and frameLength
1559       // the fragment length, and again store this info
1560       // in a RLEFramesInfo.
1561       long rleSegmentLength[15];
1562       // skipping (not reading) RLE Segments
1563       if ( nbRleSegments > 1)
1564       {
1565          for(unsigned int k = 1; k <= nbRleSegments-1; k++)
1566          {
1567              rleSegmentLength[k] =  rleSegmentOffsetTable[k+1]
1568                                   - rleSegmentOffsetTable[k];
1569              SkipBytes(rleSegmentLength[k]);
1570           }
1571        }
1572
1573        rleSegmentLength[nbRleSegments] = frameLength
1574                                       - rleSegmentOffsetTable[nbRleSegments];
1575        SkipBytes(rleSegmentLength[nbRleSegments]);
1576
1577        // Store the collected info
1578        RLEFrame *newFrame = new RLEFrame;
1579        newFrame->SetNumberOfFragments(nbRleSegments);
1580        for( unsigned int uk = 1; uk <= nbRleSegments; uk++ )
1581        {
1582           newFrame->SetOffset(uk,frameOffset + rleSegmentOffsetTable[uk]);
1583           newFrame->SetLength(uk,rleSegmentLength[uk]);
1584        }
1585        RLEInfo->AddFrame(newFrame);
1586    }
1587
1588    // Make sure that  we encounter a 'Sequence Delimiter Item'
1589    // at the end of the item :
1590    if ( !ReadTag(0xfffe, 0xe0dd) )
1591    {
1592       gdcmWarningMacro( "No sequence delimiter item at end of RLE item sequence");
1593    }
1594 }
1595
1596 /**
1597  * \brief Parse pixel data from disk of [multi-]fragment Jpeg encoding.
1598  *        Compute the jpeg extra information (fragment[s] offset[s] and
1599  *        length) and store it[them] in \ref JPEGInfo for later pixel
1600  *        retrieval usage.
1601  */
1602 void File::ComputeJPEGFragmentInfo()
1603 {
1604    // If you need to, look for comments of ComputeRLEInfo().
1605    std::string ts = GetTransferSyntax();
1606    if ( ! Global::GetTS()->IsJPEG(ts) )
1607    {
1608       return;
1609    }
1610
1611    ReadAndSkipEncapsulatedBasicOffsetTable();
1612
1613    // Loop on the fragments[s] and store the parsed information in a
1614    // JPEGInfo.
1615    long fragmentLength;
1616    while ( (fragmentLength = ReadTagLength(0xfffe, 0xe000)) != 0 )
1617    {
1618       long fragmentOffset = Fp->tellg();
1619
1620        // Store the collected info
1621        JPEGFragment *newFragment = new JPEGFragment;
1622        newFragment->SetOffset(fragmentOffset);
1623        newFragment->SetLength(fragmentLength);
1624        JPEGInfo->AddFragment(newFragment);
1625
1626        SkipBytes(fragmentLength);
1627    }
1628
1629    // Make sure that  we encounter a 'Sequence Delimiter Item'
1630    // at the end of the item :
1631    if ( !ReadTag(0xfffe, 0xe0dd) )
1632    {
1633       gdcmWarningMacro( "No sequence delimiter item at end of JPEG item sequence");
1634    }
1635 }
1636
1637 /**
1638  * \brief   Assuming the internal file pointer \ref Document::Fp
1639  *          is placed at the beginning of a tag check whether this
1640  *          tag is (TestGroup, TestElem).
1641  * \warning On success the internal file pointer \ref Document::Fp
1642  *          is modified to point after the tag.
1643  *          On failure (i.e. when the tag wasn't the expected tag
1644  *          (TestGroup, TestElem) the internal file pointer
1645  *          \ref Document::Fp is restored to it's original position.
1646  * @param   testGroup The expected group   of the tag.
1647  * @param   testElem  The expected Element of the tag.
1648  * @return  True on success, false otherwise.
1649  */
1650 bool File::ReadTag(uint16_t testGroup, uint16_t testElem)
1651 {
1652    long positionOnEntry = Fp->tellg();
1653    long currentPosition = Fp->tellg();          // On debugging purposes
1654
1655    // Read the Item Tag group and element, and make
1656    // sure they are what we expected:
1657    uint16_t itemTagGroup;
1658    uint16_t itemTagElem;
1659    try
1660    {
1661       itemTagGroup = ReadInt16();
1662       itemTagElem  = ReadInt16();
1663    }
1664    catch ( FormatError e )
1665    {
1666       //std::cerr << e << std::endl;
1667       return false;
1668    }
1669    if ( itemTagGroup != testGroup || itemTagElem != testElem )
1670    {
1671       gdcmWarningMacro( "Wrong Item Tag found:"
1672        << "   We should have found tag ("
1673        << std::hex << testGroup << "," << testElem << ")" << std::endl
1674        << "   but instead we encountered tag ("
1675        << std::hex << itemTagGroup << "," << itemTagElem << ")"
1676        << "  at address: " << "  0x(" << (unsigned int)currentPosition  << ")"
1677        ) ;
1678       Fp->seekg(positionOnEntry, std::ios::beg);
1679
1680       return false;
1681    }
1682    return true;
1683 }
1684
1685 /**
1686  * \brief   Assuming the internal file pointer \ref Document::Fp
1687  *          is placed at the beginning of a tag (TestGroup, TestElement),
1688  *          read the length associated to the Tag.
1689  * \warning On success the internal file pointer \ref Document::Fp
1690  *          is modified to point after the tag and it's length.
1691  *          On failure (i.e. when the tag wasn't the expected tag
1692  *          (TestGroup, TestElement) the internal file pointer
1693  *          \ref Document::Fp is restored to it's original position.
1694  * @param   testGroup The expected Group   of the tag.
1695  * @param   testElem  The expected Element of the tag.
1696  * @return  On success returns the length associated to the tag. On failure
1697  *          returns 0.
1698  */
1699 uint32_t File::ReadTagLength(uint16_t testGroup, uint16_t testElem)
1700 {
1701
1702    if ( !ReadTag(testGroup, testElem) )
1703    {
1704       return 0;
1705    }
1706
1707    //// Then read the associated Item Length
1708    long currentPosition = Fp->tellg();
1709    uint32_t itemLength  = ReadInt32();
1710    {
1711       gdcmWarningMacro( "Basic Item Length is: "
1712         << itemLength << std::endl
1713         << "  at address: " << std::hex << (unsigned int)currentPosition);
1714    }
1715    return itemLength;
1716 }
1717
1718 /**
1719  * \brief When parsing the Pixel Data of an encapsulated file, read
1720  *        the basic offset table (when present, and BTW dump it).
1721  */
1722 void File::ReadAndSkipEncapsulatedBasicOffsetTable()
1723 {
1724    //// Read the Basic Offset Table Item Tag length...
1725    uint32_t itemLength = ReadTagLength(0xfffe, 0xe000);
1726
1727    // When present, read the basic offset table itself.
1728    // Notes: - since the presence of this basic offset table is optional
1729    //          we can't rely on it for the implementation, and we will simply
1730    //          trash it's content (when present).
1731    //        - still, when present, we could add some further checks on the
1732    //          lengths, but we won't bother with such fuses for the time being.
1733    if ( itemLength != 0 )
1734    {
1735       char *basicOffsetTableItemValue = new char[itemLength + 1];
1736       Fp->read(basicOffsetTableItemValue, itemLength);
1737
1738 #ifdef GDCM_DEBUG
1739       for (unsigned int i=0; i < itemLength; i += 4 )
1740       {
1741          uint32_t individualLength = str2num( &basicOffsetTableItemValue[i],
1742                                               uint32_t);
1743          gdcmWarningMacro( "Read one length: " <<
1744                           std::hex << individualLength );
1745       }
1746 #endif //GDCM_DEBUG
1747
1748       delete[] basicOffsetTableItemValue;
1749    }
1750 }
1751
1752 // These are the deprecated method that one day should be removed (after the next release)
1753
1754 #ifndef GDCM_LEGACY_REMOVE
1755 /**
1756  * \brief  Constructor (DEPRECATED : temporaryly kept not to break the API)
1757  * @param  filename name of the file whose header we want to analyze
1758  * @deprecated do not use any longer
1759  */
1760 File::File( std::string const &filename )
1761      :Document( )
1762 {
1763    RLEInfo  = new RLEFramesInfo;
1764    JPEGInfo = new JPEGFragmentsInfo;
1765
1766    Load( filename ); // gdcm::Document is first Loaded, then the 'File part'
1767 }
1768
1769 /**
1770  * \brief   Loader. (DEPRECATED :  temporaryly kept not to break the API)
1771  * @param   fileName file to be open for parsing
1772  * @return false if file cannot be open or no swap info was found,
1773  *         or no tag was found.
1774  * @deprecated Use the Load() [ + SetLoadMode() ] + SetFileName() functions instead
1775  */
1776 bool File::Load( std::string const &fileName )
1777 {
1778    GDCM_LEGACY_REPLACED_BODY(File::Load(std::string), "1.2",
1779                              File::Load());
1780    SetFileName( fileName );
1781    if ( ! this->Document::Load( ) )
1782       return false;
1783
1784    return DoTheLoadingJob( );
1785 }
1786 #endif
1787
1788 // -----------------------------------------------------------------------------------------
1789 //             THERALYS Algorithm to determine the most similar basic orientation
1790 //
1791 //  Transliterated from original Python code.
1792 //  Kept as close as possible to the original code
1793 //  in order to speed up any further modif of Python code :-(
1794 // ------------------------------------------------------------------------------------------
1795
1796 /**
1797  * \brief  THERALYS' Algorithm to determine the most similar basic orientation
1798  *           (Axial, Coronal, Sagital) of the image
1799  * \note Should be run on the first gdcm::File of a 'coherent' Serie
1800  * @return orientation code
1801  * @return orientation code
1802  *   #   0 :   Not Applicable (neither 0020,0037 Image Orientation Patient
1803  *   #                         nor     0020,0032Image Position     found )
1804  *   #   1 :   Axial
1805  *   #  -1 :   Axial invert
1806  *   #   2 :   Coronal
1807  *   #  -2 :   Coronal invert
1808  *   #   3 :   Sagital
1809  *   #  -3 :   Sagital invert
1810  *   #   4 :   Heart Axial
1811  *   #  -4 :   Heart Axial invert
1812  *   #   5 :   Heart Coronal
1813  *   #  -5 :   Heart Coronal invert
1814  *   #   6 :   Heart Sagital
1815  *   #  -6 :   Heart Sagital invert
1816  */
1817 float File::TypeOrientation( )
1818 {
1819    float *iop = new float[6];
1820    bool succ = GetImageOrientationPatient( iop );
1821    if ( !succ )
1822    {
1823       delete iop;
1824       return 0.;
1825    }
1826
1827    vector3D ori1;
1828    vector3D ori2;
1829
1830    ori1.x = iop[0]; ori1.y = iop[1]; ori1.z = iop[2];
1831    ori1.x = iop[3]; ori2.y = iop[4]; ori2.z = iop[5];
1832
1833    // two perpendicular vectors describe one plane
1834    float dicPlane[6][2][3] =
1835    { {  {1,    0,    0   },{0,       1,     0     }  },       // Axial
1836      {  {1,    0,    0   },{0,       0,    -1     }  },       // Coronal
1837      {  {0,    1,    0   },{0,       0,    -1     }  },       // Sagittal
1838      {  { 0.8, 0.5,  0.0 },{-0.1,    0.1 , -0.95  }  },       // Axial - HEART
1839      {  { 0.8, 0.5,  0.0 },{-0.6674, 0.687, 0.1794}  },       // Coronal - HEART
1840      {  {-0.1, 0.1, -0.95},{-0.6674, 0.687, 0.1794}  }        // Sagittal - HEART
1841    };
1842
1843    vector3D refA;
1844    vector3D refB;
1845    int i = 0;
1846    Res res;   // [ <result> , <memory of the last succes calcule> ]
1847    res.first = 0;
1848    res.second = 99999;
1849    for (int numDicPlane=0; numDicPlane<6; numDicPlane++)
1850    {
1851        i = i + 1;
1852        // refA=plane[0]
1853        refA.x = dicPlane[numDicPlane][0][0];
1854        refA.y = dicPlane[numDicPlane][0][1];
1855        refA.z = dicPlane[numDicPlane][0][2];
1856        // refB=plane[1]
1857        refB.x = dicPlane[numDicPlane][1][0];
1858        refB.y = dicPlane[numDicPlane][1][1];
1859        refB.z = dicPlane[numDicPlane][1][2];
1860        res=VerfCriterion(  i, CalculLikelyhood2Vec(refA,refB,ori1,ori2), res );
1861        res=VerfCriterion( -i, CalculLikelyhood2Vec(refB,refA,ori1,ori2), res );
1862    }
1863    delete iop;
1864    return res.first;
1865 /*
1866 //             i=0
1867 //             res=[0,99999]  ## [ <result> , <memory of the last succes calculus> ]
1868 //             for plane in dicPlane:
1869 //                 i=i+1
1870 //                 refA=plane[0]
1871 //                 refB=plane[1]
1872 //                 res=self.VerfCriterion(  i , self.CalculLikelyhood2Vec(refA,refB,ori1,ori2) , res )
1873 //                 res=self.VerfCriterion( -i , self.CalculLikelyhood2Vec(refB,refA,ori1,ori2) , res )
1874 //             return res[0]
1875 */
1876
1877 }
1878
1879 Res File::VerfCriterion(int typeCriterion, float criterionNew, Res res)
1880 {
1881    float type      = res.first;
1882    float criterion = res.second;
1883    if (criterionNew < criterion)
1884    {
1885       res.first  = criterionNew;
1886       res.second = typeCriterion;
1887    }
1888 /*
1889 //   type = res[0]
1890 //   criterion = res[1]
1891 // #     if criterionNew<0.1 and criterionNew<criterion:
1892 //   if criterionNew<criterion:
1893 //      criterion=criterionNew
1894 //      type=typeCriterion
1895 //   return [ type , criterion ]
1896 */
1897    return res;
1898 }
1899
1900 float File::CalculLikelyhood2Vec(vector3D refA, vector3D refB,
1901                                  vector3D ori1, vector3D ori2)
1902 {
1903 //   # ------------------------- Purpose : -----------------------------------
1904 //   # - This function determines the orientation similarity of two planes.
1905 //   #   Each plane is described by two vectors.
1906 //   # ------------------------- Parameters : --------------------------------
1907 //   # - <refA>  : - type : vector 3D (float)
1908 //   # - <refB>  : - type : vector 3D (float)
1909 //   #             - Description of the first plane
1910 //   # - <ori1>  : - type : vector 3D (float)
1911 //   # - <ori2>  : - type : vector 3D (float)
1912 //   #             - Description of the second plane
1913 //   # ------------------------- Return : ------------------------------------
1914 //   #  float :   0 if the planes are perpendicular. While the difference of
1915 //   #             the orientation between the planes are big more enlarge is
1916 //   #             the criterion.
1917 //   # ------------------------- Other : -------------------------------------
1918 //   #  The calculus is based with vectors normalice
1919
1920    vector3D ori3 = ProductVectorial(ori1,ori2);
1921    vector3D refC = ProductVectorial(refA,refB);
1922    float    res  = powf(refC.x-ori3.x, 2.) +
1923                    powf(refC.y-ori3.y, 2.) +
1924                    powf(refC.z-ori3.z, 2.);
1925
1926 /*
1927 //   ori3=self.ProductVectorial(ori1,ori2)
1928 //   refC=self.ProductVectorial(refA,refB)
1929 //   res=math.pow(refC[0]-ori3[0],2) + math.pow(refC[1]-ori3[1],2) + math.pow(refC[2]-ori3[2],2)
1930 //   return math.sqrt(res)
1931 */
1932    return sqrt(res);
1933 }
1934
1935 vector3D File::ProductVectorial(vector3D vec1, vector3D vec2)
1936 {
1937
1938 //   # ------------------------- Purpose : -----------------------------------
1939 //   # - Calculus of the poduct vectorial between two vectors 3D
1940 //   # ------------------------- Parameters : --------------------------------
1941 //   # - <vec1>  : - type : vector 3D (float)
1942 //   # - <vec2>  : - type : vector 3D (float)
1943 //   # ------------------------- Return : ------------------------------------
1944 //   #  (vec) :    - Vector 3D
1945 //   # ------------------------- Other : -------------------------------------
1946
1947    vector3D vec3;
1948    vec3.x =    vec1.y*vec3.z - vec1.z*vec2.y;
1949    vec3.y = -( vec1.x*vec2.z - vec1.z*vec2.x);
1950    vec3.z =    vec1.x*vec2.y - vec1.y*vec2.x;
1951 /*
1952 //   vec3=[0,0,0]
1953 //   vec3[0]=vec1[1]*vec2[2]    - vec1[2]*vec2[1]
1954 //   vec3[1]=-( vec1[0]*vec2[2] - vec1[2]*vec2[0])
1955 //   vec3[2]=vec1[0]*vec2[1]    - vec1[1]*vec2[0]
1956 */
1957    return vec3;
1958 }
1959
1960 //-----------------------------------------------------------------------------
1961 // Print
1962
1963 //-----------------------------------------------------------------------------
1964 } // end namespace gdcm