src/gdcmPixelReadConvert.cxx

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