src/gdcmPixelConvert.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: gdcmPixelConvert.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2004/10/28 22:21:57 $
   7   Version:   $Revision: 1.22 $
   8
   9   Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
  10   l'Image). All rights reserved. See Doc/License.txt or
  11   http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
  12
  13      This software is distributed WITHOUT ANY WARRANTY; without even
  14      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  15      PURPOSE.  See the above copyright notices for more information.
  16
  17 =========================================================================*/
  18
  19 //////////////////   TEMPORARY NOTE
  20 // look for "fixMem" and convert that to a member of this class
  21 // Removing the prefix fixMem and dealing with allocations should do the trick
  22 //
  23 // grep PIXELCONVERT everywhere and clean up !
  24
  25 #include "gdcmDebug.h"
  26 #include "gdcmPixelConvert.h"
  27 #include <fstream>
  28
  29 namespace gdcm
  30 {
  31 #define str2num(str, typeNum) *((typeNum *)(str))
  32
  33 // For JPEG 2000, body in file gdcmJpeg2000.cxx
  34 bool gdcm_read_JPEG2000_file (std::ifstream* fp, void* image_buffer);
  35
  36 // For JPEG 8 Bits, body in file gdcmJpeg8.cxx
  37 bool gdcm_read_JPEG_file8    (std::ifstream* fp, void* image_buffer);
  38
  39 // For JPEG 12 Bits, body in file gdcmJpeg12.cxx
  40 bool gdcm_read_JPEG_file12   (std::ifstream* fp, void* image_buffer);
  41
  42 // For JPEG 16 Bits, body in file gdcmJpeg16.cxx
  43 // Beware this is misleading there is no 16bits DCT algorithm, only
  44 // jpeg lossless compression exist in 16bits.
  45 bool gdcm_read_JPEG_file16   (std::ifstream* fp, void* image_buffer);
  46
  47
  48 //-----------------------------------------------------------------------------
  49 // Constructor / Destructor
  50 PixelConvert::PixelConvert()
  51 {
  52    RGB = 0;
  53    RGBSize = 0;
  54    Decompressed = 0;
  55    DecompressedSize = 0;
  56    LutRGBA = 0;
  57    LutRedData = 0;
  58    LutGreenData = 0;
  59    LutBlueData =0;
  60 }
  61
  62 void PixelConvert::Squeeze()
  63 {
  64    if ( RGB )
  65    {
  66       delete [] RGB;
  67    }
  68    if ( Decompressed )
  69    {
  70       delete [] Decompressed;
  71    }
  72    if ( LutRGBA )
  73    {
  74       delete [] LutRGBA;
  75    }
  76 }
  77
  78 PixelConvert::~PixelConvert()
  79 {
  80    Squeeze();
  81 }
  82
  83 void PixelConvert::AllocateRGB()
  84 {
  85   if ( RGB ) {
  86      delete [] RGB;
  87   }
  88   RGB = new uint8_t[ RGBSize ];
  89 }
  90
  91 void PixelConvert::AllocateDecompressed()
  92 {
  93   if ( Decompressed ) {
  94      delete [] Decompressed;
  95   }
  96   Decompressed = new uint8_t[ DecompressedSize ];
  97 }
  98
  99 /**
 100  * \brief Read from file a 12 bits per pixel image and decompress it
 101  *        into a 16 bits per pixel image.
 102  */
 103 void PixelConvert::ReadAndDecompress12BitsTo16Bits( std::ifstream* fp )
 104                throw ( FormatError )
 105 {
 106    int nbPixels = XSize * YSize;
 107    uint16_t* localDecompres = (uint16_t*)Decompressed;
 108
 109    for( int p = 0; p < nbPixels; p += 2 )
 110    {
 111       uint8_t b0, b1, b2;
 112
 113       fp->read( (char*)&b0, 1);
 114       if ( fp->fail() || fp->eof() )//Fp->gcount() == 1
 115       {
 116          throw FormatError( "PixelConvert::ReadAndDecompress12BitsTo16Bits()",
 117                                 "Unfound first block" );
 118       }
 119
 120       fp->read( (char*)&b1, 1 );
 121       if ( fp->fail() || fp->eof())//Fp->gcount() == 1
 122       {
 123          throw FormatError( "PixelConvert::ReadAndDecompress12BitsTo16Bits()",
 124                                 "Unfound second block" );
 125       }
 126
 127       fp->read( (char*)&b2, 1 );
 128       if ( fp->fail() || fp->eof())//Fp->gcount() == 1
 129       {
 130          throw FormatError( "PixelConvert::ReadAndDecompress12BitsTo16Bits()",
 131                                 "Unfound second block" );
 132       }
 133
 134       // Two steps are necessary to please VC++
 135       //
 136       // 2 pixels 12bit =     [0xABCDEF]
 137       // 2 pixels 16bit = [0x0ABD] + [0x0FCE]
 138       //                        A                     B                 D
 139       *localDecompres++ =  ((b0 >> 4) << 8) + ((b0 & 0x0f) << 4) + (b1 & 0x0f);
 140       //                        F                     C                 E
 141       *localDecompres++ =  ((b2 & 0x0f) << 8) + ((b1 >> 4) << 4) + (b2 >> 4);
 142
 143       /// \todo JPR Troubles expected on Big-Endian processors ?
 144    }
 145 }
 146
 147 /**
 148  * \brief     Try to deal with RLE 16 Bits.
 149  *            We assume the RLE has allready been parsed and loaded in
 150  *            Decompressed (through \ref ReadAndDecompressJPEGFile ).
 151  *            We here need to make 16 Bits Pixels from Low Byte and
 152  *            High Byte 'Planes'...(for what it may mean)
 153  * @return    Boolean
 154  */
 155 bool PixelConvert::DecompressRLE16BitsFromRLE8Bits( int NumberOfFrames )
 156 {
 157    size_t PixelNumber = XSize * YSize;
 158    size_t decompressedSize = XSize * YSize * NumberOfFrames;
 159
 160    // We assumed Decompressed contains the decoded RLE pixels but as
 161    // 8 bits per pixel. In order to convert those pixels to 16 bits
 162    // per pixel we cannot work in place within Decompressed and hence
 163    // we copy it in a safe place, say copyDecompressed.
 164
 165    uint8_t* copyDecompressed = new uint8_t[ decompressedSize * 2 ];
 166    memmove( copyDecompressed, Decompressed, decompressedSize * 2 );
 167
 168    uint8_t* x = Decompressed;
 169    uint8_t* a = copyDecompressed;
 170    uint8_t* b = a + PixelNumber;
 171
 172    for ( int i = 0; i < NumberOfFrames; i++ )
 173    {
 174       for ( unsigned int j = 0; j < PixelNumber; j++ )
 175       {
 176          *(x++) = *(a++);
 177          *(x++) = *(b++);
 178       }
 179    }
 180
 181    delete[] copyDecompressed;
 182
 183    /// \todo check that operator new []didn't fail, and sometimes return false
 184    return true;
 185 }
 186
 187 /**
 188  * \brief Implementation of the RLE decoding algorithm for decompressing
 189  *        a RLE fragment. [refer to PS 3.5-2003, section G.3.2 p 86]
 190  * @param subDecompressed Sub region of \ref Decompressed where the de
 191  *        decoded fragment should be placed.
 192  * @param fragmentSize The length of the binary fragment as found on the disk.
 193  * @param decompressedSegmentSize The expected length of the fragment ONCE
 194  *        decompressed.
 195  * @param fp File Pointer: on entry the position should be the one of
 196  *        the fragment to be decoded.
 197  */
 198 bool PixelConvert::ReadAndDecompressRLEFragment( uint8_t* subDecompressed,
 199                                                  long fragmentSize,
 200                                                  long decompressedSegmentSize,
 201                                                  std::ifstream* fp )
 202 {
 203    int8_t count;
 204    long numberOfOutputBytes = 0;
 205    long numberOfReadBytes = 0;
 206
 207    while( numberOfOutputBytes < decompressedSegmentSize )
 208    {
 209       fp->read( (char*)&count, 1 );
 210       numberOfReadBytes += 1;
 211       if ( count >= 0 )
 212       // Note: count <= 127 comparison is always true due to limited range
 213       //       of data type int8_t [since the maximum of an exact width
 214       //       signed integer of width N is 2^(N-1) - 1, which for int8_t
 215       //       is 127].
 216       {
 217          fp->read( (char*)subDecompressed, count + 1);
 218          numberOfReadBytes   += count + 1;
 219          subDecompressed     += count + 1;
 220          numberOfOutputBytes += count + 1;
 221       }
 222       else
 223       {
 224          if ( ( count <= -1 ) && ( count >= -127 ) )
 225          {
 226             int8_t newByte;
 227             fp->read( (char*)&newByte, 1);
 228             numberOfReadBytes += 1;
 229             for( int i = 0; i < -count + 1; i++ )
 230             {
 231                subDecompressed[i] = newByte;
 232             }
 233             subDecompressed     += -count + 1;
 234             numberOfOutputBytes += -count + 1;
 235          }
 236       }
 237       // if count = 128 output nothing
 238
 239       if ( numberOfReadBytes > fragmentSize )
 240       {
 241          dbg.Verbose(0, "PixelConvert::ReadAndDecompressRLEFragment: we "
 242                         "read more bytes than the segment size.");
 243          return false;
 244       }
 245    }
 246    return true;
 247 }
 248
 249 /**
 250  * \brief     Reads from disk the Pixel Data of 'Run Length Encoded'
 251  *            Dicom encapsulated file and decompress it.
 252  * @param     fp already open File Pointer
 253  *            at which the pixel data should be copied
 254  * @return    Boolean
 255  */
 256 bool PixelConvert::ReadAndDecompressRLEFile( std::ifstream* fp )
 257 {
 258    uint8_t* subDecompressed = Decompressed;
 259    long decompressedSegmentSize = XSize * YSize;
 260
 261    // Loop on the frame[s]
 262    for( RLEFramesInfo::RLEFrameList::iterator
 263         it  = RLEInfo->Frames.begin();
 264         it != RLEInfo->Frames.end();
 265       ++it )
 266    {
 267       // Loop on the fragments
 268       for( unsigned int k = 1; k <= (*it)->NumberFragments; k++ )
 269       {
 270          fp->seekg(  (*it)->Offset[k] , std::ios_base::beg );
 271          (void)ReadAndDecompressRLEFragment( subDecompressed,
 272                                              (*it)->Length[k],
 273                                              decompressedSegmentSize,
 274                                              fp );
 275          subDecompressed += decompressedSegmentSize;
 276       }
 277    }
 278
 279    if ( BitsAllocated == 16 )
 280    {
 281       // Try to deal with RLE 16 Bits
 282       (void)DecompressRLE16BitsFromRLE8Bits( ZSize );
 283    }
 284
 285    return true;
 286 }
 287
 288 /**
 289  * \brief Swap the bytes, according to \ref SwapCode.
 290  */
 291 void PixelConvert::ConvertSwapZone()
 292 {
 293    unsigned int i;
 294
 295    if( BitsAllocated == 16 )
 296    {
 297       uint16_t* im16 = (uint16_t*)Decompressed;
 298       switch( SwapCode )
 299       {
 300          case 0:
 301          case 12:
 302          case 1234:
 303             break;
 304          case 21:
 305          case 3412:
 306          case 2143:
 307          case 4321:
 308             for( i = 0; i < DecompressedSize / 2; i++ )
 309             {
 310                im16[i]= (im16[i] >> 8) | (im16[i] << 8 );
 311             }
 312             break;
 313          default:
 314             dbg.Verbose( 0, "PixelConvert::ConvertSwapZone: SwapCode value "
 315                             "(16 bits) not allowed." );
 316       }
 317    }
 318    else if( BitsAllocated == 32 )
 319    {
 320       uint32_t s32;
 321       uint16_t high;
 322       uint16_t low;
 323       uint32_t* im32 = (uint32_t*)Decompressed;
 324       switch ( SwapCode )
 325       {
 326          case 0:
 327          case 1234:
 328             break;
 329          case 4321:
 330             for( i = 0; i < DecompressedSize / 4; i++ )
 331             {
 332                low     = im32[i] & 0x0000ffff;  // 4321
 333                high    = im32[i] >> 16;
 334                high    = ( high >> 8 ) | ( high << 8 );
 335                low     = ( low  >> 8 ) | ( low  << 8 );
 336                s32     = low;
 337                im32[i] = ( s32 << 16 ) | high;
 338             }
 339             break;
 340          case 2143:
 341             for( i = 0; i < DecompressedSize / 4; i++ )
 342             {
 343                low     = im32[i] & 0x0000ffff;   // 2143
 344                high    = im32[i] >> 16;
 345                high    = ( high >> 8 ) | ( high << 8 );
 346                low     = ( low  >> 8 ) | ( low  << 8 );
 347                s32     = high;
 348                im32[i] = ( s32 << 16 ) | low;
 349             }
 350             break;
 351          case 3412:
 352             for( i = 0; i < DecompressedSize / 4; i++ )
 353             {
 354                low     = im32[i] & 0x0000ffff; // 3412
 355                high    = im32[i] >> 16;
 356                s32     = low;
 357                im32[i] = ( s32 << 16 ) | high;
 358             }
 359             break;
 360          default:
 361             dbg.Verbose( 0, "PixelConvert::ConvertSwapZone: SwapCode value "
 362                             "(32 bits) not allowed." );
 363       }
 364    }
 365 }
 366
 367 /**
 368  * \brief Deal with endianity i.e. re-arange bytes inside the integer
 369  */
 370 void PixelConvert::ConvertReorderEndianity()
 371 {
 372    if ( BitsAllocated != 8 )
 373    {
 374       ConvertSwapZone();
 375    }
 376
 377    // Special kludge in order to deal with xmedcon broken images:
 378    if (  ( BitsAllocated == 16 )
 379        && ( BitsStored < BitsAllocated )
 380        && ( ! PixelSign ) )
 381    {
 382       int l = (int)( DecompressedSize / ( BitsAllocated / 8 ) );
 383       uint16_t *deb = (uint16_t *)Decompressed;
 384       for(int i = 0; i<l; i++)
 385       {
 386          if( *deb == 0xffff )
 387          {
 388            *deb = 0;
 389          }
 390          deb++;
 391       }
 392    }
 393 }
 394
 395 /**
 396  * \brief     Reads from disk the Pixel Data of JPEG Dicom encapsulated
 397  &            file and decompress it.
 398  * @param     fp File Pointer
 399  * @return    Boolean
 400  */
 401 bool PixelConvert::ReadAndDecompressJPEGFile( std::ifstream* fp )
 402 {
 403    uint8_t* localDecompressed = Decompressed;
 404    // Loop on the fragment[s]
 405    for( JPEGFragmentsInfo::JPEGFragmentsList::iterator
 406         it  = JPEGInfo->Fragments.begin();
 407         it != JPEGInfo->Fragments.end();
 408       ++it )
 409    {
 410       fp->seekg( (*it)->Offset, std::ios_base::beg);
 411
 412       if ( IsJPEG2000 )
 413       {
 414          if ( ! gdcm_read_JPEG2000_file( fp,localDecompressed ) )
 415          {
 416             return false;
 417          }
 418       }
 419       else if ( BitsStored == 8)
 420       {
 421          // JPEG Lossy : call to IJG 6b
 422          if ( ! gdcm_read_JPEG_file8( fp, localDecompressed ) )
 423          {
 424             return false;
 425          }
 426       }
 427       else if ( BitsStored == 12)
 428       {
 429          // Reading Fragment pixels
 430          if ( ! gdcm_read_JPEG_file12 ( fp, localDecompressed ) )
 431          {
 432             return false;
 433          }
 434       }
 435       else if ( BitsStored == 16)
 436       {
 437          // Reading Fragment pixels
 438          if ( ! gdcm_read_JPEG_file16 ( fp, localDecompressed ) )
 439          {
 440             return false;
 441          }
 442          //assert( IsJPEGLossless );
 443       }
 444       else
 445       {
 446          // other JPEG lossy not supported
 447          dbg.Error("PixelConvert::ReadAndDecompressJPEGFile: unknown "
 448                    "jpeg lossy compression ");
 449          return false;
 450       }
 451
 452       // Advance to next free location in Decompressed
 453       // for next fragment decompression (if any)
 454       int length = XSize * YSize * SamplesPerPixel;
 455       int numberBytes = BitsAllocated / 8;
 456
 457       localDecompressed += length * numberBytes;
 458    }
 459    return true;
 460 }
 461
 462 /**
 463  * \brief  Re-arrange the bits within the bytes.
 464  * @return Boolean
 465  */
 466 bool PixelConvert::ConvertReArrangeBits() throw ( FormatError )
 467 {
 468    if ( BitsStored != BitsAllocated )
 469    {
 470       int l = (int)( DecompressedSize / ( BitsAllocated / 8 ) );
 471       if ( BitsAllocated == 16 )
 472       {
 473          uint16_t mask = 0xffff;
 474          mask = mask >> ( BitsAllocated - BitsStored );
 475          uint16_t* deb = (uint16_t*)Decompressed;
 476          for(int i = 0; i<l; i++)
 477          {
 478             *deb = (*deb >> (BitsStored - HighBitPosition - 1)) & mask;
 479             deb++;
 480          }
 481       }
 482       else if ( BitsAllocated == 32 )
 483       {
 484          uint32_t mask = 0xffffffff;
 485          mask = mask >> ( BitsAllocated - BitsStored );
 486          uint32_t* deb = (uint32_t*)Decompressed;
 487          for(int i = 0; i<l; i++)
 488          {
 489             *deb = (*deb >> (BitsStored - HighBitPosition - 1)) & mask;
 490             deb++;
 491          }
 492       }
 493       else
 494       {
 495          dbg.Verbose(0, "PixelConvert::ConvertReArrangeBits: weird image");
 496          throw FormatError( "PixelConvert::ConvertReArrangeBits()",
 497                                 "weird image !?" );
 498       }
 499    }
 500    return true;
 501 }
 502
 503 /**
 504  * \brief   Convert (Y plane, cB plane, cR plane) to RGB pixels
 505  * \warning Works on all the frames at a time
 506  */
 507 void PixelConvert::ConvertYcBcRPlanesToRGBPixels()
 508 {
 509    uint8_t* localDecompressed = Decompressed;
 510    uint8_t* copyDecompressed = new uint8_t[ DecompressedSize ];
 511    memmove( copyDecompressed, localDecompressed, DecompressedSize );
 512
 513    // to see the tricks about YBR_FULL, YBR_FULL_422,
 514    // YBR_PARTIAL_422, YBR_ICT, YBR_RCT have a look at :
 515    // ftp://medical.nema.org/medical/dicom/final/sup61_ft.pdf
 516    // and be *very* affraid
 517    //
 518    int l        = XSize * YSize;
 519    int nbFrames = ZSize;
 520
 521    uint8_t* a = copyDecompressed;
 522    uint8_t* b = copyDecompressed + l;
 523    uint8_t* c = copyDecompressed + l + l;
 524    double R, G, B;
 525
 526    /// \todo : Replace by the 'well known' integer computation
 527    ///         counterpart. Refer to
 528    ///            http://lestourtereaux.free.fr/papers/data/yuvrgb.pdf
 529    ///         for code optimisation.
 530
 531    for ( int i = 0; i < nbFrames; i++ )
 532    {
 533       for ( int j = 0; j < l; j++ )
 534       {
 535          R = 1.164 *(*a-16) + 1.596 *(*c -128) + 0.5;
 536          G = 1.164 *(*a-16) - 0.813 *(*c -128) - 0.392 *(*b -128) + 0.5;
 537          B = 1.164 *(*a-16) + 2.017 *(*b -128) + 0.5;
 538
 539          if (R < 0.0)   R = 0.0;
 540          if (G < 0.0)   G = 0.0;
 541          if (B < 0.0)   B = 0.0;
 542          if (R > 255.0) R = 255.0;
 543          if (G > 255.0) G = 255.0;
 544          if (B > 255.0) B = 255.0;
 545
 546          *(localDecompressed++) = (uint8_t)R;
 547          *(localDecompressed++) = (uint8_t)G;
 548          *(localDecompressed++) = (uint8_t)B;
 549          a++;
 550          b++;
 551          c++;
 552       }
 553    }
 554    delete[] copyDecompressed;
 555 }
 556
 557 /**
 558  * \brief   Convert (Red plane, Green plane, Blue plane) to RGB pixels
 559  * \warning Works on all the frames at a time
 560  */
 561 void PixelConvert::ConvertRGBPlanesToRGBPixels()
 562 {
 563    uint8_t* localDecompressed = Decompressed;
 564    uint8_t* copyDecompressed = new uint8_t[ DecompressedSize ];
 565    memmove( copyDecompressed, localDecompressed, DecompressedSize );
 566
 567    int l = XSize * YSize * ZSize;
 568
 569    uint8_t* a = copyDecompressed;
 570    uint8_t* b = copyDecompressed + l;
 571    uint8_t* c = copyDecompressed + l + l;
 572
 573    for (int j = 0; j < l; j++)
 574    {
 575       *(localDecompressed++) = *(a++);
 576       *(localDecompressed++) = *(b++);
 577       *(localDecompressed++) = *(c++);
 578    }
 579    delete[] copyDecompressed;
 580 }
 581
 582 bool PixelConvert::ReadAndDecompressPixelData( std::ifstream* fp )
 583 {
 584    ComputeDecompressedAndRGBSizes();
 585    AllocateDecompressed();
 586    //////////////////////////////////////////////////
 587    //// First stage: get our hands on the Pixel Data.
 588    if ( !fp )
 589    {
 590      dbg.Verbose( 0, "PixelConvert::ReadAndDecompressPixelData: "
 591                      "unavailable file pointer." );
 592       return false;
 593    }
 594
 595    fp->seekg( PixelOffset, std::ios_base::beg );
 596    if( fp->fail() || fp->eof()) //Fp->gcount() == 1
 597    {
 598      dbg.Verbose( 0, "PixelConvert::ReadAndDecompressPixelData: "
 599                      "unable to find PixelOffset in file." );
 600       return false;
 601    }
 602
 603    //////////////////////////////////////////////////
 604    //// Second stage: read from disk dans decompress.
 605    if ( BitsAllocated == 12 )
 606    {
 607       ReadAndDecompress12BitsTo16Bits( fp);
 608    }
 609    else if ( IsDecompressed )
 610    {
 611       fp->read( (char*)Decompressed, PixelDataLength);
 612       if ( fp->fail() || fp->eof())//Fp->gcount() == 1
 613       {
 614          dbg.Verbose( 0, "PixelConvert::ReadAndDecompressPixelData: "
 615                          "reading of decompressed pixel data failed." );
 616          return false;
 617       }
 618    }
 619    else if ( IsRLELossless )
 620    {
 621       if ( ! ReadAndDecompressRLEFile( fp ) )
 622       {
 623          dbg.Verbose( 0, "PixelConvert::ReadAndDecompressPixelData: "
 624                          "RLE decompressor failed." );
 625          return false;
 626       }
 627    }
 628    else
 629    {
 630       // Default case concerns JPEG family
 631       if ( ! ReadAndDecompressJPEGFile( fp ) )
 632       {
 633          dbg.Verbose( 0, "PixelConvert::ReadAndDecompressPixelData: "
 634                          "JPEG decompressor failed." );
 635          return false;
 636       }
 637    }
 638
 639    ////////////////////////////////////////////
 640    //// Third stage: twigle the bytes and bits.
 641    ConvertReorderEndianity();
 642    ConvertReArrangeBits();
 643    ConvertHandleColor();
 644
 645    return true;
 646 }
 647
 648 void PixelConvert::ConvertHandleColor()
 649 {
 650    //////////////////////////////////
 651    // Deal with the color decoding i.e. handle:
 652    //   - R, G, B planes (as opposed to RGB pixels)
 653    //   - YBR (various) encodings.
 654    //   - LUT[s] (or "PALETTE COLOR").
 655    //
 656    // The classification in the color decoding schema is based on the blending
 657    // of two Dicom tags values:
 658    // * "Photometric Interpretation" for which we have the cases:
 659    //  - [Photo A] MONOCHROME[1|2] pictures,
 660    //  - [Photo B] RGB or YBR_FULL_422 (which acts as RGB),
 661    //  - [Photo C] YBR_* (with the above exception of YBR_FULL_422)
 662    //  - [Photo D] "PALETTE COLOR" which indicates the presence of LUT[s].
 663    // * "Planar Configuration" for which we have the cases:
 664    //  - [Planar 0] 0 then Pixels are already RGB
 665    //  - [Planar 1] 1 then we have 3 planes : R, G, B,
 666    //  - [Planar 2] 2 then we have 1 gray Plane and 3 LUTs
 667    //
 668    // Now in theory, one could expect some coherence when blending the above
 669    // cases. For example we should not encounter files belonging at the
 670    // time to case [Planar 0] and case [Photo D].
 671    // Alas, this was only theory ! Because in practice some odd (read ill
 672    // formated Dicom) files (e.g. gdcmData/US-PAL-8-10x-echo.dcm) we encounter:
 673    //     - "Planar Configuration" = 0,
 674    //     - "Photometric Interpretation" = "PALETTE COLOR".
 675    // Hence gdcm shall use the folowing "heuristic" in order to be tolerant
 676    // towards Dicom-non-conformance files:
 677    //   << whatever the "Planar Configuration" value might be, a
 678    //      "Photometric Interpretation" set to "PALETTE COLOR" forces
 679    //      a LUT intervention >>
 680    //
 681    // Now we are left with the following handling of the cases:
 682    // - [Planar 0] OR  [Photo A] no color decoding (since respectively
 683    //       Pixels are already RGB and monochrome pictures have no color :),
 684    // - [Planar 1] AND [Photo B] handled with ConvertRGBPlanesToRGBPixels()
 685    // - [Planar 1] AND [Photo C] handled with ConvertYcBcRPlanesToRGBPixels()
 686    // - [Planar 2] OR  [Photo D] requires LUT intervention.
 687
 688    if ( ! IsDecompressedRGB() )
 689    {
 690       // [Planar 2] OR  [Photo D]: LUT intervention done outside
 691       return;
 692    }
 693
 694    if ( PlanarConfiguration == 1 )
 695    {
 696       if ( IsYBRFull )
 697       {
 698          // [Planar 1] AND [Photo C] (remember YBR_FULL_422 acts as RGB)
 699          ConvertYcBcRPlanesToRGBPixels();
 700       }
 701       else
 702       {
 703          // [Planar 1] AND [Photo C]
 704          ConvertRGBPlanesToRGBPixels();
 705       }
 706       return;
 707    }
 708
 709    // When planarConf is 0, and RLELossless (forbidden by Dicom norm)
 710    // pixels need to be RGB-fied anyway
 711    if (IsRLELossless)
 712    {
 713      ConvertRGBPlanesToRGBPixels();
 714    }
 715    // In *normal *case, when planarConf is 0, pixels are already in RGB
 716 }
 717
 718 /**
 719  * \brief Predicate to know wether the image[s] (once decompressed) is RGB.
 720  * \note See comments of \ref ConvertHandleColor
 721  */
 722 bool PixelConvert::IsDecompressedRGB()
 723 {
 724    if (   IsMonochrome
 725        || PlanarConfiguration == 2
 726        || IsPaletteColor )
 727    {
 728       return false;
 729    }
 730    return true;
 731 }
 732
 733 void PixelConvert::ComputeDecompressedAndRGBSizes()
 734 {
 735    int bitsAllocated = BitsAllocated;
 736    // Number of "Bits Allocated" is fixed to 16 when it's 12, since
 737    // in this case we will expand the image to 16 bits (see
 738    //    \ref ReadAndDecompress12BitsTo16Bits() )
 739    if (  BitsAllocated == 12 )
 740    {
 741       bitsAllocated = 16;
 742    }
 743
 744    DecompressedSize =  XSize * YSize * ZSize
 745                      * ( bitsAllocated / 8 )
 746                      * SamplesPerPixel;
 747    if ( HasLUT )
 748    {
 749       RGBSize = 3 * DecompressedSize;
 750    }
 751
 752 }
 753
 754 void PixelConvert::GrabInformationsFromHeader( Header* header )
 755 {
 756    // Just in case some access to a Header element requires disk access.
 757    // Note: gdcmDocument::Fp is leaved open after OpenFile.
 758   std::ifstream* fp = header->OpenFile();
 759    // Number of Bits Allocated for storing a Pixel is defaulted to 16
 760    // when absent from the header.
 761    BitsAllocated = header->GetBitsAllocated();
 762    if ( BitsAllocated == 0 )
 763    {
 764       BitsAllocated = 16;
 765    }
 766
 767    // Number of "Bits Stored" defaulted to number of "Bits Allocated"
 768    // when absent from the header.
 769    BitsStored = header->GetBitsStored();
 770    if ( BitsStored == 0 )
 771    {
 772       BitsStored = BitsAllocated;
 773    }
 774
 775    // High Bit Position
 776    HighBitPosition = header->GetHighBitPosition();
 777    if ( HighBitPosition == 0 )
 778    {
 779       HighBitPosition = BitsAllocated - 1;
 780    }
 781
 782    XSize = header->GetXSize();
 783    YSize = header->GetYSize();
 784    ZSize = header->GetZSize();
 785    SamplesPerPixel = header->GetSamplesPerPixel();
 786    PixelSize = header->GetPixelSize();
 787    PixelSign = header->IsSignedPixelData();
 788    SwapCode  = header->GetSwapCode();
 789    TransferSyntaxType ts = header->GetTransferSyntax();
 790    IsDecompressed =
 791         ( ! header->IsDicomV3() )
 792      || ts == ImplicitVRLittleEndian
 793      || ts == ExplicitVRLittleEndian
 794      || ts == ExplicitVRBigEndian
 795      || ts == DeflatedExplicitVRLittleEndian;
 796    IsJPEG2000     = header->IsJPEG2000();
 797    IsJPEGLossless = header->IsJPEGLossless();
 798    IsRLELossless  =  ( ts == RLELossless );
 799    PixelOffset     = header->GetPixelOffset();
 800    PixelDataLength = header->GetPixelAreaLength();
 801    RLEInfo  = header->GetRLEInfo();
 802    JPEGInfo = header->GetJPEGInfo();
 803
 804    PlanarConfiguration = header->GetPlanarConfiguration();
 805    IsMonochrome = header->IsMonochrome();
 806    IsPaletteColor = header->IsPaletteColor();
 807    IsYBRFull = header->IsYBRFull();
 808
 809    /////////////////////////////////////////////////////////////////
 810    // LUT section:
 811    HasLUT = header->HasLUT();
 812    if ( HasLUT )
 813    {
 814       LutRedDescriptor   = header->GetEntryByNumber( 0x0028, 0x1101 );
 815       LutGreenDescriptor = header->GetEntryByNumber( 0x0028, 0x1102 );
 816       LutBlueDescriptor  = header->GetEntryByNumber( 0x0028, 0x1103 );
 817
 818       // Depending on the value of Document::MAX_SIZE_LOAD_ELEMENT_VALUE
 819       // [ refer to invocation of Document::SetMaxSizeLoadEntry() in
 820       // Document::Document() ], the loading of the value (content) of a
 821       // [Bin|Val]Entry occurence migth have been hindered (read simply NOT
 822       // loaded). Hence, we first try to obtain the LUTs data from the header
 823       // and when this fails we read the LUTs data directely from disk.
 824       /// \todo Reading a [Bin|Val]Entry directly from disk is a kludge.
 825       ///       We should NOT bypass the [Bin|Val]Entry class. Instead
 826       ///       an access to an UNLOADED content of a [Bin|Val]Entry occurence
 827       ///       (e.g. BinEntry::GetBinArea()) should force disk access from
 828       ///       within the [Bin|Val]Entry class itself. The only problem
 829       ///       is that the [Bin|Val]Entry is unaware of the FILE* is was
 830       ///       parsed from. Fix that. FIXME.
 831
 832       ////// Red round:
 833       LutRedData = (uint8_t*)header->GetEntryBinAreaByNumber( 0x0028, 0x1201 );
 834       if ( ! LutRedData )
 835       {
 836          // Read the Lut Data from disk
 837          DocEntry* lutRedDataEntry = header->GetDocEntryByNumber( 0x0028,
 838                                                                   0x1201 );
 839          LutRedData = new uint8_t[ lutRedDataEntry->GetLength() ];
 840          fp->seekg(  lutRedDataEntry->GetOffset() ,std::ios_base::beg );
 841          fp->read( (char*)LutRedData, (size_t)lutRedDataEntry->GetLength());
 842          if ( fp->fail() || fp->eof())//Fp->gcount() == 1
 843          {
 844             dbg.Verbose(0, "PixelConvert::GrabInformationsFromHeader: "
 845                             "unable to read red LUT data" );
 846             return;
 847          }
 848       }
 849
 850       ////// Green round:
 851       LutGreenData = (uint8_t*)header->GetEntryBinAreaByNumber(0x0028, 0x1202 );
 852       if ( ! LutGreenData)
 853       {
 854          // Read the Lut Data from disk
 855          DocEntry* lutGreenDataEntry = header->GetDocEntryByNumber( 0x0028,
 856                                                                     0x1202 );
 857          LutGreenData = new uint8_t[ lutGreenDataEntry->GetLength() ];
 858          fp->seekg( lutGreenDataEntry->GetOffset() , std::ios_base::beg );
 859          fp->read( (char*)LutGreenData, (size_t)lutGreenDataEntry->GetLength() );
 860          if ( fp->fail() || fp->eof())//Fp->gcount() == 1
 861          {
 862             dbg.Verbose(0, "PixelConvert::GrabInformationsFromHeader: "
 863                            "unable to read green LUT data" );
 864             return;
 865          }
 866       }
 867
 868       ////// Blue round:
 869       LutBlueData = (uint8_t*)header->GetEntryBinAreaByNumber( 0x0028, 0x1203 );
 870       if ( ! LutBlueData )
 871       {
 872          // Read the Lut Data from disk
 873          DocEntry* lutBlueDataEntry  = header->GetDocEntryByNumber( 0x0028,
 874                                                                     0x1203 );
 875          LutBlueData = new uint8_t[ lutBlueDataEntry->GetLength() ];
 876          fp->seekg(  lutBlueDataEntry->GetOffset() , std::ios_base::beg );
 877          fp->read( (char*)LutBlueData, (size_t)lutBlueDataEntry->GetLength() );
 878          if ( fp->fail() || fp->eof())//Fp->gcount() == 1
 879          {
 880             dbg.Verbose(0, "PixelConvert::GrabInformationsFromHeader: "
 881                            "unable to read blue LUT data" );
 882             return;
 883          }
 884       }
 885    }
 886
 887    header->CloseFile();
 888 }
 889
 890 /**
 891  * \brief Build Red/Green/Blue/Alpha LUT from Header
 892  *         when (0028,0004),Photometric Interpretation = [PALETTE COLOR ]
 893  *          and (0028,1101),(0028,1102),(0028,1102)
 894  *            - xxx Palette Color Lookup Table Descriptor - are found
 895  *          and (0028,1201),(0028,1202),(0028,1202)
 896  *            - xxx Palette Color Lookup Table Data - are found
 897  * \warning does NOT deal with :
 898  *   0028 1100 Gray Lookup Table Descriptor (Retired)
 899  *   0028 1221 Segmented Red Palette Color Lookup Table Data
 900  *   0028 1222 Segmented Green Palette Color Lookup Table Data
 901  *   0028 1223 Segmented Blue Palette Color Lookup Table Data
 902  *   no known Dicom reader deals with them :-(
 903  * @return a RGBA Lookup Table
 904  */
 905 void PixelConvert::BuildLUTRGBA()
 906 {
 907    if ( LutRGBA )
 908    {
 909       return;
 910    }
 911    // Not so easy : see
 912    // http://www.barre.nom.fr/medical/dicom2/limitations.html#Color%20Lookup%20Tables
 913
 914    if ( ! IsPaletteColor )
 915    {
 916       return;
 917    }
 918
 919    if (   LutRedDescriptor   == GDCM_UNFOUND
 920        || LutGreenDescriptor == GDCM_UNFOUND
 921        || LutBlueDescriptor  == GDCM_UNFOUND )
 922    {
 923       return;
 924    }
 925
 926    ////////////////////////////////////////////
 927    // Extract the info from the LUT descriptors
 928    int lengthR;   // Red LUT length in Bytes
 929    int debR;      // Subscript of the first Lut Value
 930    int nbitsR;    // Lut item size (in Bits)
 931    int nbRead = sscanf( LutRedDescriptor.c_str(),
 932                         "%d\\%d\\%d",
 933                         &lengthR, &debR, &nbitsR );
 934    if( nbRead != 3 )
 935    {
 936       dbg.Verbose(0, "PixelConvert::BuildLUTRGBA: wrong red LUT descriptor");
 937    }
 938
 939    int lengthG;  // Green LUT length in Bytes
 940    int debG;     // Subscript of the first Lut Value
 941    int nbitsG;   // Lut item size (in Bits)
 942    nbRead = sscanf( LutGreenDescriptor.c_str(),
 943                     "%d\\%d\\%d",
 944                     &lengthG, &debG, &nbitsG );
 945    if( nbRead != 3 )
 946    {
 947       dbg.Verbose(0, "PixelConvert::BuildLUTRGBA: wrong green LUT descriptor");
 948    }
 949
 950    int lengthB;  // Blue LUT length in Bytes
 951    int debB;     // Subscript of the first Lut Value
 952    int nbitsB;   // Lut item size (in Bits)
 953    nbRead = sscanf( LutRedDescriptor.c_str(),
 954                     "%d\\%d\\%d",
 955                     &lengthB, &debB, &nbitsB );
 956    if( nbRead != 3 )
 957    {
 958       dbg.Verbose(0, "PixelConvert::BuildLUTRGBA: wrong blue LUT descriptor");
 959    }
 960
 961    ////////////////////////////////////////////////////////
 962    if ( ( ! LutRedData ) || ( ! LutGreenData ) || ( ! LutBlueData ) )
 963    {
 964       return;
 965    }
 966
 967    ////////////////////////////////////////////////
 968    // forge the 4 * 8 Bits Red/Green/Blue/Alpha LUT
 969    LutRGBA = new uint8_t[ 1024 ]; // 256 * 4 (R, G, B, Alpha)
 970    if ( !LutRGBA )
 971    {
 972       return;
 973    }
 974    memset( LutRGBA, 0, 1024 );
 975
 976    int mult;
 977    if ( ( nbitsR == 16 ) && ( BitsAllocated == 8 ) )
 978    {
 979       // when LUT item size is different than pixel size
 980       mult = 2; // high byte must be = low byte
 981    }
 982    else
 983    {
 984       // See PS 3.3-2003 C.11.1.1.2 p 619
 985       mult = 1;
 986    }
 987
 988    // if we get a black image, let's just remove the '+1'
 989    // from 'i*mult+1' and check again
 990    // if it works, we shall have to check the 3 Palettes
 991    // to see which byte is ==0 (first one, or second one)
 992    // and fix the code
 993    // We give up the checking to avoid some (useless ?)overhead
 994    // (optimistic asumption)
 995    int i;
 996    uint8_t* a = LutRGBA + 0;
 997    for( i=0; i < lengthR; ++i )
 998    {
 999       *a = LutRedData[i*mult+1];
1000       a += 4;
1001    }
1002
1003    a = LutRGBA + 1;
1004    for( i=0; i < lengthG; ++i)
1005    {
1006       *a = LutGreenData[i*mult+1];
1007       a += 4;
1008    }
1009
1010    a = LutRGBA + 2;
1011    for(i=0; i < lengthB; ++i)
1012    {
1013       *a = LutBlueData[i*mult+1];
1014       a += 4;
1015    }
1016
1017    a = LutRGBA + 3;
1018    for(i=0; i < 256; ++i)
1019    {
1020       *a = 1; // Alpha component
1021       a += 4;
1022    }
1023 }
1024
1025 /**
1026  * \brief Build the RGB image from the Decompressed imagage and the LUTs.
1027  */
1028 bool PixelConvert::BuildRGBImage()
1029 {
1030    if ( RGB )
1031    {
1032       // The job is already done.
1033       return true;
1034    }
1035
1036    if ( ! Decompressed )
1037    {
1038       // The job can't be done
1039       return false;
1040    }
1041
1042    BuildLUTRGBA();
1043    if ( ! LutRGBA )
1044    {
1045       // The job can't be done
1046       return false;
1047    }
1048
1049    // Build RGB Pixels
1050    AllocateRGB();
1051    uint8_t* localRGB = RGB;
1052    for (size_t i = 0; i < DecompressedSize; ++i )
1053    {
1054       int j  = Decompressed[i] * 4;
1055       *localRGB++ = LutRGBA[j];
1056       *localRGB++ = LutRGBA[j+1];
1057       *localRGB++ = LutRGBA[j+2];
1058    }
1059    return true;
1060 }
1061
1062 /**
1063  * \brief        Print self.
1064  * @param indent Indentation string to be prepended during printing.
1065  * @param os     Stream to print to.
1066  */
1067 void PixelConvert::Print( std::string indent, std::ostream &os )
1068 {
1069    os << indent
1070       << "--- Pixel information -------------------------"
1071       << std::endl;
1072    os << indent
1073       << "Pixel Data: offset " << PixelOffset
1074       << " x" << std::hex << PixelOffset << std::dec
1075       << "   length " << PixelDataLength
1076       << " x" << std::hex << PixelDataLength << std::dec
1077       << std::endl;
1078
1079    if ( IsRLELossless )
1080    {
1081       if ( RLEInfo )
1082       {
1083          RLEInfo->Print( indent, os );
1084       }
1085       else
1086       {
1087          dbg.Verbose(0, "PixelConvert::Print: set as RLE file "
1088                         "but NO RLEinfo present.");
1089       }
1090    }
1091
1092    if ( IsJPEG2000 || IsJPEGLossless )
1093    {
1094       if ( JPEGInfo )
1095       {
1096          JPEGInfo->Print( indent, os );
1097       }
1098       else
1099       {
1100          dbg.Verbose(0, "PixelConvert::Print: set as JPEG file "
1101                         "but NO JPEGinfo present.");
1102       }
1103    }
1104 }
1105
1106 } // end namespace gdcm
1107
1108 // NOTES on File internal calls
1109 // User
1110 // ---> GetImageData
1111 //     ---> GetImageDataIntoVector
1112 //        |---> GetImageDataIntoVectorRaw
1113 //        | lut intervention
1114 // User
1115 // ---> GetImageDataRaw
1116 //     ---> GetImageDataIntoVectorRaw