src/gdcmPixelReadConvert.cxx

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