1 /*=========================================================================
4 Module: $RCSfile: gdcmPixelReadConvert.cxx,v $
6 Date: $Date: 2005/11/09 11:22:51 $
7 Version: $Revision: 1.100 $
9 Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
10 l'Image). All rights reserved. See Doc/License.txt or
11 http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
13 This software is distributed WITHOUT ANY WARRANTY; without even
14 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
15 PURPOSE. See the above copyright notices for more information.
17 =========================================================================*/
19 #include "gdcmPixelReadConvert.h"
20 #include "gdcmDebug.h"
22 #include "gdcmGlobal.h"
24 #include "gdcmDocEntry.h"
25 #include "gdcmRLEFramesInfo.h"
26 #include "gdcmJPEGFragmentsInfo.h"
29 #include <stdio.h> //for sscanf
34 bool ReadMPEGFile (std::ifstream *fp, char *inputdata, size_t lenght);
35 bool gdcm_read_JPEG2000_file (void* raw,
36 char *inputdata, size_t inputlength);
37 //-----------------------------------------------------------------------------
38 #define str2num(str, typeNum) *((typeNum *)(str))
40 //-----------------------------------------------------------------------------
41 // Constructor / Destructor
43 PixelReadConvert::PixelReadConvert()
59 /// Canonical Destructor
60 PixelReadConvert::~PixelReadConvert()
65 //-----------------------------------------------------------------------------
68 * \brief Predicate to know whether the image[s] (once Raw) is RGB.
69 * \note See comments of \ref ConvertHandleColor
71 bool PixelReadConvert::IsRawRGB()
74 || PlanarConfiguration == 2
82 * \brief Gets various usefull informations from the file header
83 * @param file gdcm::File pointer
85 void PixelReadConvert::GrabInformationsFromFile( File *file )
87 // Number of Bits Allocated for storing a Pixel is defaulted to 16
88 // when absent from the file.
89 BitsAllocated = file->GetBitsAllocated();
90 if ( BitsAllocated == 0 )
95 // Number of "Bits Stored", defaulted to number of "Bits Allocated"
96 // when absent from the file.
97 BitsStored = file->GetBitsStored();
98 if ( BitsStored == 0 )
100 BitsStored = BitsAllocated;
103 // High Bit Position, defaulted to "Bits Allocated" - 1
104 HighBitPosition = file->GetHighBitPosition();
105 if ( HighBitPosition == 0 )
107 HighBitPosition = BitsAllocated - 1;
110 XSize = file->GetXSize();
111 YSize = file->GetYSize();
112 ZSize = file->GetZSize();
113 SamplesPerPixel = file->GetSamplesPerPixel();
114 //PixelSize = file->GetPixelSize(); Useless
115 PixelSign = file->IsSignedPixelData();
116 SwapCode = file->GetSwapCode();
118 if (! file->IsDicomV3() ) // Should be ACR-NEMA file
121 // Don't loose time checking unexistant Transfer Syntax !
122 IsPrivateGETransferSyntax = IsMPEG
123 = IsJPEG2000 = IsJPEGLS = IsJPEGLossy
124 = IsJPEGLossless = IsRLELossless
129 std::string ts = file->GetTransferSyntax();
134 // mind the order : check the most usual first.
135 if( IsRaw = Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ExplicitVRLittleEndian) break;
136 if( IsRaw = Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ImplicitVRLittleEndian ) break;
137 if( IsRaw = Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ExplicitVRBigEndian) break;
138 if( IsRaw = Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ImplicitVRBigEndianPrivateGE) break;
139 if( IsRaw = Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::DeflatedExplicitVRLittleEndian) break;
143 IsPrivateGETransferSyntax =
144 ( Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ImplicitVRBigEndianPrivateGE );
146 IsMPEG = IsJPEG2000 = IsJPEGLS = IsJPEGLossy = IsJPEGLossless = IsRLELossless = false;
152 // mind the order : check the most usual first.
153 if( IsJPEGLossy = Global::GetTS()->IsJPEGLossy(ts) ) break;
154 if( IsJPEGLossless = Global::GetTS()->IsJPEGLossless(ts) ) break;
155 if( IsRLELossless = Global::GetTS()->IsRLELossless(ts) ) break;
156 if( IsJPEG2000 = Global::GetTS()->IsJPEG2000(ts) ) break;
157 if( IsMPEG = Global::GetTS()->IsMPEG(ts) ) break;
158 if( IsJPEGLS = Global::GetTS()->IsJPEGLS(ts) ) break;
159 gdcmWarningMacro("Unexepected Transfer Syntax :[" << ts << "]");
165 PixelOffset = file->GetPixelOffset();
166 PixelDataLength = file->GetPixelAreaLength();
167 RLEInfo = file->GetRLEInfo();
168 JPEGInfo = file->GetJPEGInfo();
170 IsMonochrome = file->IsMonochrome();
171 IsMonochrome1 = file->IsMonochrome1();
172 IsPaletteColor = file->IsPaletteColor();
173 IsYBRFull = file->IsYBRFull();
175 PlanarConfiguration = file->GetPlanarConfiguration();
177 /////////////////////////////////////////////////////////////////
179 HasLUT = file->HasLUT();
182 // Just in case some access to a File element requires disk access.
183 LutRedDescriptor = file->GetEntryString( 0x0028, 0x1101 );
184 LutGreenDescriptor = file->GetEntryString( 0x0028, 0x1102 );
185 LutBlueDescriptor = file->GetEntryString( 0x0028, 0x1103 );
187 // FIXME : The following comment is probabely meaningless, since LUT are *always*
188 // loaded at parsing time, whatever their length is.
190 // Depending on the value of Document::MAX_SIZE_LOAD_ELEMENT_VALUE
191 // [ refer to invocation of Document::SetMaxSizeLoadEntry() in
192 // Document::Document() ], the loading of the value (content) of a
193 // [Bin|Val]Entry occurence migth have been hindered (read simply NOT
194 // loaded). Hence, we first try to obtain the LUTs data from the file
195 // and when this fails we read the LUTs data directly from disk.
196 // \TODO Reading a [Bin|Val]Entry directly from disk is a kludge.
197 // We should NOT bypass the [Bin|Val]Entry class. Instead
198 // an access to an UNLOADED content of a [Bin|Val]Entry occurence
199 // (e.g. DataEntry::GetBinArea()) should force disk access from
200 // within the [Bin|Val]Entry class itself. The only problem
201 // is that the [Bin|Val]Entry is unaware of the FILE* is was
202 // parsed from. Fix that. FIXME.
205 file->LoadEntryBinArea(0x0028, 0x1201);
206 LutRedData = (uint8_t*)file->GetEntryBinArea( 0x0028, 0x1201 );
209 gdcmWarningMacro("Unable to read Red Palette Color Lookup Table data");
213 file->LoadEntryBinArea(0x0028, 0x1202);
214 LutGreenData = (uint8_t*)file->GetEntryBinArea(0x0028, 0x1202 );
217 gdcmWarningMacro("Unable to read Green Palette Color Lookup Table data");
221 file->LoadEntryBinArea(0x0028, 0x1203);
222 LutBlueData = (uint8_t*)file->GetEntryBinArea( 0x0028, 0x1203 );
225 gdcmWarningMacro("Unable to read Blue Palette Color Lookup Table data");
230 ComputeRawAndRGBSizes();
233 /// \brief Reads from disk and decompresses Pixels
234 bool PixelReadConvert::ReadAndDecompressPixelData( std::ifstream *fp )
236 // ComputeRawAndRGBSizes is already made by
237 // ::GrabInformationsFromfile. So, the structure sizes are
241 //////////////////////////////////////////////////
242 //// First stage: get our hands on the Pixel Data.
245 gdcmWarningMacro( "Unavailable file pointer." );
249 fp->seekg( PixelOffset, std::ios::beg );
250 if ( fp->fail() || fp->eof() )
252 gdcmWarningMacro( "Unable to find PixelOffset in file." );
258 //////////////////////////////////////////////////
259 //// Second stage: read from disk and decompress.
260 if ( BitsAllocated == 12 )
262 ReadAndDecompress12BitsTo16Bits( fp);
266 // This problem can be found when some obvious informations are found
267 // after the field containing the image data. In this case, these
268 // bad data are added to the size of the image (in the PixelDataLength
269 // variable). But RawSize is the right size of the image !
270 if ( PixelDataLength != RawSize )
272 gdcmWarningMacro( "Mismatch between PixelReadConvert : "
273 << PixelDataLength << " and RawSize : " << RawSize );
275 if ( PixelDataLength > RawSize )
277 fp->read( (char*)Raw, RawSize);
281 fp->read( (char*)Raw, PixelDataLength);
284 if ( fp->fail() || fp->eof())
286 gdcmWarningMacro( "Reading of Raw pixel data failed." );
290 else if ( IsRLELossless )
292 if ( ! RLEInfo->DecompressRLEFile
293 ( fp, Raw, XSize, YSize, ZSize, BitsAllocated ) )
295 gdcmWarningMacro( "RLE decompressor failed." );
301 //gdcmWarningMacro( "Sorry, MPEG not yet taken into account" );
303 // fp has already been seek to start of mpeg
304 ReadMPEGFile(fp, (char*)Raw, PixelDataLength);
309 // Default case concerns JPEG family
310 if ( ! ReadAndDecompressJPEGFile( fp ) )
312 gdcmWarningMacro( "JPEG decompressor ( ReadAndDecompressJPEGFile()"
313 << " method ) failed." );
318 ////////////////////////////////////////////
319 //// Third stage: twigle the bytes and bits.
320 ConvertReorderEndianity();
321 ConvertReArrangeBits();
322 ConvertFixGreyLevels();
323 if (UserFunction) // user is allowed to Mirror, TopDown, Rotate,...the image
324 UserFunction( Raw, FileInternal);
325 ConvertHandleColor();
330 /// Deletes Pixels Area
331 void PixelReadConvert::Squeeze()
347 * \brief Build the RGB image from the Raw image and the LUTs.
349 bool PixelReadConvert::BuildRGBImage()
353 // The job is already done.
359 // The job can't be done
366 // The job can't be done
370 gdcmDebugMacro( "--> BuildRGBImage" );
376 if ( BitsAllocated <= 8 )
378 uint8_t *localRGB = RGB;
379 for (size_t i = 0; i < RawSize; ++i )
382 *localRGB++ = LutRGBA[j];
383 *localRGB++ = LutRGBA[j+1];
384 *localRGB++ = LutRGBA[j+2];
388 else // deal with 16 bits pixels and 16 bits Palette color
390 uint16_t *localRGB = (uint16_t *)RGB;
391 for (size_t i = 0; i < RawSize/2; ++i )
393 j = ((uint16_t *)Raw)[i] * 4;
394 *localRGB++ = ((uint16_t *)LutRGBA)[j];
395 *localRGB++ = ((uint16_t *)LutRGBA)[j+1];
396 *localRGB++ = ((uint16_t *)LutRGBA)[j+2];
403 //-----------------------------------------------------------------------------
406 //-----------------------------------------------------------------------------
409 * \brief Read from file a 12 bits per pixel image and decompress it
410 * into a 16 bits per pixel image.
412 void PixelReadConvert::ReadAndDecompress12BitsTo16Bits( std::ifstream *fp )
413 throw ( FormatError )
415 int nbPixels = XSize * YSize;
416 uint16_t *localDecompres = (uint16_t*)Raw;
418 for( int p = 0; p < nbPixels; p += 2 )
422 fp->read( (char*)&b0, 1);
423 if ( fp->fail() || fp->eof() )
425 throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
426 "Unfound first block" );
429 fp->read( (char*)&b1, 1 );
430 if ( fp->fail() || fp->eof())
432 throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
433 "Unfound second block" );
436 fp->read( (char*)&b2, 1 );
437 if ( fp->fail() || fp->eof())
439 throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
440 "Unfound second block" );
443 // Two steps are necessary to please VC++
445 // 2 pixels 12bit = [0xABCDEF]
446 // 2 pixels 16bit = [0x0ABD] + [0x0FCE]
448 *localDecompres++ = ((b0 >> 4) << 8) + ((b0 & 0x0f) << 4) + (b1 & 0x0f);
450 *localDecompres++ = ((b2 & 0x0f) << 8) + ((b1 >> 4) << 4) + (b2 >> 4);
452 /// \todo JPR Troubles expected on Big-Endian processors ?
457 * \brief Reads from disk the Pixel Data of JPEG Dicom encapsulated
458 * file and decompress it.
459 * @param fp File Pointer
462 bool PixelReadConvert::ReadAndDecompressJPEGFile( std::ifstream *fp )
466 // make sure this is the right JPEG compression
467 assert( !IsJPEGLossless || !IsJPEGLossy || !IsJPEGLS );
468 // FIXME this is really ugly but it seems I have to load the complete
469 // jpeg2000 stream to use jasper:
470 // I don't think we'll ever be able to deal with multiple fragments properly
472 unsigned long inputlength = 0;
473 JPEGFragment *jpegfrag = JPEGInfo->GetFirstFragment();
476 inputlength += jpegfrag->GetLength();
477 jpegfrag = JPEGInfo->GetNextFragment();
479 gdcmAssertMacro( inputlength != 0);
480 uint8_t *inputdata = new uint8_t[inputlength];
481 char *pinputdata = (char*)inputdata;
482 jpegfrag = JPEGInfo->GetFirstFragment();
485 fp->seekg( jpegfrag->GetOffset(), std::ios::beg);
486 fp->read(pinputdata, jpegfrag->GetLength());
487 pinputdata += jpegfrag->GetLength();
488 jpegfrag = JPEGInfo->GetNextFragment();
490 // Warning the inputdata buffer is delete in the function
491 if ( ! gdcm_read_JPEG2000_file( Raw,
492 (char*)inputdata, inputlength ) )
496 // wow what happen, must be an error
497 gdcmWarningMacro( "gdcm_read_JPEG2000_file() failed ");
502 // make sure this is the right JPEG compression
503 assert( !IsJPEGLossless || !IsJPEGLossy || !IsJPEG2000 );
504 // WARNING : JPEG-LS is NOT the 'classical' Jpeg Lossless :
505 // [JPEG-LS is the basis for new lossless/near-lossless compression
506 // standard for continuous-tone images intended for JPEG2000. The standard
507 // is based on the LOCO-I algorithm (LOw COmplexity LOssless COmpression
508 // for Images) developed at Hewlett-Packard Laboratories]
510 // see http://datacompression.info/JPEGLS.shtml
513 std::cerr << "count:" << JPEGInfo->GetFragmentCount() << std::endl;
514 unsigned long inputlength = 0;
515 JPEGFragment *jpegfrag = JPEGInfo->GetFirstFragment();
518 inputlength += jpegfrag->GetLength();
519 jpegfrag = JPEGInfo->GetNextFragment();
521 gdcmAssertMacro( inputlength != 0);
522 uint8_t *inputdata = new uint8_t[inputlength];
523 char *pinputdata = (char*)inputdata;
524 jpegfrag = JPEGInfo->GetFirstFragment();
527 fp->seekg( jpegfrag->GetOffset(), std::ios::beg);
528 fp->read(pinputdata, jpegfrag->GetLength());
529 pinputdata += jpegfrag->GetLength();
530 jpegfrag = JPEGInfo->GetNextFragment();
533 //fp->read((char*)Raw, PixelDataLength);
535 std::ofstream out("/tmp/jpegls.jpg");
536 out.write((char*)inputdata, inputlength);
541 gdcmWarningMacro( "Sorry, JPEG-LS not yet taken into account" );
542 fp->seekg( JPEGInfo->GetFirstFragment()->GetOffset(), std::ios::beg);
543 // if ( ! gdcm_read_JPEGLS_file( fp,Raw ) )
548 // make sure this is the right JPEG compression
549 assert( !IsJPEGLS || !IsJPEG2000 );
550 // Precompute the offset localRaw will be shifted with
551 int length = XSize * YSize * SamplesPerPixel;
552 int numberBytes = BitsAllocated / 8;
554 JPEGInfo->DecompressFromFile(fp, Raw, BitsStored, numberBytes, length );
560 * \brief Build Red/Green/Blue/Alpha LUT from File when :
561 * - (0028,0004) : Photometric Interpretation == [PALETTE COLOR ]
563 * - (0028,1101),(0028,1102),(0028,1102)
564 * xxx Palette Color Lookup Table Descriptor are found
566 * - (0028,1201),(0028,1202),(0028,1202)
567 * xxx Palette Color Lookup Table Data - are found
568 * \warning does NOT deal with :
569 * - 0028 1100 Gray Lookup Table Descriptor (Retired)
570 * - 0028 1221 Segmented Red Palette Color Lookup Table Data
571 * - 0028 1222 Segmented Green Palette Color Lookup Table Data
572 * - 0028 1223 Segmented Blue Palette Color Lookup Table Data
573 * no known Dicom reader deals with them :-(
574 * @return a RGBA Lookup Table
576 void PixelReadConvert::BuildLUTRGBA()
579 // Note to code reviewers :
580 // The problem is *much more* complicated, since a lot of manufacturers
581 // Don't follow the norm :
582 // have a look at David Clunie's remark at the end of this .cxx file.
589 // http://www.barre.nom.fr/medical/dicom2/limitations.html#Color%20Lookup%20Tables
591 if ( ! IsPaletteColor )
596 if ( LutRedDescriptor == GDCM_UNFOUND
597 || LutGreenDescriptor == GDCM_UNFOUND
598 || LutBlueDescriptor == GDCM_UNFOUND )
600 gdcmWarningMacro( "(At least) a LUT Descriptor is missing" );
604 ////////////////////////////////////////////
605 // Extract the info from the LUT descriptors
606 int lengthR; // Red LUT length in Bytes
607 int debR; // Subscript of the first Lut Value
608 int nbitsR; // Lut item size (in Bits)
609 int nbRead; // nb of items in LUT descriptor (must be = 3)
611 nbRead = sscanf( LutRedDescriptor.c_str(),
613 &lengthR, &debR, &nbitsR );
616 gdcmWarningMacro( "Wrong Red LUT descriptor" );
618 int lengthG; // Green LUT length in Bytes
619 int debG; // Subscript of the first Lut Value
620 int nbitsG; // Lut item size (in Bits)
622 nbRead = sscanf( LutGreenDescriptor.c_str(),
624 &lengthG, &debG, &nbitsG );
627 gdcmWarningMacro( "Wrong Green LUT descriptor" );
630 int lengthB; // Blue LUT length in Bytes
631 int debB; // Subscript of the first Lut Value
632 int nbitsB; // Lut item size (in Bits)
633 nbRead = sscanf( LutRedDescriptor.c_str(),
635 &lengthB, &debB, &nbitsB );
638 gdcmWarningMacro( "Wrong Blue LUT descriptor" );
641 gdcmDebugMacro(" lengthR " << lengthR << " debR "
642 << debR << " nbitsR " << nbitsR);
643 gdcmDebugMacro(" lengthG " << lengthG << " debG "
644 << debG << " nbitsG " << nbitsG);
645 gdcmDebugMacro(" lengthB " << lengthB << " debB "
646 << debB << " nbitsB " << nbitsB);
648 if ( !lengthR ) // if = 2^16, this shall be 0 see : CP-143
650 if ( !lengthG ) // if = 2^16, this shall be 0
652 if ( !lengthB ) // if = 2^16, this shall be 0
655 ////////////////////////////////////////////////////////
657 if ( ( ! LutRedData ) || ( ! LutGreenData ) || ( ! LutBlueData ) )
659 gdcmWarningMacro( "(At least) a LUT is missing" );
663 // -------------------------------------------------------------
665 if ( BitsAllocated <= 8 )
667 // forge the 4 * 8 Bits Red/Green/Blue/Alpha LUT
668 LutRGBA = new uint8_t[ 1024 ]; // 256 * 4 (R, G, B, Alpha)
673 memset( LutRGBA, 0, 1024 );
676 if ( ( nbitsR == 16 ) && ( BitsAllocated == 8 ) )
678 // when LUT item size is different than pixel size
679 mult = 2; // high byte must be = low byte
683 // See PS 3.3-2003 C.11.1.1.2 p 619
687 // if we get a black image, let's just remove the '+1'
688 // from 'i*mult+1' and check again
689 // if it works, we shall have to check the 3 Palettes
690 // to see which byte is ==0 (first one, or second one)
692 // We give up the checking to avoid some (useless ?) overhead
693 // (optimistic asumption)
697 //take "Subscript of the first Lut Value" (debR,debG,debB) into account!
699 //FIXME : +1 : to get 'low value' byte
700 // Trouble expected on Big Endian Processors ?
701 // 16 BIts Per Pixel Palette Color to be swapped?
703 a = LutRGBA + 0 + debR;
704 for( i=0; i < lengthR; ++i )
706 *a = LutRedData[i*mult+1];
710 a = LutRGBA + 1 + debG;
711 for( i=0; i < lengthG; ++i)
713 *a = LutGreenData[i*mult+1];
717 a = LutRGBA + 2 + debB;
718 for(i=0; i < lengthB; ++i)
720 *a = LutBlueData[i*mult+1];
725 for(i=0; i < 256; ++i)
727 *a = 1; // Alpha component
733 // Probabely the same stuff is to be done for 16 Bits Pixels
734 // with 65536 entries LUT ?!?
735 // Still looking for accurate info on the web :-(
737 gdcmWarningMacro( "Sorry Palette Color Lookup Tables not yet dealt with"
738 << " for 16 Bits Per Pixel images" );
740 // forge the 4 * 16 Bits Red/Green/Blue/Alpha LUT
742 LutRGBA = (uint8_t *)new uint16_t[ 65536*4 ]; // 2^16 * 4 (R, G, B, Alpha)
745 memset( LutRGBA, 0, 65536*4*2 ); // 16 bits = 2 bytes ;-)
747 LutItemNumber = 65536;
753 //take "Subscript of the first Lut Value" (debR,debG,debB) into account!
755 a16 = (uint16_t*)LutRGBA + 0 + debR;
756 for( i=0; i < lengthR; ++i )
758 *a16 = ((uint16_t*)LutRedData)[i];
762 a16 = (uint16_t*)LutRGBA + 1 + debG;
763 for( i=0; i < lengthG; ++i)
765 *a16 = ((uint16_t*)LutGreenData)[i];
769 a16 = (uint16_t*)LutRGBA + 2 + debB;
770 for(i=0; i < lengthB; ++i)
772 *a16 = ((uint16_t*)LutBlueData)[i];
776 a16 = (uint16_t*)LutRGBA + 3 ;
777 for(i=0; i < 65536; ++i)
779 *a16 = 1; // Alpha component
782 /* Just to 'see' the LUT, at debug time
783 // Don't remove this commented out code.
785 a16=(uint16_t*)LutRGBA;
786 for (int j=0;j<65536;j++)
788 std::cout << *a16 << " " << *(a16+1) << " "
789 << *(a16+2) << " " << *(a16+3) << std::endl;
797 * \brief Swap the bytes, according to \ref SwapCode.
799 void PixelReadConvert::ConvertSwapZone()
802 uint16_t localSwapCode = SwapCode;
804 // If this file is 'ImplicitVR BigEndian PrivateGE Transfer Syntax',
805 // then the header is in little endian format and the pixel data is in
806 // big endian format. When reading the header, GDCM has already established
807 // a byte swapping code suitable for this machine to read the
808 // header. In TS::ImplicitVRBigEndianPrivateGE, this code will need
809 // to be switched in order to read the pixel data. This must be
810 // done REGARDLESS of the processor endianess!
812 // Example: Assume we are on a little endian machine. When
813 // GDCM reads the header, the header will match the machine
814 // endianess and the swap code will be established as a no-op.
815 // When GDCM reaches the pixel data, it will need to switch the
816 // swap code to do big endian to little endian conversion.
818 // Now, assume we are on a big endian machine. When GDCM reads the
819 // header, the header will be recognized as a different endianess
820 // than the machine endianess, and a swap code will be established
821 // to convert from little endian to big endian. When GDCM readers
822 // the pixel data, the pixel data endianess will now match the
823 // machine endianess. But we currently have a swap code that
824 // converts from little endian to big endian. In this case, we
825 // need to switch the swap code to a no-op.
827 // Therefore, in either case, if the file is in
828 // 'ImplicitVR BigEndian PrivateGE Transfer Syntax', then GDCM needs to switch
829 // the byte swapping code when entering the pixel data.
831 if ( IsPrivateGETransferSyntax )
833 // PrivateGETransferSyntax only exists for 'true' Dicom images
834 // we assume there is no 'exotic' 32 bits endianess!
835 switch (localSwapCode)
838 localSwapCode = 4321;
841 localSwapCode = 1234;
846 if ( BitsAllocated == 16 )
848 uint16_t *im16 = (uint16_t*)Raw;
849 switch( localSwapCode )
856 for( i = 0; i < RawSize / 2; i++ )
858 im16[i]= (im16[i] >> 8) | (im16[i] << 8 );
862 gdcmWarningMacro("SwapCode value (16 bits) not allowed.");
865 else if ( BitsAllocated == 32 )
870 uint32_t *im32 = (uint32_t*)Raw;
871 switch ( localSwapCode )
876 for( i = 0; i < RawSize / 4; i++ )
878 low = im32[i] & 0x0000ffff; // 4321
879 high = im32[i] >> 16;
880 high = ( high >> 8 ) | ( high << 8 );
881 low = ( low >> 8 ) | ( low << 8 );
883 im32[i] = ( s32 << 16 ) | high;
887 for( i = 0; i < RawSize / 4; i++ )
889 low = im32[i] & 0x0000ffff; // 2143
890 high = im32[i] >> 16;
891 high = ( high >> 8 ) | ( high << 8 );
892 low = ( low >> 8 ) | ( low << 8 );
894 im32[i] = ( s32 << 16 ) | low;
898 for( i = 0; i < RawSize / 4; i++ )
900 low = im32[i] & 0x0000ffff; // 3412
901 high = im32[i] >> 16;
903 im32[i] = ( s32 << 16 ) | high;
907 gdcmWarningMacro("SwapCode value (32 bits) not allowed." );
913 * \brief Deal with endianness i.e. re-arange bytes inside the integer
915 void PixelReadConvert::ConvertReorderEndianity()
917 if ( BitsAllocated != 8 )
922 // Special kludge in order to deal with xmedcon broken images:
923 if ( BitsAllocated == 16
924 && BitsStored < BitsAllocated
927 int l = (int)( RawSize / ( BitsAllocated / 8 ) );
928 uint16_t *deb = (uint16_t *)Raw;
929 for(int i = 0; i<l; i++)
931 if ( *deb == 0xffff )
941 * \brief Deal with Grey levels i.e. re-arange them
942 * to have low values = dark, high values = bright
944 void PixelReadConvert::ConvertFixGreyLevels()
949 uint32_t i; // to please M$VC6
954 if ( BitsAllocated == 8 )
956 uint8_t *deb = (uint8_t *)Raw;
957 for (i=0; i<RawSize; i++)
965 if ( BitsAllocated == 16 )
968 for (j=0; j<BitsStored-1; j++)
970 mask = (mask << 1) +1; // will be fff when BitsStored=12
973 uint16_t *deb = (uint16_t *)Raw;
974 for (i=0; i<RawSize/2; i++)
984 if ( BitsAllocated == 8 )
986 uint8_t smask8 = 255;
987 uint8_t *deb = (uint8_t *)Raw;
988 for (i=0; i<RawSize; i++)
990 *deb = smask8 - *deb;
995 if ( BitsAllocated == 16 )
997 uint16_t smask16 = 65535;
998 uint16_t *deb = (uint16_t *)Raw;
999 for (i=0; i<RawSize/2; i++)
1001 *deb = smask16 - *deb;
1010 * \brief Re-arrange the bits within the bytes.
1011 * @return Boolean always true
1013 bool PixelReadConvert::ConvertReArrangeBits() throw ( FormatError )
1016 if ( BitsStored != BitsAllocated )
1018 int l = (int)( RawSize / ( BitsAllocated / 8 ) );
1019 if ( BitsAllocated == 16 )
1021 // pmask : to mask the 'unused bits' (may contain overlays)
1022 uint16_t pmask = 0xffff;
1023 pmask = pmask >> ( BitsAllocated - BitsStored );
1025 uint16_t *deb = (uint16_t*)Raw;
1027 if ( !PixelSign ) // Pixels are unsigned
1029 for(int i = 0; i<l; i++)
1031 *deb = (*deb >> (BitsStored - HighBitPosition - 1)) & pmask;
1035 else // Pixels are signed
1037 // smask : to check the 'sign' when BitsStored != BitsAllocated
1038 uint16_t smask = 0x0001;
1039 smask = smask << ( 16 - (BitsAllocated - BitsStored + 1) );
1040 // nmask : to propagate sign bit on negative values
1041 int16_t nmask = (int16_t)0x8000;
1042 nmask = nmask >> ( BitsAllocated - BitsStored - 1 );
1044 for(int i = 0; i<l; i++)
1046 *deb = *deb >> (BitsStored - HighBitPosition - 1);
1049 *deb = *deb | nmask;
1053 *deb = *deb & pmask;
1059 else if ( BitsAllocated == 32 )
1061 // pmask : to mask the 'unused bits' (may contain overlays)
1062 uint32_t pmask = 0xffffffff;
1063 pmask = pmask >> ( BitsAllocated - BitsStored );
1065 uint32_t *deb = (uint32_t*)Raw;
1069 for(int i = 0; i<l; i++)
1071 *deb = (*deb >> (BitsStored - HighBitPosition - 1)) & pmask;
1077 // smask : to check the 'sign' when BitsStored != BitsAllocated
1078 uint32_t smask = 0x00000001;
1079 smask = smask >> ( 32 - (BitsAllocated - BitsStored +1 ));
1080 // nmask : to propagate sign bit on negative values
1081 int32_t nmask = 0x80000000;
1082 nmask = nmask >> ( BitsAllocated - BitsStored -1 );
1084 for(int i = 0; i<l; i++)
1086 *deb = *deb >> (BitsStored - HighBitPosition - 1);
1088 *deb = *deb | nmask;
1090 *deb = *deb & pmask;
1097 gdcmWarningMacro("Weird image (BitsAllocated !=8, 12, 16, 32)");
1098 throw FormatError( "Weird image !?" );
1105 * \brief Convert (Red plane, Green plane, Blue plane) to RGB pixels
1106 * \warning Works on all the frames at a time
1108 void PixelReadConvert::ConvertRGBPlanesToRGBPixels()
1110 gdcmWarningMacro("--> ConvertRGBPlanesToRGBPixels");
1112 uint8_t *localRaw = Raw;
1113 uint8_t *copyRaw = new uint8_t[ RawSize ];
1114 memmove( copyRaw, localRaw, RawSize );
1116 int l = XSize * YSize * ZSize;
1118 uint8_t *a = copyRaw;
1119 uint8_t *b = copyRaw + l;
1120 uint8_t *c = copyRaw + l + l;
1122 for (int j = 0; j < l; j++)
1124 *(localRaw++) = *(a++);
1125 *(localRaw++) = *(b++);
1126 *(localRaw++) = *(c++);
1132 * \brief Convert (cY plane, cB plane, cR plane) to RGB pixels
1133 * \warning Works on all the frames at a time
1135 void PixelReadConvert::ConvertYcBcRPlanesToRGBPixels()
1137 // Remarks for YBR newbees :
1138 // YBR_FULL works very much like RGB, i.e. three samples per pixel,
1139 // just the color space is YCbCr instead of RGB. This is particularly useful
1140 // for doppler ultrasound where most of the image is grayscale
1141 // (i.e. only populates the Y components) and Cb and Cr are mostly zero,
1142 // except for the few patches of color on the image.
1143 // On such images, RLE achieves a compression ratio that is much better
1144 // than the compression ratio on an equivalent RGB image.
1146 gdcmWarningMacro("--> ConvertYcBcRPlanesToRGBPixels");
1148 uint8_t *localRaw = Raw;
1149 uint8_t *copyRaw = new uint8_t[ RawSize ];
1150 memmove( copyRaw, localRaw, RawSize );
1152 // to see the tricks about YBR_FULL, YBR_FULL_422,
1153 // YBR_PARTIAL_422, YBR_ICT, YBR_RCT have a look at :
1154 // ftp://medical.nema.org/medical/dicom/final/sup61_ft.pdf
1155 // and be *very* affraid
1157 int l = XSize * YSize;
1158 int nbFrames = ZSize;
1160 uint8_t *a = copyRaw + 0;
1161 uint8_t *b = copyRaw + l;
1162 uint8_t *c = copyRaw + l+ l;
1165 /// We replaced easy to understand but time consuming floating point
1166 /// computations by the 'well known' integer computation counterpart
1168 /// http://lestourtereaux.free.fr/papers/data/yuvrgb.pdf
1169 /// for code optimisation.
1171 for ( int i = 0; i < nbFrames; i++ )
1173 for ( int j = 0; j < l; j++ )
1175 R = 38142 *(*a-16) + 52298 *(*c -128);
1176 G = 38142 *(*a-16) - 26640 *(*c -128) - 12845 *(*b -128);
1177 B = 38142 *(*a-16) + 66093 *(*b -128);
1186 if (R > 255) R = 255;
1187 if (G > 255) G = 255;
1188 if (B > 255) B = 255;
1190 *(localRaw++) = (uint8_t)R;
1191 *(localRaw++) = (uint8_t)G;
1192 *(localRaw++) = (uint8_t)B;
1201 /// \brief Deals with the color decoding i.e. handle:
1202 /// - R, G, B planes (as opposed to RGB pixels)
1203 /// - YBR (various) encodings.
1204 /// - LUT[s] (or "PALETTE COLOR").
1206 void PixelReadConvert::ConvertHandleColor()
1208 //////////////////////////////////
1209 // Deal with the color decoding i.e. handle:
1210 // - R, G, B planes (as opposed to RGB pixels)
1211 // - YBR (various) encodings.
1212 // - LUT[s] (or "PALETTE COLOR").
1214 // The classification in the color decoding schema is based on the blending
1215 // of two Dicom tags values:
1216 // * "Photometric Interpretation" for which we have the cases:
1217 // - [Photo A] MONOCHROME[1|2] pictures,
1218 // - [Photo B] RGB or YBR_FULL_422 (which acts as RGB),
1219 // - [Photo C] YBR_* (with the above exception of YBR_FULL_422)
1220 // - [Photo D] "PALETTE COLOR" which indicates the presence of LUT[s].
1221 // * "Planar Configuration" for which we have the cases:
1222 // - [Planar 0] 0 then Pixels are already RGB
1223 // - [Planar 1] 1 then we have 3 planes : R, G, B,
1224 // - [Planar 2] 2 then we have 1 gray Plane and 3 LUTs
1226 // Now in theory, one could expect some coherence when blending the above
1227 // cases. For example we should not encounter files belonging at the
1228 // time to case [Planar 0] and case [Photo D].
1229 // Alas, this was only theory ! Because in practice some odd (read ill
1230 // formated Dicom) files (e.g. gdcmData/US-PAL-8-10x-echo.dcm) we encounter:
1231 // - "Planar Configuration" = 0,
1232 // - "Photometric Interpretation" = "PALETTE COLOR".
1233 // Hence gdcm will use the folowing "heuristic" in order to be tolerant
1234 // towards Dicom-non-conformant files:
1235 // << whatever the "Planar Configuration" value might be, a
1236 // "Photometric Interpretation" set to "PALETTE COLOR" forces
1237 // a LUT intervention >>
1239 // Now we are left with the following handling of the cases:
1240 // - [Planar 0] OR [Photo A] no color decoding (since respectively
1241 // Pixels are already RGB and monochrome pictures have no color :),
1242 // - [Planar 1] AND [Photo B] handled with ConvertRGBPlanesToRGBPixels()
1243 // - [Planar 1] AND [Photo C] handled with ConvertYcBcRPlanesToRGBPixels()
1244 // - [Planar 2] OR [Photo D] requires LUT intervention.
1246 gdcmDebugMacro("--> ConvertHandleColor "
1247 << "Planar Configuration " << PlanarConfiguration );
1251 // [Planar 2] OR [Photo D]: LUT intervention done outside
1252 gdcmDebugMacro("--> RawRGB : LUT intervention done outside");
1256 if ( PlanarConfiguration == 1 )
1260 // [Planar 1] AND [Photo C] (remember YBR_FULL_422 acts as RGB)
1261 gdcmDebugMacro("--> YBRFull");
1262 ConvertYcBcRPlanesToRGBPixels();
1266 // [Planar 1] AND [Photo C]
1267 gdcmDebugMacro("--> YBRFull");
1268 ConvertRGBPlanesToRGBPixels();
1273 // When planarConf is 0, and RLELossless (forbidden by Dicom norm)
1274 // pixels need to be RGB-fyied anyway
1278 gdcmDebugMacro("--> RLE Lossless");
1279 ConvertRGBPlanesToRGBPixels();
1282 // In *normal *case, when planarConf is 0, pixels are already in RGB
1285 /// Computes the Pixels Size
1286 void PixelReadConvert::ComputeRawAndRGBSizes()
1288 int bitsAllocated = BitsAllocated;
1289 // Number of "Bits Allocated" is fixed to 16 when it's 12, since
1290 // in this case we will expand the image to 16 bits (see
1291 // \ref ReadAndDecompress12BitsTo16Bits() )
1292 if ( BitsAllocated == 12 )
1297 RawSize = XSize * YSize * ZSize
1298 * ( bitsAllocated / 8 )
1302 RGBSize = 3 * RawSize; // works for 8 and 16 bits per Pixel
1310 /// Allocates room for RGB Pixels
1311 void PixelReadConvert::AllocateRGB()
1315 RGB = new uint8_t[RGBSize];
1318 /// Allocates room for RAW Pixels
1319 void PixelReadConvert::AllocateRaw()
1323 Raw = new uint8_t[RawSize];
1326 //-----------------------------------------------------------------------------
1329 * \brief Print self.
1330 * @param indent Indentation string to be prepended during printing.
1331 * @param os Stream to print to.
1333 void PixelReadConvert::Print( std::ostream &os, std::string const &indent )
1336 << "--- Pixel information -------------------------"
1339 << "Pixel Data: offset " << PixelOffset
1340 << " x(" << std::hex << PixelOffset << std::dec
1341 << ") length " << PixelDataLength
1342 << " x(" << std::hex << PixelDataLength << std::dec
1343 << ")" << std::endl;
1345 if ( IsRLELossless )
1349 RLEInfo->Print( os, indent );
1353 gdcmWarningMacro("Set as RLE file but NO RLEinfo present.");
1357 if ( IsJPEG2000 || IsJPEGLossless || IsJPEGLossy || IsJPEGLS )
1361 JPEGInfo->Print( os, indent );
1365 gdcmWarningMacro("Set as JPEG file but NO JPEGinfo present.");
1370 //-----------------------------------------------------------------------------
1371 } // end namespace gdcm
1373 // Note to developpers :
1374 // Here is a very detailled post from David Clunie, on the troubles caused
1375 // 'non standard' LUT and LUT description
1376 // We shall have to take it into accound in our code.
1381 Subject: Problem with VOI LUTs in Agfa and Fuji CR and GE DX images, was Re: VOI LUT issues
1382 Date: Sun, 06 Feb 2005 17:13:40 GMT
1383 From: David Clunie <dclunie@dclunie.com>
1384 Reply-To: dclunie@dclunie.com
1385 Newsgroups: comp.protocols.dicom
1386 References: <1107553502.040221.189550@o13g2000cwo.googlegroups.com>
1388 > THE LUT that comes with [my] image claims to be 16-bit, but none of the
1389 > values goes higher than 4095. That being said, though, none of my
1390 > original pixel values goes higher than that, either. I have read
1391 > elsewhere on this group that when that happens you are supposed to
1392 > adjust the LUT. Can someone be more specific? There was a thread from
1393 > 2002 where Marco and David were mentioning doing precisely that.
1400 You have encountered the well known "we know what the standard says but
1401 we are going to ignore it and do what we have been doing for almost
1402 a decade regardless" CR vendor bug. Agfa started this, but they are not
1403 the only vendor doing this now; GE and Fuji may have joined the club.
1405 Sadly, one needs to look at the LUT Data, figure out what the maximum
1406 value actually encoded is, and find the next highest power of 2 (e.g.
1407 212 in this case), to figure out what the range of the data is
1408 supposed to be. I have assumed that if the maximum value in the LUT
1409 data is less than a power of 2 minus 1 (e.g. 0xebc) then the intent
1410 of the vendor was not to use the maximum available grayscale range
1411 of the display (e.g. the maximum is 0xfff in this case). An alternative
1412 would be to scale to the actual maximum rather than a power of two.
1414 Very irritating, and in theory not totally reliable if one really
1415 intended the full 16 bits and only used, say 15, but that is extremely
1416 unlikely since everything would be too dark, and this heuristic
1419 There has never been anything in the standard that describes having
1420 to go through these convolutions. Since the only value in the
1421 standard that describes the bit depth of the LUT values is LUT
1422 Descriptor value 3 and that is (usually) always required to be
1423 either 8 or 16, it mystifies me how the creators' of these images
1424 imagine that the receiver is going to divine the range that is intended. Further, the standard is quite explicit that this 3rd
1425 value defines the range of LUT values, but as far as I am aware, all
1426 the vendors are ignoring the standard and indeed sending a third value
1429 This problem is not confined to CR, and is also seen with DX products.
1431 Typically I have seen:
1433 - Agfa CR, which usually (always ?) sends LUTs, values up to 0x0fff
1434 - Fuji CR, which occasionally send LUTs, values up to 0x03ff
1435 - GE DX, for presentation, which always have LUTs, up to 0x3fff
1437 Swissray, Siemens, Philips, Canon and Kodak never seem to send VOI LUTs
1438 at this point (which is a whole other problem). Note that the presence
1439 or absence of a VOI LUT as opposed to window values may be configurable
1440 on the modality in some cases, and I have just looked at what I happen
1441 to have received from a myriad of sites over whose configuration I have
1442 no control. This may be why the majority of Fuji images have no VOI LUTs,
1443 but a few do (or it may be the Siemens system that these Fuji images went
1444 through that perhaps added it). I do have some test Hologic DX images that
1445 are not from a clinical site that do actually get this right (a value
1446 of 12 for the third value and a max of 0xfff).
1448 Since almost every vendor that I have encountered that encodes LUTs
1449 makes this mistake, perhaps it is time to amend the standard to warn
1450 implementor's of receivers and/or sanction this bad behavior. We have
1451 talked about this in the past in WG 6 but so far everyone has been
1452 reluctant to write into the standard such a comment. Maybe it is time
1453 to try again, since if one is not aware of this problem, one cannot
1454 effectively implement display using VOI LUTs, and there is a vast
1455 installed base to contend with.
1457 I did not check presentation states, in which VOI LUTs could also be
1458 encountered, for the prevalence of this mistake, nor did I look at the
1459 encoding of Modality LUT's, which are unusual. Nor did I check digital
1460 mammography images. I would be interested to hear from anyone who has.
1464 PS. The following older thread in this newsgroup discusses this:
1466 "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"
1468 PPS. From a historical perspective, the following may be of interest.
1470 In the original standard in 1993, all that was said about this was a
1471 reference to the corresponding such where Modality LUTs are described
1474 "The third value specifies the number of bits for each entry in the
1475 LUT Data. It shall take the value 8 or 16. The LUT Data shall be stored
1476 in a format equivalent to 8 or 16 bits allocated and high bit equal
1479 Since the high bit hint was not apparently explicit enough, a very
1480 early CP, CP 15 (submitted by Agfa as it happens), replaced this with:
1482 "The third value conveys the range of LUT entry values. It shall take
1483 the value 8 or 16, corresponding with the LUT entry value range of
1486 Note: The third value is not required for describing the
1487 LUT data and is only included for informational usage
1488 and for maintaining compatibility with ACRNEMA 2.0.
1490 The LUT Data contains the LUT entry values."
1492 That is how it read in the 1996, 1998 and 1999 editions.
1494 By the 2000 edition, Supplement 33 that introduced presentation states
1495 extensively reworked this entire section and tried to explain this in
1498 "The output range is from 0 to 2^n-1 where n is the third value of LUT
1499 Descriptor. This range is always unsigned."
1501 and also added a note to spell out what the output range meant in the
1504 "9. The output of the Window Center/Width or VOI LUT transformation
1505 is either implicitly scaled to the full range of the display device
1506 if there is no succeeding transformation defined, or implicitly scaled
1507 to the full input range of the succeeding transformation step (such as
1508 the Presentation LUT), if present. See C.11.6.1."
1510 It still reads this way in the 2004 edition.
1512 Note that LUTs in other applications than the general VOI LUT allow for
1513 values other than 8 or 16 in the third value of LUT descriptor to permit
1514 ranges other than 0 to 255 or 65535.
1516 In addition, the DX Image Module specializes the VOI LUT
1517 attributes as follows, in PS 3.3 section C.8.11.3.1.5 (added in Sup 32):
1519 "The third value specifies the number of bits for each entry in the LUT
1520 Data (analogous to ìbits storedî). It shall be between 10-16. The LUT
1521 Data shall be stored in a format equivalent to 16 ìbits allocatedî and
1522 ìhigh bitî equal to ìbits storedî - 1. The third value conveys the range
1523 of LUT entry values. These unsigned LUT entry values shall range between
1524 0 and 2^n-1, where n is the third value of the LUT Descriptor."
1526 So in the case of the GE DX for presentation images, the third value of
1527 LUT descriptor is allowed to be and probably should be 14 rather than 16.