src/gdcmPixelReadConvert.cxx

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