src/gdcmPixelReadConvert.cxx

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