src/gdcmFile.cxx

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