src/gdcmPixelReadConvert.cxx

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