src/gdcmPixelReadConvert.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: gdcmPixelReadConvert.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2005/10/26 06:57:02 $
   7   Version:   $Revision: 1.89 $
   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 "gdcmPixelReadConvert.h"
  20 #include "gdcmDebug.h"
  21 #include "gdcmFile.h"
  22 #include "gdcmGlobal.h"
  23 #include "gdcmTS.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 //bool ReadMPEGFile (std::ifstream *fp, void *image_buffer, size_t lenght);
  35 bool gdcm_read_JPEG2000_file (void* raw,
  36                               char *inputdata, size_t inputlength);
  37 //-----------------------------------------------------------------------------
  38 #define str2num(str, typeNum) *((typeNum *)(str))
  39
  40 //-----------------------------------------------------------------------------
  41 // Constructor / Destructor
  42 /// Constructor
  43 PixelReadConvert::PixelReadConvert()
  44 {
  45    RGB          = 0;
  46    RGBSize      = 0;
  47    Raw          = 0;
  48    RawSize      = 0;
  49    LutRGBA      = 0;
  50    LutRedData   = 0;
  51    LutGreenData = 0;
  52    LutBlueData  = 0;
  53    RLEInfo      = 0;
  54    JPEGInfo     = 0;
  55    UserFunction = 0;
  56    FileInternal = 0;
  57 }
  58
  59 /// Canonical Destructor
  60 PixelReadConvert::~PixelReadConvert()
  61 {
  62    Squeeze();
  63 }
  64
  65 //-----------------------------------------------------------------------------
  66 // Public
  67 /**
  68  * \brief Predicate to know whether the image[s] (once Raw) is RGB.
  69  * \note See comments of \ref ConvertHandleColor
  70  */
  71 bool PixelReadConvert::IsRawRGB()
  72 {
  73    if (   IsMonochrome
  74        || PlanarConfiguration == 2
  75        || IsPaletteColor )
  76    {
  77       return false;
  78    }
  79    return true;
  80 }
  81 /**
  82  * \brief Gets various usefull informations from the file header
  83  * @param file gdcm::File pointer
  84  */
  85 void PixelReadConvert::GrabInformationsFromFile( File *file )
  86 {
  87    // Number of Bits Allocated for storing a Pixel is defaulted to 16
  88    // when absent from the file.
  89    BitsAllocated = file->GetBitsAllocated();
  90    if ( BitsAllocated == 0 )
  91    {
  92       BitsAllocated = 16;
  93    }
  94
  95    // Number of "Bits Stored", defaulted to number of "Bits Allocated"
  96    // when absent from the file.
  97    BitsStored = file->GetBitsStored();
  98    if ( BitsStored == 0 )
  99    {
 100       BitsStored = BitsAllocated;
 101    }
 102
 103    // High Bit Position, defaulted to "Bits Allocated" - 1
 104    HighBitPosition = file->GetHighBitPosition();
 105    if ( HighBitPosition == 0 )
 106    {
 107       HighBitPosition = BitsAllocated - 1;
 108    }
 109
 110    XSize           = file->GetXSize();
 111    YSize           = file->GetYSize();
 112    ZSize           = file->GetZSize();
 113    SamplesPerPixel = file->GetSamplesPerPixel();
 114    //PixelSize       = file->GetPixelSize();  Useless
 115    PixelSign       = file->IsSignedPixelData();
 116    SwapCode        = file->GetSwapCode();
 117    std::string ts  = file->GetTransferSyntax();
 118
 119
 120 //   if ( ts == GDCM_UNKNOWN )
 121 //     {
 122 //     gdcmErrorMacro( "Could someone tell me how in the world could this happen !" );
 123 // -->
 124 //--> on ALL acr-nema images ! JPRx
 125 //-->
 126 //     abort(); // DO NOT REMOVE.  WE SHOULD NEVER READ SUCH IMAGE EVER (only gdcm can write such broekn dicom file)
 127 //     }
 128
 129    IsRaw =
 130         ( ! file->IsDicomV3() )  // Should be ACR-NEMA file
 131      || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ImplicitVRLittleEndian
 132      || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ImplicitVRBigEndianPrivateGE
 133      || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ExplicitVRLittleEndian
 134      || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ExplicitVRBigEndian
 135      || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::DeflatedExplicitVRLittleEndian;
 136
 137    IsPrivateGETransferSyntax =
 138                 ( Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ImplicitVRBigEndianPrivateGE );
 139
 140    IsMPEG          = Global::GetTS()->IsMPEG(ts);
 141    IsJPEG2000      = Global::GetTS()->IsJPEG2000(ts);
 142    IsJPEGLS        = Global::GetTS()->IsJPEGLS(ts);
 143    IsJPEGLossy     = Global::GetTS()->IsJPEGLossy(ts);
 144    IsJPEGLossless  = Global::GetTS()->IsJPEGLossless(ts);
 145    IsRLELossless   = Global::GetTS()->IsRLELossless(ts);
 146
 147    PixelOffset     = file->GetPixelOffset();
 148    PixelDataLength = file->GetPixelAreaLength();
 149    RLEInfo         = file->GetRLEInfo();
 150    JPEGInfo        = file->GetJPEGInfo();
 151
 152    IsMonochrome    = file->IsMonochrome();
 153    IsMonochrome1   = file->IsMonochrome1();
 154    IsPaletteColor  = file->IsPaletteColor();
 155    IsYBRFull       = file->IsYBRFull();
 156
 157    PlanarConfiguration = file->GetPlanarConfiguration();
 158
 159    /////////////////////////////////////////////////////////////////
 160    // LUT section:
 161    HasLUT = file->HasLUT();
 162    if ( HasLUT )
 163    {
 164       // Just in case some access to a File element requires disk access.
 165       LutRedDescriptor   = file->GetEntryString( 0x0028, 0x1101 );
 166       LutGreenDescriptor = file->GetEntryString( 0x0028, 0x1102 );
 167       LutBlueDescriptor  = file->GetEntryString( 0x0028, 0x1103 );
 168
 169       // The following comment is probabely meaningless, since LUT are *always*
 170       // loaded at parsing time, whatever their length is.
 171
 172       // Depending on the value of Document::MAX_SIZE_LOAD_ELEMENT_VALUE
 173       // [ refer to invocation of Document::SetMaxSizeLoadEntry() in
 174       // Document::Document() ], the loading of the value (content) of a
 175       // [Bin|Val]Entry occurence migth have been hindered (read simply NOT
 176       // loaded). Hence, we first try to obtain the LUTs data from the file
 177       // and when this fails we read the LUTs data directly from disk.
 178       // \TODO Reading a [Bin|Val]Entry directly from disk is a kludge.
 179       //       We should NOT bypass the [Bin|Val]Entry class. Instead
 180       //       an access to an UNLOADED content of a [Bin|Val]Entry occurence
 181       //       (e.g. DataEntry::GetBinArea()) should force disk access from
 182       //       within the [Bin|Val]Entry class itself. The only problem
 183       //       is that the [Bin|Val]Entry is unaware of the FILE* is was
 184       //       parsed from. Fix that. FIXME.
 185
 186       // //// Red round
 187       file->LoadEntryBinArea(0x0028, 0x1201);
 188       LutRedData = (uint8_t*)file->GetEntryBinArea( 0x0028, 0x1201 );
 189       if ( ! LutRedData )
 190       {
 191          gdcmWarningMacro("Unable to read Red Palette Color Lookup Table data");
 192       }
 193
 194       // //// Green round:
 195       file->LoadEntryBinArea(0x0028, 0x1202);
 196       LutGreenData = (uint8_t*)file->GetEntryBinArea(0x0028, 0x1202 );
 197       if ( ! LutGreenData)
 198       {
 199          gdcmWarningMacro("Unable to read Green Palette Color Lookup Table data");
 200       }
 201
 202       // //// Blue round:
 203       file->LoadEntryBinArea(0x0028, 0x1203);
 204       LutBlueData = (uint8_t*)file->GetEntryBinArea( 0x0028, 0x1203 );
 205       if ( ! LutBlueData )
 206       {
 207          gdcmWarningMacro("Unable to read Blue Palette Color Lookup Table data");
 208       }
 209    }
 210    FileInternal = file;
 211
 212    ComputeRawAndRGBSizes();
 213 }
 214
 215 /// \brief Reads from disk and decompresses Pixels
 216 bool PixelReadConvert::ReadAndDecompressPixelData( std::ifstream *fp )
 217 {
 218    // ComputeRawAndRGBSizes is already made by
 219    // ::GrabInformationsFromfile. So, the structure sizes are
 220    // correct
 221    Squeeze();
 222
 223    //////////////////////////////////////////////////
 224    //// First stage: get our hands on the Pixel Data.
 225    if ( !fp )
 226    {
 227       gdcmWarningMacro( "Unavailable file pointer." );
 228       return false;
 229    }
 230
 231    fp->seekg( PixelOffset, std::ios::beg );
 232    if ( fp->fail() || fp->eof() )
 233    {
 234       gdcmWarningMacro( "Unable to find PixelOffset in file." );
 235       return false;
 236    }
 237
 238    AllocateRaw();
 239
 240    //////////////////////////////////////////////////
 241    //// Second stage: read from disk and decompress.
 242    if ( BitsAllocated == 12 )
 243    {
 244       ReadAndDecompress12BitsTo16Bits( fp);
 245    }
 246    else if ( IsRaw )
 247    {
 248       // This problem can be found when some obvious informations are found
 249       // after the field containing the image data. In this case, these
 250       // bad data are added to the size of the image (in the PixelDataLength
 251       // variable). But RawSize is the right size of the image !
 252       if ( PixelDataLength != RawSize )
 253       {
 254          gdcmWarningMacro( "Mismatch between PixelReadConvert : "
 255                           << PixelDataLength << " and RawSize : " << RawSize );
 256       }
 257       if ( PixelDataLength > RawSize )
 258       {
 259          fp->read( (char*)Raw, RawSize);
 260       }
 261       else
 262       {
 263          fp->read( (char*)Raw, PixelDataLength);
 264       }
 265
 266       if ( fp->fail() || fp->eof())
 267       {
 268          gdcmWarningMacro( "Reading of Raw pixel data failed." );
 269          return false;
 270       }
 271    }
 272    else if ( IsRLELossless )
 273    {
 274       if ( ! RLEInfo->DecompressRLEFile
 275                                ( fp, Raw, XSize, YSize, ZSize, BitsAllocated ) )
 276       {
 277          gdcmWarningMacro( "RLE decompressor failed." );
 278          return false;
 279       }
 280    }
 281    else if ( IsMPEG )
 282    {
 283       //gdcmWarningMacro( "Sorry, MPEG not yet taken into account" );
 284       //return false;
 285       // fp has already been seek to start of mpeg
 286       //ReadMPEGFile(fp, Raw, PixelDataLength);
 287       return true;
 288    }
 289    else
 290    {
 291       // Default case concerns JPEG family
 292       if ( ! ReadAndDecompressJPEGFile( fp ) )
 293       {
 294          gdcmWarningMacro( "JPEG decompressor failed." );
 295          return false;
 296       }
 297    }
 298
 299    ////////////////////////////////////////////
 300    //// Third stage: twigle the bytes and bits.
 301    ConvertReorderEndianity();
 302    ConvertReArrangeBits();
 303    ConvertFixGreyLevels();
 304    if (UserFunction) // user is allowed to Mirror, TopDown, Rotate,...the image
 305       UserFunction( Raw, FileInternal);
 306    ConvertHandleColor();
 307
 308    return true;
 309 }
 310
 311 /// Deletes Pixels Area
 312 void PixelReadConvert::Squeeze()
 313 {
 314    if ( RGB )
 315       delete [] RGB;
 316    RGB = 0;
 317
 318    if ( Raw )
 319       delete [] Raw;
 320    Raw = 0;
 321
 322    if ( LutRGBA )
 323       delete [] LutRGBA;
 324    LutRGBA = 0;
 325 }
 326
 327 /**
 328  * \brief Build the RGB image from the Raw image and the LUTs.
 329  */
 330 bool PixelReadConvert::BuildRGBImage()
 331 {
 332    if ( RGB )
 333    {
 334       // The job is already done.
 335       return true;
 336    }
 337
 338    if ( ! Raw )
 339    {
 340       // The job can't be done
 341       return false;
 342    }
 343
 344    BuildLUTRGBA();
 345    if ( ! LutRGBA )
 346    {
 347       // The job can't be done
 348       return false;
 349    }
 350
 351    gdcmWarningMacro( "--> BuildRGBImage" );
 352
 353    // Build RGB Pixels
 354    AllocateRGB();
 355
 356    int j;
 357    if ( BitsAllocated <= 8 )
 358    {
 359       uint8_t *localRGB = RGB;
 360       for (size_t i = 0; i < RawSize; ++i )
 361       {
 362          j  = Raw[i] * 4;
 363          *localRGB++ = LutRGBA[j];
 364          *localRGB++ = LutRGBA[j+1];
 365          *localRGB++ = LutRGBA[j+2];
 366       }
 367     }
 368
 369     else  // deal with 16 bits pixels and 16 bits Palette color
 370     {
 371       uint16_t *localRGB = (uint16_t *)RGB;
 372       for (size_t i = 0; i < RawSize/2; ++i )
 373       {
 374          j  = ((uint16_t *)Raw)[i] * 4;
 375          *localRGB++ = ((uint16_t *)LutRGBA)[j];
 376          *localRGB++ = ((uint16_t *)LutRGBA)[j+1];
 377          *localRGB++ = ((uint16_t *)LutRGBA)[j+2];
 378       }
 379     }
 380
 381    return true;
 382 }
 383
 384 //-----------------------------------------------------------------------------
 385 // Protected
 386
 387 //-----------------------------------------------------------------------------
 388 // Private
 389 /**
 390  * \brief Read from file a 12 bits per pixel image and decompress it
 391  *        into a 16 bits per pixel image.
 392  */
 393 void PixelReadConvert::ReadAndDecompress12BitsTo16Bits( std::ifstream *fp )
 394                throw ( FormatError )
 395 {
 396    int nbPixels = XSize * YSize;
 397    uint16_t *localDecompres = (uint16_t*)Raw;
 398
 399    for( int p = 0; p < nbPixels; p += 2 )
 400    {
 401       uint8_t b0, b1, b2;
 402
 403       fp->read( (char*)&b0, 1);
 404       if ( fp->fail() || fp->eof() )
 405       {
 406          throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
 407                                 "Unfound first block" );
 408       }
 409
 410       fp->read( (char*)&b1, 1 );
 411       if ( fp->fail() || fp->eof())
 412       {
 413          throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
 414                                 "Unfound second block" );
 415       }
 416
 417       fp->read( (char*)&b2, 1 );
 418       if ( fp->fail() || fp->eof())
 419       {
 420          throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
 421                                 "Unfound second block" );
 422       }
 423
 424       // Two steps are necessary to please VC++
 425       //
 426       // 2 pixels 12bit =     [0xABCDEF]
 427       // 2 pixels 16bit = [0x0ABD] + [0x0FCE]
 428       //                        A                     B                 D
 429       *localDecompres++ =  ((b0 >> 4) << 8) + ((b0 & 0x0f) << 4) + (b1 & 0x0f);
 430       //                        F                     C                 E
 431       *localDecompres++ =  ((b2 & 0x0f) << 8) + ((b1 >> 4) << 4) + (b2 >> 4);
 432
 433       /// \todo JPR Troubles expected on Big-Endian processors ?
 434    }
 435 }
 436
 437 /**
 438  * \brief     Reads from disk the Pixel Data of JPEG Dicom encapsulated
 439  *            file and decompress it.
 440  * @param     fp File Pointer
 441  * @return    Boolean
 442  */
 443 bool PixelReadConvert::ReadAndDecompressJPEGFile( std::ifstream *fp )
 444 {
 445    if ( IsJPEG2000 )
 446    {
 447      // make sure this is the right JPEG compression
 448      assert( !IsJPEGLossless || !IsJPEGLossy || !IsJPEGLS );
 449      // FIXME this is really ugly but it seems I have to load the complete
 450      // jpeg2000 stream to use jasper:
 451      // I don't think we'll ever be able to deal with multiple fragments properly
 452
 453       unsigned long inputlength = 0;
 454       JPEGFragment *jpegfrag = JPEGInfo->GetFirstFragment();
 455       while( jpegfrag )
 456       {
 457          inputlength += jpegfrag->GetLength();
 458          jpegfrag = JPEGInfo->GetNextFragment();
 459       }
 460       gdcmAssertMacro( inputlength != 0);
 461       uint8_t *inputdata = new uint8_t[inputlength];
 462       char *pinputdata = (char*)inputdata;
 463       jpegfrag = JPEGInfo->GetFirstFragment();
 464       while( jpegfrag )
 465       {
 466          fp->seekg( jpegfrag->GetOffset(), std::ios::beg);
 467          fp->read(pinputdata, jpegfrag->GetLength());
 468          pinputdata += jpegfrag->GetLength();
 469          jpegfrag = JPEGInfo->GetNextFragment();
 470       }
 471       // Warning the inputdata buffer is delete in the function
 472       if ( ! gdcm_read_JPEG2000_file( Raw,
 473           (char*)inputdata, inputlength ) )
 474       {
 475          return true;
 476       }
 477       // wow what happen, must be an error
 478       return false;
 479    }
 480    else if ( IsJPEGLS )
 481    {
 482      // make sure this is the right JPEG compression
 483      assert( !IsJPEGLossless || !IsJPEGLossy || !IsJPEG2000 );
 484    // WARNING : JPEG-LS is NOT the 'classical' Jpeg Lossless :
 485    // [JPEG-LS is the basis for new lossless/near-lossless compression
 486    // standard for continuous-tone images intended for JPEG2000. The standard
 487    // is based on the LOCO-I algorithm (LOw COmplexity LOssless COmpression
 488    // for Images) developed at Hewlett-Packard Laboratories]
 489    //
 490    // see http://datacompression.info/JPEGLS.shtml
 491    //
 492 #if 0
 493    std::cerr << "count:" << JPEGInfo->GetFragmentCount() << std::endl;
 494       unsigned long inputlength = 0;
 495       JPEGFragment *jpegfrag = JPEGInfo->GetFirstFragment();
 496       while( jpegfrag )
 497       {
 498          inputlength += jpegfrag->GetLength();
 499          jpegfrag = JPEGInfo->GetNextFragment();
 500       }
 501       gdcmAssertMacro( inputlength != 0);
 502       uint8_t *inputdata = new uint8_t[inputlength];
 503       char *pinputdata = (char*)inputdata;
 504       jpegfrag = JPEGInfo->GetFirstFragment();
 505       while( jpegfrag )
 506       {
 507          fp->seekg( jpegfrag->GetOffset(), std::ios::beg);
 508          fp->read(pinputdata, jpegfrag->GetLength());
 509          pinputdata += jpegfrag->GetLength();
 510          jpegfrag = JPEGInfo->GetNextFragment();
 511       }
 512
 513   //fp->read((char*)Raw, PixelDataLength);
 514
 515   std::ofstream out("/tmp/jpegls.jpg");
 516   out.write((char*)inputdata, inputlength);
 517   out.close();
 518   delete[] inputdata;
 519 #endif
 520
 521       gdcmWarningMacro( "Sorry, JPEG-LS not yet taken into account" );
 522       fp->seekg( JPEGInfo->GetFirstFragment()->GetOffset(), std::ios::beg);
 523 //    if ( ! gdcm_read_JPEGLS_file( fp,Raw ) )
 524          return false;
 525    }
 526    else
 527    {
 528      // make sure this is the right JPEG compression
 529      assert( !IsJPEGLS || !IsJPEG2000 );
 530      // Precompute the offset localRaw will be shifted with
 531      int length = XSize * YSize * SamplesPerPixel;
 532      int numberBytes = BitsAllocated / 8;
 533
 534      JPEGInfo->DecompressFromFile(fp, Raw, BitsStored, numberBytes, length );
 535      return true;
 536    }
 537 }
 538
 539 /**
 540  * \brief Build Red/Green/Blue/Alpha LUT from File when :
 541  *         - (0028,0004) : Photometric Interpretation == [PALETTE COLOR ]
 542  *         and
 543  *         - (0028,1101),(0028,1102),(0028,1102)
 544  *            xxx Palette Color Lookup Table Descriptor are found
 545  *          and
 546  *         - (0028,1201),(0028,1202),(0028,1202)
 547  *           xxx Palette Color Lookup Table Data - are found
 548  * \warning does NOT deal with :
 549  *   - 0028 1100 Gray Lookup Table Descriptor (Retired)
 550  *   - 0028 1221 Segmented Red Palette Color Lookup Table Data
 551  *   - 0028 1222 Segmented Green Palette Color Lookup Table Data
 552  *   - 0028 1223 Segmented Blue Palette Color Lookup Table Data
 553  *   no known Dicom reader deals with them :-(
 554  * @return a RGBA Lookup Table
 555  */
 556 void PixelReadConvert::BuildLUTRGBA()
 557 {
 558
 559    // Note to code reviewers :
 560    // The problem is *much more* complicated, since a lot of manufacturers
 561    // Don't follow the norm :
 562    // have a look at David Clunie's remark at the end of this .cxx file.
 563    if ( LutRGBA )
 564
 565    {
 566       return;
 567    }
 568    // Not so easy : see
 569    // http://www.barre.nom.fr/medical/dicom2/limitations.html#Color%20Lookup%20Tables
 570
 571    if ( ! IsPaletteColor )
 572    {
 573       return;
 574    }
 575
 576    if (   LutRedDescriptor   == GDCM_UNFOUND
 577        || LutGreenDescriptor == GDCM_UNFOUND
 578        || LutBlueDescriptor  == GDCM_UNFOUND )
 579    {
 580       gdcmWarningMacro( "(At least) a LUT Descriptor is missing" );
 581       return;
 582    }
 583
 584    ////////////////////////////////////////////
 585    // Extract the info from the LUT descriptors
 586    int lengthR;   // Red LUT length in Bytes
 587    int debR;      // Subscript of the first Lut Value
 588    int nbitsR;    // Lut item size (in Bits)
 589    int nbRead;    // nb of items in LUT descriptor (must be = 3)
 590
 591    nbRead = sscanf( LutRedDescriptor.c_str(),
 592                         "%d\\%d\\%d",
 593                         &lengthR, &debR, &nbitsR );
 594    if ( nbRead != 3 )
 595    {
 596       gdcmWarningMacro( "Wrong Red LUT descriptor" );
 597    }
 598    int lengthG;  // Green LUT length in Bytes
 599    int debG;     // Subscript of the first Lut Value
 600    int nbitsG;   // Lut item size (in Bits)
 601
 602    nbRead = sscanf( LutGreenDescriptor.c_str(),
 603                     "%d\\%d\\%d",
 604                     &lengthG, &debG, &nbitsG );
 605    if ( nbRead != 3 )
 606    {
 607       gdcmWarningMacro( "Wrong Green LUT descriptor" );
 608    }
 609
 610    int lengthB;  // Blue LUT length in Bytes
 611    int debB;     // Subscript of the first Lut Value
 612    int nbitsB;   // Lut item size (in Bits)
 613    nbRead = sscanf( LutRedDescriptor.c_str(),
 614                     "%d\\%d\\%d",
 615                     &lengthB, &debB, &nbitsB );
 616    if ( nbRead != 3 )
 617    {
 618       gdcmWarningMacro( "Wrong Blue LUT descriptor" );
 619    }
 620
 621    gdcmWarningMacro(" lengthR " << lengthR << " debR "
 622                  << debR << " nbitsR " << nbitsR);
 623    gdcmWarningMacro(" lengthG " << lengthG << " debG "
 624                  << debG << " nbitsG " << nbitsG);
 625    gdcmWarningMacro(" lengthB " << lengthB << " debB "
 626                  << debB << " nbitsB " << nbitsB);
 627
 628    if ( !lengthR ) // if = 2^16, this shall be 0 see : CP-143
 629       lengthR=65536;
 630    if ( !lengthG ) // if = 2^16, this shall be 0
 631       lengthG=65536;
 632    if ( !lengthB ) // if = 2^16, this shall be 0
 633       lengthB=65536;
 634
 635    ////////////////////////////////////////////////////////
 636
 637    if ( ( ! LutRedData ) || ( ! LutGreenData ) || ( ! LutBlueData ) )
 638    {
 639       gdcmWarningMacro( "(At least) a LUT is missing" );
 640       return;
 641    }
 642
 643    // -------------------------------------------------------------
 644
 645    if ( BitsAllocated <= 8 )
 646    {
 647       // forge the 4 * 8 Bits Red/Green/Blue/Alpha LUT
 648       LutRGBA = new uint8_t[ 1024 ]; // 256 * 4 (R, G, B, Alpha)
 649       if ( !LutRGBA )
 650          return;
 651       LutItemNumber = 256;
 652       LutItemSize   = 8;
 653       memset( LutRGBA, 0, 1024 );
 654
 655       int mult;
 656       if ( ( nbitsR == 16 ) && ( BitsAllocated == 8 ) )
 657       {
 658          // when LUT item size is different than pixel size
 659          mult = 2; // high byte must be = low byte
 660       }
 661       else
 662       {
 663          // See PS 3.3-2003 C.11.1.1.2 p 619
 664          mult = 1;
 665       }
 666
 667       // if we get a black image, let's just remove the '+1'
 668       // from 'i*mult+1' and check again
 669       // if it works, we shall have to check the 3 Palettes
 670       // to see which byte is ==0 (first one, or second one)
 671       // and fix the code
 672       // We give up the checking to avoid some (useless ?) overhead
 673       // (optimistic asumption)
 674       int i;
 675       uint8_t *a;
 676
 677       //take "Subscript of the first Lut Value" (debR,debG,debB) into account!
 678
 679       //FIXME :  +1 : to get 'low value' byte
 680       //         Trouble expected on Big Endian Processors ?
 681       //         16 BIts Per Pixel Palette Color to be swapped?
 682
 683       a = LutRGBA + 0 + debR;
 684       for( i=0; i < lengthR; ++i )
 685       {
 686          *a = LutRedData[i*mult+1];
 687          a += 4;
 688       }
 689
 690       a = LutRGBA + 1 + debG;
 691       for( i=0; i < lengthG; ++i)
 692       {
 693          *a = LutGreenData[i*mult+1];
 694          a += 4;
 695       }
 696
 697       a = LutRGBA + 2 + debB;
 698       for(i=0; i < lengthB; ++i)
 699       {
 700          *a = LutBlueData[i*mult+1];
 701          a += 4;
 702       }
 703
 704       a = LutRGBA + 3 ;
 705       for(i=0; i < 256; ++i)
 706       {
 707          *a = 1; // Alpha component
 708          a += 4;
 709       }
 710    }
 711    else
 712    {
 713       // Probabely the same stuff is to be done for 16 Bits Pixels
 714       // with 65536 entries LUT ?!?
 715       // Still looking for accurate info on the web :-(
 716
 717       gdcmWarningMacro( "Sorry Palette Color Lookup Tables not yet dealt with"
 718                          << " for 16 Bits Per Pixel images" );
 719
 720       // forge the 4 * 16 Bits Red/Green/Blue/Alpha LUT
 721
 722       LutRGBA = (uint8_t *)new uint16_t[ 65536*4 ]; // 2^16 * 4 (R, G, B, Alpha)
 723       if ( !LutRGBA )
 724          return;
 725       memset( LutRGBA, 0, 65536*4*2 );  // 16 bits = 2 bytes ;-)
 726
 727       LutItemNumber = 65536;
 728       LutItemSize   = 16;
 729
 730       int i;
 731       uint16_t *a16;
 732
 733       //take "Subscript of the first Lut Value" (debR,debG,debB) into account!
 734
 735       a16 = (uint16_t*)LutRGBA + 0 + debR;
 736       for( i=0; i < lengthR; ++i )
 737       {
 738          *a16 = ((uint16_t*)LutRedData)[i];
 739          a16 += 4;
 740       }
 741
 742       a16 = (uint16_t*)LutRGBA + 1 + debG;
 743       for( i=0; i < lengthG; ++i)
 744       {
 745          *a16 = ((uint16_t*)LutGreenData)[i];
 746          a16 += 4;
 747       }
 748
 749       a16 = (uint16_t*)LutRGBA + 2 + debB;
 750       for(i=0; i < lengthB; ++i)
 751       {
 752          *a16 = ((uint16_t*)LutBlueData)[i];
 753          a16 += 4;
 754       }
 755
 756       a16 = (uint16_t*)LutRGBA + 3 ;
 757       for(i=0; i < 65536; ++i)
 758       {
 759          *a16 = 1; // Alpha component
 760          a16 += 4;
 761       }
 762 /* Just to 'see' the LUT, at debug time
 763 // Don't remove this commented out code.
 764
 765       a16=(uint16_t*)LutRGBA;
 766       for (int j=0;j<65536;j++)
 767       {
 768          std::cout << *a16     << " " << *(a16+1) << " "
 769                    << *(a16+2) << " " << *(a16+3) << std::endl;
 770          a16+=4;
 771       }
 772 */
 773    }
 774 }
 775
 776 /**
 777  * \brief Swap the bytes, according to \ref SwapCode.
 778  */
 779 void PixelReadConvert::ConvertSwapZone()
 780 {
 781    unsigned int i;
 782    uint16_t localSwapCode = SwapCode;
 783
 784    // If this file is 'ImplicitVR BigEndian PrivateGE Transfer Syntax',
 785    // then the header is in little endian format and the pixel data is in
 786    // big endian format.  When reading the header, GDCM has already established
 787    // a byte swapping code suitable for this machine to read the
 788    // header. In TS::ImplicitVRBigEndianPrivateGE, this code will need
 789    // to be switched in order to read the pixel data.  This must be
 790    // done REGARDLESS of the processor endianess!
 791    //
 792    // Example:  Assume we are on a little endian machine.  When
 793    // GDCM reads the header, the header will match the machine
 794    // endianess and the swap code will be established as a no-op.
 795    // When GDCM reaches the pixel data, it will need to switch the
 796    // swap code to do big endian to little endian conversion.
 797    //
 798    // Now, assume we are on a big endian machine.  When GDCM reads the
 799    // header, the header will be recognized as a different endianess
 800    // than the machine endianess, and a swap code will be established
 801    // to convert from little endian to big endian.  When GDCM readers
 802    // the pixel data, the pixel data endianess will now match the
 803    // machine endianess.  But we currently have a swap code that
 804    // converts from little endian to big endian.  In this case, we
 805    // need to switch the swap code to a no-op.
 806    //
 807    // Therefore, in either case, if the file is in
 808    // 'ImplicitVR BigEndian PrivateGE Transfer Syntax', then GDCM needs to switch
 809    // the byte swapping code when entering the pixel data.
 810
 811 /* //Let me check something.
 812    //I wait for the Dashboard !
 813    if ( IsPrivateGETransferSyntax )
 814    {
 815       // PrivateGETransferSyntax only exists for 'true' Dicom images
 816       // we assume there is no 'exotic' 32 bits endianess!
 817       switch (localSwapCode)
 818       {
 819          case 1234:
 820             localSwapCode = 4321;
 821             break;
 822          case 4321:
 823             localSwapCode = 1234;
 824             break;
 825       }
 826    }
 827 */
 828    if ( BitsAllocated == 16 )
 829    {
 830       uint16_t *im16 = (uint16_t*)Raw;
 831       switch( localSwapCode )
 832       {
 833          case 1234:
 834             break;
 835          case 3412:
 836          case 2143:
 837          case 4321:
 838             for( i = 0; i < RawSize / 2; i++ )
 839             {
 840                im16[i]= (im16[i] >> 8) | (im16[i] << 8 );
 841             }
 842             break;
 843          default:
 844             gdcmWarningMacro("SwapCode value (16 bits) not allowed.");
 845       }
 846    }
 847    else if ( BitsAllocated == 32 )
 848    {
 849       uint32_t s32;
 850       uint16_t high;
 851       uint16_t low;
 852       uint32_t *im32 = (uint32_t*)Raw;
 853       switch ( localSwapCode )
 854       {
 855          case 1234:
 856             break;
 857          case 4321:
 858             for( i = 0; i < RawSize / 4; i++ )
 859             {
 860                low     = im32[i] & 0x0000ffff;  // 4321
 861                high    = im32[i] >> 16;
 862                high    = ( high >> 8 ) | ( high << 8 );
 863                low     = ( low  >> 8 ) | ( low  << 8 );
 864                s32     = low;
 865                im32[i] = ( s32 << 16 ) | high;
 866             }
 867             break;
 868          case 2143:
 869             for( i = 0; i < RawSize / 4; i++ )
 870             {
 871                low     = im32[i] & 0x0000ffff;   // 2143
 872                high    = im32[i] >> 16;
 873                high    = ( high >> 8 ) | ( high << 8 );
 874                low     = ( low  >> 8 ) | ( low  << 8 );
 875                s32     = high;
 876                im32[i] = ( s32 << 16 ) | low;
 877             }
 878             break;
 879          case 3412:
 880             for( i = 0; i < RawSize / 4; i++ )
 881             {
 882                low     = im32[i] & 0x0000ffff; // 3412
 883                high    = im32[i] >> 16;
 884                s32     = low;
 885                im32[i] = ( s32 << 16 ) | high;
 886             }
 887             break;
 888          default:
 889             gdcmWarningMacro("SwapCode value (32 bits) not allowed." );
 890       }
 891    }
 892 }
 893
 894 /**
 895  * \brief Deal with endianness i.e. re-arange bytes inside the integer
 896  */
 897 void PixelReadConvert::ConvertReorderEndianity()
 898 {
 899    if ( BitsAllocated != 8 )
 900    {
 901       ConvertSwapZone();
 902    }
 903
 904    // Special kludge in order to deal with xmedcon broken images:
 905    if ( BitsAllocated == 16
 906      && BitsStored < BitsAllocated
 907      && !PixelSign )
 908    {
 909       int l = (int)( RawSize / ( BitsAllocated / 8 ) );
 910       uint16_t *deb = (uint16_t *)Raw;
 911       for(int i = 0; i<l; i++)
 912       {
 913          if ( *deb == 0xffff )
 914          {
 915            *deb = 0;
 916          }
 917          deb++;
 918       }
 919    }
 920 }
 921
 922 /**
 923  * \brief Deal with Grey levels i.e. re-arange them
 924  *        to have low values = dark, high values = bright
 925  */
 926 void PixelReadConvert::ConvertFixGreyLevels()
 927 {
 928    if (!IsMonochrome1)
 929       return;
 930
 931    uint32_t i; // to please M$VC6
 932    int16_t j;
 933
 934    if (!PixelSign)
 935    {
 936       if ( BitsAllocated == 8 )
 937       {
 938          uint8_t *deb = (uint8_t *)Raw;
 939          for (i=0; i<RawSize; i++)
 940          {
 941             *deb = 255 - *deb;
 942             deb++;
 943          }
 944          return;
 945       }
 946
 947       if ( BitsAllocated == 16 )
 948       {
 949          uint16_t mask =1;
 950          for (j=0; j<BitsStored-1; j++)
 951          {
 952             mask = (mask << 1) +1; // will be fff when BitsStored=12
 953          }
 954
 955          uint16_t *deb = (uint16_t *)Raw;
 956          for (i=0; i<RawSize/2; i++)
 957          {
 958             *deb = mask - *deb;
 959             deb++;
 960          }
 961          return;
 962        }
 963    }
 964    else
 965    {
 966       if ( BitsAllocated == 8 )
 967       {
 968          uint8_t smask8 = 255;
 969          uint8_t *deb = (uint8_t *)Raw;
 970          for (i=0; i<RawSize; i++)
 971          {
 972             *deb = smask8 - *deb;
 973             deb++;
 974          }
 975          return;
 976       }
 977       if ( BitsAllocated == 16 )
 978       {
 979          uint16_t smask16 = 65535;
 980          uint16_t *deb = (uint16_t *)Raw;
 981          for (i=0; i<RawSize/2; i++)
 982          {
 983             *deb = smask16 - *deb;
 984             deb++;
 985          }
 986          return;
 987       }
 988    }
 989 }
 990
 991 /**
 992  * \brief  Re-arrange the bits within the bytes.
 993  * @return Boolean always true
 994  */
 995 bool PixelReadConvert::ConvertReArrangeBits() throw ( FormatError )
 996 {
 997
 998    if ( BitsStored != BitsAllocated )
 999    {
1000       int l = (int)( RawSize / ( BitsAllocated / 8 ) );
1001       if ( BitsAllocated == 16 )
1002       {
1003          // pmask : to mask the 'unused bits' (may contain overlays)
1004          uint16_t pmask = 0xffff;
1005          pmask = pmask >> ( BitsAllocated - BitsStored );
1006
1007          uint16_t *deb = (uint16_t*)Raw;
1008
1009          if ( !PixelSign )  // Pixels are unsigned
1010          {
1011             for(int i = 0; i<l; i++)
1012             {
1013                *deb = (*deb >> (BitsStored - HighBitPosition - 1)) & pmask;
1014                deb++;
1015             }
1016          }
1017          else // Pixels are signed
1018          {
1019             // smask : to check the 'sign' when BitsStored != BitsAllocated
1020             uint16_t smask = 0x0001;
1021             smask = smask << ( 16 - (BitsAllocated - BitsStored + 1) );
1022             // nmask : to propagate sign bit on negative values
1023             int16_t nmask = (int16_t)0x8000;
1024             nmask = nmask >> ( BitsAllocated - BitsStored - 1 );
1025
1026             for(int i = 0; i<l; i++)
1027             {
1028                *deb = *deb >> (BitsStored - HighBitPosition - 1);
1029                if ( *deb & smask )
1030                {
1031                   *deb = *deb | nmask;
1032                }
1033                else
1034                {
1035                   *deb = *deb & pmask;
1036                }
1037                deb++;
1038             }
1039          }
1040       }
1041       else if ( BitsAllocated == 32 )
1042       {
1043          // pmask : to mask the 'unused bits' (may contain overlays)
1044          uint32_t pmask = 0xffffffff;
1045          pmask = pmask >> ( BitsAllocated - BitsStored );
1046
1047          uint32_t *deb = (uint32_t*)Raw;
1048
1049          if ( !PixelSign )
1050          {
1051             for(int i = 0; i<l; i++)
1052             {
1053                *deb = (*deb >> (BitsStored - HighBitPosition - 1)) & pmask;
1054                deb++;
1055             }
1056          }
1057          else
1058          {
1059             // smask : to check the 'sign' when BitsStored != BitsAllocated
1060             uint32_t smask = 0x00000001;
1061             smask = smask >> ( 32 - (BitsAllocated - BitsStored +1 ));
1062             // nmask : to propagate sign bit on negative values
1063             int32_t nmask = 0x80000000;
1064             nmask = nmask >> ( BitsAllocated - BitsStored -1 );
1065
1066             for(int i = 0; i<l; i++)
1067             {
1068                *deb = *deb >> (BitsStored - HighBitPosition - 1);
1069                if ( *deb & smask )
1070                   *deb = *deb | nmask;
1071                else
1072                   *deb = *deb & pmask;
1073                deb++;
1074             }
1075          }
1076       }
1077       else
1078       {
1079          gdcmWarningMacro("Weird image (BitsAllocated !=8, 12, 16, 32)");
1080          throw FormatError( "Weird image !?" );
1081       }
1082    }
1083    return true;
1084 }
1085
1086 /**
1087  * \brief   Convert (Red plane, Green plane, Blue plane) to RGB pixels
1088  * \warning Works on all the frames at a time
1089  */
1090 void PixelReadConvert::ConvertRGBPlanesToRGBPixels()
1091 {
1092    gdcmWarningMacro("--> ConvertRGBPlanesToRGBPixels");
1093
1094    uint8_t *localRaw = Raw;
1095    uint8_t *copyRaw = new uint8_t[ RawSize ];
1096    memmove( copyRaw, localRaw, RawSize );
1097
1098    int l = XSize * YSize * ZSize;
1099
1100    uint8_t *a = copyRaw;
1101    uint8_t *b = copyRaw + l;
1102    uint8_t *c = copyRaw + l + l;
1103
1104    for (int j = 0; j < l; j++)
1105    {
1106       *(localRaw++) = *(a++);
1107       *(localRaw++) = *(b++);
1108       *(localRaw++) = *(c++);
1109    }
1110    delete[] copyRaw;
1111 }
1112
1113 /**
1114  * \brief   Convert (cY plane, cB plane, cR plane) to RGB pixels
1115  * \warning Works on all the frames at a time
1116  */
1117 void PixelReadConvert::ConvertYcBcRPlanesToRGBPixels()
1118 {
1119   // Remarks for YBR newbees :
1120   // YBR_FULL works very much like RGB, i.e. three samples per pixel,
1121   // just the color space is YCbCr instead of RGB. This is particularly useful
1122   // for doppler ultrasound where most of the image is grayscale
1123   // (i.e. only populates the Y components) and Cb and Cr are mostly zero,
1124   // except for the few patches of color on the image.
1125   // On such images, RLE achieves a compression ratio that is much better
1126   // than the compression ratio on an equivalent RGB image.
1127
1128    gdcmWarningMacro("--> ConvertYcBcRPlanesToRGBPixels");
1129
1130    uint8_t *localRaw = Raw;
1131    uint8_t *copyRaw = new uint8_t[ RawSize ];
1132    memmove( copyRaw, localRaw, RawSize );
1133
1134    // to see the tricks about YBR_FULL, YBR_FULL_422,
1135    // YBR_PARTIAL_422, YBR_ICT, YBR_RCT have a look at :
1136    // ftp://medical.nema.org/medical/dicom/final/sup61_ft.pdf
1137    // and be *very* affraid
1138    //
1139    int l        = XSize * YSize;
1140    int nbFrames = ZSize;
1141
1142    uint8_t *a = copyRaw + 0;
1143    uint8_t *b = copyRaw + l;
1144    uint8_t *c = copyRaw + l+ l;
1145    int32_t R, G, B;
1146
1147    ///  We replaced easy to understand but time consuming floating point
1148    ///  computations by the 'well known' integer computation counterpart
1149    ///  Refer to :
1150    ///            http://lestourtereaux.free.fr/papers/data/yuvrgb.pdf
1151    ///  for code optimisation.
1152
1153    for ( int i = 0; i < nbFrames; i++ )
1154    {
1155       for ( int j = 0; j < l; j++ )
1156       {
1157          R = 38142 *(*a-16) + 52298 *(*c -128);
1158          G = 38142 *(*a-16) - 26640 *(*c -128) - 12845 *(*b -128);
1159          B = 38142 *(*a-16) + 66093 *(*b -128);
1160
1161          R = (R+16384)>>15;
1162          G = (G+16384)>>15;
1163          B = (B+16384)>>15;
1164
1165          if (R < 0)   R = 0;
1166          if (G < 0)   G = 0;
1167          if (B < 0)   B = 0;
1168          if (R > 255) R = 255;
1169          if (G > 255) G = 255;
1170          if (B > 255) B = 255;
1171
1172          *(localRaw++) = (uint8_t)R;
1173          *(localRaw++) = (uint8_t)G;
1174          *(localRaw++) = (uint8_t)B;
1175          a++;
1176          b++;
1177          c++;
1178       }
1179    }
1180    delete[] copyRaw;
1181 }
1182
1183 /// \brief Deals with the color decoding i.e. handle:
1184 ///   - R, G, B planes (as opposed to RGB pixels)
1185 ///   - YBR (various) encodings.
1186 ///   - LUT[s] (or "PALETTE COLOR").
1187
1188 void PixelReadConvert::ConvertHandleColor()
1189 {
1190    //////////////////////////////////
1191    // Deal with the color decoding i.e. handle:
1192    //   - R, G, B planes (as opposed to RGB pixels)
1193    //   - YBR (various) encodings.
1194    //   - LUT[s] (or "PALETTE COLOR").
1195    //
1196    // The classification in the color decoding schema is based on the blending
1197    // of two Dicom tags values:
1198    // * "Photometric Interpretation" for which we have the cases:
1199    //  - [Photo A] MONOCHROME[1|2] pictures,
1200    //  - [Photo B] RGB or YBR_FULL_422 (which acts as RGB),
1201    //  - [Photo C] YBR_* (with the above exception of YBR_FULL_422)
1202    //  - [Photo D] "PALETTE COLOR" which indicates the presence of LUT[s].
1203    // * "Planar Configuration" for which we have the cases:
1204    //  - [Planar 0] 0 then Pixels are already RGB
1205    //  - [Planar 1] 1 then we have 3 planes : R, G, B,
1206    //  - [Planar 2] 2 then we have 1 gray Plane and 3 LUTs
1207    //
1208    // Now in theory, one could expect some coherence when blending the above
1209    // cases. For example we should not encounter files belonging at the
1210    // time to case [Planar 0] and case [Photo D].
1211    // Alas, this was only theory ! Because in practice some odd (read ill
1212    // formated Dicom) files (e.g. gdcmData/US-PAL-8-10x-echo.dcm) we encounter:
1213    //     - "Planar Configuration" = 0,
1214    //     - "Photometric Interpretation" = "PALETTE COLOR".
1215    // Hence gdcm will use the folowing "heuristic" in order to be tolerant
1216    // towards Dicom-non-conformant files:
1217    //   << whatever the "Planar Configuration" value might be, a
1218    //      "Photometric Interpretation" set to "PALETTE COLOR" forces
1219    //      a LUT intervention >>
1220    //
1221    // Now we are left with the following handling of the cases:
1222    // - [Planar 0] OR  [Photo A] no color decoding (since respectively
1223    //       Pixels are already RGB and monochrome pictures have no color :),
1224    // - [Planar 1] AND [Photo B] handled with ConvertRGBPlanesToRGBPixels()
1225    // - [Planar 1] AND [Photo C] handled with ConvertYcBcRPlanesToRGBPixels()
1226    // - [Planar 2] OR  [Photo D] requires LUT intervention.
1227
1228    gdcmWarningMacro("--> ConvertHandleColor"
1229                     << "Planar Configuration " << PlanarConfiguration );
1230
1231    if ( ! IsRawRGB() )
1232    {
1233       // [Planar 2] OR  [Photo D]: LUT intervention done outside
1234       gdcmWarningMacro("--> RawRGB : LUT intervention done outside");
1235       return;
1236    }
1237
1238    if ( PlanarConfiguration == 1 )
1239    {
1240       if ( IsYBRFull )
1241       {
1242          // [Planar 1] AND [Photo C] (remember YBR_FULL_422 acts as RGB)
1243          gdcmWarningMacro("--> YBRFull");
1244          ConvertYcBcRPlanesToRGBPixels();
1245       }
1246       else
1247       {
1248          // [Planar 1] AND [Photo C]
1249          gdcmWarningMacro("--> YBRFull");
1250          ConvertRGBPlanesToRGBPixels();
1251       }
1252       return;
1253    }
1254
1255    // When planarConf is 0, and RLELossless (forbidden by Dicom norm)
1256    // pixels need to be RGB-fyied anyway
1257
1258    if (IsRLELossless)
1259    {
1260      gdcmWarningMacro("--> RLE Lossless");
1261      ConvertRGBPlanesToRGBPixels();
1262    }
1263
1264    // In *normal *case, when planarConf is 0, pixels are already in RGB
1265 }
1266
1267 /// Computes the Pixels Size
1268 void PixelReadConvert::ComputeRawAndRGBSizes()
1269 {
1270    int bitsAllocated = BitsAllocated;
1271    // Number of "Bits Allocated" is fixed to 16 when it's 12, since
1272    // in this case we will expand the image to 16 bits (see
1273    //    \ref ReadAndDecompress12BitsTo16Bits() )
1274    if (  BitsAllocated == 12 )
1275    {
1276       bitsAllocated = 16;
1277    }
1278
1279    RawSize =  XSize * YSize * ZSize
1280                      * ( bitsAllocated / 8 )
1281                      * SamplesPerPixel;
1282    if ( HasLUT )
1283    {
1284       RGBSize = 3 * RawSize; // works for 8 and 16 bits per Pixel
1285    }
1286    else
1287    {
1288       RGBSize = RawSize;
1289    }
1290 }
1291
1292 /// Allocates room for RGB Pixels
1293 void PixelReadConvert::AllocateRGB()
1294 {
1295   if ( RGB )
1296      delete [] RGB;
1297   RGB = new uint8_t[RGBSize];
1298 }
1299
1300 /// Allocates room for RAW Pixels
1301 void PixelReadConvert::AllocateRaw()
1302 {
1303   if ( Raw )
1304      delete [] Raw;
1305   Raw = new uint8_t[RawSize];
1306 }
1307
1308 //-----------------------------------------------------------------------------
1309 // Print
1310 /**
1311  * \brief        Print self.
1312  * @param indent Indentation string to be prepended during printing.
1313  * @param os     Stream to print to.
1314  */
1315 void PixelReadConvert::Print( std::ostream &os, std::string const &indent )
1316 {
1317    os << indent
1318       << "--- Pixel information -------------------------"
1319       << std::endl;
1320    os << indent
1321       << "Pixel Data: offset " << PixelOffset
1322       << " x(" << std::hex << PixelOffset << std::dec
1323       << ")   length " << PixelDataLength
1324       << " x(" << std::hex << PixelDataLength << std::dec
1325       << ")" << std::endl;
1326
1327    if ( IsRLELossless )
1328    {
1329       if ( RLEInfo )
1330       {
1331          RLEInfo->Print( os, indent );
1332       }
1333       else
1334       {
1335          gdcmWarningMacro("Set as RLE file but NO RLEinfo present.");
1336       }
1337    }
1338
1339    if ( IsJPEG2000 || IsJPEGLossless || IsJPEGLossy || IsJPEGLS )
1340    {
1341       if ( JPEGInfo )
1342       {
1343          JPEGInfo->Print( os, indent );
1344       }
1345       else
1346       {
1347          gdcmWarningMacro("Set as JPEG file but NO JPEGinfo present.");
1348       }
1349    }
1350 }
1351
1352 //-----------------------------------------------------------------------------
1353 } // end namespace gdcm
1354
1355 // Note to developpers :
1356 // Here is a very detailled post from David Clunie, on the troubles caused
1357 // 'non standard' LUT and LUT description
1358 // We shall have to take it into accound in our code.
1359 // Some day ...
1360
1361
1362 /*
1363 Subject: Problem with VOI LUTs in Agfa and Fuji CR and GE DX images, was Re: VOI LUT issues
1364 Date: Sun, 06 Feb 2005 17:13:40 GMT
1365 From: David Clunie <dclunie@dclunie.com>
1366 Reply-To: dclunie@dclunie.com
1367 Newsgroups: comp.protocols.dicom
1368 References: <1107553502.040221.189550@o13g2000cwo.googlegroups.com>
1369
1370 > THE LUT that comes with [my] image claims to be 16-bit, but none of the
1371 > values goes higher than 4095.  That being said, though, none of my
1372 > original pixel values goes higher than that, either.  I have read
1373 > elsewhere on this group that when that happens you are supposed to
1374 > adjust the LUT.  Can someone be more specific?  There was a thread from
1375 > 2002 where Marco and David were mentioning doing precisely that.
1376 >
1377 > Thanks
1378 >
1379 > -carlos rodriguez
1380
1381
1382 You have encountered the well known "we know what the standard says but
1383 we are going to ignore it and do what we have been doing for almost
1384 a decade regardless" CR vendor bug. Agfa started this, but they are not
1385 the only vendor doing this now; GE and Fuji may have joined the club.
1386
1387 Sadly, one needs to look at the LUT Data, figure out what the maximum
1388 value actually encoded is, and find the next highest power of 2 (e.g.
1389 212 in this case), to figure out what the range of the data is
1390 supposed to be. I have assumed that if the maximum value in the LUT
1391 data is less than a power of 2 minus 1 (e.g. 0xebc) then the intent
1392 of the vendor was not to use the maximum available grayscale range
1393 of the display (e.g. the maximum is 0xfff in this case). An alternative
1394 would be to scale to the actual maximum rather than a power of two.
1395
1396 Very irritating, and in theory not totally reliable if one really
1397 intended the full 16 bits and only used, say 15, but that is extremely
1398 unlikely since everything would be too dark, and this heuristic
1399 seems to work OK.
1400
1401 There has never been anything in the standard that describes having
1402 to go through these convolutions. Since the only value in the
1403 standard that describes the bit depth of the LUT values is LUT
1404 Descriptor value 3 and that is (usually) always required to be
1405 either 8 or 16, it mystifies me how the creators' of these images
1406 imagine that the receiver is going to divine the range that is intended. Further, the standard is quite explicit that this 3rd
1407 value defines the range of LUT values, but as far as I am aware, all
1408 the vendors are ignoring the standard and indeed sending a third value
1409 of 16 in all cases.
1410
1411 This problem is not confined to CR, and is also seen with DX products.
1412
1413 Typically I have seen:
1414
1415 - Agfa CR, which usually (always ?) sends LUTs, values up to 0x0fff
1416 - Fuji CR, which occasionally send LUTs, values up to 0x03ff
1417 - GE DX, for presentation, which always have LUTs, up to 0x3fff
1418
1419 Swissray, Siemens, Philips, Canon and Kodak never seem to send VOI LUTs
1420 at this point (which is a whole other problem). Note that the presence
1421 or absence of a VOI LUT as opposed to window values may be configurable
1422 on the modality in some cases, and I have just looked at what I happen
1423 to have received from a myriad of sites over whose configuration I have
1424 no control. This may be why the majority of Fuji images have no VOI LUTs,
1425 but a few do (or it may be the Siemens system that these Fuji images went
1426 through that perhaps added it). I do have some test Hologic DX images that
1427 are not from a clinical site that do actually get this right (a value
1428 of 12 for the third value and a max of 0xfff).
1429
1430 Since almost every vendor that I have encountered that encodes LUTs
1431 makes this mistake, perhaps it is time to amend the standard to warn
1432 implementor's of receivers and/or sanction this bad behavior. We have
1433 talked about this in the past in WG 6 but so far everyone has been
1434 reluctant to write into the standard such a comment. Maybe it is time
1435 to try again, since if one is not aware of this problem, one cannot
1436 effectively implement display using VOI LUTs, and there is a vast
1437 installed base to contend with.
1438
1439 I did not check presentation states, in which VOI LUTs could also be
1440 encountered, for the prevalence of this mistake, nor did I look at the
1441 encoding of Modality LUT's, which are unusual. Nor did I check digital
1442 mammography images. I would be interested to hear from anyone who has.
1443
1444 David
1445
1446 PS. The following older thread in this newsgroup discusses this:
1447
1448 "http://groups-beta.google.com/group/comp.protocols.dicom/browse_frm/t hread/6a033444802a35fc/0f0a9a1e35c1468e?q=voi+lut&_done=%2Fgroup%2Fcom p.protocols.dicom%2Fsearch%3Fgroup%3Dcomp.protocols.dicom%26q%3Dvoi+lu t%26qt_g%3D1%26searchnow%3DSearch+this+group%26&_doneTitle=Back+to+Sea rch&&d#0f0a9a1e35c1468e"
1449
1450 PPS. From a historical perspective, the following may be of interest.
1451
1452 In the original standard in 1993, all that was said about this was a
1453 reference to the corresponding such where Modality LUTs are described
1454 that said:
1455
1456 "The third value specifies the number of bits for each entry in the
1457 LUT Data. It shall take the value 8 or 16. The LUT Data shall be stored
1458 in a format equivalent to 8 or 16 bits allocated and high bit equal
1459 1-bits allocated."
1460
1461 Since the high bit hint was not apparently explicit enough, a very
1462 early CP, CP 15 (submitted by Agfa as it happens), replaced this with:
1463
1464 "The third value conveys the range of LUT entry values. It shall take
1465 the value 8 or 16, corresponding with the LUT entry value range of
1466 256 or 65536.
1467
1468 Note:    The third value is not required for describing the
1469     LUT data and is only included for informational usage
1470     and for maintaining compatibility with ACRNEMA 2.0.
1471
1472 The LUT Data contains the LUT entry values."
1473
1474 That is how it read in the 1996, 1998 and 1999 editions.
1475
1476 By the 2000 edition, Supplement 33 that introduced presentation states
1477 extensively reworked this entire section and tried to explain this in
1478 different words:
1479
1480 "The output range is from 0 to 2^n-1 where n is the third value of LUT
1481 Descriptor. This range is always unsigned."
1482
1483 and also added a note to spell out what the output range meant in the
1484 VOI LUT section:
1485
1486 "9. The output of the Window Center/Width or VOI LUT transformation
1487 is either implicitly scaled to the full range of the display device
1488 if there is no succeeding transformation defined, or implicitly scaled
1489 to the full input range of the succeeding transformation step (such as
1490 the Presentation LUT), if present. See C.11.6.1."
1491
1492 It still reads this way in the 2004 edition.
1493
1494 Note that LUTs in other applications than the general VOI LUT allow for
1495 values other than 8 or 16 in the third value of LUT descriptor to permit
1496 ranges other than 0 to 255 or 65535.
1497
1498 In addition, the DX Image Module specializes the VOI LUT
1499 attributes as follows, in PS 3.3 section C.8.11.3.1.5 (added in Sup 32):
1500
1501 "The third value specifies the number of bits for each entry in the LUT
1502 Data (analogous to ìbits storedî). It shall be between 10-16. The LUT
1503 Data shall be stored in a format equivalent to 16 ìbits allocatedî and
1504 ìhigh bitî equal to ìbits storedî - 1. The third value conveys the range
1505 of LUT entry values. These unsigned LUT entry values shall range between
1506 0 and 2^n-1, where n is the third value of the LUT Descriptor."
1507
1508 So in the case of the GE DX for presentation images, the third value of
1509 LUT descriptor is allowed to be and probably should be 14 rather than 16.
1510
1511 */