src/gdcmFile.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: gdcmFile.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2004/09/13 07:49:36 $
   7   Version:   $Revision: 1.125 $
   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.htm for details.
  12
  13      This software is distributed WITHOUT ANY WARRANTY; without even
  14      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  15      PURPOSE.  See the above copyright notices for more information.
  16
  17 =========================================================================*/
  18
  19 #include "gdcmFile.h"
  20 #include "gdcmDebug.h"
  21 #include "jpeg/ljpg/jpegless.h"
  22
  23 typedef std::pair<TagDocEntryHT::iterator,TagDocEntryHT::iterator> IterHT;
  24
  25 //-------------------------------------------------------------------------
  26 // Constructor / Destructor
  27 /**
  28  * \ingroup   gdcmFile
  29  * \brief Constructor dedicated to deal with the *pixels* area of a ACR/DICOMV3
  30  *        file (gdcmHeader only deals with the ... header)
  31  *        Opens (in read only and when possible) an existing file and checks
  32  *        for DICOM compliance. Returns NULL on failure.
  33  *        It will be up to the user to load the pixels into memory
  34  *        (see GetImageData, GetImageDataRaw)
  35  * \note  the in-memory representation of all available tags found in
  36  *        the DICOM header is post-poned to first header information access.
  37  *        This avoid a double parsing of public part of the header when
  38  *        user sets an a posteriori shadow dictionary (efficiency can be
  39  *        seen as a side effect).
  40  * @param header already built gdcmHeader
  41  */
  42 gdcmFile::gdcmFile(gdcmHeader *header)
  43 {
  44    Header     = header;
  45    SelfHeader = false;
  46    SetInitialValues();
  47 }
  48
  49 /**
  50  * \ingroup   gdcmFile
  51  * \brief Constructor dedicated to deal with the *pixels* area of a ACR/DICOMV3
  52  *        file (gdcmHeader only deals with the ... header)
  53  *        Opens (in read only and when possible) an existing file and checks
  54  *        for DICOM compliance. Returns NULL on failure.
  55  *        It will be up to the user to load the pixels into memory
  56  *        (see GetImageData, GetImageDataRaw)
  57  * \note  the in-memory representation of all available tags found in
  58  *        the DICOM header is post-poned to first header information access.
  59  *        This avoid a double parsing of public part of the header when
  60  *        one sets an a posteriori shadow dictionary (efficiency can be
  61  *        seen as a side effect).
  62  * @param filename file to be opened for parsing
  63  */
  64 gdcmFile::gdcmFile(std::string const & filename )
  65 {
  66    Header = new gdcmHeader( filename );
  67    SelfHeader = true;
  68    SetInitialValues();
  69 }
  70
  71 /**
  72  * \ingroup   gdcmFile
  73  * \brief Sets some initial for the Constructor
  74  */
  75 void gdcmFile::SetInitialValues()
  76 {
  77    PixelRead  = -1; // no ImageData read yet.
  78    LastAllocatedPixelDataLength = 0;
  79    Pixel_Data = 0;
  80
  81    InitialSpp = "";
  82    InitialPhotInt = "";
  83    InitialPlanConfig = "";
  84    InitialBitsAllocated = "";
  85
  86    InitialRedLUTDescr   = 0;
  87    InitialGreenLUTDescr = 0;
  88    InitialBlueLUTDescr  = 0;
  89    InitialRedLUTData    = 0;
  90    InitialGreenLUTData  = 0;
  91    InitialBlueLUTData   = 0;
  92
  93    if ( Header->IsReadable() )
  94    {
  95       SetPixelDataSizeFromHeader();
  96
  97       // the following values *may* be modified
  98       // by gdcmFile::GetImageDataIntoVectorRaw
  99       // we save their initial value.
 100       InitialSpp           = Header->GetEntryByNumber(0x0028,0x0002);
 101       InitialPhotInt       = Header->GetEntryByNumber(0x0028,0x0004);
 102       InitialPlanConfig    = Header->GetEntryByNumber(0x0028,0x0006);
 103       InitialBitsAllocated = Header->GetEntryByNumber(0x0028,0x0100);
 104
 105       // the following entries *may* be removed from the H table
 106       // (NOT deleted ...) by gdcmFile::GetImageDataIntoVectorRaw
 107       // we keep a pointer on them.
 108       InitialRedLUTDescr   = Header->GetDocEntryByNumber(0x0028,0x1101);
 109       InitialGreenLUTDescr = Header->GetDocEntryByNumber(0x0028,0x1102);
 110       InitialBlueLUTDescr  = Header->GetDocEntryByNumber(0x0028,0x1103);
 111       InitialRedLUTData    = Header->GetDocEntryByNumber(0x0028,0x1201);
 112       InitialGreenLUTData  = Header->GetDocEntryByNumber(0x0028,0x1202);
 113       InitialBlueLUTData   = Header->GetDocEntryByNumber(0x0028,0x1203);
 114    }
 115 }
 116
 117 /**
 118  * \ingroup   gdcmFile
 119  * \brief canonical destructor
 120  * \note  If the gdcmHeader was created by the gdcmFile constructor,
 121  *        it is destroyed by the gdcmFile
 122  */
 123 gdcmFile::~gdcmFile()
 124 {
 125    if( SelfHeader )
 126    {
 127       delete Header;
 128    }
 129    Header = 0;
 130
 131 // InitialLutDescriptors and InitialLutData
 132 // will have to be deleted if the don't belong any longer
 133 // to the Header H table ...
 134 }
 135
 136 //-----------------------------------------------------------------------------
 137 // Print
 138
 139 //-----------------------------------------------------------------------------
 140 // Public
 141
 142 /**
 143  * \ingroup   gdcmFile
 144  * \brief     computes the length (in bytes) we must ALLOCATE to receive the
 145  *            image(s) pixels (multiframes taken into account)
 146  * \warning : it is NOT the group 7FE0 length
 147  *          (no interest for compressed images).
 148  * \warning : not end user intended ?
 149  */
 150 void gdcmFile::SetPixelDataSizeFromHeader()
 151 {
 152    // see PS 3.3-2003 : C.7.6.3.2.1
 153    //
 154    //   MONOCHROME1
 155    //   MONOCHROME2
 156    //   PALETTE COLOR
 157    //   RGB
 158    //   HSV  (Retired)
 159    //   ARGB (Retired)
 160    //   CMYK (Retired)
 161    //   YBR_FULL
 162    //   YBR_FULL_422 (no LUT, no Palette)
 163    //   YBR_PARTIAL_422
 164    //   YBR_ICT
 165    //   YBR_RCT
 166
 167    // LUT's
 168    // ex : gdcm-US-ALOKA-16.dcm
 169    // 0028|1221 [OW]   [Segmented Red Palette Color Lookup Table Data]
 170    // 0028|1222 [OW]   [Segmented Green Palette Color Lookup Table Data]
 171    // 0028|1223 [OW]   [Segmented Blue Palette Color Lookup Table Data]
 172
 173    // ex  : OT-PAL-8-face.dcm
 174    // 0028|1201 [US]   [Red Palette Color Lookup Table Data]
 175    // 0028|1202 [US]   [Green Palette Color Lookup Table Data]
 176    // 0028|1203 [US]   [Blue Palette Color Lookup Table Data]
 177
 178    // Number of "Bits Allocated"
 179    int nb;
 180    std::string str_nb = Header->GetEntryByNumber(0x0028,0x0100);
 181
 182    if ( str_nb == GDCM_UNFOUND )
 183    {
 184       nb = 16;
 185    }
 186    else
 187    {
 188       nb = atoi( str_nb.c_str() );
 189       if (nb == 12)
 190       {
 191          nb =16;
 192       }
 193    }
 194    ImageDataSize =
 195    ImageDataSizeRaw = Header->GetXSize() * Header->GetYSize()
 196                 * Header->GetZSize() * (nb/8) * Header->GetSamplesPerPixel();
 197    std::string str_PhotometricInterpretation =
 198                              Header->GetEntryByNumber(0x0028,0x0004);
 199
 200    // if ( str_PhotometricInterpretation == "PALETTE COLOR " ),
 201    // pb when undealt Segmented Palette Color
 202
 203    if ( Header->HasLUT() )
 204    {
 205       ImageDataSize *= 3;
 206    }
 207 }
 208
 209 /**
 210  * \ingroup gdcmFile
 211  * \brief   - Allocates necessary memory,
 212  *          - Reads the pixels from disk (uncompress if necessary),
 213  *          - Transforms YBR pixels, if any, into RGB pixels
 214  *          - Transforms 3 planes R, G, B, if any, into a single RGB Plane
 215  *          - Transforms single Grey plane + 3 Palettes into a RGB Plane
 216  *          - Copies the pixel data (image[s]/volume[s]) to newly allocated zone.
 217  * @return  Pointer to newly allocated pixel data.
 218  *          NULL if alloc fails
 219  */
 220 void *gdcmFile::GetImageData()
 221 {
 222    // FIXME (Mathieu)
 223    // I need to deallocate Pixel_Data before doing any allocation:
 224
 225    if ( Pixel_Data )
 226      if ( LastAllocatedPixelDataLength != ImageDataSize )
 227         free(Pixel_Data);
 228    if ( !Pixel_Data )
 229       Pixel_Data = new uint8_t[ImageDataSize];
 230
 231    if ( Pixel_Data )
 232    {
 233       LastAllocatedPixelDataLength = ImageDataSize;
 234
 235       // we load the pixels (and transform grey level + LUT into RGB)
 236       GetImageDataIntoVector(Pixel_Data, ImageDataSize);
 237
 238       // We say the value *is* loaded.
 239       GetHeader()->SetEntryByNumber( GDCM_BINLOADED,
 240          GetHeader()->GetGrPixel(), GetHeader()->GetNumPixel());
 241
 242       // Will be 7fe0, 0010 in standard case
 243       GetHeader()->SetEntryVoidAreaByNumber( Pixel_Data,
 244          GetHeader()->GetGrPixel(), GetHeader()->GetNumPixel());
 245    }
 246    PixelRead = 0; // no PixelRaw
 247
 248    return Pixel_Data;
 249 }
 250
 251 /**
 252  * \ingroup gdcmFile
 253  * \brief
 254  *          Read the pixels from disk (uncompress if necessary),
 255  *          Transforms YBR pixels, if any, into RGB pixels
 256  *          Transforms 3 planes R, G, B, if any, into a single RGB Plane
 257  *          Transforms single Grey plane + 3 Palettes into a RGB Plane
 258  *          Copies at most MaxSize bytes of pixel data to caller allocated
 259  *          memory space.
 260  * \warning This function allows people that want to build a volume
 261  *          from an image stack *not to* have, first to get the image pixels,
 262  *          and then move them to the volume area.
 263  *          It's absolutely useless for any VTK user since vtk chooses
 264  *          to invert the lines of an image, that is the last line comes first
 265  *          (for some axis related reasons?). Hence he will have
 266  *          to load the image line by line, starting from the end.
 267  *          VTK users have to call GetImageData
 268  *
 269  * @param   destination Address (in caller's memory space) at which the
 270  *          pixel data should be copied
 271  * @param   maxSize Maximum number of bytes to be copied. When MaxSize
 272  *          is not sufficient to hold the pixel data the copy is not
 273  *          executed (i.e. no partial copy).
 274  * @return  On success, the number of bytes actually copied. Zero on
 275  *          failure e.g. MaxSize is lower than necessary.
 276  */
 277 size_t gdcmFile::GetImageDataIntoVector (void* destination, size_t maxSize)
 278 {
 279    GetImageDataIntoVectorRaw (destination, maxSize);
 280    PixelRead = 0 ; // =0 : no ImageDataRaw
 281    if ( !Header->HasLUT() )
 282    {
 283       return ImageDataSize;
 284    }
 285
 286    // from Lut R + Lut G + Lut B
 287    uint8_t *newDest = new uint8_t[ImageDataSize];
 288    uint8_t *a       = (uint8_t *)destination;
 289    uint8_t *lutRGBA = Header->GetLUTRGBA();
 290
 291    if ( lutRGBA )
 292    {
 293       int j;
 294       // move Gray pixels to temp area
 295       memmove(newDest, destination, ImageDataSizeRaw);
 296       for (size_t i=0; i<ImageDataSizeRaw; ++i)
 297       {
 298          // Build RGB Pixels
 299          j    = newDest[i]*4;
 300          *a++ = lutRGBA[j];
 301          *a++ = lutRGBA[j+1];
 302          *a++ = lutRGBA[j+2];
 303       }
 304       delete[] newDest;
 305
 306    // now, it's an RGB image
 307    // Lets's write it in the Header
 308
 309    // FIXME : Better use CreateOrReplaceIfExist ?
 310
 311    std::string spp = "3";        // Samples Per Pixel
 312    Header->SetEntryByNumber(spp,0x0028,0x0002);
 313    std::string rgb = "RGB ";     // Photometric Interpretation
 314    Header->SetEntryByNumber(rgb,0x0028,0x0004);
 315    std::string planConfig = "0"; // Planar Configuration
 316    Header->SetEntryByNumber(planConfig,0x0028,0x0006);
 317
 318    }
 319    else  // GetLUTRGBA() failed
 320    {
 321       // (gdcm-US-ALOKA-16.dcm), contains Segmented xxx Palette Color
 322       // that are *more* than 65535 long ?!?
 323       // No idea how to manage such an image !
 324       // Need to make RGB Pixels (?) from grey Pixels (?!) and Gray Lut  (!?!)
 325       // It seems that *no Dicom Viewer* has any idea :-(
 326
 327       std::string photomInterp = "MONOCHROME1 ";  // Photometric Interpretation
 328       Header->SetEntryByNumber(photomInterp,0x0028,0x0004);
 329    }
 330
 331    /// \todo Drop Palette Color out of the Header?
 332    return ImageDataSize;
 333 }
 334
 335 /**
 336  * \ingroup gdcmFile
 337  * \brief   Allocates necessary memory,
 338  *          Transforms YBR pixels (if any) into RGB pixels
 339  *          Transforms 3 planes R, G, B  (if any) into a single RGB Plane
 340  *          Copies the pixel data (image[s]/volume[s]) to newly allocated zone.
 341  *          DOES NOT transform Grey plane + 3 Palettes into a RGB Plane
 342  * @return  Pointer to newly allocated pixel data.
 343  * \        NULL if alloc fails
 344  */
 345 void * gdcmFile::GetImageDataRaw ()
 346 {
 347    size_t imgDataSize;
 348    if ( Header->HasLUT() )
 349       /// \todo Let gdcmHeader user a chance to get the right value
 350       imgDataSize = ImageDataSizeRaw;
 351    else
 352       imgDataSize = ImageDataSize;
 353
 354    // FIXME (Mathieu)
 355    // I need to deallocate Pixel_Data before doing any allocation:
 356
 357    if ( Pixel_Data )
 358       if ( LastAllocatedPixelDataLength != imgDataSize )
 359          free(Pixel_Data);
 360    if ( !Pixel_Data )
 361       Pixel_Data = new uint8_t[imgDataSize];
 362
 363    if ( Pixel_Data )
 364    {
 365       LastAllocatedPixelDataLength = imgDataSize;
 366
 367       // we load the pixels ( grey level or RGB, but NO transformation)
 368        GetImageDataIntoVectorRaw(Pixel_Data, imgDataSize);
 369
 370       // We say the value *is* loaded.
 371       GetHeader()->SetEntryByNumber( GDCM_BINLOADED,
 372          GetHeader()->GetGrPixel(), GetHeader()->GetNumPixel());
 373
 374       // will be 7fe0, 0010 in standard cases
 375       GetHeader()->SetEntryVoidAreaByNumber(Pixel_Data,
 376          GetHeader()->GetGrPixel(), GetHeader()->GetNumPixel());
 377    }
 378    PixelRead = 1; // PixelRaw
 379
 380    return Pixel_Data;
 381 }
 382
 383 /**
 384  * \ingroup gdcmFile
 385  * \brief   Copies at most MaxSize bytes of pixel data to caller's
 386  *          memory space.
 387  * \warning This function was designed to avoid people that want to build
 388  *          a volume from an image stack to need first to get the image pixels
 389  *          and then move them to the volume area.
 390  *          It's absolutely useless for any VTK user since vtk chooses
 391  *          to invert the lines of an image, that is the last line comes first
 392  *          (for some axis related reasons?). Hence he will have
 393  *          to load the image line by line, starting from the end.
 394  *          VTK users hace to call GetImageData
 395  * \warning DOES NOT transform the Grey Plane + Palette Color (if any)
 396  *                   into a single RGB Pixels Plane
 397  *          the (VTK) user will manage the palettes
 398  *
 399  * @param   destination Address (in caller's memory space) at which the
 400  *          pixel data should be copied
 401  * @param   maxSize Maximum number of bytes to be copied. When MaxSize
 402  *          is not sufficient to hold the pixel data the copy is not
 403  *          executed (i.e. no partial copy).
 404  * @return  On success, the number of bytes actually copied. Zero on
 405  *          failure e.g. MaxSize is lower than necessary.
 406  */
 407 size_t gdcmFile::GetImageDataIntoVectorRaw (void *destination, size_t maxSize)
 408 {
 409    int nb, nbu, highBit, sign;
 410
 411   // we save the initial values of the following
 412   // in order to be able to restore the header in a disk-consistent state
 413   // (if user asks twice to get the pixels from disk)
 414
 415    if ( PixelRead == -1 ) // File was never "read" before
 416    {
 417       InitialSpp           = Header->GetEntryByNumber(0x0028,0x0002);
 418       InitialPhotInt       = Header->GetEntryByNumber(0x0028,0x0004);
 419       InitialPlanConfig    = Header->GetEntryByNumber(0x0028,0x0006);
 420       InitialBitsAllocated = Header->GetEntryByNumber(0x0028,0x0100);
 421    }
 422    else // File was already "read", the following *may* have been modified
 423         // we restore them to be in a disk-consistent state
 424    {
 425        // FIXME : What happened with the LUTs ?
 426        Header->SetEntryByNumber(InitialSpp,0x0028,0x0002);
 427        Header->SetEntryByNumber(InitialPhotInt,0x0028,0x0004);
 428        Header->SetEntryByNumber(InitialPlanConfig,0x0028,0x0006);
 429        Header->SetEntryByNumber(InitialBitsAllocated,0x0028,0x0100);
 430    }
 431
 432    PixelRead = 1 ; // PixelRaw
 433
 434    if ( ImageDataSize > maxSize )
 435    {
 436       dbg.Verbose(0, "gdcmFile::GetImageDataIntoVector: pixel data bigger"
 437                      "than caller's expected MaxSize");
 438       return (size_t)0;
 439    }
 440
 441    ReadPixelData( destination );
 442
 443    // Number of Bits Allocated for storing a Pixel
 444    std::string str_nb = Header->GetEntryByNumber(0x0028,0x0100);
 445    if ( str_nb == GDCM_UNFOUND )
 446    {
 447       nb = 16;
 448    }
 449    else
 450    {
 451       nb = atoi( str_nb.c_str() );
 452       // FIXME
 453       // From reading SetPixelDataSizeFromHeader, it seems 12 should be treated
 454       // separately, correct ?
 455    }
 456
 457    // Number of Bits actually used
 458    std::string str_nbu = Header->GetEntryByNumber(0x0028,0x0101);
 459    if ( str_nbu == GDCM_UNFOUND )
 460    {
 461       nbu = nb;
 462    }
 463    else
 464    {
 465       nbu = atoi( str_nbu.c_str() );
 466    }
 467
 468    // High Bit Position
 469    std::string str_highBit = Header->GetEntryByNumber(0x0028,0x0102);
 470    if ( str_highBit == GDCM_UNFOUND )
 471    {
 472       highBit = nb - 1;
 473    }
 474    else
 475    {
 476       highBit = atoi( str_highBit.c_str() );
 477    }
 478
 479    // Pixel sign
 480    // 0 = Unsigned
 481    // 1 = Signed
 482    std::string str_sign = Header->GetEntryByNumber(0x0028,0x0103);
 483    if ( str_sign == GDCM_UNFOUND )
 484    {
 485       sign = 0;  // default is unsigned
 486    }
 487    else
 488    {
 489       sign = atoi( str_sign.c_str() );
 490    }
 491
 492    // re arange bytes inside the integer (processor endianity)
 493    if ( nb != 8 )
 494    {
 495       SwapZone(destination, Header->GetSwapCode(), ImageDataSize, nb);
 496    }
 497
 498    // to avoid pb with some xmedcon breakers images
 499    if ( nb == 16 && nbu < nb && sign == 0)
 500    {
 501       int l = (int)(ImageDataSize / (nb/8));
 502       uint16_t *deb = (uint16_t *)destination;
 503       for(int i = 0; i<l; i++)
 504       {
 505          if( *deb == 0xffff )
 506          {
 507            *deb = 0;
 508          }
 509          deb++;
 510       }
 511    }
 512
 513    // re arange bits inside the bytes
 514    if ( nbu != nb )
 515    {
 516       int l = (int)(ImageDataSize / (nb/8));
 517       if ( nb == 16 )
 518       {
 519          uint16_t mask = 0xffff;
 520          mask = mask >> (nb-nbu);
 521          uint16_t *deb = (uint16_t *)destination;
 522          for(int i = 0; i<l; i++)
 523          {
 524             *deb = (*deb >> (nbu - highBit - 1)) & mask;
 525             deb++;
 526          }
 527       }
 528       else if ( nb == 32 )
 529       {
 530          uint32_t mask = 0xffffffff;
 531          mask         = mask >> (nb - nbu);
 532          uint32_t *deb = (uint32_t *)destination;
 533          for(int i = 0; i<l; i++)
 534          {
 535             *deb = (*deb >> (nbu - highBit - 1)) & mask;
 536             deb++;
 537          }
 538       }
 539       else
 540       {
 541          dbg.Verbose(0, "gdcmFile::GetImageDataIntoVector: weird image");
 542          return 0;
 543       }
 544    }
 545 // DO NOT remove this code commented out.
 546 // Nobody knows what's expecting you ...
 547 // Just to 'see' what was actually read on disk :-(
 548
 549 //   FILE * f2;
 550 //   f2 = fopen("SpuriousFile.RAW","wb");
 551 //   fwrite(destination,ImageDataSize,1,f2);
 552 //   fclose(f2);
 553
 554    // Deal with the color
 555    // -------------------
 556
 557    std::string str_PhotometricInterpretation =
 558       Header->GetEntryByNumber(0x0028,0x0004);
 559
 560    if ( str_PhotometricInterpretation == "MONOCHROME1 "
 561      || str_PhotometricInterpretation == "MONOCHROME2 " )
 562    {
 563       return ImageDataSize;
 564    }
 565
 566    // Planar configuration = 0 : Pixels are already RGB
 567    // Planar configuration = 1 : 3 planes : R, G, B
 568    // Planar configuration = 2 : 1 gray Plane + 3 LUT
 569
 570    // Well ... supposed to be !
 571    // See US-PAL-8-10x-echo.dcm: PlanarConfiguration=0,
 572    //                            PhotometricInterpretation=PALETTE COLOR
 573    // and heuristic has to be found :-(
 574
 575    int planConf = Header->GetPlanarConfiguration();  // 0028,0006
 576
 577    // Whatever Planar Configuration is,
 578    // "PALETTE COLOR " implies that we deal with the palette.
 579    if ( str_PhotometricInterpretation == "PALETTE COLOR ")
 580    {
 581       planConf = 2;
 582    }
 583
 584    switch ( planConf )
 585    {
 586       case 0:
 587          // Pixels are already RGB
 588          break;
 589       case 1:
 590          if (str_PhotometricInterpretation == "YBR_FULL")
 591          {
 592             // Warning : YBR_FULL_422 acts as RGB
 593             //         : we need to make RGB Pixels from Planes Y,cB,cR
 594
 595             // to see the tricks about YBR_FULL, YBR_FULL_422,
 596             // YBR_PARTIAL_422, YBR_ICT, YBR_RCT have a look at :
 597             // ftp://medical.nema.org/medical/dicom/final/sup61_ft.pdf
 598             // and be *very* affraid
 599             //
 600             int l        = Header->GetXSize() * Header->GetYSize();
 601             int nbFrames = Header->GetZSize();
 602
 603             uint8_t* newDest = new uint8_t[ImageDataSize];
 604             uint8_t* x       = newDest;
 605             uint8_t* a       = (uint8_t*)destination;
 606             uint8_t* b       = a + l;
 607             uint8_t* c       = b + l;
 608             double R,G,B;
 609
 610             /// \todo : Replace by the 'well known' integer computation
 611             /// counterpart
 612             /// see http://lestourtereaux.free.fr/papers/data/yuvrgb.pdf
 613             /// for code optimisation
 614
 615             for (int i = 0; i < nbFrames; i++)
 616             {
 617                for (int j = 0; j < l; j++)
 618                {
 619                   R = 1.164 *(*a-16) + 1.596 *(*c -128) + 0.5;
 620                   G = 1.164 *(*a-16) - 0.813 *(*c -128) - 0.392 *(*b -128) + 0.5;
 621                   B = 1.164 *(*a-16) + 2.017 *(*b -128) + 0.5;
 622
 623                   if (R < 0.0)   R = 0.0;
 624                   if (G < 0.0)   G = 0.0;
 625                   if (B < 0.0)   B = 0.0;
 626                   if (R > 255.0) R = 255.0;
 627                   if (G > 255.0) G = 255.0;
 628                   if (B > 255.0) B = 255.0;
 629
 630                   *(x++) = (uint8_t)R;
 631                   *(x++) = (uint8_t)G;
 632                   *(x++) = (uint8_t)B;
 633                   a++; b++; c++;
 634                }
 635             }
 636             memmove(destination, newDest, ImageDataSize);
 637             delete[] newDest;
 638          }
 639          else
 640          {
 641             // need to make RGB Pixels from R,G,B Planes
 642             // (all the Frames at a time)
 643
 644             int l = Header->GetXSize() * Header->GetYSize() * Header->GetZSize();
 645
 646             uint8_t *newDest = new uint8_t[ImageDataSize];
 647             uint8_t *x       = newDest;
 648             uint8_t *a       = (uint8_t *)destination;
 649             uint8_t *b       = a + l;
 650             uint8_t *c       = b + l;
 651
 652             for (int j = 0; j < l; j++)
 653             {
 654                *(x++) = *(a++);
 655                *(x++) = *(b++);
 656                *(x++) = *(c++);
 657             }
 658             memmove(destination, newDest, ImageDataSize);
 659             delete[] newDest;
 660          }
 661          break;
 662       case 2:
 663          // Palettes were found
 664          // Let the user deal with them !
 665          return ImageDataSize;
 666    }
 667    // now, it's an RGB image
 668    // Lets's write it in the Header
 669
 670    // CreateOrReplaceIfExist ?
 671
 672    std::string spp = "3";            // Samples Per Pixel
 673    std::string photInt = "RGB ";     // Photometric Interpretation
 674    std::string planConfig = "0";     // Planar Configuration
 675
 676
 677
 678    Header->SetEntryByNumber(spp,0x0028,0x0002);
 679    Header->SetEntryByNumber(photInt,0x0028,0x0004);
 680    Header->SetEntryByNumber(planConfig,0x0028,0x0006);
 681
 682    /// \todo Drop Palette Color out of the Header?
 683    return ImageDataSize;
 684 }
 685
 686 /**
 687  * \ingroup   gdcmFile
 688  * \brief performs a shalow copy (not a deep copy) of the user given
 689  *        pixel area.
 690  *        'image' Pixels are presented as C-like 2D arrays : line per line.
 691  *        'volume'Pixels are presented as C-like 3D arrays : plane per plane
 692  * \warning user is kindly requested NOT TO 'free' the Pixel area
 693  * @param inData user supplied pixel area
 694  * @param expectedSize total image size, in Bytes
 695  *
 696  * @return boolean
 697  */
 698 bool gdcmFile::SetImageData(void *inData, size_t expectedSize)
 699 {
 700    Header->SetImageDataSize( expectedSize );
 701 // FIXME : if already allocated, memory leak !
 702    Pixel_Data     = inData;
 703    ImageDataSize = ImageDataSizeRaw = expectedSize;
 704    PixelRead     = 1;
 705 // FIXME : 7fe0, 0010 IS NOT set ...
 706    return true;
 707 }
 708
 709 /**
 710  * \ingroup   gdcmFile
 711  * \brief Writes on disk A SINGLE Dicom file
 712  *        NO test is performed on  processor "Endiannity".
 713  *        It's up to the user to call his Reader properly
 714  * @param fileName name of the file to be created
 715  *                 (any already existing file is over written)
 716  * @return false if write fails
 717  */
 718
 719 bool gdcmFile::WriteRawData(std::string const & fileName)
 720 {
 721    FILE *fp1 = fopen(fileName.c_str(), "wb");
 722    if (fp1 == NULL)
 723    {
 724       printf("Fail to open (write) file [%s] \n", fileName.c_str());
 725       return false;
 726    }
 727    fwrite (Pixel_Data, ImageDataSize, 1, fp1);
 728    fclose (fp1);
 729
 730    return true;
 731 }
 732
 733 /**
 734  * \ingroup   gdcmFile
 735  * \brief Writes on disk A SINGLE Dicom file,
 736  *        using the Implicit Value Representation convention
 737  *        NO test is performed on  processor "Endiannity".
 738  * @param fileName name of the file to be created
 739  *                 (any already existing file is overwritten)
 740  * @return false if write fails
 741  */
 742
 743 bool gdcmFile::WriteDcmImplVR (std::string const & fileName)
 744 {
 745    return WriteBase(fileName, gdcmImplicitVR);
 746 }
 747
 748 /**
 749  * \ingroup   gdcmFile
 750 * \brief Writes on disk A SINGLE Dicom file,
 751  *        using the Explicit Value Representation convention
 752  *        NO test is performed on  processor "Endiannity". * @param fileName name of the file to be created
 753  *                 (any already existing file is overwritten)
 754  * @return false if write fails
 755  */
 756
 757 bool gdcmFile::WriteDcmExplVR (std::string const & fileName)
 758 {
 759    return WriteBase(fileName, gdcmExplicitVR);
 760 }
 761
 762 /**
 763  * \ingroup   gdcmFile
 764  * \brief Writes on disk A SINGLE Dicom file,
 765  *        using the ACR-NEMA convention
 766  *        NO test is performed on  processor "Endiannity".
 767  *        (a l'attention des logiciels cliniques
 768  *        qui ne prennent en entrée QUE des images ACR ...
 769  * \warning if a DICOM_V3 header is supplied,
 770  *         groups < 0x0008 and shadow groups are ignored
 771  * \warning NO TEST is performed on processor "Endiannity".
 772  * @param fileName name of the file to be created
 773  *                 (any already existing file is overwritten)
 774  * @return false if write fails
 775  */
 776
 777 bool gdcmFile::WriteAcr (std::string const & fileName)
 778 {
 779    return WriteBase(fileName, gdcmACR);
 780 }
 781
 782 //-----------------------------------------------------------------------------
 783 // Protected
 784 /**
 785  * \ingroup   gdcmFile
 786  * \brief NOT a end user inteded function
 787  *        (used by WriteDcmExplVR, WriteDcmImplVR, WriteAcr, etc)
 788  * @param fileName name of the file to be created
 789  *                 (any already existing file is overwritten)
 790  * @param  type file type (ExplicitVR, ImplicitVR, ...)
 791  * @return false if write fails
 792  */
 793 bool gdcmFile::WriteBase (std::string const & fileName, FileType type)
 794 {
 795    if ( PixelRead == -1 && type != gdcmExplicitVR)
 796    {
 797       return false;
 798    }
 799
 800    FILE *fp1 = fopen(fileName.c_str(), "wb");
 801    if (fp1 == NULL)
 802    {
 803       printf("Failed to open (write) File [%s] \n", fileName.c_str());
 804       return false;
 805    }
 806
 807    if ( type == gdcmImplicitVR || type == gdcmExplicitVR )
 808    {
 809       // writing Dicom File Preamble
 810       uint8_t* filePreamble = new uint8_t[128];
 811       memset(filePreamble, 0, 128);
 812       fwrite(filePreamble, 128, 1, fp1);
 813       fwrite("DICM", 4, 1, fp1);
 814
 815       delete[] filePreamble;
 816    }
 817
 818    // --------------------------------------------------------------
 819    // Special Patch to allow gdcm to re-write ACR-LibIDO formated images
 820    //
 821    // if recognition code tells us we dealt with a LibIDO image
 822    // we reproduce on disk the switch between lineNumber and columnNumber
 823    // just before writting ...
 824
 825    /// \todo the best trick would be *change* the recognition code
 826    ///       but pb expected if user deals with, e.g. COMPLEX images
 827
 828    std::string rows, columns;
 829    if ( Header->GetFileType() == gdcmACR_LIBIDO)
 830    {
 831       rows    = Header->GetEntryByNumber(0x0028, 0x0010);
 832       columns = Header->GetEntryByNumber(0x0028, 0x0011);
 833
 834       Header->SetEntryByNumber(columns,  0x0028, 0x0010);
 835       Header->SetEntryByNumber(rows   ,  0x0028, 0x0011);
 836    }
 837    // ----------------- End of Special Patch ----------------
 838
 839    uint16_t grPixel  = Header->GetGrPixel();
 840    uint16_t numPixel = Header->GetNumPixel();;
 841
 842    gdcmDocEntry* PixelElement =
 843       GetHeader()->GetDocEntryByNumber(grPixel, numPixel);
 844
 845    if ( PixelRead == 1 )
 846    {
 847       // we read pixel 'as is' (no tranformation LUT -> RGB)
 848       PixelElement->SetLength( ImageDataSizeRaw );
 849    }
 850    else if ( PixelRead == 0 )
 851    {
 852       // we tranformed GrayLevel pixels + LUT into RGB Pixel
 853       PixelElement->SetLength( ImageDataSize );
 854    }
 855
 856    Header->Write(fp1, type);
 857
 858    // --------------------------------------------------------------
 859    // Special Patch to allow gdcm to re-write ACR-LibIDO formated images
 860    //
 861    // ...and we restore the Header to be Dicom Compliant again
 862    // just after writting
 863
 864    if ( Header->GetFileType() == gdcmACR_LIBIDO )
 865    {
 866       Header->SetEntryByNumber(rows   , 0x0028, 0x0010);
 867       Header->SetEntryByNumber(columns, 0x0028, 0x0011);
 868    }
 869    // ----------------- End of Special Patch ----------------
 870
 871    // fwrite(Pixel_Data, ImageDataSize, 1, fp1);  // should be useless, now
 872    fclose (fp1);
 873
 874    return true;
 875 }
 876
 877 //-----------------------------------------------------------------------------
 878 // Private
 879 /**
 880  * \ingroup gdcmFile
 881  * \brief   Swap the bytes, according to swap code.
 882  * \warning not end user intended
 883  * @param   im area to deal with
 884  * @param   swap swap code
 885  * @param   lgr Area Length
 886  * @param   nb Pixels Bit number
 887  */
 888 void gdcmFile::SwapZone(void *im, int swap, int lgr, int nb)
 889 {
 890    int i;
 891
 892    if( nb == 16 )
 893    {
 894       uint16_t* im16 = (uint16_t*)im;
 895       switch( swap )
 896       {
 897          case 0:
 898          case 12:
 899          case 1234:
 900             break;
 901          case 21:
 902          case 3412:
 903          case 2143:
 904          case 4321:
 905             for(i=0; i < lgr/2; i++)
 906             {
 907                im16[i]= (im16[i] >> 8) | (im16[i] << 8 );
 908             }
 909             break;
 910          default:
 911             std::cout << "SWAP value (16 bits) not allowed :i" << swap <<
 912             std::endl;
 913       }
 914    }
 915    else if( nb == 32 )
 916    {
 917       uint32_t s32;
 918       uint16_t fort, faible;
 919       uint32_t* im32 = (uint32_t*)im;
 920       switch ( swap )
 921       {
 922          case 0:
 923          case 1234:
 924             break;
 925          case 4321:
 926             for(i = 0; i < lgr/4; i++)
 927             {
 928                faible  = im32[i] & 0x0000ffff;  // 4321
 929                fort    = im32[i] >> 16;
 930                fort    = ( fort >> 8   ) | ( fort << 8 );
 931                faible  = ( faible >> 8 ) | ( faible << 8);
 932                s32     = faible;
 933                im32[i] = ( s32 << 16 ) | fort;
 934             }
 935             break;
 936          case 2143:
 937             for(i = 0; i < lgr/4; i++)
 938             {
 939                faible  = im32[i] & 0x0000ffff;   // 2143
 940                fort    = im32[i] >> 16;
 941                fort    = ( fort >> 8 ) | ( fort << 8 );
 942                faible  = ( faible >> 8) | ( faible << 8);
 943                s32     = fort;
 944                im32[i] = ( s32 << 16 ) | faible;
 945             }
 946             break;
 947          case 3412:
 948             for(i = 0; i < lgr/4; i++)
 949             {
 950                faible  = im32[i] & 0x0000ffff; // 3412
 951                fort    = im32[i] >> 16;
 952                s32     = faible;
 953                im32[i] = ( s32 << 16 ) | fort;
 954             }
 955             break;
 956          default:
 957             std::cout << "SWAP value (32 bits) not allowed : " << swap <<
 958             std::endl;
 959       }
 960    }
 961 }
 962
 963 /**
 964  * \ingroup gdcmFile
 965  * \brief   Read pixel data from disk (optionaly decompressing) into the
 966  *          caller specified memory location.
 967  * @param   destination where the pixel data should be stored.
 968  *
 969  */
 970 bool gdcmFile::ReadPixelData(void *destination)
 971 {
 972    FILE *fp = Header->OpenFile();
 973
 974    if ( !fp )
 975    {
 976       return false;
 977    }
 978    if ( fseek(fp, Header->GetPixelOffset(), SEEK_SET) == -1 )
 979    {
 980       Header->CloseFile();
 981       return false;
 982    }
 983
 984    // ----------------------  Compacted File (12 Bits Per Pixel)
 985    // unpack 12 Bits pixels into 16 Bits pixels
 986    // 2 pixels 12bit =     [0xABCDEF]
 987    // 2 pixels 16bit = [0x0ABD] + [0x0FCE]
 988
 989    if ( Header->GetBitsAllocated() == 12 )
 990    {
 991       int nbPixels = Header->GetXSize() * Header->GetYSize();
 992       uint8_t b0, b1, b2;
 993
 994       uint16_t* pdestination = (uint16_t*)destination;
 995       for(int p = 0; p < nbPixels; p += 2 )
 996       {
 997          fread(&b0,1,1,fp);
 998          fread(&b1,1,1,fp);
 999          fread(&b2,1,1,fp);
1000
1001          //Two steps is necessary to please VC++
1002          *pdestination++ =  ((b0 >> 4) << 8) + ((b0 & 0x0f) << 4) + (b1 & 0x0f);
1003          //                     A                     B                 D
1004          *pdestination++ =  ((b2 & 0x0f) << 8) + ((b1 >> 4) << 4) + (b2 >> 4);
1005          //                     F                     C                 E
1006
1007          // Troubles expected on Big-Endian processors ?
1008       }
1009
1010       Header->CloseFile();
1011       return true;
1012    }
1013
1014    // ----------------------  Uncompressed File
1015    if ( !Header->IsDicomV3()                             ||
1016         Header->IsImplicitVRLittleEndianTransferSyntax() ||
1017         Header->IsExplicitVRLittleEndianTransferSyntax() ||
1018         Header->IsExplicitVRBigEndianTransferSyntax()    ||
1019         Header->IsDeflatedExplicitVRLittleEndianTransferSyntax() )
1020    {
1021       size_t ItemRead = fread(destination, Header->GetPixelAreaLength(), 1, fp);
1022       Header->CloseFile();
1023       if ( ItemRead != 1 )
1024       {
1025          return false;
1026       }
1027       else
1028       {
1029          return true;
1030       }
1031    }
1032
1033    // ---------------------- Run Length Encoding
1034    if ( Header->IsRLELossLessTransferSyntax() )
1035    {
1036       bool res = gdcm_read_RLE_file (fp,destination);
1037       Header->CloseFile();
1038       return res;
1039    }
1040
1041    // --------------- SingleFrame/Multiframe JPEG Lossless/Lossy/2000
1042    int nb;
1043    std::string str_nb = Header->GetEntryByNumber(0x0028,0x0100);
1044    if ( str_nb == GDCM_UNFOUND )
1045    {
1046       nb = 16;
1047    }
1048    else
1049    {
1050       nb = atoi( str_nb.c_str() );
1051       if ( nb == 12 )
1052       {
1053          nb = 16;  // ?? 12 should be ACR-NEMA only
1054       }
1055    }
1056
1057    int nBytes= nb/8;
1058    int taille = Header->GetXSize() * Header->GetYSize()
1059                 * Header->GetSamplesPerPixel();
1060    long fragmentBegining; // for ftell, fseek
1061
1062    bool jpg2000     = Header->IsJPEG2000();
1063    bool jpgLossless = Header->IsJPEGLossless();
1064
1065    bool res = true;
1066    uint16_t ItemTagGr, ItemTagEl;
1067    int ln;
1068
1069    //  Position on begining of Jpeg Pixels
1070
1071    fread(&ItemTagGr,2,1,fp);  // Reading (fffe) : Item Tag Gr
1072    fread(&ItemTagEl,2,1,fp);  // Reading (e000) : Item Tag El
1073    if(Header->GetSwapCode())
1074    {
1075       ItemTagGr = Header->SwapShort(ItemTagGr);
1076       ItemTagEl = Header->SwapShort(ItemTagEl);
1077    }
1078
1079    fread(&ln,4,1,fp);
1080    if( Header->GetSwapCode() )
1081    {
1082       ln = Header->SwapLong( ln );    // Basic Offset Table Item length
1083    }
1084
1085    if ( ln != 0 )
1086    {
1087       // What is it used for ?!?
1088       uint8_t* BasicOffsetTableItemValue = new uint8_t[ln+1];
1089       fread(BasicOffsetTableItemValue,ln,1,fp);
1090       //delete[] BasicOffsetTableItemValue;
1091    }
1092
1093    // first Fragment initialisation
1094    fread(&ItemTagGr,2,1,fp);  // Reading (fffe) : Item Tag Gr
1095    fread(&ItemTagEl,2,1,fp);  // Reading (e000) : Item Tag El
1096    if( Header->GetSwapCode() )
1097    {
1098       ItemTagGr = Header->SwapShort( ItemTagGr );
1099       ItemTagEl = Header->SwapShort( ItemTagEl );
1100    }
1101
1102    // parsing fragments until Sequence Delim. Tag found
1103    while ( ItemTagGr == 0xfffe && ItemTagEl != 0xe0dd )
1104    {
1105       // --- for each Fragment
1106       fread(&ln,4,1,fp);
1107       if( Header->GetSwapCode() )
1108       {
1109          ln = Header->SwapLong(ln);    // Fragment Item length
1110       }
1111       fragmentBegining = ftell( fp );
1112
1113       if ( jpg2000 )
1114       {
1115       // JPEG 2000 :    call to ???
1116       res = gdcm_read_JPEG2000_file (fp,destination);  // Not Yet written
1117       // ------------------------------------- endif (JPEG2000)
1118       }
1119       else if (jpgLossless)
1120       {
1121          // JPEG LossLess : call to xmedcom Lossless JPEG
1122          // Reading Fragment pixels
1123          JPEGLosslessDecodeImage (fp, (uint16_t*)destination,
1124                Header->GetPixelSize() * 8 * Header->GetSamplesPerPixel(), ln);
1125          res = 1; // in order not to break the loop
1126
1127       } // ------------------------------------- endif (JPEGLossless)
1128       else
1129       {
1130          // JPEG Lossy : call to IJG 6b
1131          if ( Header->GetBitsStored() == 8)
1132          {
1133             // Reading Fragment pixels
1134             res = gdcm_read_JPEG_file (fp,destination);
1135          }
1136          else
1137          {
1138             // Reading Fragment pixels
1139             res = gdcm_read_JPEG_file12 (fp,destination);
1140          }
1141          // ------------------------------------- endif (JPEGLossy)
1142       }
1143
1144       if ( !res )
1145       {
1146          break;
1147       }
1148
1149       // location in user's memory
1150       // for next fragment (if any)
1151       destination = (uint8_t*)destination + taille * nBytes;
1152
1153       fseek(fp,fragmentBegining, SEEK_SET); // To be sure we start
1154       fseek(fp,ln,SEEK_CUR);                // at the begining of next fragment
1155
1156       ItemTagGr = ItemTagEl = 0;
1157       fread(&ItemTagGr,2,1,fp);  // Reading (fffe) : Item Tag Gr
1158       fread(&ItemTagEl,2,1,fp);  // Reading (e000) : Item Tag El
1159       if( Header->GetSwapCode() )
1160       {
1161          ItemTagGr = Header->SwapShort( ItemTagGr );
1162          ItemTagEl = Header->SwapShort( ItemTagEl );
1163       }
1164    }
1165    // endWhile parsing fragments until Sequence Delim. Tag found
1166
1167    Header->CloseFile();
1168    return res;
1169 }
1170 //-----------------------------------------------------------------------------