src/gdcmFile.cxx

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