src/gdcmPixelReadConvert.cxx

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