X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=src%2FgdcmPixelReadConvert.cxx;h=45aaf812dd75df17a02bc981bebf4232b6f29c06;hb=56ad78cdf3b5de32f0a94a38cd32ddc9f403ffc3;hp=35d29222c50ae044182c4c986f82939e233d3959;hpb=25531920a24515b24ddaec0aff11b976d1a4a7cf;p=gdcm.git diff --git a/src/gdcmPixelReadConvert.cxx b/src/gdcmPixelReadConvert.cxx index 35d29222..45aaf812 100644 --- a/src/gdcmPixelReadConvert.cxx +++ b/src/gdcmPixelReadConvert.cxx @@ -1,10 +1,10 @@ /*========================================================================= - + Program: gdcm Module: $RCSfile: gdcmPixelReadConvert.cxx,v $ Language: C++ - Date: $Date: 2005/02/02 16:34:55 $ - Version: $Revision: 1.45 $ + Date: $Date: 2008/05/20 09:22:03 $ + Version: $Revision: 1.128 $ Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de l'Image). All rights reserved. See Doc/License.txt or @@ -16,26 +16,36 @@ =========================================================================*/ +#include "gdcmPixelReadConvert.h" #include "gdcmDebug.h" #include "gdcmFile.h" #include "gdcmGlobal.h" #include "gdcmTS.h" -#include "gdcmPixelReadConvert.h" #include "gdcmDocEntry.h" #include "gdcmRLEFramesInfo.h" #include "gdcmJPEGFragmentsInfo.h" +#include "gdcmSegmentedPalette.h" #include #include //for sscanf -namespace gdcm +#if defined(__BORLANDC__) + #include // for memset +#endif + +namespace GDCM_NAME_SPACE { + +//bool ReadMPEGFile (std::ifstream *fp, char *inputdata, size_t lenght); +bool gdcm_read_JPEG2000_file (void* raw, + char *inputdata, size_t inputlength); //----------------------------------------------------------------------------- #define str2num(str, typeNum) *((typeNum *)(str)) //----------------------------------------------------------------------------- // Constructor / Destructor -PixelReadConvert::PixelReadConvert() +/// Constructor +PixelReadConvert::PixelReadConvert() { RGB = 0; RGBSize = 0; @@ -45,8 +55,13 @@ PixelReadConvert::PixelReadConvert() LutRedData = 0; LutGreenData = 0; LutBlueData = 0; + RLEInfo = 0; + JPEGInfo = 0; + UserFunction = 0; + FileInternal = 0; } +/// Canonical Destructor PixelReadConvert::~PixelReadConvert() { Squeeze(); @@ -55,8 +70,8 @@ PixelReadConvert::~PixelReadConvert() //----------------------------------------------------------------------------- // Public /** - * \brief Predicate to know wether the image[s] (once Raw) is RGB. - * \note See comments of \ref ConvertHandleColor + * \brief Predicate to know whether the image[s] (once Raw) is RGB. + * \note See comments of ConvertHandleColor */ bool PixelReadConvert::IsRawRGB() { @@ -68,8 +83,13 @@ bool PixelReadConvert::IsRawRGB() } return true; } - -void PixelReadConvert::GrabInformationsFromFile( File *file ) +/** + * \brief Gets various usefull informations from the file header + * @param file gdcm::File pointer + * @param fileHelper gdcm::FileHelper pointer + */ +void PixelReadConvert::GrabInformationsFromFile( File *file, + FileHelper *fileHelper ) { // Number of Bits Allocated for storing a Pixel is defaulted to 16 // when absent from the file. @@ -78,6 +98,10 @@ void PixelReadConvert::GrabInformationsFromFile( File *file ) { BitsAllocated = 16; } + else if ( BitsAllocated > 8 && BitsAllocated < 16 && BitsAllocated != 12 ) + { + BitsAllocated = 16; + } // Number of "Bits Stored", defaulted to number of "Bits Allocated" // when absent from the file. @@ -94,90 +118,199 @@ void PixelReadConvert::GrabInformationsFromFile( File *file ) HighBitPosition = BitsAllocated - 1; } - XSize = file->GetXSize(); - YSize = file->GetYSize(); - ZSize = file->GetZSize(); + XSize = file->GetXSize(); + YSize = file->GetYSize(); + ZSize = file->GetZSize(); + TSize = file->GetTSize(); SamplesPerPixel = file->GetSamplesPerPixel(); - PixelSize = file->GetPixelSize(); - PixelSign = file->IsSignedPixelData(); - SwapCode = file->GetSwapCode(); - std::string ts = file->GetTransferSyntax(); - IsRaw = - ( ! file->IsDicomV3() ) - || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ImplicitVRLittleEndian - || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ImplicitVRLittleEndianDLXGE - || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ExplicitVRLittleEndian - || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ExplicitVRBigEndian - || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::DeflatedExplicitVRLittleEndian; - - IsJPEG2000 = Global::GetTS()->IsJPEG2000(ts); - IsJPEGLS = Global::GetTS()->IsJPEGLS(ts); - IsJPEGLossy = Global::GetTS()->IsJPEGLossy(ts); - IsJPEGLossless = Global::GetTS()->IsJPEGLossless(ts); - IsRLELossless = Global::GetTS()->IsRLELossless(ts); + //PixelSize = file->GetPixelSize(); Useless + PixelSign = file->IsSignedPixelData(); + SwapCode = file->GetSwapCode(); + + IsPrivateGETransferSyntax = IsMPEG + = IsJPEG2000 = IsJPEGLS = IsJPEGLossy + = IsJPEGLossless = IsRLELossless + = false; + + if (! file->IsDicomV3() ) // Should be ACR-NEMA file + { + IsRaw = true; + } + else + { + std::string ts = file->GetTransferSyntax(); + + IsRaw = false; + while (true) // shorter to write than 'if elseif elseif elseif' ... + { + // mind the order : check the most usual first. + if( IsRaw = (Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ExplicitVRLittleEndian)) break; + if( IsRaw = (Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ImplicitVRLittleEndian)) break; + if( IsRaw = (Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ExplicitVRBigEndian)) break; + if( IsRaw = (Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ImplicitVRBigEndianPrivateGE)) break; + // DeflatedExplicitVRLittleEndian syntax means the whole Dataset (Header + Pixels) is compressed ! + // Not dealt with ! (Parser hangs) + //if( IsRaw = Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::DeflatedExplicitVRLittleEndian) break; + break; + } + // cache whether this is a strange GE transfer syntax (which uses + // a little endian transfer syntax for the header and a big endian + // transfer syntax for the pixel data). + IsPrivateGETransferSyntax = + ( Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ImplicitVRBigEndianPrivateGE ); + + IsMPEG = IsJPEG2000 = IsJPEGLS = IsJPEGLossy = IsJPEGLossless = IsRLELossless = false; + if (!IsRaw) + { + while(true) + { + // mind the order : check the most usual first. + if( IsJPEGLossy = (Global::GetTS()->IsJPEGLossy(ts))) break; + if( IsJPEGLossless = (Global::GetTS()->IsJPEGLossless(ts))) break; + if( IsRLELossless = (Global::GetTS()->IsRLELossless(ts))) break; + if( IsJPEG2000 = (Global::GetTS()->IsJPEG2000(ts))) break; + if( IsMPEG = (Global::GetTS()->IsMPEG(ts))) break; + if( IsJPEGLS = (Global::GetTS()->IsJPEGLS(ts))) break; + // DeflatedExplicitVRLittleEndian is considered as 'Unexpected' + // (we don't know yet how to process !) + gdcmWarningMacro("Unexpected Transfer Syntax :[" << ts << "]"); + break; + } + } + } PixelOffset = file->GetPixelOffset(); PixelDataLength = file->GetPixelAreaLength(); - RLEInfo = file->GetRLEInfo(); - JPEGInfo = file->GetJPEGInfo(); + RLEInfo = file->GetRLEInfo(); + JPEGInfo = file->GetJPEGInfo(); + + IsMonochrome = file->IsMonochrome(); + IsMonochrome1 = file->IsMonochrome1(); + IsPaletteColor = file->IsPaletteColor(); + IsYBRFull = file->IsYBRFull(); PlanarConfiguration = file->GetPlanarConfiguration(); - IsMonochrome = file->IsMonochrome(); - IsPaletteColor = file->IsPaletteColor(); - IsYBRFull = file->IsYBRFull(); ///////////////////////////////////////////////////////////////// // LUT section: HasLUT = file->HasLUT(); if ( HasLUT ) { - // Just in case some access to a File element requires disk access. - LutRedDescriptor = file->GetEntryValue( 0x0028, 0x1101 ); - LutGreenDescriptor = file->GetEntryValue( 0x0028, 0x1102 ); - LutBlueDescriptor = file->GetEntryValue( 0x0028, 0x1103 ); - - // Depending on the value of Document::MAX_SIZE_LOAD_ELEMENT_VALUE - // [ refer to invocation of Document::SetMaxSizeLoadEntry() in - // Document::Document() ], the loading of the value (content) of a - // [Bin|Val]Entry occurence migth have been hindered (read simply NOT - // loaded). Hence, we first try to obtain the LUTs data from the file - // and when this fails we read the LUTs data directly from disk. - /// \TODO Reading a [Bin|Val]Entry directly from disk is a kludge. - /// We should NOT bypass the [Bin|Val]Entry class. Instead - /// an access to an UNLOADED content of a [Bin|Val]Entry occurence - /// (e.g. BinEntry::GetBinArea()) should force disk access from - /// within the [Bin|Val]Entry class itself. The only problem - /// is that the [Bin|Val]Entry is unaware of the FILE* is was - /// parsed from. Fix that. FIXME. - - ////// Red round - file->LoadEntryBinArea(0x0028, 0x1201); - LutRedData = (uint8_t*)file->GetEntryBinArea( 0x0028, 0x1201 ); - if ( ! LutRedData ) - { - gdcmVerboseMacro( "Unable to read Red LUT data" ); - } - - ////// Green round: - file->LoadEntryBinArea(0x0028, 0x1202); - LutGreenData = (uint8_t*)file->GetEntryBinArea(0x0028, 0x1202 ); - if ( ! LutGreenData) - { - gdcmVerboseMacro( "Unable to read Green LUT data" ); - } - - ////// Blue round: - file->LoadEntryBinArea(0x0028, 0x1203); - LutBlueData = (uint8_t*)file->GetEntryBinArea( 0x0028, 0x1203 ); - if ( ! LutBlueData ) - { - gdcmVerboseMacro( "Unable to read Blue LUT data" ); - } +/* + C.7.6.3.1.5 + The three values of Palette Color Lookup Table Descriptor (0028,1101-1103) + describe the format of the Lookup Table Data in the corresponding + Data Element (0028,1201-1203) or (0028,1221-1223). + + The first value is the number of entries in the lookup table. + When the number of table entries is equal to 2**16 then this value shall be 0. + + The second value is the first stored pixel value mapped. + This pixel value is mapped to the first entry in the Lookup Table Data. + All image pixel values less than the first entry value mapped are also + mapped to the first entry in the Lookup Table Data. + An image pixel value one greater than the first entry value mapped is + mapped to the second entry in the Lookup Table Data. + Subsequent image pixel values are mapped to the subsequent entries in + the Lookup Table Data up to an image pixel value equal to number of + entries + first entry value mapped - 1 which is mapped to the last entry + in the Lookup Table Data. + Image pixel values greater than or equal to number of entries + first entry + value mapped are also mapped to the last entry in the Lookup Table Data. + + The third value specifies the number of bits for each entry in the Lookup + Table Data. It shall take the value of 8 or 16. + The LUT Data shall be stored in a format equivalent to 8 or 16 bits + allocated where the high bit is equal to bits allocated-1. + + When the Palette Color Lookup Table Descriptor (0028,1101-1103) are used as + part of the Palette Color Lookup Table Module, the third value shall be + equal to 16. + + Note: A value of 16 indicates the Lookup Table Data will range from (0,0,0) + minimum intensity to (65535,65535,65535) maximum intensity. + +*/ + + // Just in case some access to a File element requires disk access. + LutRedDescriptor = file->GetEntryString( 0x0028, 0x1101 ); + LutGreenDescriptor = file->GetEntryString( 0x0028, 0x1102 ); + LutBlueDescriptor = file->GetEntryString( 0x0028, 0x1103 ); + // Is it a Segmented Palette ? Check if we find the red one: + if( file->GetDocEntry(0x0028,0x1221) ) // no need to check for blue & green + { + GDCM_NAME_SPACE::TagKey DCM_RedPaletteColorLookupTableDescriptor (0x0028, 0x1101); + GDCM_NAME_SPACE::TagKey DCM_GreenPaletteColorLookupTableDescriptor (0x0028, 0x1102); + GDCM_NAME_SPACE::TagKey DCM_BluePaletteColorLookupTableDescriptor (0x0028, 0x1103); + + GDCM_NAME_SPACE::TagKey DCM_SegmentedRedPaletteColorLookupTableData (0x0028, 0x1221); + GDCM_NAME_SPACE::TagKey DCM_SegmentedGreenPaletteColorLookupTableData (0x0028, 0x1222); + GDCM_NAME_SPACE::TagKey DCM_SegmentedBluePaletteColorLookupTableData (0x0028, 0x1223); + + + LutRedData = new uint8_t[65535*2]; // FIXME: leak + LutGreenData = new uint8_t[65535*2]; + LutBlueData = new uint8_t[65535*2]; + // TODO need to check file is indeed PALETTE COLOR: + ReadPaletteInto(file, DCM_RedPaletteColorLookupTableDescriptor, + DCM_SegmentedRedPaletteColorLookupTableData,LutRedData); + ReadPaletteInto(file, DCM_GreenPaletteColorLookupTableDescriptor, + DCM_SegmentedGreenPaletteColorLookupTableData,LutGreenData); + ReadPaletteInto(file, DCM_BluePaletteColorLookupTableDescriptor, + DCM_SegmentedBluePaletteColorLookupTableData,LutBlueData); + + } + else + { + + // FIXME : The following comment is probabely meaningless, since LUT are *always* + // loaded at parsing time, whatever their length is. + + // Depending on the value of Document::MAX_SIZE_LOAD_ELEMENT_VALUE + // [ refer to invocation of Document::SetMaxSizeLoadEntry() in + // Document::Document() ], the loading of the value (content) of a + // [Bin|Val]Entry occurence migth have been hindered (read simply NOT + // loaded). Hence, we first try to obtain the LUTs data from the file + // and when this fails we read the LUTs data directly from disk. + // \TODO Reading a [Bin|Val]Entry directly from disk is a kludge. + // We should NOT bypass the [Bin|Val]Entry class. Instead + // an access to an UNLOADED content of a [Bin|Val]Entry occurence + // (e.g. DataEntry::GetBinArea()) should force disk access from + // within the [Bin|Val]Entry class itself. The only problem + // is that the [Bin|Val]Entry is unaware of the FILE* is was + // parsed from. Fix that. FIXME. + + // //// Red round + file->LoadEntryBinArea(0x0028, 0x1201); + LutRedData = (uint8_t*)file->GetEntryBinArea( 0x0028, 0x1201 ); + if ( ! LutRedData ) + { + gdcmWarningMacro("Unable to read Red Palette Color Lookup Table data"); + } + + // //// Green round: + file->LoadEntryBinArea(0x0028, 0x1202); + LutGreenData = (uint8_t*)file->GetEntryBinArea(0x0028, 0x1202 ); + if ( ! LutGreenData) + { + gdcmWarningMacro("Unable to read Green Palette Color Lookup Table data"); + } + + // //// Blue round: + file->LoadEntryBinArea(0x0028, 0x1203); + LutBlueData = (uint8_t*)file->GetEntryBinArea( 0x0028, 0x1203 ); + if ( ! LutBlueData ) + { + gdcmWarningMacro("Unable to read Blue Palette Color Lookup Table data"); + } + } } - + FileInternal = file; + FH = fileHelper; ComputeRawAndRGBSizes(); } +/// \brief Reads from disk and decompresses Pixels bool PixelReadConvert::ReadAndDecompressPixelData( std::ifstream *fp ) { // ComputeRawAndRGBSizes is already made by @@ -189,22 +322,33 @@ bool PixelReadConvert::ReadAndDecompressPixelData( std::ifstream *fp ) //// First stage: get our hands on the Pixel Data. if ( !fp ) { - gdcmVerboseMacro( "Unavailable file pointer." ); + gdcmWarningMacro( "Unavailable file pointer." ); return false; } fp->seekg( PixelOffset, std::ios::beg ); - if( fp->fail() || fp->eof()) + if ( fp->fail() || fp->eof() ) { - gdcmVerboseMacro( "Unable to find PixelOffset in file." ); + gdcmWarningMacro( "Unable to find PixelOffset in file." ); return false; } AllocateRaw(); ////////////////////////////////////////////////// - //// Second stage: read from disk dans decompress. - if ( BitsAllocated == 12 ) + + CallStartMethod(); // for progress bar + unsigned int count = 0; + unsigned int frameSize; + unsigned int bitsAllocated = BitsAllocated; + //if(bitsAllocated == 12) + if(bitsAllocated > 8 && bitsAllocated < 16) + bitsAllocated = 16; + frameSize = XSize*YSize*SamplesPerPixel*bitsAllocated/8; + + //// Second stage: read from disk and decompress. + + if ( BitsAllocated == 12 ) // We suppose 'BitsAllocated' = 12 only exist for uncompressed files { ReadAndDecompress12BitsTo16Bits( fp); } @@ -214,39 +358,67 @@ bool PixelReadConvert::ReadAndDecompressPixelData( std::ifstream *fp ) // after the field containing the image data. In this case, these // bad data are added to the size of the image (in the PixelDataLength // variable). But RawSize is the right size of the image ! - if( PixelDataLength != RawSize) + if ( PixelDataLength != RawSize ) { - gdcmVerboseMacro( "Mismatch between PixelReadConvert and RawSize." ); - } - if( PixelDataLength > RawSize) - { - fp->read( (char*)Raw, RawSize); + gdcmWarningMacro( "Mismatch between PixelReadConvert : " + << PixelDataLength << " and RawSize : " << RawSize ); } + + //todo : is it the right patch? + char *raw = (char*)Raw; + uint32_t remainingLength; + unsigned int i; + unsigned int lengthToRead; + + if ( PixelDataLength > RawSize ) + lengthToRead = RawSize; else + lengthToRead = PixelDataLength; + + // perform a frame by frame reading + remainingLength = lengthToRead; + unsigned int nbFrames = lengthToRead / frameSize; + for (i=0;iread( (char*)Raw, PixelDataLength); + Progress = (float)(count+1)/(float)nbFrames; + fp->read( raw, frameSize); + raw += frameSize; + remainingLength -= frameSize; + count++; } + if (remainingLength !=0 ) + fp->read( raw, remainingLength); if ( fp->fail() || fp->eof()) { - gdcmVerboseMacro( "Reading of Raw pixel data failed." ); + gdcmWarningMacro( "Reading of Raw pixel data failed." ); return false; } } else if ( IsRLELossless ) { - if ( ! RLEInfo->DecompressRLEFile( fp, Raw, XSize, YSize, ZSize, BitsAllocated ) ) + if ( ! RLEInfo->DecompressRLEFile + ( fp, Raw, XSize, YSize, ZSize, TSize, BitsAllocated ) ) { - gdcmVerboseMacro( "RLE decompressor failed." ); + gdcmWarningMacro( "RLE decompressor failed." ); return false; } } + else if ( IsMPEG ) + { + //gdcmWarningMacro( "Sorry, MPEG not yet taken into account" ); + //return false; + // fp has already been seek to start of mpeg + //ReadMPEGFile(fp, (char*)Raw, PixelDataLength); + return true; + } else { // Default case concerns JPEG family if ( ! ReadAndDecompressJPEGFile( fp ) ) { - gdcmVerboseMacro( "JPEG decompressor failed." ); + gdcmWarningMacro( "JPEG decompressor ( ReadAndDecompressJPEGFile()" + << " method ) failed." ); return false; } } @@ -255,11 +427,15 @@ bool PixelReadConvert::ReadAndDecompressPixelData( std::ifstream *fp ) //// Third stage: twigle the bytes and bits. ConvertReorderEndianity(); ConvertReArrangeBits(); + ConvertFixGreyLevels(); + if (UserFunction) // user is allowed to Mirror, TopDown, Rotate,...the image + UserFunction( Raw, FileInternal); ConvertHandleColor(); return true; } +/// Deletes Pixels Area void PixelReadConvert::Squeeze() { if ( RGB ) @@ -270,13 +446,13 @@ void PixelReadConvert::Squeeze() delete [] Raw; Raw = 0; - if ( LutRGBA ) - delete [] LutRGBA; - LutRGBA = 0; + //if ( LutRGBA ) + // delete [] LutRGBA; + //LutRGBA = 0; } /** - * \brief Build the RGB image from the Raw imagage and the LUTs. + * \brief Build the RGB image from the Raw image and the LUTs. */ bool PixelReadConvert::BuildRGBImage() { @@ -298,17 +474,37 @@ bool PixelReadConvert::BuildRGBImage() // The job can't be done return false; } + + gdcmDebugMacro( "--> BuildRGBImage" ); // Build RGB Pixels AllocateRGB(); - uint8_t *localRGB = RGB; - for (size_t i = 0; i < RawSize; ++i ) + + int j; + if ( BitsAllocated <= 8 ) { - int j = Raw[i] * 4; - *localRGB++ = LutRGBA[j]; - *localRGB++ = LutRGBA[j+1]; - *localRGB++ = LutRGBA[j+2]; - } + uint8_t *localRGB = RGB; + for (size_t i = 0; i < RawSize; ++i ) + { + j = Raw[i] * 4; + *localRGB++ = LutRGBA[j]; + *localRGB++ = LutRGBA[j+1]; + *localRGB++ = LutRGBA[j+2]; + } + } + + else // deal with 16 bits pixels and 16 bits Palette color + { + uint16_t *localRGB = (uint16_t *)RGB; + for (size_t i = 0; i < RawSize/2; ++i ) + { + j = ((uint16_t *)Raw)[i] * 4; + *localRGB++ = ((uint16_t *)LutRGBA)[j]; + *localRGB++ = ((uint16_t *)LutRGBA)[j+1]; + *localRGB++ = ((uint16_t *)LutRGBA)[j+2]; + } + } + return true; } @@ -324,7 +520,8 @@ bool PixelReadConvert::BuildRGBImage() void PixelReadConvert::ReadAndDecompress12BitsTo16Bits( std::ifstream *fp ) throw ( FormatError ) { - int nbPixels = XSize * YSize; + /// \todo Fix the 3D, 4D pb + int nbPixels = XSize * YSize * TSize; uint16_t *localDecompres = (uint16_t*)Raw; for( int p = 0; p < nbPixels; p += 2 ) @@ -375,47 +572,160 @@ bool PixelReadConvert::ReadAndDecompressJPEGFile( std::ifstream *fp ) { if ( IsJPEG2000 ) { - gdcmVerboseMacro( "Sorry, JPEG2000 not yet taken into account" ); - fp->seekg( JPEGInfo->GetFirstFragment()->GetOffset(), std::ios::beg); -// if ( ! gdcm_read_JPEG2000_file( fp,Raw ) ) - return false; + // make sure this is the right JPEG compression + assert( !IsJPEGLossless || !IsJPEGLossy || !IsJPEGLS ); + // FIXME this is really ugly but it seems I have to load the complete + // jpeg2000 stream to use jasper: + // I don't think we'll ever be able to deal with multiple fragments properly + + if( ZSize == 1 ) + { + unsigned long inputlength = 0; + JPEGFragment *jpegfrag = JPEGInfo->GetFirstFragment(); + while( jpegfrag ) + { + inputlength += jpegfrag->GetLength(); + jpegfrag = JPEGInfo->GetNextFragment(); + } + gdcmAssertMacro( inputlength != 0); + uint8_t *inputdata = new uint8_t[inputlength]; + char *pinputdata = (char*)inputdata; + jpegfrag = JPEGInfo->GetFirstFragment(); + while( jpegfrag ) + { + fp->seekg( jpegfrag->GetOffset(), std::ios::beg); + fp->read(pinputdata, jpegfrag->GetLength()); + pinputdata += jpegfrag->GetLength(); + jpegfrag = JPEGInfo->GetNextFragment(); + } + // Warning the inputdata buffer is deleted in the function + if ( gdcm_read_JPEG2000_file( Raw, + (char*)inputdata, inputlength ) ) + { + return true; + } + // wow what happen, must be an error + gdcmWarningMacro( "gdcm_read_JPEG2000_file() failed "); + return false; + } + else + { + if( (unsigned int)ZSize != JPEGInfo->GetFragmentCount() ) + { + gdcmErrorMacro( "Sorry GDCM does not handle this type of fragments" ); + return false; + } + // Hopefully every dicom fragment is *exactly* the j2k stream + JPEGFragment *jpegfrag = JPEGInfo->GetFirstFragment(); + char *praw = (char*)Raw; + while( jpegfrag ) + { + unsigned long inputlength = jpegfrag->GetLength(); + char *inputdata = new char[inputlength]; + fp->seekg( jpegfrag->GetOffset(), std::ios::beg); + fp->read(inputdata, jpegfrag->GetLength()); + // Warning the inputdata buffer is deleted in the function + gdcm_read_JPEG2000_file( praw, + inputdata, inputlength) ; + praw += XSize*YSize*SamplesPerPixel*(BitsAllocated/8); + jpegfrag = JPEGInfo->GetNextFragment(); + } + return true; + } } - - if ( IsJPEGLS ) + else if ( IsJPEGLS ) { - gdcmVerboseMacro( "Sorry, JPEG-LS not yet taken into account" ); + // make sure this is the right JPEG compression + assert( !IsJPEGLossless || !IsJPEGLossy || !IsJPEG2000 ); + // WARNING : JPEG-LS is NOT the 'classical' Jpeg Lossless : + // [JPEG-LS is the basis for new lossless/near-lossless compression + // standard for continuous-tone images intended for JPEG2000. The standard + // is based on the LOCO-I algorithm (LOw COmplexity LOssless COmpression + // for Images) developed at Hewlett-Packard Laboratories] + // + // see http://datacompression.info/JPEGLS.shtml + // +#if 0 + std::cerr << "count:" << JPEGInfo->GetFragmentCount() << std::endl; + unsigned long inputlength = 0; + JPEGFragment *jpegfrag = JPEGInfo->GetFirstFragment(); + while( jpegfrag ) + { + inputlength += jpegfrag->GetLength(); + jpegfrag = JPEGInfo->GetNextFragment(); + } + gdcmAssertMacro( inputlength != 0); + uint8_t *inputdata = new uint8_t[inputlength]; + char *pinputdata = (char*)inputdata; + jpegfrag = JPEGInfo->GetFirstFragment(); + while( jpegfrag ) + { + fp->seekg( jpegfrag->GetOffset(), std::ios::beg); + fp->read(pinputdata, jpegfrag->GetLength()); + pinputdata += jpegfrag->GetLength(); + jpegfrag = JPEGInfo->GetNextFragment(); + } + + //fp->read((char*)Raw, PixelDataLength); + + std::ofstream out("/tmp/jpegls.jpg"); + out.write((char*)inputdata, inputlength); + out.close(); + delete[] inputdata; +#endif + + gdcmWarningMacro( "Sorry, JPEG-LS not yet taken into account" ); fp->seekg( JPEGInfo->GetFirstFragment()->GetOffset(), std::ios::beg); // if ( ! gdcm_read_JPEGLS_file( fp,Raw ) ) return false; } - - // else ?? - // Precompute the offset localRaw will be shifted with - int length = XSize * YSize * SamplesPerPixel; - int numberBytes = BitsAllocated / 8; - - JPEGInfo->DecompressFromFile(fp, Raw, BitsStored, numberBytes, length ); - return true; + else if( JPEGInfo ) + { + // make sure this is the right JPEG compression + assert( !IsJPEGLS || !IsJPEG2000 ); + // Precompute the offset localRaw will be shifted with + int length = XSize * YSize * ZSize * SamplesPerPixel; + int numberBytes = BitsAllocated / 8; + + // to avoid major troubles when BitsStored == 8 && BitsAllocated==16 ! + int dummy; + if (BitsStored == 8 && BitsAllocated==16) + dummy = 16; + else + dummy = BitsStored; + JPEGInfo->DecompressFromFile(fp, Raw, dummy, numberBytes, length ); + return true; + } + //else (not sure how get there...), must be one of those crazy DICOM file + return false; } /** - * \brief Build Red/Green/Blue/Alpha LUT from File - * when (0028,0004),Photometric Interpretation = [PALETTE COLOR ] - * and (0028,1101),(0028,1102),(0028,1102) - * - xxx Palette Color Lookup Table Descriptor - are found - * and (0028,1201),(0028,1202),(0028,1202) - * - xxx Palette Color Lookup Table Data - are found + * \brief Build Red/Green/Blue/Alpha LUT from File when : + * - (0028,0004) : Photometric Interpretation == [PALETTE COLOR ] + * and + * - (0028,1101),(0028,1102),(0028,1102) + * xxx Palette Color Lookup Table Descriptor are found + * and + * - (0028,1201),(0028,1202),(0028,1202) + * xxx Palette Color Lookup Table Data - are found * \warning does NOT deal with : - * 0028 1100 Gray Lookup Table Descriptor (Retired) - * 0028 1221 Segmented Red Palette Color Lookup Table Data - * 0028 1222 Segmented Green Palette Color Lookup Table Data - * 0028 1223 Segmented Blue Palette Color Lookup Table Data + * - 0028 1100 Gray Lookup Table Descriptor (Retired) + * - 0028 1221 Segmented Red Palette Color Lookup Table Data + * - 0028 1222 Segmented Green Palette Color Lookup Table Data + * - 0028 1223 Segmented Blue Palette Color Lookup Table Data * no known Dicom reader deals with them :-( * @return a RGBA Lookup Table */ void PixelReadConvert::BuildLUTRGBA() { + + // Note to code reviewers : + // The problem is *much more* complicated, since a lot of manufacturers + // Don't follow the norm : + // have a look at David Clunie's remark at the end of this .cxx file. if ( LutRGBA ) + { return; } @@ -431,6 +741,7 @@ void PixelReadConvert::BuildLUTRGBA() || LutGreenDescriptor == GDCM_UNFOUND || LutBlueDescriptor == GDCM_UNFOUND ) { + gdcmWarningMacro( "(At least) a LUT Descriptor is missing" ); return; } @@ -439,23 +750,25 @@ void PixelReadConvert::BuildLUTRGBA() int lengthR; // Red LUT length in Bytes int debR; // Subscript of the first Lut Value int nbitsR; // Lut item size (in Bits) - int nbRead = sscanf( LutRedDescriptor.c_str(), + int nbRead; // nb of items in LUT descriptor (must be = 3) + + nbRead = sscanf( LutRedDescriptor.c_str(), "%d\\%d\\%d", &lengthR, &debR, &nbitsR ); - if( nbRead != 3 ) + if ( nbRead != 3 ) { - gdcmVerboseMacro( "Wrong Red LUT descriptor" ); - } - + gdcmWarningMacro( "Wrong Red LUT descriptor" ); + } int lengthG; // Green LUT length in Bytes int debG; // Subscript of the first Lut Value int nbitsG; // Lut item size (in Bits) + nbRead = sscanf( LutGreenDescriptor.c_str(), "%d\\%d\\%d", - &lengthG, &debG, &nbitsG ); - if( nbRead != 3 ) + &lengthG, &debG, &nbitsG ); + if ( nbRead != 3 ) { - gdcmVerboseMacro( "Wrong Green LUT descriptor" ); + gdcmWarningMacro( "Wrong Green LUT descriptor" ); } int lengthB; // Blue LUT length in Bytes @@ -464,85 +777,220 @@ void PixelReadConvert::BuildLUTRGBA() nbRead = sscanf( LutRedDescriptor.c_str(), "%d\\%d\\%d", &lengthB, &debB, &nbitsB ); - if( nbRead != 3 ) + if ( nbRead != 3 ) { - gdcmVerboseMacro( "Wrong Blue LUT descriptor" ); + gdcmWarningMacro( "Wrong Blue LUT descriptor" ); } + + gdcmDebugMacro(" lengthR " << lengthR << " debR " + << debR << " nbitsR " << nbitsR); + gdcmDebugMacro(" lengthG " << lengthG << " debG " + << debG << " nbitsG " << nbitsG); + gdcmDebugMacro(" lengthB " << lengthB << " debB " + << debB << " nbitsB " << nbitsB); + + if ( !lengthR ) // if = 2^16, this shall be 0 see : CP-143 + lengthR=65536; + if ( !lengthG ) // if = 2^16, this shall be 0 + lengthG=65536; + if ( !lengthB ) // if = 2^16, this shall be 0 + lengthB=65536; //////////////////////////////////////////////////////// + if ( ( ! LutRedData ) || ( ! LutGreenData ) || ( ! LutBlueData ) ) { + gdcmWarningMacro( "(At least) a LUT is missing" ); return; } - //////////////////////////////////////////////// - // forge the 4 * 8 Bits Red/Green/Blue/Alpha LUT - LutRGBA = new uint8_t[ 1024 ]; // 256 * 4 (R, G, B, Alpha) - if ( !LutRGBA ) - return; - - memset( LutRGBA, 0, 1024 ); - - int mult; - if ( ( nbitsR == 16 ) && ( BitsAllocated == 8 ) ) - { - // when LUT item size is different than pixel size - mult = 2; // high byte must be = low byte - } - else + // ------------------------------------------------------------- + + if ( BitsAllocated <= 8 ) { - // See PS 3.3-2003 C.11.1.1.2 p 619 - mult = 1; - } + // forge the 4 * 8 Bits Red/Green/Blue/Alpha LUT + LutRGBA = new uint8_t[ 1024 ]; // 256 * 4 (R, G, B, Alpha) + if ( !LutRGBA ) + return; + LutItemNumber = 256; + LutItemSize = 8; + memset( LutRGBA, 0, 1024 ); - // if we get a black image, let's just remove the '+1' - // from 'i*mult+1' and check again - // if it works, we shall have to check the 3 Palettes - // to see which byte is ==0 (first one, or second one) - // and fix the code - // We give up the checking to avoid some (useless ?) overhead - // (optimistic asumption) - int i; - uint8_t *a = LutRGBA + 0; - for( i=0; i < lengthR; ++i ) - { - *a = LutRedData[i*mult+1]; - a += 4; - } + int mult; + if ( ( nbitsR == 16 ) && ( BitsAllocated == 8 ) ) + { + // when LUT item size is different than pixel size + mult = 2; // high byte must be = low byte + } + else + { + // See PS 3.3-2003 C.11.1.1.2 p 619 + mult = 1; + } - a = LutRGBA + 1; - for( i=0; i < lengthG; ++i) - { - *a = LutGreenData[i*mult+1]; - a += 4; - } + // if we get a black image, let's just remove the '+1' + // from 'i*mult+1' and check again + // if it works, we shall have to check the 3 Palettes + // to see which byte is ==0 (first one, or second one) + // and fix the code + // We give up the checking to avoid some (useless ?) overhead + // (optimistic asumption) + int i; + uint8_t *a; + + //take "Subscript of the first Lut Value" (debR,debG,debB) into account! + + //FIXME : +1 : to get 'low value' byte + // Trouble expected on Big Endian Processors ? + // 16 BIts Per Pixel Palette Color to be swapped? + + a = LutRGBA + 0 + debR; + for( i=0; i < lengthR; ++i ) + { + *a = LutRedData[i*mult+1]; + a += 4; + } - a = LutRGBA + 2; - for(i=0; i < lengthB; ++i) - { - *a = LutBlueData[i*mult+1]; - a += 4; - } + a = LutRGBA + 1 + debG; + for( i=0; i < lengthG; ++i) + { + *a = LutGreenData[i*mult+1]; + a += 4; + } - a = LutRGBA + 3; - for(i=0; i < 256; ++i) + a = LutRGBA + 2 + debB; + for(i=0; i < lengthB; ++i) + { + *a = LutBlueData[i*mult+1]; + a += 4; + } + + a = LutRGBA + 3 ; + for(i=0; i < 256; ++i) + { + *a = 1; // Alpha component + a += 4; + } + } + else { - *a = 1; // Alpha component - a += 4; + // Probabely the same stuff is to be done for 16 Bits Pixels + // with 65536 entries LUT ?!? + // Still looking for accurate info on the web :-( + + gdcmWarningMacro( "Sorry Palette Color Lookup Tables not yet dealt with" + << " for 16 Bits Per Pixel images" ); + + // forge the 4 * 16 Bits Red/Green/Blue/Alpha LUT + + LutRGBA = (uint8_t *)new uint16_t[ 65536*4 ]; // 2^16 * 4 (R, G, B, Alpha) + if ( !LutRGBA ) + return; + memset( LutRGBA, 0, 65536*4*2 ); // 16 bits = 2 bytes ;-) + + LutItemNumber = 65536; + LutItemSize = 16; + + int i; + uint16_t *a16; + + //take "Subscript of the first Lut Value" (debR,debG,debB) into account! + + a16 = (uint16_t*)LutRGBA + 0 + debR; + for( i=0; i < lengthR; ++i ) + { + *a16 = ((uint16_t*)LutRedData)[i]; + a16 += 4; + } + + a16 = (uint16_t*)LutRGBA + 1 + debG; + for( i=0; i < lengthG; ++i) + { + *a16 = ((uint16_t*)LutGreenData)[i]; + a16 += 4; + } + + a16 = (uint16_t*)LutRGBA + 2 + debB; + for(i=0; i < lengthB; ++i) + { + *a16 = ((uint16_t*)LutBlueData)[i]; + a16 += 4; + } + + a16 = (uint16_t*)LutRGBA + 3 ; + for(i=0; i < 65536; ++i) + { + *a16 = 1; // Alpha component + a16 += 4; + } +// Just to 'see' the LUT, at debug time +// Don't remove this commented out code. + + a16=(uint16_t*)LutRGBA; + for (int j=0;j<65536;j++) + { + std::cout << *a16 << " " << *(a16+1) << " " + << *(a16+2) << " " << *(a16+3) << std::endl; + a16+=4; + } + } } /** - * \brief Swap the bytes, according to \ref SwapCode. + * \brief Swap the bytes, according to SwapCode. */ void PixelReadConvert::ConvertSwapZone() { unsigned int i; - if( BitsAllocated == 16 ) + // If this file is 'ImplicitVR BigEndian PrivateGE Transfer Syntax', + // then the header is in little endian format and the pixel data is in + // big endian format. When reading the header, GDCM has already established + // a byte swapping code suitable for this machine to read the + // header. In TS::ImplicitVRBigEndianPrivateGE, this code will need + // to be switched in order to read the pixel data. This must be + // done REGARDLESS of the processor endianess! + // + // Example: Assume we are on a little endian machine. When + // GDCM reads the header, the header will match the machine + // endianess and the swap code will be established as a no-op. + // When GDCM reaches the pixel data, it will need to switch the + // swap code to do big endian to little endian conversion. + // + // Now, assume we are on a big endian machine. When GDCM reads the + // header, the header will be recognized as a different endianess + // than the machine endianess, and a swap code will be established + // to convert from little endian to big endian. When GDCM readers + // the pixel data, the pixel data endianess will now match the + // machine endianess. But we currently have a swap code that + // converts from little endian to big endian. In this case, we + // need to switch the swap code to a no-op. + // + // Therefore, in either case, if the file is in + // 'ImplicitVR BigEndian PrivateGE Transfer Syntax', then GDCM needs to switch + // the byte swapping code when entering the pixel data. + + int tempSwapCode = SwapCode; + if ( IsPrivateGETransferSyntax ) + { + gdcmWarningMacro(" IsPrivateGETransferSyntax found; turn the SwapCode"); + // PrivateGETransferSyntax only exists for 'true' Dicom images + // we assume there is no 'exotic' 32 bits endianess! + if (SwapCode == 1234) + { + tempSwapCode = 4321; + } + else if (SwapCode == 4321) + { + tempSwapCode = 1234; + } + } + + if ( BitsAllocated == 16 ) { uint16_t *im16 = (uint16_t*)Raw; - switch( SwapCode ) + switch( tempSwapCode ) { case 1234: break; @@ -555,16 +1003,17 @@ void PixelReadConvert::ConvertSwapZone() } break; default: - gdcmVerboseMacro("SwapCode value (16 bits) not allowed."); + gdcmWarningMacro("SwapCode value (16 bits) not allowed." + << tempSwapCode); } } - else if( BitsAllocated == 32 ) + else if ( BitsAllocated == 32 ) { uint32_t s32; uint16_t high; uint16_t low; uint32_t *im32 = (uint32_t*)Raw; - switch ( SwapCode ) + switch ( tempSwapCode ) { case 1234: break; @@ -600,7 +1049,7 @@ void PixelReadConvert::ConvertSwapZone() } break; default: - gdcmVerboseMacro("SwapCode value (32 bits) not allowed." ); + gdcmWarningMacro("SwapCode value (32 bits) not allowed." << tempSwapCode ); } } } @@ -624,7 +1073,7 @@ void PixelReadConvert::ConvertReorderEndianity() uint16_t *deb = (uint16_t *)Raw; for(int i = 0; i> ( BitsAllocated - BitsStored ); + // pmask : to mask the 'unused bits' (may contain overlays) + uint16_t pmask = 0xffff; + + // It's up to the user to decide if he wants to ignore overlays (if any), + // not to gdcm, without asking. + // default is NOT TO LOAD, in order not to confuse ITK users (and others!). + + if ( !FH->GetKeepOverlays() ) // mask spurious bits ! (overlay are NOT loaded!) + { + pmask = pmask >> ( BitsAllocated - BitsStored ); + } + // else : it's up to the user to manage the 'pixels + overlays' he just loaded! + uint16_t *deb = (uint16_t*)Raw; - for(int i = 0; i> (BitsStored - HighBitPosition - 1)) & mask; - deb++; + for(int i = 0; i> (BitsStored - HighBitPosition - 1)) & pmask; + deb++; + } + } + else // Pixels are signed + { + // Hope there is never ACR-NEMA-like overlays within signed pixels (?!?) + + // smask : to check the 'sign' when BitsStored != BitsAllocated + uint16_t smask = 0x0001; + smask = smask << ( 16 - (BitsAllocated - BitsStored + 1) ); + // nmask : to propagate sign bit on negative values + int16_t nmask = (int16_t)0x8000; + nmask = nmask >> ( BitsAllocated - BitsStored - 1 ); + + for(int i = 0; i> (BitsStored - HighBitPosition - 1); + if ( *deb & smask ) + { + *deb = *deb | nmask; + } + else + { + *deb = *deb & pmask; + } + deb++; + } } } else if ( BitsAllocated == 32 ) { - uint32_t mask = 0xffffffff; - mask = mask >> ( BitsAllocated - BitsStored ); + // pmask : to mask the 'unused bits' (may contain overlays) + uint32_t pmask = 0xffffffff; + pmask = pmask >> ( BitsAllocated - BitsStored ); + uint32_t *deb = (uint32_t*)Raw; - for(int i = 0; i> (BitsStored - HighBitPosition - 1)) & mask; - deb++; + for(int i = 0; i> (BitsStored - HighBitPosition - 1)) & pmask; + deb++; + } + } + else + { + // smask : to check the 'sign' when BitsStored != BitsAllocated + uint32_t smask = 0x00000001; + smask = smask >> ( 32 - (BitsAllocated - BitsStored +1 )); + // nmask : to propagate sign bit on negative values + int32_t nmask = 0x80000000; + nmask = nmask >> ( BitsAllocated - BitsStored -1 ); + + for(int i = 0; i> (BitsStored - HighBitPosition - 1); + if ( *deb & smask ) + *deb = *deb | nmask; + else + *deb = *deb & pmask; + deb++; + } } } else { - gdcmVerboseMacro("Weird image"); + gdcmWarningMacro("Weird image (BitsAllocated !=8, 12, 16, 32)"); throw FormatError( "Weird image !?" ); } } @@ -679,6 +1263,8 @@ bool PixelReadConvert::ConvertReArrangeBits() throw ( FormatError ) */ void PixelReadConvert::ConvertRGBPlanesToRGBPixels() { + gdcmWarningMacro("--> ConvertRGBPlanesToRGBPixels"); + uint8_t *localRaw = Raw; uint8_t *copyRaw = new uint8_t[ RawSize ]; memmove( copyRaw, localRaw, RawSize ); @@ -704,6 +1290,17 @@ void PixelReadConvert::ConvertRGBPlanesToRGBPixels() */ void PixelReadConvert::ConvertYcBcRPlanesToRGBPixels() { + // Remarks for YBR newbees : + // YBR_FULL works very much like RGB, i.e. three samples per pixel, + // just the color space is YCbCr instead of RGB. This is particularly useful + // for doppler ultrasound where most of the image is grayscale + // (i.e. only populates the Y components) and Cb and Cr are mostly zero, + // except for the few patches of color on the image. + // On such images, RLE achieves a compression ratio that is much better + // than the compression ratio on an equivalent RGB image. + + gdcmWarningMacro("--> ConvertYcBcRPlanesToRGBPixels"); + uint8_t *localRaw = Raw; uint8_t *copyRaw = new uint8_t[ RawSize ]; memmove( copyRaw, localRaw, RawSize ); @@ -713,33 +1310,40 @@ void PixelReadConvert::ConvertYcBcRPlanesToRGBPixels() // ftp://medical.nema.org/medical/dicom/final/sup61_ft.pdf // and be *very* affraid // - int l = XSize * YSize; + + /// \todo : find an example to see how 3rd dim and 4th dim work together + int l = XSize * YSize * TSize; int nbFrames = ZSize; - uint8_t *a = copyRaw; + uint8_t *a = copyRaw + 0; uint8_t *b = copyRaw + l; - uint8_t *c = copyRaw + l + l; - double R, G, B; + uint8_t *c = copyRaw + l+ l; + int32_t R, G, B; - /// \todo : Replace by the 'well known' integer computation - /// counterpart. Refer to + /// We replaced easy to understand but time consuming floating point + /// computations by the 'well known' integer computation counterpart + /// Refer to : /// http://lestourtereaux.free.fr/papers/data/yuvrgb.pdf - /// for code optimisation. + /// for code optimisation. for ( int i = 0; i < nbFrames; i++ ) { for ( int j = 0; j < l; j++ ) { - R = 1.164 *(*a-16) + 1.596 *(*c -128) + 0.5; - G = 1.164 *(*a-16) - 0.813 *(*c -128) - 0.392 *(*b -128) + 0.5; - B = 1.164 *(*a-16) + 2.017 *(*b -128) + 0.5; + R = 38142 *(*a-16) + 52298 *(*c -128); + G = 38142 *(*a-16) - 26640 *(*c -128) - 12845 *(*b -128); + B = 38142 *(*a-16) + 66093 *(*b -128); - if (R < 0.0) R = 0.0; - if (G < 0.0) G = 0.0; - if (B < 0.0) B = 0.0; - if (R > 255.0) R = 255.0; - if (G > 255.0) G = 255.0; - if (B > 255.0) B = 255.0; + R = (R+16384)>>15; + G = (G+16384)>>15; + B = (B+16384)>>15; + + if (R < 0) R = 0; + if (G < 0) G = 0; + if (B < 0) B = 0; + if (R > 255) R = 255; + if (G > 255) G = 255; + if (B > 255) B = 255; *(localRaw++) = (uint8_t)R; *(localRaw++) = (uint8_t)G; @@ -752,6 +1356,11 @@ void PixelReadConvert::ConvertYcBcRPlanesToRGBPixels() delete[] copyRaw; } +/// \brief Deals with the color decoding i.e. handle: +/// - R, G, B planes (as opposed to RGB pixels) +/// - YBR (various) encodings. +/// - LUT[s] (or "PALETTE COLOR"). + void PixelReadConvert::ConvertHandleColor() { ////////////////////////////////// @@ -780,7 +1389,7 @@ void PixelReadConvert::ConvertHandleColor() // - "Planar Configuration" = 0, // - "Photometric Interpretation" = "PALETTE COLOR". // Hence gdcm will use the folowing "heuristic" in order to be tolerant - // towards Dicom-non-conformance files: + // towards Dicom-non-conformant files: // << whatever the "Planar Configuration" value might be, a // "Photometric Interpretation" set to "PALETTE COLOR" forces // a LUT intervention >> @@ -792,9 +1401,13 @@ void PixelReadConvert::ConvertHandleColor() // - [Planar 1] AND [Photo C] handled with ConvertYcBcRPlanesToRGBPixels() // - [Planar 2] OR [Photo D] requires LUT intervention. + gdcmDebugMacro("--> ConvertHandleColor " + << "Planar Configuration " << PlanarConfiguration ); + if ( ! IsRawRGB() ) { // [Planar 2] OR [Photo D]: LUT intervention done outside + gdcmDebugMacro("--> RawRGB : LUT intervention done outside"); return; } @@ -803,49 +1416,58 @@ void PixelReadConvert::ConvertHandleColor() if ( IsYBRFull ) { // [Planar 1] AND [Photo C] (remember YBR_FULL_422 acts as RGB) + gdcmDebugMacro("--> YBRFull"); ConvertYcBcRPlanesToRGBPixels(); } else { // [Planar 1] AND [Photo C] + gdcmDebugMacro("--> YBRFull"); ConvertRGBPlanesToRGBPixels(); } return; } // When planarConf is 0, and RLELossless (forbidden by Dicom norm) - // pixels need to be RGB-fied anyway + // pixels need to be RGB-fyied anyway + if (IsRLELossless) - { + { + gdcmDebugMacro("--> RLE Lossless"); ConvertRGBPlanesToRGBPixels(); } + // In *normal *case, when planarConf is 0, pixels are already in RGB } +/// Computes the Pixels Size void PixelReadConvert::ComputeRawAndRGBSizes() { int bitsAllocated = BitsAllocated; // Number of "Bits Allocated" is fixed to 16 when it's 12, since // in this case we will expand the image to 16 bits (see - // \ref ReadAndDecompress12BitsTo16Bits() ) + // ReadAndDecompress12BitsTo16Bits() ) if ( BitsAllocated == 12 ) { bitsAllocated = 16; } - RawSize = XSize * YSize * ZSize + RawSize = XSize * YSize * ZSize * TSize * ( bitsAllocated / 8 ) * SamplesPerPixel; if ( HasLUT ) { - RGBSize = 3 * RawSize; + RGBSize = 3 * RawSize; // works for 8 and 16 bits per Pixel } else { RGBSize = RawSize; } + RawSize += RawSize%2; + RGBSize += RGBSize%2; } +/// Allocates room for RGB Pixels void PixelReadConvert::AllocateRGB() { if ( RGB ) @@ -853,6 +1475,7 @@ void PixelReadConvert::AllocateRGB() RGB = new uint8_t[RGBSize]; } +/// Allocates room for RAW Pixels void PixelReadConvert::AllocateRaw() { if ( Raw ) @@ -887,7 +1510,7 @@ void PixelReadConvert::Print( std::ostream &os, std::string const &indent ) } else { - gdcmVerboseMacro("Set as RLE file but NO RLEinfo present."); + gdcmWarningMacro("Set as RLE file but NO RLEinfo present."); } } @@ -899,21 +1522,195 @@ void PixelReadConvert::Print( std::ostream &os, std::string const &indent ) } else { - gdcmVerboseMacro("Set as JPEG file but NO JPEGinfo present."); + gdcmWarningMacro("Set as JPEG file but NO JPEGinfo present."); } } } +/** + * \brief CallStartMethod + */ +void PixelReadConvert::CallStartMethod() +{ + Progress = 0.0f; + Abort = false; + CommandManager::ExecuteCommand(FH,CMD_STARTPROGRESS); +} + +/** + * \brief CallProgressMethod + */ +void PixelReadConvert::CallProgressMethod() +{ + CommandManager::ExecuteCommand(FH,CMD_PROGRESS); +} + +/** + * \brief CallEndMethod + */ +void PixelReadConvert::CallEndMethod() +{ + Progress = 1.0f; + CommandManager::ExecuteCommand(FH,CMD_ENDPROGRESS); +} + //----------------------------------------------------------------------------- } // end namespace gdcm -// NOTES on File internal calls -// User -// ---> GetImageData -// ---> GetImageDataIntoVector -// |---> GetImageDataIntoVectorRaw -// | lut intervention -// User -// ---> GetImageDataRaw -// ---> GetImageDataIntoVectorRaw +// Note to developpers : +// Here is a very detailled post from David Clunie, on the troubles caused +// 'non standard' LUT and LUT description +// We shall have to take it into accound in our code. +// Some day ... + + +/* +Subject: Problem with VOI LUTs in Agfa and Fuji CR and GE DX images, was Re: VOI LUT issues +Date: Sun, 06 Feb 2005 17:13:40 GMT +From: David Clunie +Reply-To: dclunie@dclunie.com +Newsgroups: comp.protocols.dicom +References: <1107553502.040221.189550@o13g2000cwo.googlegroups.com> + +> THE LUT that comes with [my] image claims to be 16-bit, but none of the +> values goes higher than 4095. That being said, though, none of my +> original pixel values goes higher than that, either. I have read +> elsewhere on this group that when that happens you are supposed to +> adjust the LUT. Can someone be more specific? There was a thread from +> 2002 where Marco and David were mentioning doing precisely that. +> +> Thanks +> +> -carlos rodriguez + + +You have encountered the well known "we know what the standard says but +we are going to ignore it and do what we have been doing for almost +a decade regardless" CR vendor bug. Agfa started this, but they are not +the only vendor doing this now; GE and Fuji may have joined the club. + +Sadly, one needs to look at the LUT Data, figure out what the maximum +value actually encoded is, and find the next highest power of 2 (e.g. +212 in this case), to figure out what the range of the data is +supposed to be. I have assumed that if the maximum value in the LUT +data is less than a power of 2 minus 1 (e.g. 0xebc) then the intent +of the vendor was not to use the maximum available grayscale range +of the display (e.g. the maximum is 0xfff in this case). An alternative +would be to scale to the actual maximum rather than a power of two. + +Very irritating, and in theory not totally reliable if one really +intended the full 16 bits and only used, say 15, but that is extremely +unlikely since everything would be too dark, and this heuristic +seems to work OK. + +There has never been anything in the standard that describes having +to go through these convolutions. Since the only value in the +standard that describes the bit depth of the LUT values is LUT +Descriptor value 3 and that is (usually) always required to be +either 8 or 16, it mystifies me how the creators' of these images +imagine that the receiver is going to divine the range that is intended. Further, the standard is quite explicit that this 3rd +value defines the range of LUT values, but as far as I am aware, all +the vendors are ignoring the standard and indeed sending a third value +of 16 in all cases. + +This problem is not confined to CR, and is also seen with DX products. + +Typically I have seen: + +- Agfa CR, which usually (always ?) sends LUTs, values up to 0x0fff +- Fuji CR, which occasionally send LUTs, values up to 0x03ff +- GE DX, for presentation, which always have LUTs, up to 0x3fff + +Swissray, Siemens, Philips, Canon and Kodak never seem to send VOI LUTs +at this point (which is a whole other problem). Note that the presence +or absence of a VOI LUT as opposed to window values may be configurable +on the modality in some cases, and I have just looked at what I happen +to have received from a myriad of sites over whose configuration I have +no control. This may be why the majority of Fuji images have no VOI LUTs, +but a few do (or it may be the Siemens system that these Fuji images went +through that perhaps added it). I do have some test Hologic DX images that +are not from a clinical site that do actually get this right (a value +of 12 for the third value and a max of 0xfff). + +Since almost every vendor that I have encountered that encodes LUTs +makes this mistake, perhaps it is time to amend the standard to warn +implementor's of receivers and/or sanction this bad behavior. We have +talked about this in the past in WG 6 but so far everyone has been +reluctant to write into the standard such a comment. Maybe it is time +to try again, since if one is not aware of this problem, one cannot +effectively implement display using VOI LUTs, and there is a vast +installed base to contend with. + +I did not check presentation states, in which VOI LUTs could also be +encountered, for the prevalence of this mistake, nor did I look at the +encoding of Modality LUT's, which are unusual. Nor did I check digital +mammography images. I would be interested to hear from anyone who has. + +David + +PS. The following older thread in this newsgroup discusses this: + +"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" + +PPS. From a historical perspective, the following may be of interest. + +In the original standard in 1993, all that was said about this was a +reference to the corresponding such where Modality LUTs are described +that said: + +"The third value specifies the number of bits for each entry in the +LUT Data. It shall take the value 8 or 16. The LUT Data shall be stored +in a format equivalent to 8 or 16 bits allocated and high bit equal +1-bits allocated." + +Since the high bit hint was not apparently explicit enough, a very +early CP, CP 15 (submitted by Agfa as it happens), replaced this with: + +"The third value conveys the range of LUT entry values. It shall take +the value 8 or 16, corresponding with the LUT entry value range of +256 or 65536. + +Note: The third value is not required for describing the + LUT data and is only included for informational usage + and for maintaining compatibility with ACRNEMA 2.0. + +The LUT Data contains the LUT entry values." + +That is how it read in the 1996, 1998 and 1999 editions. + +By the 2000 edition, Supplement 33 that introduced presentation states +extensively reworked this entire section and tried to explain this in +different words: + +"The output range is from 0 to 2^n-1 where n is the third value of LUT +Descriptor. This range is always unsigned." + +and also added a note to spell out what the output range meant in the +VOI LUT section: + +"9. The output of the Window Center/Width or VOI LUT transformation +is either implicitly scaled to the full range of the display device +if there is no succeeding transformation defined, or implicitly scaled +to the full input range of the succeeding transformation step (such as +the Presentation LUT), if present. See C.11.6.1." + +It still reads this way in the 2004 edition. + +Note that LUTs in other applications than the general VOI LUT allow for +values other than 8 or 16 in the third value of LUT descriptor to permit +ranges other than 0 to 255 or 65535. + +In addition, the DX Image Module specializes the VOI LUT +attributes as follows, in PS 3.3 section C.8.11.3.1.5 (added in Sup 32): + +"The third value specifies the number of bits for each entry in the LUT +Data (analogous to ìbits storedî). It shall be between 10-16. The LUT +Data shall be stored in a format equivalent to 16 ìbits allocatedî and +ìhigh bitî equal to ìbits storedî - 1. The third value conveys the range +of LUT entry values. These unsigned LUT entry values shall range between +0 and 2^n-1, where n is the third value of the LUT Descriptor." + +So in the case of the GE DX for presentation images, the third value of +LUT descriptor is allowed to be and probably should be 14 rather than 16. +*/