src/gdcmPixelReadConvert.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: gdcmPixelReadConvert.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2005/07/01 11:25:51 $
   7   Version:   $Revision: 1.74 $
   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 }
  56
  57 /// Canonical Destructor
  58 PixelReadConvert::~PixelReadConvert()
  59 {
  60    Squeeze();
  61 }
  62
  63 //-----------------------------------------------------------------------------
  64 // Public
  65 /**
  66  * \brief Predicate to know whether the image[s] (once Raw) is RGB.
  67  * \note See comments of \ref ConvertHandleColor
  68  */
  69 bool PixelReadConvert::IsRawRGB()
  70 {
  71    if (   IsMonochrome
  72        || PlanarConfiguration == 2
  73        || IsPaletteColor )
  74    {
  75       return false;
  76    }
  77    return true;
  78 }
  79 /**
  80  * \brief Gets various usefull informations from the file header
  81  * @param file gdcm::File pointer
  82  */
  83 void PixelReadConvert::GrabInformationsFromFile( File *file )
  84 {
  85    // Number of Bits Allocated for storing a Pixel is defaulted to 16
  86    // when absent from the file.
  87    BitsAllocated = file->GetBitsAllocated();
  88    if ( BitsAllocated == 0 )
  89    {
  90       BitsAllocated = 16;
  91    }
  92
  93    // Number of "Bits Stored", defaulted to number of "Bits Allocated"
  94    // when absent from the file.
  95    BitsStored = file->GetBitsStored();
  96    if ( BitsStored == 0 )
  97    {
  98       BitsStored = BitsAllocated;
  99    }
 100
 101    // High Bit Position, defaulted to "Bits Allocated" - 1
 102    HighBitPosition = file->GetHighBitPosition();
 103    if ( HighBitPosition == 0 )
 104    {
 105       HighBitPosition = BitsAllocated - 1;
 106    }
 107
 108    XSize           = file->GetXSize();
 109    YSize           = file->GetYSize();
 110    ZSize           = file->GetZSize();
 111    SamplesPerPixel = file->GetSamplesPerPixel();
 112    PixelSize       = file->GetPixelSize();
 113    PixelSign       = file->IsSignedPixelData();
 114    SwapCode        = file->GetSwapCode();
 115    std::string ts  = file->GetTransferSyntax();
 116    IsRaw =
 117         ( ! file->IsDicomV3() )
 118      || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ImplicitVRLittleEndian
 119      || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ImplicitVRLittleEndianDLXGE
 120      || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ExplicitVRLittleEndian
 121      || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ExplicitVRBigEndian
 122      || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::DeflatedExplicitVRLittleEndian;
 123
 124    IsMPEG          = Global::GetTS()->IsMPEG(ts);
 125    IsJPEG2000      = Global::GetTS()->IsJPEG2000(ts);
 126    IsJPEGLS        = Global::GetTS()->IsJPEGLS(ts);
 127    IsJPEGLossy     = Global::GetTS()->IsJPEGLossy(ts);
 128    IsJPEGLossless  = Global::GetTS()->IsJPEGLossless(ts);
 129    IsRLELossless   = Global::GetTS()->IsRLELossless(ts);
 130
 131    PixelOffset     = file->GetPixelOffset();
 132    PixelDataLength = file->GetPixelAreaLength();
 133    RLEInfo         = file->GetRLEInfo();
 134    JPEGInfo        = file->GetJPEGInfo();
 135
 136    IsMonochrome    = file->IsMonochrome();
 137    IsMonochrome1   = file->IsMonochrome1();
 138    IsPaletteColor  = file->IsPaletteColor();
 139    IsYBRFull       = file->IsYBRFull();
 140
 141    PlanarConfiguration = file->GetPlanarConfiguration();
 142
 143    /////////////////////////////////////////////////////////////////
 144    // LUT section:
 145    HasLUT = file->HasLUT();
 146    if ( HasLUT )
 147    {
 148       // Just in case some access to a File element requires disk access.
 149       LutRedDescriptor   = file->GetEntryValue( 0x0028, 0x1101 );
 150       LutGreenDescriptor = file->GetEntryValue( 0x0028, 0x1102 );
 151       LutBlueDescriptor  = file->GetEntryValue( 0x0028, 0x1103 );
 152
 153       // Depending on the value of Document::MAX_SIZE_LOAD_ELEMENT_VALUE
 154       // [ refer to invocation of Document::SetMaxSizeLoadEntry() in
 155       // Document::Document() ], the loading of the value (content) of a
 156       // [Bin|Val]Entry occurence migth have been hindered (read simply NOT
 157       // loaded). Hence, we first try to obtain the LUTs data from the file
 158       // and when this fails we read the LUTs data directly from disk.
 159       // \TODO Reading a [Bin|Val]Entry directly from disk is a kludge.
 160       //       We should NOT bypass the [Bin|Val]Entry class. Instead
 161       //       an access to an UNLOADED content of a [Bin|Val]Entry occurence
 162       //       (e.g. BinEntry::GetBinArea()) should force disk access from
 163       //       within the [Bin|Val]Entry class itself. The only problem
 164       //       is that the [Bin|Val]Entry is unaware of the FILE* is was
 165       //       parsed from. Fix that. FIXME.
 166
 167       // //// Red round
 168       file->LoadEntryBinArea(0x0028, 0x1201);
 169       LutRedData = (uint8_t*)file->GetEntryBinArea( 0x0028, 0x1201 );
 170       if ( ! LutRedData )
 171       {
 172          gdcmWarningMacro( "Unable to read Red Palette Color Lookup Table data" );
 173       }
 174
 175       // //// Green round:
 176       file->LoadEntryBinArea(0x0028, 0x1202);
 177       LutGreenData = (uint8_t*)file->GetEntryBinArea(0x0028, 0x1202 );
 178       if ( ! LutGreenData)
 179       {
 180          gdcmWarningMacro( "Unable to read Green Palette Color Lookup Table data" );
 181       }
 182
 183       // //// Blue round:
 184       file->LoadEntryBinArea(0x0028, 0x1203);
 185       LutBlueData = (uint8_t*)file->GetEntryBinArea( 0x0028, 0x1203 );
 186       if ( ! LutBlueData )
 187       {
 188          gdcmWarningMacro( "Unable to read Blue Palette Color Lookup Table data" );
 189       }
 190    }
 191
 192    ComputeRawAndRGBSizes();
 193 }
 194
 195 /// \brief Reads from disk and decompresses Pixels
 196 bool PixelReadConvert::ReadAndDecompressPixelData( std::ifstream *fp )
 197 {
 198    // ComputeRawAndRGBSizes is already made by
 199    // ::GrabInformationsFromfile. So, the structure sizes are
 200    // correct
 201    Squeeze();
 202
 203    //////////////////////////////////////////////////
 204    //// First stage: get our hands on the Pixel Data.
 205    if ( !fp )
 206    {
 207       gdcmWarningMacro( "Unavailable file pointer." );
 208       return false;
 209    }
 210
 211    fp->seekg( PixelOffset, std::ios::beg );
 212    if ( fp->fail() || fp->eof() )
 213    {
 214       gdcmWarningMacro( "Unable to find PixelOffset in file." );
 215       return false;
 216    }
 217
 218    AllocateRaw();
 219
 220    //////////////////////////////////////////////////
 221    //// Second stage: read from disk and decompress.
 222    if ( BitsAllocated == 12 )
 223    {
 224       ReadAndDecompress12BitsTo16Bits( fp);
 225    }
 226    else if ( IsRaw )
 227    {
 228       // This problem can be found when some obvious informations are found
 229       // after the field containing the image data. In this case, these
 230       // bad data are added to the size of the image (in the PixelDataLength
 231       // variable). But RawSize is the right size of the image !
 232       if ( PixelDataLength != RawSize )
 233       {
 234          gdcmWarningMacro( "Mismatch between PixelReadConvert : "
 235                             << PixelDataLength << " and RawSize : " << RawSize );
 236       }
 237       if ( PixelDataLength > RawSize )
 238       {
 239          fp->read( (char*)Raw, RawSize);
 240       }
 241       else
 242       {
 243          fp->read( (char*)Raw, PixelDataLength);
 244       }
 245
 246       if ( fp->fail() || fp->eof())
 247       {
 248          gdcmWarningMacro( "Reading of Raw pixel data failed." );
 249          return false;
 250       }
 251    }
 252    else if ( IsRLELossless )
 253    {
 254       if ( ! RLEInfo->DecompressRLEFile( fp, Raw, XSize, YSize, ZSize, BitsAllocated ) )
 255       {
 256          gdcmWarningMacro( "RLE decompressor failed." );
 257          return false;
 258       }
 259    }
 260    else if ( IsMPEG )
 261    {
 262       //gdcmWarningMacro( "Sorry, MPEG not yet taken into account" );
 263       //return false;
 264       //ReadMPEGFile(fp, Raw, PixelDataLength); // fp has already been seek to start of mpeg
 265       return true;
 266    }
 267    else
 268    {
 269       // Default case concerns JPEG family
 270       if ( ! ReadAndDecompressJPEGFile( fp ) )
 271       {
 272          gdcmWarningMacro( "JPEG decompressor failed." );
 273          return false;
 274       }
 275    }
 276
 277    ////////////////////////////////////////////
 278    //// Third stage: twigle the bytes and bits.
 279    ConvertReorderEndianity();
 280    ConvertReArrangeBits();
 281    ConvertFixGreyLevels();
 282    ConvertHandleColor();
 283
 284    return true;
 285 }
 286
 287 /// Deletes Pixels Area
 288 void PixelReadConvert::Squeeze()
 289 {
 290    if ( RGB )
 291       delete [] RGB;
 292    RGB = 0;
 293
 294    if ( Raw )
 295       delete [] Raw;
 296    Raw = 0;
 297
 298    if ( LutRGBA )
 299       delete [] LutRGBA;
 300    LutRGBA = 0;
 301 }
 302
 303 /**
 304  * \brief Build the RGB image from the Raw image and the LUTs.
 305  */
 306 bool PixelReadConvert::BuildRGBImage()
 307 {
 308    if ( RGB )
 309    {
 310       // The job is already done.
 311       return true;
 312    }
 313
 314    if ( ! Raw )
 315    {
 316       // The job can't be done
 317       return false;
 318    }
 319
 320    BuildLUTRGBA();
 321    if ( ! LutRGBA )
 322    {
 323       // The job can't be done
 324       return false;
 325    }
 326
 327    gdcmWarningMacro( "--> BuildRGBImage" );
 328
 329    // Build RGB Pixels
 330    AllocateRGB();
 331
 332    int j;
 333    if ( BitsAllocated <= 8 )
 334    {
 335       uint8_t *localRGB = RGB;
 336       for (size_t i = 0; i < RawSize; ++i )
 337       {
 338          j  = Raw[i] * 4;
 339          *localRGB++ = LutRGBA[j];
 340          *localRGB++ = LutRGBA[j+1];
 341          *localRGB++ = LutRGBA[j+2];
 342       }
 343     }
 344
 345     else  // deal with 16 bits pixels and 16 bits Palette color
 346     {
 347       uint16_t *localRGB = (uint16_t *)RGB;
 348       for (size_t i = 0; i < RawSize/2; ++i )
 349       {
 350          j  = ((uint16_t *)Raw)[i] * 4;
 351          *localRGB++ = ((uint16_t *)LutRGBA)[j];
 352          *localRGB++ = ((uint16_t *)LutRGBA)[j+1];
 353          *localRGB++ = ((uint16_t *)LutRGBA)[j+2];
 354       }
 355     }
 356
 357    return true;
 358 }
 359
 360 //-----------------------------------------------------------------------------
 361 // Protected
 362
 363 //-----------------------------------------------------------------------------
 364 // Private
 365 /**
 366  * \brief Read from file a 12 bits per pixel image and decompress it
 367  *        into a 16 bits per pixel image.
 368  */
 369 void PixelReadConvert::ReadAndDecompress12BitsTo16Bits( std::ifstream *fp )
 370                throw ( FormatError )
 371 {
 372    int nbPixels = XSize * YSize;
 373    uint16_t *localDecompres = (uint16_t*)Raw;
 374
 375    for( int p = 0; p < nbPixels; p += 2 )
 376    {
 377       uint8_t b0, b1, b2;
 378
 379       fp->read( (char*)&b0, 1);
 380       if ( fp->fail() || fp->eof() )
 381       {
 382          throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
 383                                 "Unfound first block" );
 384       }
 385
 386       fp->read( (char*)&b1, 1 );
 387       if ( fp->fail() || fp->eof())
 388       {
 389          throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
 390                                 "Unfound second block" );
 391       }
 392
 393       fp->read( (char*)&b2, 1 );
 394       if ( fp->fail() || fp->eof())
 395       {
 396          throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
 397                                 "Unfound second block" );
 398       }
 399
 400       // Two steps are necessary to please VC++
 401       //
 402       // 2 pixels 12bit =     [0xABCDEF]
 403       // 2 pixels 16bit = [0x0ABD] + [0x0FCE]
 404       //                        A                     B                 D
 405       *localDecompres++ =  ((b0 >> 4) << 8) + ((b0 & 0x0f) << 4) + (b1 & 0x0f);
 406       //                        F                     C                 E
 407       *localDecompres++ =  ((b2 & 0x0f) << 8) + ((b1 >> 4) << 4) + (b2 >> 4);
 408
 409       /// \todo JPR Troubles expected on Big-Endian processors ?
 410    }
 411 }
 412
 413 /**
 414  * \brief     Reads from disk the Pixel Data of JPEG Dicom encapsulated
 415  *            file and decompress it.
 416  * @param     fp File Pointer
 417  * @return    Boolean
 418  */
 419 bool PixelReadConvert::ReadAndDecompressJPEGFile( std::ifstream *fp )
 420 {
 421    if ( IsJPEG2000 )
 422    {
 423      // make sure this is the right JPEG compression
 424      assert( !IsJPEGLossless || !IsJPEGLossy || !IsJPEGLS );
 425      // FIXME this is really ugly but it seems I have to load the complete
 426      // jpeg2000 stream to use jasper:
 427      // I don't think we'll ever be able to deal with multiple fragments properly
 428
 429       unsigned long inputlength = 0;
 430       JPEGFragment *jpegfrag = JPEGInfo->GetFirstFragment();
 431       while( jpegfrag )
 432       {
 433          inputlength += jpegfrag->GetLength();
 434          jpegfrag = JPEGInfo->GetNextFragment();
 435       }
 436       gdcmAssertMacro( inputlength != 0);
 437       uint8_t *inputdata = new uint8_t[inputlength];
 438       char *pinputdata = (char*)inputdata;
 439       jpegfrag = JPEGInfo->GetFirstFragment();
 440       while( jpegfrag )
 441       {
 442          fp->seekg( jpegfrag->GetOffset(), std::ios::beg);
 443          fp->read(pinputdata, jpegfrag->GetLength());
 444          pinputdata += jpegfrag->GetLength();
 445          jpegfrag = JPEGInfo->GetNextFragment();
 446       }
 447       // Warning the inputdata buffer is delete in the function
 448       if ( ! gdcm_read_JPEG2000_file( Raw,
 449           (char*)inputdata, inputlength ) )
 450       {
 451          return true;
 452       }
 453       // wow what happen, must be an error
 454       return false;
 455    }
 456    else if ( IsJPEGLS )
 457    {
 458      // make sure this is the right JPEG compression
 459      assert( !IsJPEGLossless || !IsJPEGLossy || !IsJPEG2000 );
 460    // WARNING : JPEG-LS is NOT the 'classical' Jpeg Lossless :
 461    // [JPEG-LS is the basis for new lossless/near-lossless compression
 462    // standard for continuous-tone images intended for JPEG2000. The standard
 463    // is based on the LOCO-I algorithm (LOw COmplexity LOssless COmpression
 464    // for Images) developed at Hewlett-Packard Laboratories]
 465    //
 466    // see http://datacompression.info/JPEGLS.shtml
 467    //
 468 #if 0
 469    std::cerr << "count:" << JPEGInfo->GetFragmentCount() << std::endl;
 470       unsigned long inputlength = 0;
 471       JPEGFragment *jpegfrag = JPEGInfo->GetFirstFragment();
 472       while( jpegfrag )
 473       {
 474          inputlength += jpegfrag->GetLength();
 475          jpegfrag = JPEGInfo->GetNextFragment();
 476       }
 477       gdcmAssertMacro( inputlength != 0);
 478       uint8_t *inputdata = new uint8_t[inputlength];
 479       char *pinputdata = (char*)inputdata;
 480       jpegfrag = JPEGInfo->GetFirstFragment();
 481       while( jpegfrag )
 482       {
 483          fp->seekg( jpegfrag->GetOffset(), std::ios::beg);
 484          fp->read(pinputdata, jpegfrag->GetLength());
 485          pinputdata += jpegfrag->GetLength();
 486          jpegfrag = JPEGInfo->GetNextFragment();
 487       }
 488
 489   //fp->read((char*)Raw, PixelDataLength);
 490
 491   std::ofstream out("/tmp/jpegls.jpg");
 492   out.write((char*)inputdata, inputlength);
 493   out.close();
 494   delete[] inputdata;
 495 #endif
 496
 497       gdcmWarningMacro( "Sorry, JPEG-LS not yet taken into account" );
 498       fp->seekg( JPEGInfo->GetFirstFragment()->GetOffset(), std::ios::beg);
 499 //    if ( ! gdcm_read_JPEGLS_file( fp,Raw ) )
 500          return false;
 501    }
 502    else
 503    {
 504      // make sure this is the right JPEG compression
 505      assert( !IsJPEGLS || !IsJPEG2000 );
 506      // Precompute the offset localRaw will be shifted with
 507      int length = XSize * YSize * SamplesPerPixel;
 508      int numberBytes = BitsAllocated / 8;
 509
 510      JPEGInfo->DecompressFromFile(fp, Raw, BitsStored, numberBytes, length );
 511      return true;
 512    }
 513 }
 514
 515 /**
 516  * \brief Build Red/Green/Blue/Alpha LUT from File
 517  *         when (0028,0004),Photometric Interpretation = [PALETTE COLOR ]
 518  *          and (0028,1101),(0028,1102),(0028,1102)
 519  *            - xxx Palette Color Lookup Table Descriptor - are found
 520  *          and (0028,1201),(0028,1202),(0028,1202)
 521  *            - xxx Palette Color Lookup Table Data - are found
 522  * \warning does NOT deal with :
 523  *   0028 1100 Gray Lookup Table Descriptor (Retired)
 524  *   0028 1221 Segmented Red Palette Color Lookup Table Data
 525  *   0028 1222 Segmented Green Palette Color Lookup Table Data
 526  *   0028 1223 Segmented Blue Palette Color Lookup Table Data
 527  *   no known Dicom reader deals with them :-(
 528  * @return a RGBA Lookup Table
 529  */
 530 void PixelReadConvert::BuildLUTRGBA()
 531 {
 532    if ( LutRGBA )
 533    {
 534       return;
 535    }
 536    // Not so easy : see
 537    // http://www.barre.nom.fr/medical/dicom2/limitations.html#Color%20Lookup%20Tables
 538
 539    if ( ! IsPaletteColor )
 540    {
 541       return;
 542    }
 543
 544    if (   LutRedDescriptor   == GDCM_UNFOUND
 545        || LutGreenDescriptor == GDCM_UNFOUND
 546        || LutBlueDescriptor  == GDCM_UNFOUND )
 547    {
 548       gdcmWarningMacro( "(At least) a LUT Descriptor is missing" );
 549       return;
 550    }
 551
 552    ////////////////////////////////////////////
 553    // Extract the info from the LUT descriptors
 554    int lengthR;   // Red LUT length in Bytes
 555    int debR;      // Subscript of the first Lut Value
 556    int nbitsR;    // Lut item size (in Bits)
 557    int nbRead;    // nb of items in LUT descriptor (must be = 3)
 558
 559    nbRead = sscanf( LutRedDescriptor.c_str(),
 560                         "%d\\%d\\%d",
 561                         &lengthR, &debR, &nbitsR );
 562    if ( nbRead != 3 )
 563    {
 564       gdcmWarningMacro( "Wrong Red LUT descriptor" );
 565    }
 566    int lengthG;  // Green LUT length in Bytes
 567    int debG;     // Subscript of the first Lut Value
 568    int nbitsG;   // Lut item size (in Bits)
 569
 570    nbRead = sscanf( LutGreenDescriptor.c_str(),
 571                     "%d\\%d\\%d",
 572                     &lengthG, &debG, &nbitsG );
 573    if ( nbRead != 3 )
 574    {
 575       gdcmWarningMacro( "Wrong Green LUT descriptor" );
 576    }
 577
 578    int lengthB;  // Blue LUT length in Bytes
 579    int debB;     // Subscript of the first Lut Value
 580    int nbitsB;   // Lut item size (in Bits)
 581    nbRead = sscanf( LutRedDescriptor.c_str(),
 582                     "%d\\%d\\%d",
 583                     &lengthB, &debB, &nbitsB );
 584    if ( nbRead != 3 )
 585    {
 586       gdcmWarningMacro( "Wrong Blue LUT descriptor" );
 587    }
 588
 589    gdcmWarningMacro(" lengthR " << lengthR << " debR "
 590                  << debR << " nbitsR " << nbitsR);
 591    gdcmWarningMacro(" lengthG " << lengthG << " debG "
 592                  << debG << " nbitsG " << nbitsG);
 593    gdcmWarningMacro(" lengthB " << lengthB << " debB "
 594                  << debB << " nbitsB " << nbitsB);
 595
 596    if ( !lengthR ) // if = 2^16, this shall be 0 see : CP-143
 597       lengthR=65536;
 598    if ( !lengthG ) // if = 2^16, this shall be 0
 599       lengthG=65536;
 600    if ( !lengthB ) // if = 2^16, this shall be 0
 601       lengthB=65536;
 602
 603    ////////////////////////////////////////////////////////
 604
 605    if ( ( ! LutRedData ) || ( ! LutGreenData ) || ( ! LutBlueData ) )
 606    {
 607       gdcmWarningMacro( "(At least) a LUT is missing" );
 608       return;
 609    }
 610
 611    // -------------------------------------------------------------
 612
 613    if ( BitsAllocated <= 8 )
 614    {
 615       // forge the 4 * 8 Bits Red/Green/Blue/Alpha LUT
 616       LutRGBA = new uint8_t[ 1024 ]; // 256 * 4 (R, G, B, Alpha)
 617       if ( !LutRGBA )
 618          return;
 619       LutItemNumber = 256;
 620       LutItemSize   = 8;
 621       memset( LutRGBA, 0, 1024 );
 622
 623       int mult;
 624       if ( ( nbitsR == 16 ) && ( BitsAllocated == 8 ) )
 625       {
 626          // when LUT item size is different than pixel size
 627          mult = 2; // high byte must be = low byte
 628       }
 629       else
 630       {
 631          // See PS 3.3-2003 C.11.1.1.2 p 619
 632          mult = 1;
 633       }
 634
 635       // if we get a black image, let's just remove the '+1'
 636       // from 'i*mult+1' and check again
 637       // if it works, we shall have to check the 3 Palettes
 638       // to see which byte is ==0 (first one, or second one)
 639       // and fix the code
 640       // We give up the checking to avoid some (useless ?) overhead
 641       // (optimistic asumption)
 642       int i;
 643       uint8_t *a;
 644
 645       //take "Subscript of the first Lut Value" (debR,debG,debB) into account!
 646
 647       //FIXME :  +1 : to get 'low value' byte
 648       //         Trouble expected on Big Endian Processors ?
 649       //         16 BIts Per Pixel Palette Color to be swapped?
 650
 651       a = LutRGBA + 0 + debR;
 652       for( i=0; i < lengthR; ++i )
 653       {
 654          *a = LutRedData[i*mult+1];
 655          a += 4;
 656       }
 657
 658       a = LutRGBA + 1 + debG;
 659       for( i=0; i < lengthG; ++i)
 660       {
 661          *a = LutGreenData[i*mult+1];
 662          a += 4;
 663       }
 664
 665       a = LutRGBA + 2 + debB;
 666       for(i=0; i < lengthB; ++i)
 667       {
 668          *a = LutBlueData[i*mult+1];
 669          a += 4;
 670       }
 671
 672       a = LutRGBA + 3 ;
 673       for(i=0; i < 256; ++i)
 674       {
 675          *a = 1; // Alpha component
 676          a += 4;
 677       }
 678    }
 679    else
 680    {
 681       // Probabely the same stuff is to be done for 16 Bits Pixels
 682       // with 65536 entries LUT ?!?
 683       // Still looking for accurate info on the web :-(
 684
 685       gdcmWarningMacro( "Sorry Palette Color Lookup Tables not yet dealt with"
 686                          << " for 16 Bits Per Pixel images" );
 687
 688       // forge the 4 * 16 Bits Red/Green/Blue/Alpha LUT
 689
 690       LutRGBA = (uint8_t *)new uint16_t[ 65536*4 ]; // 2^16 * 4 (R, G, B, Alpha)
 691       if ( !LutRGBA )
 692          return;
 693       memset( LutRGBA, 0, 65536*4*2 );  // 16 bits = 2 bytes ;-)
 694
 695       LutItemNumber = 65536;
 696       LutItemSize   = 16;
 697
 698       int i;
 699       uint16_t *a16;
 700
 701       //take "Subscript of the first Lut Value" (debR,debG,debB) into account!
 702
 703       a16 = (uint16_t*)LutRGBA + 0 + debR;
 704       for( i=0; i < lengthR; ++i )
 705       {
 706          *a16 = ((uint16_t*)LutRedData)[i];
 707          a16 += 4;
 708       }
 709
 710       a16 = (uint16_t*)LutRGBA + 1 + debG;
 711       for( i=0; i < lengthG; ++i)
 712       {
 713          *a16 = ((uint16_t*)LutGreenData)[i];
 714          a16 += 4;
 715       }
 716
 717       a16 = (uint16_t*)LutRGBA + 2 + debB;
 718       for(i=0; i < lengthB; ++i)
 719       {
 720          *a16 = ((uint16_t*)LutBlueData)[i];
 721          a16 += 4;
 722       }
 723
 724       a16 = (uint16_t*)LutRGBA + 3 ;
 725       for(i=0; i < 65536; ++i)
 726       {
 727          *a16 = 1; // Alpha component
 728          a16 += 4;
 729       }
 730 /* Just to 'see' the LUT, at debug time
 731
 732       a16=(uint16_t*)LutRGBA;
 733       for (int j=0;j<65536;j++)
 734       {
 735          std::cout << *a16     << " " << *(a16+1) << " "
 736                    << *(a16+2) << " " << *(a16+3) << std::endl;
 737          a16+=4;
 738       }
 739 */
 740    }
 741 }
 742
 743 /**
 744  * \brief Swap the bytes, according to \ref SwapCode.
 745  */
 746 void PixelReadConvert::ConvertSwapZone()
 747 {
 748    unsigned int i;
 749
 750    if ( BitsAllocated == 16 )
 751    {
 752       uint16_t *im16 = (uint16_t*)Raw;
 753       switch( SwapCode )
 754       {
 755          case 1234:
 756             break;
 757          case 3412:
 758          case 2143:
 759          case 4321:
 760             for( i = 0; i < RawSize / 2; i++ )
 761             {
 762                im16[i]= (im16[i] >> 8) | (im16[i] << 8 );
 763             }
 764             break;
 765          default:
 766             gdcmWarningMacro("SwapCode value (16 bits) not allowed.");
 767       }
 768    }
 769    else if ( BitsAllocated == 32 )
 770    {
 771       uint32_t s32;
 772       uint16_t high;
 773       uint16_t low;
 774       uint32_t *im32 = (uint32_t*)Raw;
 775       switch ( SwapCode )
 776       {
 777          case 1234:
 778             break;
 779          case 4321:
 780             for( i = 0; i < RawSize / 4; i++ )
 781             {
 782                low     = im32[i] & 0x0000ffff;  // 4321
 783                high    = im32[i] >> 16;
 784                high    = ( high >> 8 ) | ( high << 8 );
 785                low     = ( low  >> 8 ) | ( low  << 8 );
 786                s32     = low;
 787                im32[i] = ( s32 << 16 ) | high;
 788             }
 789             break;
 790          case 2143:
 791             for( i = 0; i < RawSize / 4; i++ )
 792             {
 793                low     = im32[i] & 0x0000ffff;   // 2143
 794                high    = im32[i] >> 16;
 795                high    = ( high >> 8 ) | ( high << 8 );
 796                low     = ( low  >> 8 ) | ( low  << 8 );
 797                s32     = high;
 798                im32[i] = ( s32 << 16 ) | low;
 799             }
 800             break;
 801          case 3412:
 802             for( i = 0; i < RawSize / 4; i++ )
 803             {
 804                low     = im32[i] & 0x0000ffff; // 3412
 805                high    = im32[i] >> 16;
 806                s32     = low;
 807                im32[i] = ( s32 << 16 ) | high;
 808             }
 809             break;
 810          default:
 811             gdcmWarningMacro("SwapCode value (32 bits) not allowed." );
 812       }
 813    }
 814 }
 815
 816 /**
 817  * \brief Deal with endianness i.e. re-arange bytes inside the integer
 818  */
 819 void PixelReadConvert::ConvertReorderEndianity()
 820 {
 821    if ( BitsAllocated != 8 )
 822    {
 823       ConvertSwapZone();
 824    }
 825
 826    // Special kludge in order to deal with xmedcon broken images:
 827    if ( BitsAllocated == 16
 828      && BitsStored < BitsAllocated
 829      && !PixelSign )
 830    {
 831       int l = (int)( RawSize / ( BitsAllocated / 8 ) );
 832       uint16_t *deb = (uint16_t *)Raw;
 833       for(int i = 0; i<l; i++)
 834       {
 835          if ( *deb == 0xffff )
 836          {
 837            *deb = 0;
 838          }
 839          deb++;
 840       }
 841    }
 842 }
 843
 844 /**
 845  * \brief Deal with Grey levels i.e. re-arange them
 846  *        to have low values = dark, high values = bright
 847  */
 848 void PixelReadConvert::ConvertFixGreyLevels()
 849 {
 850    if (!IsMonochrome1)
 851       return;
 852
 853    uint32_t i; // to please M$VC6
 854    int16_t j;
 855
 856    if (!PixelSign)
 857    {
 858       if ( BitsAllocated == 8 )
 859       {
 860          uint8_t *deb = (uint8_t *)Raw;
 861          for (i=0; i<RawSize; i++)
 862          {
 863             *deb = 255 - *deb;
 864             deb++;
 865          }
 866          return;
 867       }
 868
 869       if ( BitsAllocated == 16 )
 870       {
 871          uint16_t mask =1;
 872          for (j=0; j<BitsStored-1; j++)
 873          {
 874             mask = (mask << 1) +1; // will be fff when BitsStored=12
 875          }
 876
 877          uint16_t *deb = (uint16_t *)Raw;
 878          for (i=0; i<RawSize/2; i++)
 879          {
 880             *deb = mask - *deb;
 881             deb++;
 882          }
 883          return;
 884        }
 885    }
 886    else
 887    {
 888       if ( BitsAllocated == 8 )
 889       {
 890          uint8_t smask8 = 255;
 891          uint8_t *deb = (uint8_t *)Raw;
 892          for (i=0; i<RawSize; i++)
 893          {
 894             *deb = smask8 - *deb;
 895             deb++;
 896          }
 897          return;
 898       }
 899       if ( BitsAllocated == 16 )
 900       {
 901          uint16_t smask16 = 65535;
 902          uint16_t *deb = (uint16_t *)Raw;
 903          for (i=0; i<RawSize/2; i++)
 904          {
 905             *deb = smask16 - *deb;
 906             deb++;
 907          }
 908          return;
 909       }
 910    }
 911 }
 912
 913 /**
 914  * \brief  Re-arrange the bits within the bytes.
 915  * @return Boolean always true
 916  */
 917 bool PixelReadConvert::ConvertReArrangeBits() throw ( FormatError )
 918 {
 919
 920    if ( BitsStored != BitsAllocated )
 921    {
 922       int l = (int)( RawSize / ( BitsAllocated / 8 ) );
 923       if ( BitsAllocated == 16 )
 924       {
 925          // pmask : to mask the 'unused bits' (may contain overlays)
 926          uint16_t pmask = 0xffff;
 927          pmask = pmask >> ( BitsAllocated - BitsStored );
 928
 929          uint16_t *deb = (uint16_t*)Raw;
 930
 931          if ( !PixelSign )  // Pixels are unsigned
 932          {
 933             for(int i = 0; i<l; i++)
 934             {
 935                *deb = (*deb >> (BitsStored - HighBitPosition - 1)) & pmask;
 936                deb++;
 937             }
 938          }
 939          else // Pixels are signed
 940          {
 941             // smask : to check the 'sign' when BitsStored != BitsAllocated
 942             uint16_t smask = 0x0001;
 943             smask = smask << ( 16 - (BitsAllocated - BitsStored + 1) );
 944             // nmask : to propagate sign bit on negative values
 945             int16_t nmask = (int16_t)0x8000;
 946             nmask = nmask >> ( BitsAllocated - BitsStored - 1 );
 947 /*
 948 std::cout  << "BitsStored "     << BitsStored
 949            << " BitsAllocated " << BitsAllocated
 950            << std::endl;
 951 std::cout  << std::hex  << "pmask " << pmask
 952            << " smask " << smask
 953            << " nmask " << nmask
 954            << std::endl;
 955 */
 956             for(int i = 0; i<l; i++)
 957             {
 958                *deb = *deb >> (BitsStored - HighBitPosition - 1);
 959                if ( *deb & smask )
 960                {
 961                   *deb = *deb | nmask;
 962                }
 963                else
 964                {
 965                   *deb = *deb & pmask;
 966                }
 967                deb++;
 968             }
 969          }
 970       }
 971       else if ( BitsAllocated == 32 )
 972       {
 973          // pmask : to mask the 'unused bits' (may contain overlays)
 974          uint32_t pmask = 0xffffffff;
 975          pmask = pmask >> ( BitsAllocated - BitsStored );
 976
 977          uint32_t *deb = (uint32_t*)Raw;
 978
 979          if ( !PixelSign )
 980          {
 981             for(int i = 0; i<l; i++)
 982             {
 983                *deb = (*deb >> (BitsStored - HighBitPosition - 1)) & pmask;
 984                deb++;
 985             }
 986          }
 987          else
 988          {
 989             // smask : to check the 'sign' when BitsStored != BitsAllocated
 990             uint32_t smask = 0x00000001;
 991             smask = smask >> ( 32 - (BitsAllocated - BitsStored +1 ));
 992             // nmask : to propagate sign bit on negative values
 993             int32_t nmask = 0x80000000;
 994             nmask = nmask >> ( BitsAllocated - BitsStored -1 );
 995
 996             for(int i = 0; i<l; i++)
 997             {
 998                *deb = *deb >> (BitsStored - HighBitPosition - 1);
 999                if ( *deb & smask )
1000                   *deb = *deb | nmask;
1001                else
1002                   *deb = *deb & pmask;
1003                deb++;
1004             }
1005          }
1006       }
1007       else
1008       {
1009          gdcmWarningMacro("Weird image (BitsAllocated !=8, 12, 16, 32)");
1010          throw FormatError( "Weird image !?" );
1011       }
1012    }
1013    return true;
1014 }
1015
1016 /**
1017  * \brief   Convert (Red plane, Green plane, Blue plane) to RGB pixels
1018  * \warning Works on all the frames at a time
1019  */
1020 void PixelReadConvert::ConvertRGBPlanesToRGBPixels()
1021 {
1022    gdcmWarningMacro("--> ConvertRGBPlanesToRGBPixels");
1023
1024    uint8_t *localRaw = Raw;
1025    uint8_t *copyRaw = new uint8_t[ RawSize ];
1026    memmove( copyRaw, localRaw, RawSize );
1027
1028    int l = XSize * YSize * ZSize;
1029
1030    uint8_t *a = copyRaw;
1031    uint8_t *b = copyRaw + l;
1032    uint8_t *c = copyRaw + l + l;
1033
1034    for (int j = 0; j < l; j++)
1035    {
1036       *(localRaw++) = *(a++);
1037       *(localRaw++) = *(b++);
1038       *(localRaw++) = *(c++);
1039    }
1040    delete[] copyRaw;
1041 }
1042
1043 /**
1044  * \brief   Convert (cY plane, cB plane, cR plane) to RGB pixels
1045  * \warning Works on all the frames at a time
1046  */
1047 void PixelReadConvert::ConvertYcBcRPlanesToRGBPixels()
1048 {
1049   // Remarks for YBR newbees :
1050   // YBR_FULL works very much like RGB, i.e. three samples per pixel,
1051   // just the color space is YCbCr instead of RGB. This is particularly useful
1052   // for doppler ultrasound where most of the image is grayscale
1053   // (i.e. only populates the Y components) and Cb and Cr are mostly zero,
1054   // except for the few patches of color on the image.
1055   // On such images, RLE achieves a compression ratio that is much better
1056   // than the compression ratio on an equivalent RGB image.
1057
1058    gdcmWarningMacro("--> ConvertYcBcRPlanesToRGBPixels");
1059
1060    uint8_t *localRaw = Raw;
1061    uint8_t *copyRaw = new uint8_t[ RawSize ];
1062    memmove( copyRaw, localRaw, RawSize );
1063
1064    // to see the tricks about YBR_FULL, YBR_FULL_422,
1065    // YBR_PARTIAL_422, YBR_ICT, YBR_RCT have a look at :
1066    // ftp://medical.nema.org/medical/dicom/final/sup61_ft.pdf
1067    // and be *very* affraid
1068    //
1069    int l        = XSize * YSize;
1070    int nbFrames = ZSize;
1071
1072    uint8_t *a = copyRaw + 0;
1073    uint8_t *b = copyRaw + l;
1074    uint8_t *c = copyRaw + l+ l;
1075    int32_t R, G, B;
1076
1077    ///  We replaced easy to understand but time consuming floating point
1078    ///  computations by the 'well known' integer computation counterpart
1079    ///  Refer to :
1080    ///            http://lestourtereaux.free.fr/papers/data/yuvrgb.pdf
1081    ///  for code optimisation.
1082
1083    for ( int i = 0; i < nbFrames; i++ )
1084    {
1085       for ( int j = 0; j < l; j++ )
1086       {
1087          R = 38142 *(*a-16) + 52298 *(*c -128);
1088          G = 38142 *(*a-16) - 26640 *(*c -128) - 12845 *(*b -128);
1089          B = 38142 *(*a-16) + 66093 *(*b -128);
1090
1091          R = (R+16384)>>15;
1092          G = (G+16384)>>15;
1093          B = (B+16384)>>15;
1094
1095          if (R < 0)   R = 0;
1096          if (G < 0)   G = 0;
1097          if (B < 0)   B = 0;
1098          if (R > 255) R = 255;
1099          if (G > 255) G = 255;
1100          if (B > 255) B = 255;
1101
1102          *(localRaw++) = (uint8_t)R;
1103          *(localRaw++) = (uint8_t)G;
1104          *(localRaw++) = (uint8_t)B;
1105          a++;
1106          b++;
1107          c++;
1108       }
1109    }
1110    delete[] copyRaw;
1111 }
1112
1113 /// \brief Deals with the color decoding i.e. handle:
1114 ///   - R, G, B planes (as opposed to RGB pixels)
1115 ///   - YBR (various) encodings.
1116 ///   - LUT[s] (or "PALETTE COLOR").
1117
1118 void PixelReadConvert::ConvertHandleColor()
1119 {
1120    //////////////////////////////////
1121    // Deal with the color decoding i.e. handle:
1122    //   - R, G, B planes (as opposed to RGB pixels)
1123    //   - YBR (various) encodings.
1124    //   - LUT[s] (or "PALETTE COLOR").
1125    //
1126    // The classification in the color decoding schema is based on the blending
1127    // of two Dicom tags values:
1128    // * "Photometric Interpretation" for which we have the cases:
1129    //  - [Photo A] MONOCHROME[1|2] pictures,
1130    //  - [Photo B] RGB or YBR_FULL_422 (which acts as RGB),
1131    //  - [Photo C] YBR_* (with the above exception of YBR_FULL_422)
1132    //  - [Photo D] "PALETTE COLOR" which indicates the presence of LUT[s].
1133    // * "Planar Configuration" for which we have the cases:
1134    //  - [Planar 0] 0 then Pixels are already RGB
1135    //  - [Planar 1] 1 then we have 3 planes : R, G, B,
1136    //  - [Planar 2] 2 then we have 1 gray Plane and 3 LUTs
1137    //
1138    // Now in theory, one could expect some coherence when blending the above
1139    // cases. For example we should not encounter files belonging at the
1140    // time to case [Planar 0] and case [Photo D].
1141    // Alas, this was only theory ! Because in practice some odd (read ill
1142    // formated Dicom) files (e.g. gdcmData/US-PAL-8-10x-echo.dcm) we encounter:
1143    //     - "Planar Configuration" = 0,
1144    //     - "Photometric Interpretation" = "PALETTE COLOR".
1145    // Hence gdcm will use the folowing "heuristic" in order to be tolerant
1146    // towards Dicom-non-conformant files:
1147    //   << whatever the "Planar Configuration" value might be, a
1148    //      "Photometric Interpretation" set to "PALETTE COLOR" forces
1149    //      a LUT intervention >>
1150    //
1151    // Now we are left with the following handling of the cases:
1152    // - [Planar 0] OR  [Photo A] no color decoding (since respectively
1153    //       Pixels are already RGB and monochrome pictures have no color :),
1154    // - [Planar 1] AND [Photo B] handled with ConvertRGBPlanesToRGBPixels()
1155    // - [Planar 1] AND [Photo C] handled with ConvertYcBcRPlanesToRGBPixels()
1156    // - [Planar 2] OR  [Photo D] requires LUT intervention.
1157
1158    gdcmWarningMacro("--> ConvertHandleColor"
1159                     << "Planar Configuration " << PlanarConfiguration );
1160
1161    if ( ! IsRawRGB() )
1162    {
1163       // [Planar 2] OR  [Photo D]: LUT intervention done outside
1164       gdcmWarningMacro("--> RawRGB : LUT intervention done outside");
1165       return;
1166    }
1167
1168    if ( PlanarConfiguration == 1 )
1169    {
1170       if ( IsYBRFull )
1171       {
1172          // [Planar 1] AND [Photo C] (remember YBR_FULL_422 acts as RGB)
1173          gdcmWarningMacro("--> YBRFull");
1174          ConvertYcBcRPlanesToRGBPixels();
1175       }
1176       else
1177       {
1178          // [Planar 1] AND [Photo C]
1179          gdcmWarningMacro("--> YBRFull");
1180          ConvertRGBPlanesToRGBPixels();
1181       }
1182       return;
1183    }
1184
1185    // When planarConf is 0, and RLELossless (forbidden by Dicom norm)
1186    // pixels need to be RGB-fyied anyway
1187
1188    if (IsRLELossless)
1189    {
1190      gdcmWarningMacro("--> RLE Lossless");
1191      ConvertRGBPlanesToRGBPixels();
1192    }
1193
1194    // In *normal *case, when planarConf is 0, pixels are already in RGB
1195 }
1196
1197 /// Computes the Pixels Size
1198 void PixelReadConvert::ComputeRawAndRGBSizes()
1199 {
1200    int bitsAllocated = BitsAllocated;
1201    // Number of "Bits Allocated" is fixed to 16 when it's 12, since
1202    // in this case we will expand the image to 16 bits (see
1203    //    \ref ReadAndDecompress12BitsTo16Bits() )
1204    if (  BitsAllocated == 12 )
1205    {
1206       bitsAllocated = 16;
1207    }
1208
1209    RawSize =  XSize * YSize * ZSize
1210                      * ( bitsAllocated / 8 )
1211                      * SamplesPerPixel;
1212    if ( HasLUT )
1213    {
1214       RGBSize = 3 * RawSize; // works for 8 and 16 bits per Pixel
1215    }
1216    else
1217    {
1218       RGBSize = RawSize;
1219    }
1220 }
1221
1222 /// Allocates room for RGB Pixels
1223 void PixelReadConvert::AllocateRGB()
1224 {
1225   if ( RGB )
1226      delete [] RGB;
1227   RGB = new uint8_t[RGBSize];
1228 }
1229
1230 /// Allocates room for RAW Pixels
1231 void PixelReadConvert::AllocateRaw()
1232 {
1233   if ( Raw )
1234      delete [] Raw;
1235   Raw = new uint8_t[RawSize];
1236 }
1237
1238 //-----------------------------------------------------------------------------
1239 // Print
1240 /**
1241  * \brief        Print self.
1242  * @param indent Indentation string to be prepended during printing.
1243  * @param os     Stream to print to.
1244  */
1245 void PixelReadConvert::Print( std::ostream &os, std::string const &indent )
1246 {
1247    os << indent
1248       << "--- Pixel information -------------------------"
1249       << std::endl;
1250    os << indent
1251       << "Pixel Data: offset " << PixelOffset
1252       << " x(" << std::hex << PixelOffset << std::dec
1253       << ")   length " << PixelDataLength
1254       << " x(" << std::hex << PixelDataLength << std::dec
1255       << ")" << std::endl;
1256
1257    if ( IsRLELossless )
1258    {
1259       if ( RLEInfo )
1260       {
1261          RLEInfo->Print( os, indent );
1262       }
1263       else
1264       {
1265          gdcmWarningMacro("Set as RLE file but NO RLEinfo present.");
1266       }
1267    }
1268
1269    if ( IsJPEG2000 || IsJPEGLossless || IsJPEGLossy || IsJPEGLS )
1270    {
1271       if ( JPEGInfo )
1272       {
1273          JPEGInfo->Print( os, indent );
1274       }
1275       else
1276       {
1277          gdcmWarningMacro("Set as JPEG file but NO JPEGinfo present.");
1278       }
1279    }
1280 }
1281
1282 //-----------------------------------------------------------------------------
1283 } // end namespace gdcm
1284
1285 // NOTES on File internal calls
1286 // User
1287 // ---> GetImageData
1288 //     ---> GetImageDataIntoVector
1289 //        |---> GetImageDataIntoVectorRaw
1290 //        | lut intervention
1291 // User
1292 // ---> GetImageDataRaw
1293 //     ---> GetImageDataIntoVectorRaw
1294