src/gdcmPixelReadConvert.cxx

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