src/gdcmPixelReadConvert.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: gdcmPixelReadConvert.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2005/02/03 10:03:07 $
   7   Version:   $Revision: 1.46 $
   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 //-----------------------------------------------------------------------------
  34 #define str2num(str, typeNum) *((typeNum *)(str))
  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 PixelReadConvert::~PixelReadConvert()
  51 {
  52    Squeeze();
  53 }
  54
  55 //-----------------------------------------------------------------------------
  56 // Public
  57 /**
  58  * \brief Predicate to know whether the image[s] (once Raw) is RGB.
  59  * \note See comments of \ref ConvertHandleColor
  60  */
  61 bool PixelReadConvert::IsRawRGB()
  62 {
  63    if (   IsMonochrome
  64        || PlanarConfiguration == 2
  65        || IsPaletteColor )
  66    {
  67       return false;
  68    }
  69    return true;
  70 }
  71 /**
  72  * \brief Gets various usefull informations from the file header
  73  * @param file gdcm::File pointer
  74  */
  75 void PixelReadConvert::GrabInformationsFromFile( File *file )
  76 {
  77    // Number of Bits Allocated for storing a Pixel is defaulted to 16
  78    // when absent from the file.
  79    BitsAllocated = file->GetBitsAllocated();
  80    if ( BitsAllocated == 0 )
  81    {
  82       BitsAllocated = 16;
  83    }
  84
  85    // Number of "Bits Stored", defaulted to number of "Bits Allocated"
  86    // when absent from the file.
  87    BitsStored = file->GetBitsStored();
  88    if ( BitsStored == 0 )
  89    {
  90       BitsStored = BitsAllocated;
  91    }
  92
  93    // High Bit Position, defaulted to "Bits Allocated" - 1
  94    HighBitPosition = file->GetHighBitPosition();
  95    if ( HighBitPosition == 0 )
  96    {
  97       HighBitPosition = BitsAllocated - 1;
  98    }
  99
 100    XSize           = file->GetXSize();
 101    YSize           = file->GetYSize();
 102    ZSize           = file->GetZSize();
 103    SamplesPerPixel = file->GetSamplesPerPixel();
 104    PixelSize       = file->GetPixelSize();
 105    PixelSign       = file->IsSignedPixelData();
 106    SwapCode        = file->GetSwapCode();
 107    std::string ts  = file->GetTransferSyntax();
 108    IsRaw =
 109         ( ! file->IsDicomV3() )
 110      || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ImplicitVRLittleEndian
 111      || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ImplicitVRLittleEndianDLXGE
 112      || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ExplicitVRLittleEndian
 113      || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ExplicitVRBigEndian
 114      || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::DeflatedExplicitVRLittleEndian;
 115
 116    IsJPEG2000      = Global::GetTS()->IsJPEG2000(ts);
 117    IsJPEGLS        = Global::GetTS()->IsJPEGLS(ts);
 118    IsJPEGLossy     = Global::GetTS()->IsJPEGLossy(ts);
 119    IsJPEGLossless  = Global::GetTS()->IsJPEGLossless(ts);
 120    IsRLELossless   = Global::GetTS()->IsRLELossless(ts);
 121
 122    PixelOffset     = file->GetPixelOffset();
 123    PixelDataLength = file->GetPixelAreaLength();
 124    RLEInfo  = file->GetRLEInfo();
 125    JPEGInfo = file->GetJPEGInfo();
 126
 127    PlanarConfiguration = file->GetPlanarConfiguration();
 128    IsMonochrome = file->IsMonochrome();
 129    IsPaletteColor = file->IsPaletteColor();
 130    IsYBRFull = file->IsYBRFull();
 131
 132    /////////////////////////////////////////////////////////////////
 133    // LUT section:
 134    HasLUT = file->HasLUT();
 135    if ( HasLUT )
 136    {
 137       // Just in case some access to a File element requires disk access.
 138       LutRedDescriptor   = file->GetEntryValue( 0x0028, 0x1101 );
 139       LutGreenDescriptor = file->GetEntryValue( 0x0028, 0x1102 );
 140       LutBlueDescriptor  = file->GetEntryValue( 0x0028, 0x1103 );
 141
 142       // Depending on the value of Document::MAX_SIZE_LOAD_ELEMENT_VALUE
 143       // [ refer to invocation of Document::SetMaxSizeLoadEntry() in
 144       // Document::Document() ], the loading of the value (content) of a
 145       // [Bin|Val]Entry occurence migth have been hindered (read simply NOT
 146       // loaded). Hence, we first try to obtain the LUTs data from the file
 147       // and when this fails we read the LUTs data directly from disk.
 148       /// \TODO Reading a [Bin|Val]Entry directly from disk is a kludge.
 149       ///       We should NOT bypass the [Bin|Val]Entry class. Instead
 150       ///       an access to an UNLOADED content of a [Bin|Val]Entry occurence
 151       ///       (e.g. BinEntry::GetBinArea()) should force disk access from
 152       ///       within the [Bin|Val]Entry class itself. The only problem
 153       ///       is that the [Bin|Val]Entry is unaware of the FILE* is was
 154       ///       parsed from. Fix that. FIXME.
 155
 156       ////// Red round
 157       file->LoadEntryBinArea(0x0028, 0x1201);
 158       LutRedData = (uint8_t*)file->GetEntryBinArea( 0x0028, 0x1201 );
 159       if ( ! LutRedData )
 160       {
 161          gdcmVerboseMacro( "Unable to read Red LUT data" );
 162       }
 163
 164       ////// Green round:
 165       file->LoadEntryBinArea(0x0028, 0x1202);
 166       LutGreenData = (uint8_t*)file->GetEntryBinArea(0x0028, 0x1202 );
 167       if ( ! LutGreenData)
 168       {
 169          gdcmVerboseMacro( "Unable to read Green LUT data" );
 170       }
 171
 172       ////// Blue round:
 173       file->LoadEntryBinArea(0x0028, 0x1203);
 174       LutBlueData = (uint8_t*)file->GetEntryBinArea( 0x0028, 0x1203 );
 175       if ( ! LutBlueData )
 176       {
 177          gdcmVerboseMacro( "Unable to read Blue LUT data" );
 178       }
 179    }
 180
 181    ComputeRawAndRGBSizes();
 182 }
 183
 184 bool PixelReadConvert::ReadAndDecompressPixelData( std::ifstream *fp )
 185 {
 186    // ComputeRawAndRGBSizes is already made by
 187    // ::GrabInformationsFromfile. So, the structure sizes are
 188    // correct
 189    Squeeze();
 190
 191    //////////////////////////////////////////////////
 192    //// First stage: get our hands on the Pixel Data.
 193    if ( !fp )
 194    {
 195       gdcmVerboseMacro( "Unavailable file pointer." );
 196       return false;
 197    }
 198
 199    fp->seekg( PixelOffset, std::ios::beg );
 200    if( fp->fail() || fp->eof())
 201    {
 202       gdcmVerboseMacro( "Unable to find PixelOffset in file." );
 203       return false;
 204    }
 205
 206    AllocateRaw();
 207
 208    //////////////////////////////////////////////////
 209    //// Second stage: read from disk dans decompress.
 210    if ( BitsAllocated == 12 )
 211    {
 212       ReadAndDecompress12BitsTo16Bits( fp);
 213    }
 214    else if ( IsRaw )
 215    {
 216       // This problem can be found when some obvious informations are found
 217       // after the field containing the image data. In this case, these
 218       // bad data are added to the size of the image (in the PixelDataLength
 219       // variable). But RawSize is the right size of the image !
 220       if( PixelDataLength != RawSize)
 221       {
 222          gdcmVerboseMacro( "Mismatch between PixelReadConvert : "
 223                             << PixelDataLength << " and RawSize : " << RawSize );
 224       }
 225       if( PixelDataLength > RawSize)
 226       {
 227          fp->read( (char*)Raw, RawSize);
 228       }
 229       else
 230       {
 231          fp->read( (char*)Raw, PixelDataLength);
 232       }
 233
 234       if ( fp->fail() || fp->eof())
 235       {
 236          gdcmVerboseMacro( "Reading of Raw pixel data failed." );
 237          return false;
 238       }
 239    }
 240    else if ( IsRLELossless )
 241    {
 242       if ( ! RLEInfo->DecompressRLEFile( fp, Raw, XSize, YSize, ZSize, BitsAllocated ) )
 243       {
 244          gdcmVerboseMacro( "RLE decompressor failed." );
 245          return false;
 246       }
 247    }
 248    else
 249    {
 250       // Default case concerns JPEG family
 251       if ( ! ReadAndDecompressJPEGFile( fp ) )
 252       {
 253          gdcmVerboseMacro( "JPEG decompressor failed." );
 254          return false;
 255       }
 256    }
 257
 258    ////////////////////////////////////////////
 259    //// Third stage: twigle the bytes and bits.
 260    ConvertReorderEndianity();
 261    ConvertReArrangeBits();
 262    ConvertHandleColor();
 263
 264    return true;
 265 }
 266
 267 void PixelReadConvert::Squeeze()
 268 {
 269    if ( RGB )
 270       delete [] RGB;
 271    RGB = 0;
 272
 273    if ( Raw )
 274       delete [] Raw;
 275    Raw = 0;
 276
 277    if ( LutRGBA )
 278       delete [] LutRGBA;
 279    LutRGBA = 0;
 280 }
 281
 282 /**
 283  * \brief Build the RGB image from the Raw imagage and the LUTs.
 284  */
 285 bool PixelReadConvert::BuildRGBImage()
 286 {
 287    if ( RGB )
 288    {
 289       // The job is already done.
 290       return true;
 291    }
 292
 293    if ( ! Raw )
 294    {
 295       // The job can't be done
 296       return false;
 297    }
 298
 299    BuildLUTRGBA();
 300    if ( ! LutRGBA )
 301    {
 302       // The job can't be done
 303       return false;
 304    }
 305
 306    // Build RGB Pixels
 307    AllocateRGB();
 308    uint8_t *localRGB = RGB;
 309    for (size_t i = 0; i < RawSize; ++i )
 310    {
 311       int j  = Raw[i] * 4;
 312       *localRGB++ = LutRGBA[j];
 313       *localRGB++ = LutRGBA[j+1];
 314       *localRGB++ = LutRGBA[j+2];
 315    }
 316    return true;
 317 }
 318
 319 //-----------------------------------------------------------------------------
 320 // Protected
 321
 322 //-----------------------------------------------------------------------------
 323 // Private
 324 /**
 325  * \brief Read from file a 12 bits per pixel image and decompress it
 326  *        into a 16 bits per pixel image.
 327  */
 328 void PixelReadConvert::ReadAndDecompress12BitsTo16Bits( std::ifstream *fp )
 329                throw ( FormatError )
 330 {
 331    int nbPixels = XSize * YSize;
 332    uint16_t *localDecompres = (uint16_t*)Raw;
 333
 334    for( int p = 0; p < nbPixels; p += 2 )
 335    {
 336       uint8_t b0, b1, b2;
 337
 338       fp->read( (char*)&b0, 1);
 339       if ( fp->fail() || fp->eof() )
 340       {
 341          throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
 342                                 "Unfound first block" );
 343       }
 344
 345       fp->read( (char*)&b1, 1 );
 346       if ( fp->fail() || fp->eof())
 347       {
 348          throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
 349                                 "Unfound second block" );
 350       }
 351
 352       fp->read( (char*)&b2, 1 );
 353       if ( fp->fail() || fp->eof())
 354       {
 355          throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
 356                                 "Unfound second block" );
 357       }
 358
 359       // Two steps are necessary to please VC++
 360       //
 361       // 2 pixels 12bit =     [0xABCDEF]
 362       // 2 pixels 16bit = [0x0ABD] + [0x0FCE]
 363       //                        A                     B                 D
 364       *localDecompres++ =  ((b0 >> 4) << 8) + ((b0 & 0x0f) << 4) + (b1 & 0x0f);
 365       //                        F                     C                 E
 366       *localDecompres++ =  ((b2 & 0x0f) << 8) + ((b1 >> 4) << 4) + (b2 >> 4);
 367
 368       /// \todo JPR Troubles expected on Big-Endian processors ?
 369    }
 370 }
 371
 372 /**
 373  * \brief     Reads from disk the Pixel Data of JPEG Dicom encapsulated
 374  *            file and decompress it.
 375  * @param     fp File Pointer
 376  * @return    Boolean
 377  */
 378 bool PixelReadConvert::ReadAndDecompressJPEGFile( std::ifstream *fp )
 379 {
 380    if ( IsJPEG2000 )
 381    {
 382       gdcmVerboseMacro( "Sorry, JPEG2000 not yet taken into account" );
 383       fp->seekg( JPEGInfo->GetFirstFragment()->GetOffset(), std::ios::beg);
 384 //    if ( ! gdcm_read_JPEG2000_file( fp,Raw ) )
 385           return false;
 386    }
 387
 388    if ( IsJPEGLS )
 389    {
 390       gdcmVerboseMacro( "Sorry, JPEG-LS not yet taken into account" );
 391       fp->seekg( JPEGInfo->GetFirstFragment()->GetOffset(), std::ios::beg);
 392 //    if ( ! gdcm_read_JPEGLS_file( fp,Raw ) )
 393          return false;
 394    }
 395
 396    // else ??
 397    // Precompute the offset localRaw will be shifted with
 398    int length = XSize * YSize * SamplesPerPixel;
 399    int numberBytes = BitsAllocated / 8;
 400
 401    JPEGInfo->DecompressFromFile(fp, Raw, BitsStored, numberBytes, length );
 402    return true;
 403 }
 404
 405 /**
 406  * \brief Build Red/Green/Blue/Alpha LUT from File
 407  *         when (0028,0004),Photometric Interpretation = [PALETTE COLOR ]
 408  *          and (0028,1101),(0028,1102),(0028,1102)
 409  *            - xxx Palette Color Lookup Table Descriptor - are found
 410  *          and (0028,1201),(0028,1202),(0028,1202)
 411  *            - xxx Palette Color Lookup Table Data - are found
 412  * \warning does NOT deal with :
 413  *   0028 1100 Gray Lookup Table Descriptor (Retired)
 414  *   0028 1221 Segmented Red Palette Color Lookup Table Data
 415  *   0028 1222 Segmented Green Palette Color Lookup Table Data
 416  *   0028 1223 Segmented Blue Palette Color Lookup Table Data
 417  *   no known Dicom reader deals with them :-(
 418  * @return a RGBA Lookup Table
 419  */
 420 void PixelReadConvert::BuildLUTRGBA()
 421 {
 422    if ( LutRGBA )
 423    {
 424       return;
 425    }
 426    // Not so easy : see
 427    // http://www.barre.nom.fr/medical/dicom2/limitations.html#Color%20Lookup%20Tables
 428
 429    if ( ! IsPaletteColor )
 430    {
 431       return;
 432    }
 433
 434    if (   LutRedDescriptor   == GDCM_UNFOUND
 435        || LutGreenDescriptor == GDCM_UNFOUND
 436        || LutBlueDescriptor  == GDCM_UNFOUND )
 437    {
 438       return;
 439    }
 440
 441    ////////////////////////////////////////////
 442    // Extract the info from the LUT descriptors
 443    int lengthR;   // Red LUT length in Bytes
 444    int debR;      // Subscript of the first Lut Value
 445    int nbitsR;    // Lut item size (in Bits)
 446    int nbRead = sscanf( LutRedDescriptor.c_str(),
 447                         "%d\\%d\\%d",
 448                         &lengthR, &debR, &nbitsR );
 449    if( nbRead != 3 )
 450    {
 451       gdcmVerboseMacro( "Wrong Red LUT descriptor" );
 452    }
 453
 454    int lengthG;  // Green LUT length in Bytes
 455    int debG;     // Subscript of the first Lut Value
 456    int nbitsG;   // Lut item size (in Bits)
 457    nbRead = sscanf( LutGreenDescriptor.c_str(),
 458                     "%d\\%d\\%d",
 459                     &lengthG, &debG, &nbitsG );
 460    if( nbRead != 3 )
 461    {
 462       gdcmVerboseMacro( "Wrong Green LUT descriptor" );
 463    }
 464
 465    int lengthB;  // Blue LUT length in Bytes
 466    int debB;     // Subscript of the first Lut Value
 467    int nbitsB;   // Lut item size (in Bits)
 468    nbRead = sscanf( LutRedDescriptor.c_str(),
 469                     "%d\\%d\\%d",
 470                     &lengthB, &debB, &nbitsB );
 471    if( nbRead != 3 )
 472    {
 473       gdcmVerboseMacro( "Wrong Blue LUT descriptor" );
 474    }
 475
 476    ////////////////////////////////////////////////////////
 477    if ( ( ! LutRedData ) || ( ! LutGreenData ) || ( ! LutBlueData ) )
 478    {
 479       return;
 480    }
 481
 482    ////////////////////////////////////////////////
 483    // forge the 4 * 8 Bits Red/Green/Blue/Alpha LUT
 484    LutRGBA = new uint8_t[ 1024 ]; // 256 * 4 (R, G, B, Alpha)
 485    if ( !LutRGBA )
 486       return;
 487
 488    memset( LutRGBA, 0, 1024 );
 489
 490    int mult;
 491    if ( ( nbitsR == 16 ) && ( BitsAllocated == 8 ) )
 492    {
 493       // when LUT item size is different than pixel size
 494       mult = 2; // high byte must be = low byte
 495    }
 496    else
 497    {
 498       // See PS 3.3-2003 C.11.1.1.2 p 619
 499       mult = 1;
 500    }
 501
 502    // if we get a black image, let's just remove the '+1'
 503    // from 'i*mult+1' and check again
 504    // if it works, we shall have to check the 3 Palettes
 505    // to see which byte is ==0 (first one, or second one)
 506    // and fix the code
 507    // We give up the checking to avoid some (useless ?) overhead
 508    // (optimistic asumption)
 509    int i;
 510    uint8_t *a = LutRGBA + 0;
 511    for( i=0; i < lengthR; ++i )
 512    {
 513       *a = LutRedData[i*mult+1];
 514       a += 4;
 515    }
 516
 517    a = LutRGBA + 1;
 518    for( i=0; i < lengthG; ++i)
 519    {
 520       *a = LutGreenData[i*mult+1];
 521       a += 4;
 522    }
 523
 524    a = LutRGBA + 2;
 525    for(i=0; i < lengthB; ++i)
 526    {
 527       *a = LutBlueData[i*mult+1];
 528       a += 4;
 529    }
 530
 531    a = LutRGBA + 3;
 532    for(i=0; i < 256; ++i)
 533    {
 534       *a = 1; // Alpha component
 535       a += 4;
 536    }
 537 }
 538
 539 /**
 540  * \brief Swap the bytes, according to \ref SwapCode.
 541  */
 542 void PixelReadConvert::ConvertSwapZone()
 543 {
 544    unsigned int i;
 545
 546    if( BitsAllocated == 16 )
 547    {
 548       uint16_t *im16 = (uint16_t*)Raw;
 549       switch( SwapCode )
 550       {
 551          case 1234:
 552             break;
 553          case 3412:
 554          case 2143:
 555          case 4321:
 556             for( i = 0; i < RawSize / 2; i++ )
 557             {
 558                im16[i]= (im16[i] >> 8) | (im16[i] << 8 );
 559             }
 560             break;
 561          default:
 562             gdcmVerboseMacro("SwapCode value (16 bits) not allowed.");
 563       }
 564    }
 565    else if( BitsAllocated == 32 )
 566    {
 567       uint32_t s32;
 568       uint16_t high;
 569       uint16_t low;
 570       uint32_t *im32 = (uint32_t*)Raw;
 571       switch ( SwapCode )
 572       {
 573          case 1234:
 574             break;
 575          case 4321:
 576             for( i = 0; i < RawSize / 4; i++ )
 577             {
 578                low     = im32[i] & 0x0000ffff;  // 4321
 579                high    = im32[i] >> 16;
 580                high    = ( high >> 8 ) | ( high << 8 );
 581                low     = ( low  >> 8 ) | ( low  << 8 );
 582                s32     = low;
 583                im32[i] = ( s32 << 16 ) | high;
 584             }
 585             break;
 586          case 2143:
 587             for( i = 0; i < RawSize / 4; i++ )
 588             {
 589                low     = im32[i] & 0x0000ffff;   // 2143
 590                high    = im32[i] >> 16;
 591                high    = ( high >> 8 ) | ( high << 8 );
 592                low     = ( low  >> 8 ) | ( low  << 8 );
 593                s32     = high;
 594                im32[i] = ( s32 << 16 ) | low;
 595             }
 596             break;
 597          case 3412:
 598             for( i = 0; i < RawSize / 4; i++ )
 599             {
 600                low     = im32[i] & 0x0000ffff; // 3412
 601                high    = im32[i] >> 16;
 602                s32     = low;
 603                im32[i] = ( s32 << 16 ) | high;
 604             }
 605             break;
 606          default:
 607             gdcmVerboseMacro("SwapCode value (32 bits) not allowed." );
 608       }
 609    }
 610 }
 611
 612 /**
 613  * \brief Deal with endianness i.e. re-arange bytes inside the integer
 614  */
 615 void PixelReadConvert::ConvertReorderEndianity()
 616 {
 617    if ( BitsAllocated != 8 )
 618    {
 619       ConvertSwapZone();
 620    }
 621
 622    // Special kludge in order to deal with xmedcon broken images:
 623    if ( BitsAllocated == 16
 624      && BitsStored < BitsAllocated
 625      && !PixelSign )
 626    {
 627       int l = (int)( RawSize / ( BitsAllocated / 8 ) );
 628       uint16_t *deb = (uint16_t *)Raw;
 629       for(int i = 0; i<l; i++)
 630       {
 631          if( *deb == 0xffff )
 632          {
 633            *deb = 0;
 634          }
 635          deb++;
 636       }
 637    }
 638 }
 639
 640 /**
 641  * \brief  Re-arrange the bits within the bytes.
 642  * @return Boolean
 643  */
 644 bool PixelReadConvert::ConvertReArrangeBits() throw ( FormatError )
 645 {
 646    if ( BitsStored != BitsAllocated )
 647    {
 648       int l = (int)( RawSize / ( BitsAllocated / 8 ) );
 649       if ( BitsAllocated == 16 )
 650       {
 651          uint16_t mask = 0xffff;
 652          mask = mask >> ( BitsAllocated - BitsStored );
 653          uint16_t *deb = (uint16_t*)Raw;
 654          for(int i = 0; i<l; i++)
 655          {
 656             *deb = (*deb >> (BitsStored - HighBitPosition - 1)) & mask;
 657             deb++;
 658          }
 659       }
 660       else if ( BitsAllocated == 32 )
 661       {
 662          uint32_t mask = 0xffffffff;
 663          mask = mask >> ( BitsAllocated - BitsStored );
 664          uint32_t *deb = (uint32_t*)Raw;
 665          for(int i = 0; i<l; i++)
 666          {
 667             *deb = (*deb >> (BitsStored - HighBitPosition - 1)) & mask;
 668             deb++;
 669          }
 670       }
 671       else
 672       {
 673          gdcmVerboseMacro("Weird image");
 674          throw FormatError( "Weird image !?" );
 675       }
 676    }
 677    return true;
 678 }
 679
 680 /**
 681  * \brief   Convert (Red plane, Green plane, Blue plane) to RGB pixels
 682  * \warning Works on all the frames at a time
 683  */
 684 void PixelReadConvert::ConvertRGBPlanesToRGBPixels()
 685 {
 686    uint8_t *localRaw = Raw;
 687    uint8_t *copyRaw = new uint8_t[ RawSize ];
 688    memmove( copyRaw, localRaw, RawSize );
 689
 690    int l = XSize * YSize * ZSize;
 691
 692    uint8_t *a = copyRaw;
 693    uint8_t *b = copyRaw + l;
 694    uint8_t *c = copyRaw + l + l;
 695
 696    for (int j = 0; j < l; j++)
 697    {
 698       *(localRaw++) = *(a++);
 699       *(localRaw++) = *(b++);
 700       *(localRaw++) = *(c++);
 701    }
 702    delete[] copyRaw;
 703 }
 704
 705 /**
 706  * \brief   Convert (cY plane, cB plane, cR plane) to RGB pixels
 707  * \warning Works on all the frames at a time
 708  */
 709 void PixelReadConvert::ConvertYcBcRPlanesToRGBPixels()
 710 {
 711    uint8_t *localRaw = Raw;
 712    uint8_t *copyRaw = new uint8_t[ RawSize ];
 713    memmove( copyRaw, localRaw, RawSize );
 714
 715    // to see the tricks about YBR_FULL, YBR_FULL_422,
 716    // YBR_PARTIAL_422, YBR_ICT, YBR_RCT have a look at :
 717    // ftp://medical.nema.org/medical/dicom/final/sup61_ft.pdf
 718    // and be *very* affraid
 719    //
 720    int l        = XSize * YSize;
 721    int nbFrames = ZSize;
 722
 723    uint8_t *a = copyRaw;
 724    uint8_t *b = copyRaw + l;
 725    uint8_t *c = copyRaw + l + l;
 726    double R, G, B;
 727
 728    /// \todo : Replace by the 'well known' integer computation
 729    ///         counterpart. Refer to
 730    ///            http://lestourtereaux.free.fr/papers/data/yuvrgb.pdf
 731    ///         for code optimisation.
 732
 733    for ( int i = 0; i < nbFrames; i++ )
 734    {
 735       for ( int j = 0; j < l; j++ )
 736       {
 737          R = 1.164 *(*a-16) + 1.596 *(*c -128) + 0.5;
 738          G = 1.164 *(*a-16) - 0.813 *(*c -128) - 0.392 *(*b -128) + 0.5;
 739          B = 1.164 *(*a-16) + 2.017 *(*b -128) + 0.5;
 740
 741          if (R < 0.0)   R = 0.0;
 742          if (G < 0.0)   G = 0.0;
 743          if (B < 0.0)   B = 0.0;
 744          if (R > 255.0) R = 255.0;
 745          if (G > 255.0) G = 255.0;
 746          if (B > 255.0) B = 255.0;
 747
 748          *(localRaw++) = (uint8_t)R;
 749          *(localRaw++) = (uint8_t)G;
 750          *(localRaw++) = (uint8_t)B;
 751          a++;
 752          b++;
 753          c++;
 754       }
 755    }
 756    delete[] copyRaw;
 757 }
 758
 759 void PixelReadConvert::ConvertHandleColor()
 760 {
 761    //////////////////////////////////
 762    // Deal with the color decoding i.e. handle:
 763    //   - R, G, B planes (as opposed to RGB pixels)
 764    //   - YBR (various) encodings.
 765    //   - LUT[s] (or "PALETTE COLOR").
 766    //
 767    // The classification in the color decoding schema is based on the blending
 768    // of two Dicom tags values:
 769    // * "Photometric Interpretation" for which we have the cases:
 770    //  - [Photo A] MONOCHROME[1|2] pictures,
 771    //  - [Photo B] RGB or YBR_FULL_422 (which acts as RGB),
 772    //  - [Photo C] YBR_* (with the above exception of YBR_FULL_422)
 773    //  - [Photo D] "PALETTE COLOR" which indicates the presence of LUT[s].
 774    // * "Planar Configuration" for which we have the cases:
 775    //  - [Planar 0] 0 then Pixels are already RGB
 776    //  - [Planar 1] 1 then we have 3 planes : R, G, B,
 777    //  - [Planar 2] 2 then we have 1 gray Plane and 3 LUTs
 778    //
 779    // Now in theory, one could expect some coherence when blending the above
 780    // cases. For example we should not encounter files belonging at the
 781    // time to case [Planar 0] and case [Photo D].
 782    // Alas, this was only theory ! Because in practice some odd (read ill
 783    // formated Dicom) files (e.g. gdcmData/US-PAL-8-10x-echo.dcm) we encounter:
 784    //     - "Planar Configuration" = 0,
 785    //     - "Photometric Interpretation" = "PALETTE COLOR".
 786    // Hence gdcm will use the folowing "heuristic" in order to be tolerant
 787    // towards Dicom-non-conformance files:
 788    //   << whatever the "Planar Configuration" value might be, a
 789    //      "Photometric Interpretation" set to "PALETTE COLOR" forces
 790    //      a LUT intervention >>
 791    //
 792    // Now we are left with the following handling of the cases:
 793    // - [Planar 0] OR  [Photo A] no color decoding (since respectively
 794    //       Pixels are already RGB and monochrome pictures have no color :),
 795    // - [Planar 1] AND [Photo B] handled with ConvertRGBPlanesToRGBPixels()
 796    // - [Planar 1] AND [Photo C] handled with ConvertYcBcRPlanesToRGBPixels()
 797    // - [Planar 2] OR  [Photo D] requires LUT intervention.
 798
 799    if ( ! IsRawRGB() )
 800    {
 801       // [Planar 2] OR  [Photo D]: LUT intervention done outside
 802       return;
 803    }
 804
 805    if ( PlanarConfiguration == 1 )
 806    {
 807       if ( IsYBRFull )
 808       {
 809          // [Planar 1] AND [Photo C] (remember YBR_FULL_422 acts as RGB)
 810          ConvertYcBcRPlanesToRGBPixels();
 811       }
 812       else
 813       {
 814          // [Planar 1] AND [Photo C]
 815          ConvertRGBPlanesToRGBPixels();
 816       }
 817       return;
 818    }
 819
 820    // When planarConf is 0, and RLELossless (forbidden by Dicom norm)
 821    // pixels need to be RGB-fied anyway
 822    if (IsRLELossless)
 823    {
 824      ConvertRGBPlanesToRGBPixels();
 825    }
 826    // In *normal *case, when planarConf is 0, pixels are already in RGB
 827 }
 828
 829 void PixelReadConvert::ComputeRawAndRGBSizes()
 830 {
 831    int bitsAllocated = BitsAllocated;
 832    // Number of "Bits Allocated" is fixed to 16 when it's 12, since
 833    // in this case we will expand the image to 16 bits (see
 834    //    \ref ReadAndDecompress12BitsTo16Bits() )
 835    if (  BitsAllocated == 12 )
 836    {
 837       bitsAllocated = 16;
 838    }
 839
 840    RawSize =  XSize * YSize * ZSize
 841                      * ( bitsAllocated / 8 )
 842                      * SamplesPerPixel;
 843    if ( HasLUT )
 844    {
 845       RGBSize = 3 * RawSize;
 846    }
 847    else
 848    {
 849       RGBSize = RawSize;
 850    }
 851 }
 852
 853 void PixelReadConvert::AllocateRGB()
 854 {
 855   if ( RGB )
 856      delete [] RGB;
 857   RGB = new uint8_t[RGBSize];
 858 }
 859
 860 void PixelReadConvert::AllocateRaw()
 861 {
 862   if ( Raw )
 863      delete [] Raw;
 864   Raw = new uint8_t[RawSize];
 865 }
 866
 867 //-----------------------------------------------------------------------------
 868 // Print
 869 /**
 870  * \brief        Print self.
 871  * @param indent Indentation string to be prepended during printing.
 872  * @param os     Stream to print to.
 873  */
 874 void PixelReadConvert::Print( std::ostream &os, std::string const &indent )
 875 {
 876    os << indent
 877       << "--- Pixel information -------------------------"
 878       << std::endl;
 879    os << indent
 880       << "Pixel Data: offset " << PixelOffset
 881       << " x(" << std::hex << PixelOffset << std::dec
 882       << ")   length " << PixelDataLength
 883       << " x(" << std::hex << PixelDataLength << std::dec
 884       << ")" << std::endl;
 885
 886    if ( IsRLELossless )
 887    {
 888       if ( RLEInfo )
 889       {
 890          RLEInfo->Print( os, indent );
 891       }
 892       else
 893       {
 894          gdcmVerboseMacro("Set as RLE file but NO RLEinfo present.");
 895       }
 896    }
 897
 898    if ( IsJPEG2000 || IsJPEGLossless || IsJPEGLossy || IsJPEGLS )
 899    {
 900       if ( JPEGInfo )
 901       {
 902          JPEGInfo->Print( os, indent );
 903       }
 904       else
 905       {
 906          gdcmVerboseMacro("Set as JPEG file but NO JPEGinfo present.");
 907       }
 908    }
 909 }
 910
 911 //-----------------------------------------------------------------------------
 912 } // end namespace gdcm
 913
 914 // NOTES on File internal calls
 915 // User
 916 // ---> GetImageData
 917 //     ---> GetImageDataIntoVector
 918 //        |---> GetImageDataIntoVectorRaw
 919 //        | lut intervention
 920 // User
 921 // ---> GetImageDataRaw
 922 //     ---> GetImageDataIntoVectorRaw
 923