gdcmPixelReadConvert.cxx

Go to the documentation of this file.
00001 /*=========================================================================
00002  
00003   Program:   gdcm
00004   Module:    $RCSfile: gdcmPixelReadConvert.cxx,v $
00005   Language:  C++
00006   Date:      $Date: 2007/07/25 10:35:20 $
00007   Version:   $Revision: 1.116 $
00008                                                                                 
00009   Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
00010   l'Image). All rights reserved. See Doc/License.txt or
00011   http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
00012                                                                                 
00013      This software is distributed WITHOUT ANY WARRANTY; without even
00014      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
00015      PURPOSE.  See the above copyright notices for more information.
00016                                                                                 
00017 =========================================================================*/
00018 
00019 #include "gdcmPixelReadConvert.h"
00020 #include "gdcmDebug.h"
00021 #include "gdcmFile.h"
00022 #include "gdcmGlobal.h"
00023 #include "gdcmTS.h"
00024 #include "gdcmDocEntry.h"
00025 #include "gdcmRLEFramesInfo.h"
00026 #include "gdcmJPEGFragmentsInfo.h"
00027 
00028 #include <fstream>
00029 #include <stdio.h> //for sscanf
00030 
00031 #if defined(__BORLANDC__)
00032    #include <mem.h> // for memset
00033 #endif 
00034 
00035 namespace GDCM_NAME_SPACE
00036 {
00037 
00038 //bool ReadMPEGFile (std::ifstream *fp, char *inputdata, size_t lenght); 
00039 bool gdcm_read_JPEG2000_file (void* raw, 
00040                               char *inputdata, size_t inputlength);
00041 //-----------------------------------------------------------------------------
00042 #define str2num(str, typeNum) *((typeNum *)(str))
00043 
00044 //-----------------------------------------------------------------------------
00045 // Constructor / Destructor
00047 PixelReadConvert::PixelReadConvert() 
00048 {
00049    RGB          = 0;
00050    RGBSize      = 0;
00051    Raw          = 0;
00052    RawSize      = 0;
00053    LutRGBA      = 0;
00054    LutRedData   = 0;
00055    LutGreenData = 0;
00056    LutBlueData  = 0;
00057    RLEInfo      = 0;
00058    JPEGInfo     = 0;
00059    UserFunction = 0;
00060    FileInternal = 0;
00061 }
00062 
00064 PixelReadConvert::~PixelReadConvert() 
00065 {
00066    Squeeze();
00067 }
00068 
00069 //-----------------------------------------------------------------------------
00070 // Public
00075 bool PixelReadConvert::IsRawRGB()
00076 {
00077    if (   IsMonochrome
00078        || PlanarConfiguration == 2
00079        || IsPaletteColor )
00080    {
00081       return false;
00082    }
00083    return true;
00084 }
00090 void PixelReadConvert::GrabInformationsFromFile( File *file,
00091                                                  FileHelper *fileHelper )
00092 {
00093    // Number of Bits Allocated for storing a Pixel is defaulted to 16
00094    // when absent from the file.
00095    BitsAllocated = file->GetBitsAllocated();
00096    if ( BitsAllocated == 0 )
00097    {
00098       BitsAllocated = 16;
00099    }
00100 
00101    else if (BitsAllocated > 8 && BitsAllocated < 16 && BitsAllocated != 12)
00102    {
00103       BitsAllocated = 16;
00104    }   
00105    // Number of "Bits Stored", defaulted to number of "Bits Allocated"
00106    // when absent from the file.
00107    BitsStored = file->GetBitsStored();
00108    if ( BitsStored == 0 )
00109    {
00110       BitsStored = BitsAllocated;
00111    }
00112 
00113    // High Bit Position, defaulted to "Bits Allocated" - 1
00114    HighBitPosition = file->GetHighBitPosition();
00115    if ( HighBitPosition == 0 )
00116    {
00117       HighBitPosition = BitsAllocated - 1;
00118    }
00119 
00120    XSize           = file->GetXSize();
00121    YSize           = file->GetYSize();
00122    ZSize           = file->GetZSize();
00123    TSize           = file->GetTSize();   
00124    SamplesPerPixel = file->GetSamplesPerPixel();
00125    //PixelSize       = file->GetPixelSize();  Useless
00126    PixelSign       = file->IsSignedPixelData();
00127    SwapCode        = file->GetSwapCode();
00128 
00129    IsPrivateGETransferSyntax = IsMPEG
00130              = IsJPEG2000 = IsJPEGLS = IsJPEGLossy  
00131              = IsJPEGLossless = IsRLELossless 
00132              = false;
00133 
00134    if (! file->IsDicomV3() )  // Should be ACR-NEMA file
00135    {
00136       IsRaw = true;
00137    }
00138    else
00139    {
00140       std::string ts = file->GetTransferSyntax();
00141 
00142       IsRaw = false;
00143       while (true) // shorter to write than 'if elseif elseif elseif' ...
00144       {
00145          // mind the order : check the most usual first.
00146          if( IsRaw = Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ExplicitVRLittleEndian)         break;
00147          if( IsRaw = Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ImplicitVRLittleEndian )        break;
00148          if( IsRaw = Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ExplicitVRBigEndian)            break;
00149          if( IsRaw = Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ImplicitVRBigEndianPrivateGE)   break;
00150          // DeflatedExplicitVRLittleEndian syntax means the whole Dataset (Header + Pixels) is compressed !
00151          // Not dealt with ! (Parser hangs)
00152          //if( IsRaw = Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::DeflatedExplicitVRLittleEndian) break;
00153          break;
00154       }
00155       // cache whether this is a strange GE transfer syntax (which uses
00156       // a little endian transfer syntax for the header and a big endian
00157       // transfer syntax for the pixel data). 
00158       IsPrivateGETransferSyntax = 
00159                 ( Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ImplicitVRBigEndianPrivateGE );
00160 
00161       IsMPEG =  IsJPEG2000 =  IsJPEGLS =  IsJPEGLossy =  IsJPEGLossless = IsRLELossless = false;  
00162       if (!IsRaw)
00163       {     
00164          while(true)
00165          {
00166             // mind the order : check the most usual first.
00167             if( IsJPEGLossy     = Global::GetTS()->IsJPEGLossy(ts) )    break;
00168             if( IsJPEGLossless  = Global::GetTS()->IsJPEGLossless(ts) ) break;
00169             if( IsRLELossless   = Global::GetTS()->IsRLELossless(ts) )  break;
00170             if( IsJPEG2000      = Global::GetTS()->IsJPEG2000(ts) )     break;
00171             if( IsMPEG          = Global::GetTS()->IsMPEG(ts) )         break;
00172             if( IsJPEGLS        = Global::GetTS()->IsJPEGLS(ts) )       break;
00173             // DeflatedExplicitVRLittleEndian is considered as 'Unexpected' (we don't know yet how to process !)
00174             gdcmWarningMacro("Unexpected Transfer Syntax :[" << ts << "]");
00175             break;
00176          } 
00177       }
00178    }
00179 
00180    PixelOffset     = file->GetPixelOffset();
00181    PixelDataLength = file->GetPixelAreaLength();
00182    RLEInfo         = file->GetRLEInfo();
00183    JPEGInfo        = file->GetJPEGInfo();
00184 
00185    IsMonochrome    = file->IsMonochrome();
00186    IsMonochrome1   = file->IsMonochrome1();
00187    IsPaletteColor  = file->IsPaletteColor();
00188    IsYBRFull       = file->IsYBRFull();
00189 
00190    PlanarConfiguration = file->GetPlanarConfiguration();
00191 
00193    // LUT section:
00194    HasLUT = file->HasLUT();
00195    if ( HasLUT )
00196    {
00197       // Just in case some access to a File element requires disk access.
00198       LutRedDescriptor   = file->GetEntryString( 0x0028, 0x1101 );
00199       LutGreenDescriptor = file->GetEntryString( 0x0028, 0x1102 );
00200       LutBlueDescriptor  = file->GetEntryString( 0x0028, 0x1103 );
00201    
00202       // FIXME : The following comment is probabely meaningless, since LUT are *always*
00203       // loaded at parsing time, whatever their length is.
00204          
00205       // Depending on the value of Document::MAX_SIZE_LOAD_ELEMENT_VALUE
00206       // [ refer to invocation of Document::SetMaxSizeLoadEntry() in
00207       // Document::Document() ], the loading of the value (content) of a
00208       // [Bin|Val]Entry occurence migth have been hindered (read simply NOT
00209       // loaded). Hence, we first try to obtain the LUTs data from the file
00210       // and when this fails we read the LUTs data directly from disk.
00211       // \TODO Reading a [Bin|Val]Entry directly from disk is a kludge.
00212       //       We should NOT bypass the [Bin|Val]Entry class. Instead
00213       //       an access to an UNLOADED content of a [Bin|Val]Entry occurence
00214       //       (e.g. DataEntry::GetBinArea()) should force disk access from
00215       //       within the [Bin|Val]Entry class itself. The only problem
00216       //       is that the [Bin|Val]Entry is unaware of the FILE* is was
00217       //       parsed from. Fix that. FIXME.
00218    
00219       // //// Red round
00220       file->LoadEntryBinArea(0x0028, 0x1201);
00221       LutRedData = (uint8_t*)file->GetEntryBinArea( 0x0028, 0x1201 );
00222       if ( ! LutRedData )
00223       {
00224          gdcmWarningMacro("Unable to read Red Palette Color Lookup Table data");
00225       }
00226 
00227       // //// Green round:
00228       file->LoadEntryBinArea(0x0028, 0x1202);
00229       LutGreenData = (uint8_t*)file->GetEntryBinArea(0x0028, 0x1202 );
00230       if ( ! LutGreenData)
00231       {
00232          gdcmWarningMacro("Unable to read Green Palette Color Lookup Table data");
00233       }
00234 
00235       // //// Blue round:
00236       file->LoadEntryBinArea(0x0028, 0x1203);
00237       LutBlueData = (uint8_t*)file->GetEntryBinArea( 0x0028, 0x1203 );
00238       if ( ! LutBlueData )
00239       {
00240          gdcmWarningMacro("Unable to read Blue Palette Color Lookup Table data");
00241       }
00242    }
00243    FileInternal = file;   
00244    FH = fileHelper;
00245    ComputeRawAndRGBSizes();
00246 }
00247 
00249 bool PixelReadConvert::ReadAndDecompressPixelData( std::ifstream *fp )
00250 {
00251    // ComputeRawAndRGBSizes is already made by 
00252    // ::GrabInformationsFromfile. So, the structure sizes are
00253    // correct
00254    Squeeze();
00255 
00258    if ( !fp )
00259    {
00260       gdcmWarningMacro( "Unavailable file pointer." );
00261       return false;
00262    }
00263 
00264    fp->seekg( PixelOffset, std::ios::beg );
00265    if ( fp->fail() || fp->eof() )
00266    {
00267       gdcmWarningMacro( "Unable to find PixelOffset in file." );
00268       return false;
00269    }
00270 
00271    AllocateRaw();
00272 
00274    
00275    CallStartMethod(); // for progress bar
00276    unsigned int count = 0;
00277    unsigned int frameSize;
00278    unsigned int bitsAllocated = BitsAllocated;
00279    //if(bitsAllocated == 12)
00280    if(bitsAllocated > 8 && bitsAllocated < 16)
00281       bitsAllocated = 16;
00282    frameSize = XSize*YSize*SamplesPerPixel*bitsAllocated/8;
00283    
00285    
00286    if ( BitsAllocated == 12 ) // We suppose 'BitsAllocated' = 12 only exist for uncompressed files
00287    {
00288       ReadAndDecompress12BitsTo16Bits( fp);
00289    }
00290    else if ( IsRaw )
00291    {
00292       // This problem can be found when some obvious informations are found
00293       // after the field containing the image data. In this case, these
00294       // bad data are added to the size of the image (in the PixelDataLength
00295       // variable). But RawSize is the right size of the image !
00296       if ( PixelDataLength != RawSize )
00297       {
00298          gdcmWarningMacro( "Mismatch between PixelReadConvert : "
00299                             << PixelDataLength << " and RawSize : " << RawSize );
00300       }
00301       
00302       //todo : is it the right patch?
00303       char *raw = (char*)Raw;
00304       uint32_t remainingLength;
00305       unsigned int i; 
00306       unsigned int lengthToRead;
00307        
00308       if ( PixelDataLength > RawSize )
00309         lengthToRead =  RawSize;
00310       else
00311         lengthToRead = PixelDataLength;
00312  
00313       // perform a frame by frame reading
00314       remainingLength = lengthToRead;
00315       unsigned int nbFrames = lengthToRead / frameSize;
00316       for (i=0;i<nbFrames; i++)
00317       {
00318          Progress = (float)(count+1)/(float)nbFrames;
00319          fp->read( raw, frameSize);
00320          raw +=  frameSize;
00321          remainingLength -=  frameSize;
00322          count++;
00323       }
00324       if (remainingLength !=0 )
00325         fp->read( raw, remainingLength);
00326                  
00327       if ( fp->fail() || fp->eof())
00328       {
00329          gdcmWarningMacro( "Reading of Raw pixel data failed." );
00330          return false;
00331       }
00332    } 
00333    else if ( IsRLELossless )
00334    {
00335       if ( ! RLEInfo->DecompressRLEFile
00336                                ( fp, Raw, XSize, YSize, ZSize, TSize, BitsAllocated ) )
00337       {
00338          gdcmWarningMacro( "RLE decompressor failed." );
00339          return false;
00340       }
00341    }
00342    else if ( IsMPEG )
00343    {
00344       //gdcmWarningMacro( "Sorry, MPEG not yet taken into account" );
00345       //return false;
00346       // fp has already been seek to start of mpeg
00347       //ReadMPEGFile(fp, (char*)Raw, PixelDataLength); 
00348       return true;
00349    }
00350    else
00351    {
00352       // Default case concerns JPEG family
00353       if ( ! ReadAndDecompressJPEGFile( fp ) )
00354       {
00355          gdcmWarningMacro( "JPEG decompressor ( ReadAndDecompressJPEGFile()"
00356                               << " method ) failed." );
00357          return false;
00358       }
00359    }
00360 
00363    ConvertReorderEndianity();
00364    ConvertReArrangeBits();
00365    ConvertFixGreyLevels();
00366    if (UserFunction) // user is allowed to Mirror, TopDown, Rotate,...the image
00367       UserFunction( Raw, FileInternal);
00368    ConvertHandleColor();
00369 
00370    return true;
00371 }
00372 
00374 void PixelReadConvert::Squeeze() 
00375 {
00376    if ( RGB )
00377       delete [] RGB;
00378    RGB = 0;
00379 
00380    if ( Raw )
00381       delete [] Raw;
00382    Raw = 0;
00383 
00384    if ( LutRGBA )
00385       delete [] LutRGBA;
00386    LutRGBA = 0;
00387 }
00388 
00392 bool PixelReadConvert::BuildRGBImage()
00393 {
00394    if ( RGB )
00395    {
00396       // The job is already done.
00397       return true;
00398    }
00399 
00400    if ( ! Raw )
00401    {
00402       // The job can't be done
00403       return false;
00404    }
00405 
00406    BuildLUTRGBA();
00407    if ( ! LutRGBA )
00408    {
00409       // The job can't be done
00410       return false;
00411    }
00412 
00413    gdcmDebugMacro( "--> BuildRGBImage" );
00414                                                                                 
00415    // Build RGB Pixels
00416    AllocateRGB();
00417    
00418    int j;
00419    if ( BitsAllocated <= 8 )
00420    {
00421       uint8_t *localRGB = RGB;
00422       for (size_t i = 0; i < RawSize; ++i )
00423       {
00424          j  = Raw[i] * 4;
00425          *localRGB++ = LutRGBA[j];
00426          *localRGB++ = LutRGBA[j+1];
00427          *localRGB++ = LutRGBA[j+2];
00428       }
00429     }
00430  
00431     else  // deal with 16 bits pixels and 16 bits Palette color
00432     {
00433       uint16_t *localRGB = (uint16_t *)RGB;
00434       for (size_t i = 0; i < RawSize/2; ++i )
00435       {
00436          j  = ((uint16_t *)Raw)[i] * 4;
00437          *localRGB++ = ((uint16_t *)LutRGBA)[j];
00438          *localRGB++ = ((uint16_t *)LutRGBA)[j+1];
00439          *localRGB++ = ((uint16_t *)LutRGBA)[j+2];
00440       } 
00441     }
00442  
00443    return true;
00444 }
00445 
00446 //-----------------------------------------------------------------------------
00447 // Protected
00448 
00449 //-----------------------------------------------------------------------------
00450 // Private
00455 void PixelReadConvert::ReadAndDecompress12BitsTo16Bits( std::ifstream *fp )
00456                throw ( FormatError )
00457 {
00459    int nbPixels = XSize * YSize * TSize;
00460    uint16_t *localDecompres = (uint16_t*)Raw;
00461 
00462    for( int p = 0; p < nbPixels; p += 2 )
00463    {
00464       uint8_t b0, b1, b2;
00465 
00466       fp->read( (char*)&b0, 1);
00467       if ( fp->fail() || fp->eof() )
00468       {
00469          throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
00470                                 "Unfound first block" );
00471       }
00472 
00473       fp->read( (char*)&b1, 1 );
00474       if ( fp->fail() || fp->eof())
00475       {
00476          throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
00477                                 "Unfound second block" );
00478       }
00479 
00480       fp->read( (char*)&b2, 1 );
00481       if ( fp->fail() || fp->eof())
00482       {
00483          throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
00484                                 "Unfound second block" );
00485       }
00486 
00487       // Two steps are necessary to please VC++
00488       //
00489       // 2 pixels 12bit =     [0xABCDEF]
00490       // 2 pixels 16bit = [0x0ABD] + [0x0FCE]
00491       //                        A                     B                 D
00492       *localDecompres++ =  ((b0 >> 4) << 8) + ((b0 & 0x0f) << 4) + (b1 & 0x0f);
00493       //                        F                     C                 E
00494       *localDecompres++ =  ((b2 & 0x0f) << 8) + ((b1 >> 4) << 4) + (b2 >> 4);
00495 
00497    }
00498 }
00499 
00506 bool PixelReadConvert::ReadAndDecompressJPEGFile( std::ifstream *fp )
00507 {
00508    if ( IsJPEG2000 )
00509    {
00510      // make sure this is the right JPEG compression
00511      assert( !IsJPEGLossless || !IsJPEGLossy || !IsJPEGLS );
00512      // FIXME this is really ugly but it seems I have to load the complete
00513      // jpeg2000 stream to use jasper:
00514      // I don't think we'll ever be able to deal with multiple fragments properly
00515 
00516       unsigned long inputlength = 0;
00517       JPEGFragment *jpegfrag = JPEGInfo->GetFirstFragment();
00518       while( jpegfrag )
00519       {
00520          inputlength += jpegfrag->GetLength();
00521          jpegfrag = JPEGInfo->GetNextFragment();
00522       }
00523       gdcmAssertMacro( inputlength != 0);
00524       uint8_t *inputdata = new uint8_t[inputlength];
00525       char *pinputdata = (char*)inputdata;
00526       jpegfrag = JPEGInfo->GetFirstFragment();
00527       while( jpegfrag )
00528       {
00529          fp->seekg( jpegfrag->GetOffset(), std::ios::beg);
00530          fp->read(pinputdata, jpegfrag->GetLength());
00531          pinputdata += jpegfrag->GetLength();
00532          jpegfrag = JPEGInfo->GetNextFragment();
00533       }
00534       // Warning the inputdata buffer is delete in the function
00535       if ( ! gdcm_read_JPEG2000_file( Raw, 
00536           (char*)inputdata, inputlength ) )
00537       {
00538          return true;
00539       }
00540       // wow what happen, must be an error
00541       gdcmWarningMacro( "gdcm_read_JPEG2000_file() failed "); 
00542       return false;
00543    }
00544    else if ( IsJPEGLS )
00545    {
00546      // make sure this is the right JPEG compression
00547      assert( !IsJPEGLossless || !IsJPEGLossy || !IsJPEG2000 );
00548    // WARNING : JPEG-LS is NOT the 'classical' Jpeg Lossless : 
00549    // [JPEG-LS is the basis for new lossless/near-lossless compression
00550    // standard for continuous-tone images intended for JPEG2000. The standard
00551    // is based on the LOCO-I algorithm (LOw COmplexity LOssless COmpression
00552    // for Images) developed at Hewlett-Packard Laboratories]
00553    //
00554    // see http://datacompression.info/JPEGLS.shtml
00555    //
00556 #if 0
00557    std::cerr << "count:" << JPEGInfo->GetFragmentCount() << std::endl;
00558       unsigned long inputlength = 0;
00559       JPEGFragment *jpegfrag = JPEGInfo->GetFirstFragment();
00560       while( jpegfrag )
00561       {
00562          inputlength += jpegfrag->GetLength();
00563          jpegfrag = JPEGInfo->GetNextFragment();
00564       }
00565       gdcmAssertMacro( inputlength != 0);
00566       uint8_t *inputdata = new uint8_t[inputlength];
00567       char *pinputdata = (char*)inputdata;
00568       jpegfrag = JPEGInfo->GetFirstFragment();
00569       while( jpegfrag )
00570       {
00571          fp->seekg( jpegfrag->GetOffset(), std::ios::beg);
00572          fp->read(pinputdata, jpegfrag->GetLength());
00573          pinputdata += jpegfrag->GetLength();
00574          jpegfrag = JPEGInfo->GetNextFragment();
00575       }  
00576       
00577   //fp->read((char*)Raw, PixelDataLength);
00578 
00579   std::ofstream out("/tmp/jpegls.jpg");
00580   out.write((char*)inputdata, inputlength);
00581   out.close();
00582   delete[] inputdata;
00583 #endif
00584 
00585       gdcmWarningMacro( "Sorry, JPEG-LS not yet taken into account" );
00586       fp->seekg( JPEGInfo->GetFirstFragment()->GetOffset(), std::ios::beg);
00587 //    if ( ! gdcm_read_JPEGLS_file( fp,Raw ) )
00588          return false;
00589    }
00590    else
00591    {
00592      // make sure this is the right JPEG compression
00593      assert( !IsJPEGLS || !IsJPEG2000 );
00594      // Precompute the offset localRaw will be shifted with
00595      int length = XSize * YSize * ZSize * SamplesPerPixel;
00596      int numberBytes = BitsAllocated / 8;
00597 
00598      JPEGInfo->DecompressFromFile(fp, Raw, BitsStored, numberBytes, length );
00599      return true;
00600    }
00601 }
00602 
00620 void PixelReadConvert::BuildLUTRGBA()
00621 {
00622 
00623    // Note to code reviewers :
00624    // The problem is *much more* complicated, since a lot of manufacturers
00625    // Don't follow the norm :
00626    // have a look at David Clunie's remark at the end of this .cxx file.
00627    if ( LutRGBA )
00628    
00629    {
00630       return;
00631    }
00632    // Not so easy : see
00633    // http://www.barre.nom.fr/medical/dicom2/limitations.html#Color%20Lookup%20Tables
00634                                                                                 
00635    if ( ! IsPaletteColor )
00636    {
00637       return;
00638    }
00639                                                                                 
00640    if (   LutRedDescriptor   == GDCM_UNFOUND
00641        || LutGreenDescriptor == GDCM_UNFOUND
00642        || LutBlueDescriptor  == GDCM_UNFOUND )
00643    {
00644       gdcmWarningMacro( "(At least) a LUT Descriptor is missing" );
00645       return;
00646    }
00647 
00649    // Extract the info from the LUT descriptors
00650    int lengthR;   // Red LUT length in Bytes
00651    int debR;      // Subscript of the first Lut Value
00652    int nbitsR;    // Lut item size (in Bits)
00653    int nbRead;    // nb of items in LUT descriptor (must be = 3)
00654 
00655    nbRead = sscanf( LutRedDescriptor.c_str(),
00656                         "%d\\%d\\%d",
00657                         &lengthR, &debR, &nbitsR );
00658    if ( nbRead != 3 )
00659    {
00660       gdcmWarningMacro( "Wrong Red LUT descriptor" );
00661    }                                                                                
00662    int lengthG;  // Green LUT length in Bytes
00663    int debG;     // Subscript of the first Lut Value
00664    int nbitsG;   // Lut item size (in Bits)
00665 
00666    nbRead = sscanf( LutGreenDescriptor.c_str(),
00667                     "%d\\%d\\%d",
00668                     &lengthG, &debG, &nbitsG );  
00669    if ( nbRead != 3 )
00670    {
00671       gdcmWarningMacro( "Wrong Green LUT descriptor" );
00672    }
00673                                                                                 
00674    int lengthB;  // Blue LUT length in Bytes
00675    int debB;     // Subscript of the first Lut Value
00676    int nbitsB;   // Lut item size (in Bits)
00677    nbRead = sscanf( LutRedDescriptor.c_str(),
00678                     "%d\\%d\\%d",
00679                     &lengthB, &debB, &nbitsB );
00680    if ( nbRead != 3 )
00681    {
00682       gdcmWarningMacro( "Wrong Blue LUT descriptor" );
00683    }
00684  
00685    gdcmDebugMacro(" lengthR " << lengthR << " debR " 
00686                 << debR << " nbitsR " << nbitsR);
00687    gdcmDebugMacro(" lengthG " << lengthG << " debG " 
00688                 << debG << " nbitsG " << nbitsG);
00689    gdcmDebugMacro(" lengthB " << lengthB << " debB " 
00690                 << debB << " nbitsB " << nbitsB);
00691 
00692    if ( !lengthR ) // if = 2^16, this shall be 0 see : CP-143
00693       lengthR=65536;
00694    if ( !lengthG ) // if = 2^16, this shall be 0
00695       lengthG=65536;
00696    if ( !lengthB ) // if = 2^16, this shall be 0
00697       lengthB=65536; 
00698                                                                                 
00700 
00701    if ( ( ! LutRedData ) || ( ! LutGreenData ) || ( ! LutBlueData ) )
00702    {
00703       gdcmWarningMacro( "(At least) a LUT is missing" );
00704       return;
00705    }
00706 
00707    // -------------------------------------------------------------
00708    
00709    if ( BitsAllocated <= 8 )
00710    {
00711       // forge the 4 * 8 Bits Red/Green/Blue/Alpha LUT
00712       LutRGBA = new uint8_t[ 1024 ]; // 256 * 4 (R, G, B, Alpha)
00713       if ( !LutRGBA )
00714          return;
00715       LutItemNumber = 256;
00716       LutItemSize   = 8;
00717       memset( LutRGBA, 0, 1024 );
00718                                                                                 
00719       int mult;
00720       if ( ( nbitsR == 16 ) && ( BitsAllocated == 8 ) )
00721       {
00722          // when LUT item size is different than pixel size
00723          mult = 2; // high byte must be = low byte
00724       }
00725       else
00726       {
00727          // See PS 3.3-2003 C.11.1.1.2 p 619
00728          mult = 1;
00729       }
00730                                                                                 
00731       // if we get a black image, let's just remove the '+1'
00732       // from 'i*mult+1' and check again
00733       // if it works, we shall have to check the 3 Palettes
00734       // to see which byte is ==0 (first one, or second one)
00735       // and fix the code
00736       // We give up the checking to avoid some (useless ?) overhead
00737       // (optimistic asumption)
00738       int i;
00739       uint8_t *a;
00740 
00741       //take "Subscript of the first Lut Value" (debR,debG,debB) into account!
00742 
00743       //FIXME :  +1 : to get 'low value' byte
00744       //         Trouble expected on Big Endian Processors ?
00745       //         16 BIts Per Pixel Palette Color to be swapped?
00746 
00747       a = LutRGBA + 0 + debR;
00748       for( i=0; i < lengthR; ++i )
00749       {
00750          *a = LutRedData[i*mult+1]; 
00751          a += 4;
00752       }
00753                                                                                 
00754       a = LutRGBA + 1 + debG;
00755       for( i=0; i < lengthG; ++i)
00756       {
00757          *a = LutGreenData[i*mult+1];
00758          a += 4;
00759       }
00760                                                                                 
00761       a = LutRGBA + 2 + debB;
00762       for(i=0; i < lengthB; ++i)
00763       {
00764          *a = LutBlueData[i*mult+1];
00765          a += 4;
00766       }
00767                                     
00768       a = LutRGBA + 3 ;
00769       for(i=0; i < 256; ++i)
00770       {
00771          *a = 1; // Alpha component
00772          a += 4;
00773       }
00774    }
00775    else
00776    {
00777       // Probabely the same stuff is to be done for 16 Bits Pixels
00778       // with 65536 entries LUT ?!?
00779       // Still looking for accurate info on the web :-(
00780 
00781       gdcmWarningMacro( "Sorry Palette Color Lookup Tables not yet dealt with"
00782                          << " for 16 Bits Per Pixel images" );
00783 
00784       // forge the 4 * 16 Bits Red/Green/Blue/Alpha LUT
00785 
00786       LutRGBA = (uint8_t *)new uint16_t[ 65536*4 ]; // 2^16 * 4 (R, G, B, Alpha)
00787       if ( !LutRGBA )
00788          return;
00789       memset( LutRGBA, 0, 65536*4*2 );  // 16 bits = 2 bytes ;-)
00790 
00791       LutItemNumber = 65536;
00792       LutItemSize   = 16;
00793 
00794       int i;
00795       uint16_t *a16;
00796 
00797       //take "Subscript of the first Lut Value" (debR,debG,debB) into account!
00798 
00799       a16 = (uint16_t*)LutRGBA + 0 + debR;
00800       for( i=0; i < lengthR; ++i )
00801       {
00802          *a16 = ((uint16_t*)LutRedData)[i];
00803          a16 += 4;
00804       }
00805                                                                               
00806       a16 = (uint16_t*)LutRGBA + 1 + debG;
00807       for( i=0; i < lengthG; ++i)
00808       {
00809          *a16 = ((uint16_t*)LutGreenData)[i];
00810          a16 += 4;
00811       }
00812                                                                                 
00813       a16 = (uint16_t*)LutRGBA + 2 + debB;
00814       for(i=0; i < lengthB; ++i)
00815       {
00816          *a16 = ((uint16_t*)LutBlueData)[i];
00817          a16 += 4;
00818       }
00819                                                                              
00820       a16 = (uint16_t*)LutRGBA + 3 ;
00821       for(i=0; i < 65536; ++i)
00822       {
00823          *a16 = 1; // Alpha component
00824          a16 += 4;
00825       }
00826 /* Just to 'see' the LUT, at debug time
00827 // Don't remove this commented out code.
00828 
00829       a16=(uint16_t*)LutRGBA;
00830       for (int j=0;j<65536;j++)
00831       {
00832          std::cout << *a16     << " " << *(a16+1) << " "
00833                    << *(a16+2) << " " << *(a16+3) << std::endl;
00834          a16+=4;
00835       }
00836 */
00837    }
00838 }
00839 
00843 void PixelReadConvert::ConvertSwapZone()
00844 {
00845    unsigned int i;
00846 
00847    // If this file is 'ImplicitVR BigEndian PrivateGE Transfer Syntax', 
00848    // then the header is in little endian format and the pixel data is in 
00849    // big endian format.  When reading the header, GDCM has already established
00850    // a byte swapping code suitable for this machine to read the
00851    // header. In TS::ImplicitVRBigEndianPrivateGE, this code will need
00852    // to be switched in order to read the pixel data.  This must be
00853    // done REGARDLESS of the processor endianess!
00854    //
00855    // Example:  Assume we are on a little endian machine.  When
00856    // GDCM reads the header, the header will match the machine
00857    // endianess and the swap code will be established as a no-op.
00858    // When GDCM reaches the pixel data, it will need to switch the
00859    // swap code to do big endian to little endian conversion.
00860    //
00861    // Now, assume we are on a big endian machine.  When GDCM reads the
00862    // header, the header will be recognized as a different endianess
00863    // than the machine endianess, and a swap code will be established
00864    // to convert from little endian to big endian.  When GDCM readers
00865    // the pixel data, the pixel data endianess will now match the
00866    // machine endianess.  But we currently have a swap code that
00867    // converts from little endian to big endian.  In this case, we
00868    // need to switch the swap code to a no-op.
00869    //
00870    // Therefore, in either case, if the file is in
00871    // 'ImplicitVR BigEndian PrivateGE Transfer Syntax', then GDCM needs to switch
00872    // the byte swapping code when entering the pixel data.
00873    
00874    int tempSwapCode = SwapCode;
00875    if ( IsPrivateGETransferSyntax )
00876    {
00877       gdcmWarningMacro(" IsPrivateGETransferSyntax found; turn the SwapCode"); 
00878       // PrivateGETransferSyntax only exists for 'true' Dicom images
00879       // we assume there is no 'exotic' 32 bits endianess!
00880       if (SwapCode == 1234) 
00881       {
00882          tempSwapCode = 4321;
00883       }
00884       else if (SwapCode == 4321)
00885       {
00886          tempSwapCode = 1234;
00887       }
00888    }
00889     
00890    if ( BitsAllocated == 16 )
00891    {
00892       uint16_t *im16 = (uint16_t*)Raw;
00893       switch( tempSwapCode )
00894       {
00895          case 1234:
00896             break;
00897          case 3412:
00898          case 2143:
00899          case 4321:
00900             for( i = 0; i < RawSize / 2; i++ )
00901             {
00902                im16[i]= (im16[i] >> 8) | (im16[i] << 8 );
00903             }
00904             break;
00905          default:
00906             gdcmWarningMacro("SwapCode value (16 bits) not allowed." 
00907                         << tempSwapCode);
00908       }
00909    }
00910    else if ( BitsAllocated == 32 )
00911    {
00912       uint32_t s32;
00913       uint16_t high;
00914       uint16_t low;
00915       uint32_t *im32 = (uint32_t*)Raw;
00916       switch ( tempSwapCode )
00917       {
00918          case 1234:
00919             break;
00920          case 4321:
00921             for( i = 0; i < RawSize / 4; i++ )
00922             {
00923                low     = im32[i] & 0x0000ffff;  // 4321
00924                high    = im32[i] >> 16;
00925                high    = ( high >> 8 ) | ( high << 8 );
00926                low     = ( low  >> 8 ) | ( low  << 8 );
00927                s32     = low;
00928                im32[i] = ( s32 << 16 ) | high;
00929             }
00930             break;
00931          case 2143:
00932             for( i = 0; i < RawSize / 4; i++ )
00933             {
00934                low     = im32[i] & 0x0000ffff;   // 2143
00935                high    = im32[i] >> 16;
00936                high    = ( high >> 8 ) | ( high << 8 );
00937                low     = ( low  >> 8 ) | ( low  << 8 );
00938                s32     = high;
00939                im32[i] = ( s32 << 16 ) | low;
00940             }
00941             break;
00942          case 3412:
00943             for( i = 0; i < RawSize / 4; i++ )
00944             {
00945                low     = im32[i] & 0x0000ffff; // 3412
00946                high    = im32[i] >> 16;
00947                s32     = low;
00948                im32[i] = ( s32 << 16 ) | high;
00949             }
00950             break;
00951          default:
00952             gdcmWarningMacro("SwapCode value (32 bits) not allowed." << tempSwapCode );
00953       }
00954    }
00955 }
00956 
00960 void PixelReadConvert::ConvertReorderEndianity()
00961 {
00962    if ( BitsAllocated != 8 )
00963    {
00964       ConvertSwapZone();
00965    }
00966 
00967    // Special kludge in order to deal with xmedcon broken images:
00968    if ( BitsAllocated == 16
00969      && BitsStored < BitsAllocated
00970      && !PixelSign )
00971    {
00972       int l = (int)( RawSize / ( BitsAllocated / 8 ) );
00973       uint16_t *deb = (uint16_t *)Raw;
00974       for(int i = 0; i<l; i++)
00975       {
00976          if ( *deb == 0xffff )
00977          {
00978            *deb = 0;
00979          }
00980          deb++;
00981       }
00982    }
00983 }
00984 
00989 void PixelReadConvert::ConvertFixGreyLevels()
00990 {
00991    if (!IsMonochrome1)
00992       return;
00993 
00994    uint32_t i; // to please M$VC6
00995    int16_t j;
00996 
00997    if (!PixelSign)
00998    {
00999       if ( BitsAllocated == 8 )
01000       {
01001          uint8_t *deb = (uint8_t *)Raw;
01002          for (i=0; i<RawSize; i++)      
01003          {
01004             *deb = 255 - *deb;
01005             deb++;
01006          }
01007          return;
01008       }
01009 
01010       if ( BitsAllocated == 16 )
01011       {
01012          uint16_t mask =1;
01013          for (j=0; j<BitsStored-1; j++)
01014          {
01015             mask = (mask << 1) +1; // will be fff when BitsStored=12
01016          }
01017 
01018          uint16_t *deb = (uint16_t *)Raw;
01019          for (i=0; i<RawSize/2; i++)      
01020          {
01021             *deb = mask - *deb;
01022             deb++;
01023          }
01024          return;
01025        }
01026    }
01027    else
01028    {
01029       if ( BitsAllocated == 8 )
01030       {
01031          uint8_t smask8 = 255;
01032          uint8_t *deb = (uint8_t *)Raw;
01033          for (i=0; i<RawSize; i++)      
01034          {
01035             *deb = smask8 - *deb;
01036             deb++;
01037          }
01038          return;
01039       }
01040       if ( BitsAllocated == 16 )
01041       {
01042          uint16_t smask16 = 65535;
01043          uint16_t *deb = (uint16_t *)Raw;
01044          for (i=0; i<RawSize/2; i++)      
01045          {
01046             *deb = smask16 - *deb;
01047             deb++;
01048          }
01049          return;
01050       }
01051    }
01052 }
01053 
01058 bool PixelReadConvert::ConvertReArrangeBits() throw ( FormatError )
01059 {
01060 
01061    if ( BitsStored != BitsAllocated )
01062    {
01063       int l = (int)( RawSize / ( BitsAllocated / 8 ) );
01064       if ( BitsAllocated == 16 )
01065       {
01066          // pmask : to mask the 'unused bits' (may contain overlays)
01067          uint16_t pmask = 0xffff;
01068          pmask = pmask >> ( BitsAllocated - BitsStored );
01069 
01070          uint16_t *deb = (uint16_t*)Raw;
01071 
01072          if ( !PixelSign )  // Pixels are unsigned
01073          {
01074             for(int i = 0; i<l; i++)
01075             {   
01076                *deb = (*deb >> (BitsStored - HighBitPosition - 1)) & pmask;
01077                deb++;
01078             }
01079          }
01080          else // Pixels are signed
01081          {
01082             // smask : to check the 'sign' when BitsStored != BitsAllocated
01083             uint16_t smask = 0x0001;
01084             smask = smask << ( 16 - (BitsAllocated - BitsStored + 1) );
01085             // nmask : to propagate sign bit on negative values
01086             int16_t nmask = (int16_t)0x8000;  
01087             nmask = nmask >> ( BitsAllocated - BitsStored - 1 );
01088  
01089             for(int i = 0; i<l; i++)
01090             {
01091                *deb = *deb >> (BitsStored - HighBitPosition - 1);
01092                if ( *deb & smask )
01093                {
01094                   *deb = *deb | nmask;
01095                }
01096                else
01097                {
01098                   *deb = *deb & pmask;
01099                }
01100                deb++;
01101             }
01102          }
01103       }
01104       else if ( BitsAllocated == 32 )
01105       {
01106          // pmask : to mask the 'unused bits' (may contain overlays)
01107          uint32_t pmask = 0xffffffff;
01108          pmask = pmask >> ( BitsAllocated - BitsStored );
01109 
01110          uint32_t *deb = (uint32_t*)Raw;
01111 
01112          if ( !PixelSign )
01113          {
01114             for(int i = 0; i<l; i++)
01115             {             
01116                *deb = (*deb >> (BitsStored - HighBitPosition - 1)) & pmask;
01117                deb++;
01118             }
01119          }
01120          else
01121          {
01122             // smask : to check the 'sign' when BitsStored != BitsAllocated
01123             uint32_t smask = 0x00000001;
01124             smask = smask >> ( 32 - (BitsAllocated - BitsStored +1 ));
01125             // nmask : to propagate sign bit on negative values
01126             int32_t nmask = 0x80000000;   
01127             nmask = nmask >> ( BitsAllocated - BitsStored -1 );
01128 
01129             for(int i = 0; i<l; i++)
01130             {
01131                *deb = *deb >> (BitsStored - HighBitPosition - 1);
01132                if ( *deb & smask )
01133                   *deb = *deb | nmask;
01134                else
01135                   *deb = *deb & pmask;
01136                deb++;
01137             }
01138          }
01139       }
01140       else
01141       {
01142          gdcmWarningMacro("Weird image (BitsAllocated !=8, 12, 16, 32)");
01143          throw FormatError( "Weird image !?" );
01144       }
01145    }
01146    return true;
01147 }
01148 
01153 void PixelReadConvert::ConvertRGBPlanesToRGBPixels()
01154 {
01155    gdcmWarningMacro("--> ConvertRGBPlanesToRGBPixels");
01156 
01157    uint8_t *localRaw = Raw;
01158    uint8_t *copyRaw = new uint8_t[ RawSize ];
01159    memmove( copyRaw, localRaw, RawSize );
01160 
01161    int l = XSize * YSize * ZSize;
01162 
01163    uint8_t *a = copyRaw;
01164    uint8_t *b = copyRaw + l;
01165    uint8_t *c = copyRaw + l + l;
01166 
01167    for (int j = 0; j < l; j++)
01168    {
01169       *(localRaw++) = *(a++);
01170       *(localRaw++) = *(b++);
01171       *(localRaw++) = *(c++);
01172    }
01173    delete[] copyRaw;
01174 }
01175 
01180 void PixelReadConvert::ConvertYcBcRPlanesToRGBPixels()
01181 {
01182   // Remarks for YBR newbees :
01183   // YBR_FULL works very much like RGB, i.e. three samples per pixel, 
01184   // just the color space is YCbCr instead of RGB. This is particularly useful
01185   // for doppler ultrasound where most of the image is grayscale 
01186   // (i.e. only populates the Y components) and Cb and Cr are mostly zero,
01187   // except for the few patches of color on the image.
01188   // On such images, RLE achieves a compression ratio that is much better 
01189   // than the compression ratio on an equivalent RGB image. 
01190  
01191    gdcmWarningMacro("--> ConvertYcBcRPlanesToRGBPixels");
01192    
01193    uint8_t *localRaw = Raw;
01194    uint8_t *copyRaw = new uint8_t[ RawSize ];
01195    memmove( copyRaw, localRaw, RawSize );
01196 
01197    // to see the tricks about YBR_FULL, YBR_FULL_422,
01198    // YBR_PARTIAL_422, YBR_ICT, YBR_RCT have a look at :
01199    // ftp://medical.nema.org/medical/dicom/final/sup61_ft.pdf
01200    // and be *very* affraid
01201    //
01202    
01204    int l        = XSize * YSize * TSize;
01205    int nbFrames = ZSize;
01206 
01207    uint8_t *a = copyRaw + 0;
01208    uint8_t *b = copyRaw + l;
01209    uint8_t *c = copyRaw + l+ l;
01210    int32_t R, G, B;
01211 
01217 
01218    for ( int i = 0; i < nbFrames; i++ )
01219    {
01220       for ( int j = 0; j < l; j++ )
01221       {
01222          R = 38142 *(*a-16) + 52298 *(*c -128);
01223          G = 38142 *(*a-16) - 26640 *(*c -128) - 12845 *(*b -128);
01224          B = 38142 *(*a-16) + 66093 *(*b -128);
01225 
01226          R = (R+16384)>>15;
01227          G = (G+16384)>>15;
01228          B = (B+16384)>>15;
01229 
01230          if (R < 0)   R = 0;
01231          if (G < 0)   G = 0;
01232          if (B < 0)   B = 0;
01233          if (R > 255) R = 255;
01234          if (G > 255) G = 255;
01235          if (B > 255) B = 255;
01236 
01237          *(localRaw++) = (uint8_t)R;
01238          *(localRaw++) = (uint8_t)G;
01239          *(localRaw++) = (uint8_t)B;
01240          a++;
01241          b++;
01242          c++;
01243       }
01244    }
01245    delete[] copyRaw;
01246 }
01247 
01252 
01253 void PixelReadConvert::ConvertHandleColor()
01254 {
01256    // Deal with the color decoding i.e. handle:
01257    //   - R, G, B planes (as opposed to RGB pixels)
01258    //   - YBR (various) encodings.
01259    //   - LUT[s] (or "PALETTE COLOR").
01260    //
01261    // The classification in the color decoding schema is based on the blending
01262    // of two Dicom tags values:
01263    // * "Photometric Interpretation" for which we have the cases:
01264    //  - [Photo A] MONOCHROME[1|2] pictures,
01265    //  - [Photo B] RGB or YBR_FULL_422 (which acts as RGB),
01266    //  - [Photo C] YBR_* (with the above exception of YBR_FULL_422)
01267    //  - [Photo D] "PALETTE COLOR" which indicates the presence of LUT[s].
01268    // * "Planar Configuration" for which we have the cases:
01269    //  - [Planar 0] 0 then Pixels are already RGB
01270    //  - [Planar 1] 1 then we have 3 planes : R, G, B,
01271    //  - [Planar 2] 2 then we have 1 gray Plane and 3 LUTs
01272    //
01273    // Now in theory, one could expect some coherence when blending the above
01274    // cases. For example we should not encounter files belonging at the
01275    // time to case [Planar 0] and case [Photo D].
01276    // Alas, this was only theory ! Because in practice some odd (read ill
01277    // formated Dicom) files (e.g. gdcmData/US-PAL-8-10x-echo.dcm) we encounter:
01278    //     - "Planar Configuration" = 0,
01279    //     - "Photometric Interpretation" = "PALETTE COLOR".
01280    // Hence gdcm will use the folowing "heuristic" in order to be tolerant
01281    // towards Dicom-non-conformant files:
01282    //   << whatever the "Planar Configuration" value might be, a
01283    //      "Photometric Interpretation" set to "PALETTE COLOR" forces
01284    //      a LUT intervention >>
01285    //
01286    // Now we are left with the following handling of the cases:
01287    // - [Planar 0] OR  [Photo A] no color decoding (since respectively
01288    //       Pixels are already RGB and monochrome pictures have no color :),
01289    // - [Planar 1] AND [Photo B] handled with ConvertRGBPlanesToRGBPixels()
01290    // - [Planar 1] AND [Photo C] handled with ConvertYcBcRPlanesToRGBPixels()
01291    // - [Planar 2] OR  [Photo D] requires LUT intervention.
01292 
01293    gdcmDebugMacro("--> ConvertHandleColor "
01294                      << "Planar Configuration " << PlanarConfiguration );
01295 
01296    if ( ! IsRawRGB() )
01297    {
01298       // [Planar 2] OR  [Photo D]: LUT intervention done outside
01299       gdcmDebugMacro("--> RawRGB : LUT intervention done outside");
01300       return;
01301    }
01302                                                                                 
01303    if ( PlanarConfiguration == 1 )
01304    {
01305       if ( IsYBRFull )
01306       {
01307          // [Planar 1] AND [Photo C] (remember YBR_FULL_422 acts as RGB)
01308          gdcmDebugMacro("--> YBRFull");
01309          ConvertYcBcRPlanesToRGBPixels();
01310       }
01311       else
01312       {
01313          // [Planar 1] AND [Photo C]
01314          gdcmDebugMacro("--> YBRFull");
01315          ConvertRGBPlanesToRGBPixels();
01316       }
01317       return;
01318    }
01319                                                                                 
01320    // When planarConf is 0, and RLELossless (forbidden by Dicom norm)
01321    // pixels need to be RGB-fyied anyway
01322 
01323    if (IsRLELossless)
01324    { 
01325      gdcmDebugMacro("--> RLE Lossless");
01326      ConvertRGBPlanesToRGBPixels();
01327    }
01328 
01329    // In *normal *case, when planarConf is 0, pixels are already in RGB
01330 }
01331 
01333 void PixelReadConvert::ComputeRawAndRGBSizes()
01334 {
01335    int bitsAllocated = BitsAllocated;
01336    // Number of "Bits Allocated" is fixed to 16 when it's 12, since
01337    // in this case we will expand the image to 16 bits (see
01338    //    \ref ReadAndDecompress12BitsTo16Bits() )
01339    if (  BitsAllocated == 12 )
01340    {
01341       bitsAllocated = 16;
01342    }
01343                                                                                 
01344    RawSize =  XSize * YSize * ZSize * TSize
01345                      * ( bitsAllocated / 8 )
01346                      * SamplesPerPixel;
01347    if ( HasLUT )
01348    {
01349       RGBSize = 3 * RawSize; // works for 8 and 16 bits per Pixel
01350    }
01351    else
01352    {
01353       RGBSize = RawSize;
01354       
01355    }
01356    RawSize += RawSize%2;
01357    RGBSize += RGBSize%2;
01358 }
01359 
01361 void PixelReadConvert::AllocateRGB()
01362 {
01363   if ( RGB )
01364      delete [] RGB;
01365   RGB = new uint8_t[RGBSize];
01366 }
01367 
01369 void PixelReadConvert::AllocateRaw()
01370 {
01371   if ( Raw )
01372      delete [] Raw;
01373   Raw = new uint8_t[RawSize];
01374 }
01375 
01376 //-----------------------------------------------------------------------------
01377 // Print
01383 void PixelReadConvert::Print( std::ostream &os, std::string const &indent )
01384 {
01385    os << indent
01386       << "--- Pixel information -------------------------"
01387       << std::endl;
01388    os << indent
01389       << "Pixel Data: offset " << PixelOffset
01390       << " x(" << std::hex << PixelOffset << std::dec
01391       << ")   length " << PixelDataLength
01392       << " x(" << std::hex << PixelDataLength << std::dec
01393       << ")" << std::endl;
01394 
01395    if ( IsRLELossless )
01396    {
01397       if ( RLEInfo )
01398       {
01399          RLEInfo->Print( os, indent );
01400       }
01401       else
01402       {
01403          gdcmWarningMacro("Set as RLE file but NO RLEinfo present.");
01404       }
01405    }
01406 
01407    if ( IsJPEG2000 || IsJPEGLossless || IsJPEGLossy || IsJPEGLS )
01408    {
01409       if ( JPEGInfo )
01410       {
01411          JPEGInfo->Print( os, indent );
01412       }
01413       else
01414       {
01415          gdcmWarningMacro("Set as JPEG file but NO JPEGinfo present.");
01416       }
01417    }
01418 }
01419 
01423 void PixelReadConvert::CallStartMethod()
01424 {
01425    Progress = 0.0f;
01426    Abort    = false;
01427    CommandManager::ExecuteCommand(FH,CMD_STARTPROGRESS);
01428 }
01429 
01433 void PixelReadConvert::CallProgressMethod()
01434 {
01435    CommandManager::ExecuteCommand(FH,CMD_PROGRESS);
01436 }
01437 
01441 void PixelReadConvert::CallEndMethod()
01442 {
01443    Progress = 1.0f;
01444    CommandManager::ExecuteCommand(FH,CMD_ENDPROGRESS);
01445 }
01446 
01447 
01448 //-----------------------------------------------------------------------------
01449 } // end namespace gdcm
01450 
01451 // Note to developpers :
01452 // Here is a very detailled post from David Clunie, on the troubles caused 
01453 // 'non standard' LUT and LUT description
01454 // We shall have to take it into accound in our code.
01455 // Some day ...
01456 
01457 
01458 /*
01459 Subject: Problem with VOI LUTs in Agfa and Fuji CR and GE DX images, was Re: VOI LUT issues
01460 Date: Sun, 06 Feb 2005 17:13:40 GMT
01461 From: David Clunie <dclunie@dclunie.com>
01462 Reply-To: dclunie@dclunie.com
01463 Newsgroups: comp.protocols.dicom
01464 References: <1107553502.040221.189550@o13g2000cwo.googlegroups.com>
01465 
01466 > THE LUT that comes with [my] image claims to be 16-bit, but none of the
01467 > values goes higher than 4095.  That being said, though, none of my
01468 > original pixel values goes higher than that, either.  I have read
01469 > elsewhere on this group that when that happens you are supposed to
01470 > adjust the LUT.  Can someone be more specific?  There was a thread from
01471 > 2002 where Marco and David were mentioning doing precisely that.
01472 >
01473 > Thanks
01474 >
01475 > -carlos rodriguez
01476 
01477 
01478 You have encountered the well known "we know what the standard says but
01479 we are going to ignore it and do what we have been doing for almost
01480 a decade regardless" CR vendor bug. Agfa started this, but they are not
01481 the only vendor doing this now; GE and Fuji may have joined the club.
01482 
01483 Sadly, one needs to look at the LUT Data, figure out what the maximum
01484 value actually encoded is, and find the next highest power of 2 (e.g.
01485 212 in this case), to figure out what the range of the data is
01486 supposed to be. I have assumed that if the maximum value in the LUT
01487 data is less than a power of 2 minus 1 (e.g. 0xebc) then the intent
01488 of the vendor was not to use the maximum available grayscale range
01489 of the display (e.g. the maximum is 0xfff in this case). An alternative
01490 would be to scale to the actual maximum rather than a power of two.
01491 
01492 Very irritating, and in theory not totally reliable if one really
01493 intended the full 16 bits and only used, say 15, but that is extremely
01494 unlikely since everything would be too dark, and this heuristic
01495 seems to work OK.
01496 
01497 There has never been anything in the standard that describes having
01498 to go through these convolutions. Since the only value in the
01499 standard that describes the bit depth of the LUT values is LUT
01500 Descriptor value 3 and that is (usually) always required to be
01501 either 8 or 16, it mystifies me how the creators' of these images
01502 imagine that the receiver is going to divine the range that is intended. Further, the standard is quite explicit that this 3rd
01503 value defines the range of LUT values, but as far as I am aware, all
01504 the vendors are ignoring the standard and indeed sending a third value
01505 of 16 in all cases.
01506 
01507 This problem is not confined to CR, and is also seen with DX products.
01508 
01509 Typically I have seen:
01510 
01511 - Agfa CR, which usually (always ?) sends LUTs, values up to 0x0fff
01512 - Fuji CR, which occasionally send LUTs, values up to 0x03ff
01513 - GE DX, for presentation, which always have LUTs, up to 0x3fff
01514 
01515 Swissray, Siemens, Philips, Canon and Kodak never seem to send VOI LUTs
01516 at this point (which is a whole other problem). Note that the presence
01517 or absence of a VOI LUT as opposed to window values may be configurable
01518 on the modality in some cases, and I have just looked at what I happen
01519 to have received from a myriad of sites over whose configuration I have
01520 no control. This may be why the majority of Fuji images have no VOI LUTs,
01521 but a few do (or it may be the Siemens system that these Fuji images went
01522 through that perhaps added it). I do have some test Hologic DX images that
01523 are not from a clinical site that do actually get this right (a value
01524 of 12 for the third value and a max of 0xfff).
01525 
01526 Since almost every vendor that I have encountered that encodes LUTs
01527 makes this mistake, perhaps it is time to amend the standard to warn
01528 implementor's of receivers and/or sanction this bad behavior. We have
01529 talked about this in the past in WG 6 but so far everyone has been
01530 reluctant to write into the standard such a comment. Maybe it is time
01531 to try again, since if one is not aware of this problem, one cannot
01532 effectively implement display using VOI LUTs, and there is a vast
01533 installed base to contend with.
01534 
01535 I did not check presentation states, in which VOI LUTs could also be
01536 encountered, for the prevalence of this mistake, nor did I look at the
01537 encoding of Modality LUT's, which are unusual. Nor did I check digital
01538 mammography images. I would be interested to hear from anyone who has.
01539 
01540 David
01541 
01542 PS. The following older thread in this newsgroup discusses this:
01543 
01544 "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"
01545 
01546 PPS. From a historical perspective, the following may be of interest.
01547 
01548 In the original standard in 1993, all that was said about this was a
01549 reference to the corresponding such where Modality LUTs are described
01550 that said:
01551 
01552 "The third value specifies the number of bits for each entry in the
01553 LUT Data. It shall take the value 8 or 16. The LUT Data shall be stored
01554 in a format equivalent to 8 or 16 bits allocated and high bit equal
01555 1-bits allocated."
01556 
01557 Since the high bit hint was not apparently explicit enough, a very
01558 early CP, CP 15 (submitted by Agfa as it happens), replaced this with:
01559 
01560 "The third value conveys the range of LUT entry values. It shall take
01561 the value 8 or 16, corresponding with the LUT entry value range of
01562 256 or 65536.
01563 
01564 Note:    The third value is not required for describing the
01565     LUT data and is only included for informational usage
01566     and for maintaining compatibility with ACRNEMA 2.0.
01567 
01568 The LUT Data contains the LUT entry values."
01569 
01570 That is how it read in the 1996, 1998 and 1999 editions.
01571 
01572 By the 2000 edition, Supplement 33 that introduced presentation states
01573 extensively reworked this entire section and tried to explain this in
01574 different words:
01575 
01576 "The output range is from 0 to 2^n-1 where n is the third value of LUT
01577 Descriptor. This range is always unsigned."
01578 
01579 and also added a note to spell out what the output range meant in the
01580 VOI LUT section:
01581 
01582 "9. The output of the Window Center/Width or VOI LUT transformation
01583 is either implicitly scaled to the full range of the display device
01584 if there is no succeeding transformation defined, or implicitly scaled
01585 to the full input range of the succeeding transformation step (such as
01586 the Presentation LUT), if present. See C.11.6.1."
01587 
01588 It still reads this way in the 2004 edition.
01589 
01590 Note that LUTs in other applications than the general VOI LUT allow for
01591 values other than 8 or 16 in the third value of LUT descriptor to permit
01592 ranges other than 0 to 255 or 65535.
01593 
01594 In addition, the DX Image Module specializes the VOI LUT
01595 attributes as follows, in PS 3.3 section C.8.11.3.1.5 (added in Sup 32):
01596 
01597 "The third value specifies the number of bits for each entry in the LUT
01598 Data (analogous to ìbits storedî). It shall be between 10-16. The LUT
01599 Data shall be stored in a format equivalent to 16 ìbits allocatedî and
01600 ìhigh bitî equal to ìbits storedî - 1. The third value conveys the range
01601 of LUT entry values. These unsigned LUT entry values shall range between
01602 0 and 2^n-1, where n is the third value of the LUT Descriptor."
01603 
01604 So in the case of the GE DX for presentation images, the third value of
01605 LUT descriptor is allowed to be and probably should be 14 rather than 16.
01606 
01607 */

Generated on Fri Aug 24 12:53:18 2007 for gdcm by  doxygen 1.4.6