src/gdcmFile.cxx

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