1 /*=========================================================================
4 Module: $RCSfile: gdcmPixelReadConvert.cxx,v $
6 Date: $Date: 2005/06/30 14:50:05 $
7 Version: $Revision: 1.73 $
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, void *image_buffer, size_t lenght);
35 bool gdcm_read_JPEG2000_file (void* raw,
36 char *inputdata, size_t inputlength);
37 //-----------------------------------------------------------------------------
38 #define str2num(str, typeNum) *((typeNum *)(str))
40 //-----------------------------------------------------------------------------
41 // Constructor / Destructor
43 PixelReadConvert::PixelReadConvert()
55 /// Canonical Destructor
56 PixelReadConvert::~PixelReadConvert()
61 //-----------------------------------------------------------------------------
64 * \brief Predicate to know whether the image[s] (once Raw) is RGB.
65 * \note See comments of \ref ConvertHandleColor
67 bool PixelReadConvert::IsRawRGB()
70 || PlanarConfiguration == 2
78 * \brief Gets various usefull informations from the file header
79 * @param file gdcm::File pointer
81 void PixelReadConvert::GrabInformationsFromFile( File *file )
83 // Number of Bits Allocated for storing a Pixel is defaulted to 16
84 // when absent from the file.
85 BitsAllocated = file->GetBitsAllocated();
86 if ( BitsAllocated == 0 )
91 // Number of "Bits Stored", defaulted to number of "Bits Allocated"
92 // when absent from the file.
93 BitsStored = file->GetBitsStored();
94 if ( BitsStored == 0 )
96 BitsStored = BitsAllocated;
99 // High Bit Position, defaulted to "Bits Allocated" - 1
100 HighBitPosition = file->GetHighBitPosition();
101 if ( HighBitPosition == 0 )
103 HighBitPosition = BitsAllocated - 1;
106 XSize = file->GetXSize();
107 YSize = file->GetYSize();
108 ZSize = file->GetZSize();
109 SamplesPerPixel = file->GetSamplesPerPixel();
110 PixelSize = file->GetPixelSize();
111 PixelSign = file->IsSignedPixelData();
112 SwapCode = file->GetSwapCode();
113 std::string ts = file->GetTransferSyntax();
115 ( ! file->IsDicomV3() )
116 || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ImplicitVRLittleEndian
117 || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ImplicitVRLittleEndianDLXGE
118 || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ExplicitVRLittleEndian
119 || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ExplicitVRBigEndian
120 || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::DeflatedExplicitVRLittleEndian;
122 IsMPEG = Global::GetTS()->IsMPEG(ts);
123 IsJPEG2000 = Global::GetTS()->IsJPEG2000(ts);
124 IsJPEGLS = Global::GetTS()->IsJPEGLS(ts);
125 IsJPEGLossy = Global::GetTS()->IsJPEGLossy(ts);
126 IsJPEGLossless = Global::GetTS()->IsJPEGLossless(ts);
127 IsRLELossless = Global::GetTS()->IsRLELossless(ts);
129 PixelOffset = file->GetPixelOffset();
130 PixelDataLength = file->GetPixelAreaLength();
131 RLEInfo = file->GetRLEInfo();
132 JPEGInfo = file->GetJPEGInfo();
134 IsMonochrome = file->IsMonochrome();
135 IsMonochrome1 = file->IsMonochrome1();
136 IsPaletteColor = file->IsPaletteColor();
137 IsYBRFull = file->IsYBRFull();
139 PlanarConfiguration = file->GetPlanarConfiguration();
141 /////////////////////////////////////////////////////////////////
143 HasLUT = file->HasLUT();
146 // Just in case some access to a File element requires disk access.
147 LutRedDescriptor = file->GetEntryValue( 0x0028, 0x1101 );
148 LutGreenDescriptor = file->GetEntryValue( 0x0028, 0x1102 );
149 LutBlueDescriptor = file->GetEntryValue( 0x0028, 0x1103 );
151 // Depending on the value of Document::MAX_SIZE_LOAD_ELEMENT_VALUE
152 // [ refer to invocation of Document::SetMaxSizeLoadEntry() in
153 // Document::Document() ], the loading of the value (content) of a
154 // [Bin|Val]Entry occurence migth have been hindered (read simply NOT
155 // loaded). Hence, we first try to obtain the LUTs data from the file
156 // and when this fails we read the LUTs data directly from disk.
157 // \TODO Reading a [Bin|Val]Entry directly from disk is a kludge.
158 // We should NOT bypass the [Bin|Val]Entry class. Instead
159 // an access to an UNLOADED content of a [Bin|Val]Entry occurence
160 // (e.g. BinEntry::GetBinArea()) should force disk access from
161 // within the [Bin|Val]Entry class itself. The only problem
162 // is that the [Bin|Val]Entry is unaware of the FILE* is was
163 // parsed from. Fix that. FIXME.
166 file->LoadEntryBinArea(0x0028, 0x1201);
167 LutRedData = (uint8_t*)file->GetEntryBinArea( 0x0028, 0x1201 );
170 gdcmWarningMacro( "Unable to read Red Palette Color Lookup Table data" );
174 file->LoadEntryBinArea(0x0028, 0x1202);
175 LutGreenData = (uint8_t*)file->GetEntryBinArea(0x0028, 0x1202 );
178 gdcmWarningMacro( "Unable to read Green Palette Color Lookup Table data" );
182 file->LoadEntryBinArea(0x0028, 0x1203);
183 LutBlueData = (uint8_t*)file->GetEntryBinArea( 0x0028, 0x1203 );
186 gdcmWarningMacro( "Unable to read Blue Palette Color Lookup Table data" );
190 ComputeRawAndRGBSizes();
193 /// \brief Reads from disk and decompresses Pixels
194 bool PixelReadConvert::ReadAndDecompressPixelData( std::ifstream *fp )
196 // ComputeRawAndRGBSizes is already made by
197 // ::GrabInformationsFromfile. So, the structure sizes are
201 //////////////////////////////////////////////////
202 //// First stage: get our hands on the Pixel Data.
205 gdcmWarningMacro( "Unavailable file pointer." );
209 fp->seekg( PixelOffset, std::ios::beg );
210 if ( fp->fail() || fp->eof() )
212 gdcmWarningMacro( "Unable to find PixelOffset in file." );
218 //////////////////////////////////////////////////
219 //// Second stage: read from disk and decompress.
220 if ( BitsAllocated == 12 )
222 ReadAndDecompress12BitsTo16Bits( fp);
226 // This problem can be found when some obvious informations are found
227 // after the field containing the image data. In this case, these
228 // bad data are added to the size of the image (in the PixelDataLength
229 // variable). But RawSize is the right size of the image !
230 if ( PixelDataLength != RawSize )
232 gdcmWarningMacro( "Mismatch between PixelReadConvert : "
233 << PixelDataLength << " and RawSize : " << RawSize );
235 if ( PixelDataLength > RawSize )
237 fp->read( (char*)Raw, RawSize);
241 fp->read( (char*)Raw, PixelDataLength);
244 if ( fp->fail() || fp->eof())
246 gdcmWarningMacro( "Reading of Raw pixel data failed." );
250 else if ( IsRLELossless )
252 if ( ! RLEInfo->DecompressRLEFile( fp, Raw, XSize, YSize, ZSize, BitsAllocated ) )
254 gdcmWarningMacro( "RLE decompressor failed." );
260 //gdcmWarningMacro( "Sorry, MPEG not yet taken into account" );
262 //ReadMPEGFile(fp, Raw, PixelDataLength); // fp has already been seek to start of mpeg
267 // Default case concerns JPEG family
268 if ( ! ReadAndDecompressJPEGFile( fp ) )
270 gdcmWarningMacro( "JPEG decompressor failed." );
275 ////////////////////////////////////////////
276 //// Third stage: twigle the bytes and bits.
277 ConvertReorderEndianity();
278 ConvertReArrangeBits();
279 ConvertFixGreyLevels();
280 ConvertHandleColor();
285 /// Deletes Pixels Area
286 void PixelReadConvert::Squeeze()
302 * \brief Build the RGB image from the Raw image and the LUTs.
304 bool PixelReadConvert::BuildRGBImage()
308 // The job is already done.
314 // The job can't be done
321 // The job can't be done
325 gdcmWarningMacro( "--> BuildRGBImage" );
331 if ( BitsAllocated <= 8 )
333 uint8_t *localRGB = RGB;
334 for (size_t i = 0; i < RawSize; ++i )
337 *localRGB++ = LutRGBA[j];
338 *localRGB++ = LutRGBA[j+1];
339 *localRGB++ = LutRGBA[j+2];
343 else // deal with 16 bits pixels and 16 bits Palette color
345 uint16_t *localRGB = (uint16_t *)RGB;
346 for (size_t i = 0; i < RawSize/2; ++i )
348 j = ((uint16_t *)Raw)[i] * 4;
349 *localRGB++ = ((uint16_t *)LutRGBA)[j];
350 *localRGB++ = ((uint16_t *)LutRGBA)[j+1];
351 *localRGB++ = ((uint16_t *)LutRGBA)[j+2];
358 //-----------------------------------------------------------------------------
361 //-----------------------------------------------------------------------------
364 * \brief Read from file a 12 bits per pixel image and decompress it
365 * into a 16 bits per pixel image.
367 void PixelReadConvert::ReadAndDecompress12BitsTo16Bits( std::ifstream *fp )
368 throw ( FormatError )
370 int nbPixels = XSize * YSize;
371 uint16_t *localDecompres = (uint16_t*)Raw;
373 for( int p = 0; p < nbPixels; p += 2 )
377 fp->read( (char*)&b0, 1);
378 if ( fp->fail() || fp->eof() )
380 throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
381 "Unfound first block" );
384 fp->read( (char*)&b1, 1 );
385 if ( fp->fail() || fp->eof())
387 throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
388 "Unfound second block" );
391 fp->read( (char*)&b2, 1 );
392 if ( fp->fail() || fp->eof())
394 throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
395 "Unfound second block" );
398 // Two steps are necessary to please VC++
400 // 2 pixels 12bit = [0xABCDEF]
401 // 2 pixels 16bit = [0x0ABD] + [0x0FCE]
403 *localDecompres++ = ((b0 >> 4) << 8) + ((b0 & 0x0f) << 4) + (b1 & 0x0f);
405 *localDecompres++ = ((b2 & 0x0f) << 8) + ((b1 >> 4) << 4) + (b2 >> 4);
407 /// \todo JPR Troubles expected on Big-Endian processors ?
412 * \brief Reads from disk the Pixel Data of JPEG Dicom encapsulated
413 * file and decompress it.
414 * @param fp File Pointer
417 bool PixelReadConvert::ReadAndDecompressJPEGFile( std::ifstream *fp )
421 // make sure this is the right JPEG compression
422 assert( !IsJPEGLossless || !IsJPEGLossy || !IsJPEGLS );
423 // FIXME this is really ugly but it seems I have to load the complete
424 // jpeg2000 stream to use jasper:
425 // I don't think we'll ever be able to deal with multiple fragments properly
427 unsigned long inputlength = 0;
428 JPEGFragment *jpegfrag = JPEGInfo->GetFirstFragment();
431 inputlength += jpegfrag->GetLength();
432 jpegfrag = JPEGInfo->GetNextFragment();
434 gdcmAssertMacro( inputlength != 0);
435 uint8_t *inputdata = new uint8_t[inputlength];
436 char *pinputdata = (char*)inputdata;
437 jpegfrag = JPEGInfo->GetFirstFragment();
440 fp->seekg( jpegfrag->GetOffset(), std::ios::beg);
441 fp->read(pinputdata, jpegfrag->GetLength());
442 pinputdata += jpegfrag->GetLength();
443 jpegfrag = JPEGInfo->GetNextFragment();
445 // Warning the inputdata buffer is delete in the function
446 if ( ! gdcm_read_JPEG2000_file( Raw,
447 (char*)inputdata, inputlength ) )
451 // wow what happen, must be an error
456 // make sure this is the right JPEG compression
457 assert( !IsJPEGLossless || !IsJPEGLossy || !IsJPEG2000 );
458 // WARNING : JPEG-LS is NOT the 'classical' Jpeg Lossless :
459 // [JPEG-LS is the basis for new lossless/near-lossless compression
460 // standard for continuous-tone images intended for JPEG2000. The standard
461 // is based on the LOCO-I algorithm (LOw COmplexity LOssless COmpression
462 // for Images) developed at Hewlett-Packard Laboratories]
464 // see http://datacompression.info/JPEGLS.shtml
467 std::cerr << "count:" << JPEGInfo->GetFragmentCount() << std::endl;
468 unsigned long inputlength = 0;
469 JPEGFragment *jpegfrag = JPEGInfo->GetFirstFragment();
472 inputlength += jpegfrag->GetLength();
473 jpegfrag = JPEGInfo->GetNextFragment();
475 gdcmAssertMacro( inputlength != 0);
476 uint8_t *inputdata = new uint8_t[inputlength];
477 char *pinputdata = (char*)inputdata;
478 jpegfrag = JPEGInfo->GetFirstFragment();
481 fp->seekg( jpegfrag->GetOffset(), std::ios::beg);
482 fp->read(pinputdata, jpegfrag->GetLength());
483 pinputdata += jpegfrag->GetLength();
484 jpegfrag = JPEGInfo->GetNextFragment();
487 //fp->read((char*)Raw, PixelDataLength);
489 std::ofstream out("/tmp/jpegls.jpg");
490 out.write((char*)inputdata, inputlength);
495 gdcmWarningMacro( "Sorry, JPEG-LS not yet taken into account" );
496 fp->seekg( JPEGInfo->GetFirstFragment()->GetOffset(), std::ios::beg);
497 // if ( ! gdcm_read_JPEGLS_file( fp,Raw ) )
502 // make sure this is the right JPEG compression
503 assert( !IsJPEGLS || !IsJPEG2000 );
504 // Precompute the offset localRaw will be shifted with
505 int length = XSize * YSize * SamplesPerPixel;
506 int numberBytes = BitsAllocated / 8;
508 JPEGInfo->DecompressFromFile(fp, Raw, BitsStored, numberBytes, length );
514 * \brief Build Red/Green/Blue/Alpha LUT from File
515 * when (0028,0004),Photometric Interpretation = [PALETTE COLOR ]
516 * and (0028,1101),(0028,1102),(0028,1102)
517 * - xxx Palette Color Lookup Table Descriptor - are found
518 * and (0028,1201),(0028,1202),(0028,1202)
519 * - xxx Palette Color Lookup Table Data - are found
520 * \warning does NOT deal with :
521 * 0028 1100 Gray Lookup Table Descriptor (Retired)
522 * 0028 1221 Segmented Red Palette Color Lookup Table Data
523 * 0028 1222 Segmented Green Palette Color Lookup Table Data
524 * 0028 1223 Segmented Blue Palette Color Lookup Table Data
525 * no known Dicom reader deals with them :-(
526 * @return a RGBA Lookup Table
528 void PixelReadConvert::BuildLUTRGBA()
535 // http://www.barre.nom.fr/medical/dicom2/limitations.html#Color%20Lookup%20Tables
537 if ( ! IsPaletteColor )
542 if ( LutRedDescriptor == GDCM_UNFOUND
543 || LutGreenDescriptor == GDCM_UNFOUND
544 || LutBlueDescriptor == GDCM_UNFOUND )
546 gdcmWarningMacro( "(At least) a LUT Descriptor is missing" );
550 ////////////////////////////////////////////
551 // Extract the info from the LUT descriptors
552 int lengthR; // Red LUT length in Bytes
553 int debR; // Subscript of the first Lut Value
554 int nbitsR; // Lut item size (in Bits)
555 int nbRead; // nb of items in LUT descriptor (must be = 3)
557 nbRead = sscanf( LutRedDescriptor.c_str(),
559 &lengthR, &debR, &nbitsR );
562 gdcmWarningMacro( "Wrong Red LUT descriptor" );
564 int lengthG; // Green LUT length in Bytes
565 int debG; // Subscript of the first Lut Value
566 int nbitsG; // Lut item size (in Bits)
568 nbRead = sscanf( LutGreenDescriptor.c_str(),
570 &lengthG, &debG, &nbitsG );
573 gdcmWarningMacro( "Wrong Green LUT descriptor" );
576 int lengthB; // Blue LUT length in Bytes
577 int debB; // Subscript of the first Lut Value
578 int nbitsB; // Lut item size (in Bits)
579 nbRead = sscanf( LutRedDescriptor.c_str(),
581 &lengthB, &debB, &nbitsB );
584 gdcmWarningMacro( "Wrong Blue LUT descriptor" );
587 gdcmWarningMacro(" lengthR " << lengthR << " debR "
588 << debR << " nbitsR " << nbitsR);
589 gdcmWarningMacro(" lengthG " << lengthG << " debG "
590 << debG << " nbitsG " << nbitsG);
591 gdcmWarningMacro(" lengthB " << lengthB << " debB "
592 << debB << " nbitsB " << nbitsB);
594 if ( !lengthR ) // if = 2^16, this shall be 0 see : CP-143
596 if ( !lengthG ) // if = 2^16, this shall be 0
598 if ( !lengthB ) // if = 2^16, this shall be 0
601 ////////////////////////////////////////////////////////
603 if ( ( ! LutRedData ) || ( ! LutGreenData ) || ( ! LutBlueData ) )
605 gdcmWarningMacro( "(At least) a LUT is missing" );
609 // -------------------------------------------------------------
611 if ( BitsAllocated <= 8 )
613 // forge the 4 * 8 Bits Red/Green/Blue/Alpha LUT
614 LutRGBA = new uint8_t[ 1024 ]; // 256 * 4 (R, G, B, Alpha)
619 memset( LutRGBA, 0, 1024 );
622 if ( ( nbitsR == 16 ) && ( BitsAllocated == 8 ) )
624 // when LUT item size is different than pixel size
625 mult = 2; // high byte must be = low byte
629 // See PS 3.3-2003 C.11.1.1.2 p 619
633 // if we get a black image, let's just remove the '+1'
634 // from 'i*mult+1' and check again
635 // if it works, we shall have to check the 3 Palettes
636 // to see which byte is ==0 (first one, or second one)
638 // We give up the checking to avoid some (useless ?) overhead
639 // (optimistic asumption)
643 //take "Subscript of the first Lut Value" (debR,debG,debB) into account!
645 //FIXME : +1 : to get 'low value' byte
646 // Trouble expected on Big Endian Processors ?
647 // 16 BIts Per Pixel Palette Color to be swapped?
649 a = LutRGBA + 0 + debR;
650 for( i=0; i < lengthR; ++i )
652 *a = LutRedData[i*mult+1];
656 a = LutRGBA + 1 + debG;
657 for( i=0; i < lengthG; ++i)
659 *a = LutGreenData[i*mult+1];
663 a = LutRGBA + 2 + debB;
664 for(i=0; i < lengthB; ++i)
666 *a = LutBlueData[i*mult+1];
671 for(i=0; i < 256; ++i)
673 *a = 1; // Alpha component
679 // Probabely the same stuff is to be done for 16 Bits Pixels
680 // with 65536 entries LUT ?!?
681 // Still looking for accurate info on the web :-(
683 gdcmWarningMacro( "Sorry Palette Color Lookup Tables not yet dealt with"
684 << " for 16 Bits Per Pixel images" );
686 // forge the 4 * 16 Bits Red/Green/Blue/Alpha LUT
688 LutRGBA = (uint8_t *)new uint16_t[ 65536*4 ]; // 2^16 * 4 (R, G, B, Alpha)
691 memset( LutRGBA, 0, 65536*4*2 ); // 16 bits = 2 bytes ;-)
693 LutItemNumber = 65536;
699 //take "Subscript of the first Lut Value" (debR,debG,debB) into account!
701 a16 = (uint16_t*)LutRGBA + 0 + debR;
702 for( i=0; i < lengthR; ++i )
704 *a16 = ((uint16_t*)LutRedData)[i];
708 a16 = (uint16_t*)LutRGBA + 1 + debG;
709 for( i=0; i < lengthG; ++i)
711 *a16 = ((uint16_t*)LutGreenData)[i];
715 a16 = (uint16_t*)LutRGBA + 2 + debB;
716 for(i=0; i < lengthB; ++i)
718 *a16 = ((uint16_t*)LutBlueData)[i];
722 a16 = (uint16_t*)LutRGBA + 3 ;
723 for(i=0; i < 65536; ++i)
725 *a16 = 1; // Alpha component
728 /* Just to 'see' the LUT, at debug time
730 a16=(uint16_t*)LutRGBA;
731 for (int j=0;j<65536;j++)
733 std::cout << *a16 << " " << *(a16+1) << " "
734 << *(a16+2) << " " << *(a16+3) << std::endl;
742 * \brief Swap the bytes, according to \ref SwapCode.
744 void PixelReadConvert::ConvertSwapZone()
748 if ( BitsAllocated == 16 )
750 uint16_t *im16 = (uint16_t*)Raw;
758 for( i = 0; i < RawSize / 2; i++ )
760 im16[i]= (im16[i] >> 8) | (im16[i] << 8 );
764 gdcmWarningMacro("SwapCode value (16 bits) not allowed.");
767 else if ( BitsAllocated == 32 )
772 uint32_t *im32 = (uint32_t*)Raw;
778 for( i = 0; i < RawSize / 4; i++ )
780 low = im32[i] & 0x0000ffff; // 4321
781 high = im32[i] >> 16;
782 high = ( high >> 8 ) | ( high << 8 );
783 low = ( low >> 8 ) | ( low << 8 );
785 im32[i] = ( s32 << 16 ) | high;
789 for( i = 0; i < RawSize / 4; i++ )
791 low = im32[i] & 0x0000ffff; // 2143
792 high = im32[i] >> 16;
793 high = ( high >> 8 ) | ( high << 8 );
794 low = ( low >> 8 ) | ( low << 8 );
796 im32[i] = ( s32 << 16 ) | low;
800 for( i = 0; i < RawSize / 4; i++ )
802 low = im32[i] & 0x0000ffff; // 3412
803 high = im32[i] >> 16;
805 im32[i] = ( s32 << 16 ) | high;
809 gdcmWarningMacro("SwapCode value (32 bits) not allowed." );
815 * \brief Deal with endianness i.e. re-arange bytes inside the integer
817 void PixelReadConvert::ConvertReorderEndianity()
819 if ( BitsAllocated != 8 )
824 // Special kludge in order to deal with xmedcon broken images:
825 if ( BitsAllocated == 16
826 && BitsStored < BitsAllocated
829 int l = (int)( RawSize / ( BitsAllocated / 8 ) );
830 uint16_t *deb = (uint16_t *)Raw;
831 for(int i = 0; i<l; i++)
833 if ( *deb == 0xffff )
843 * \brief Deal with Grey levels i.e. re-arange them
844 * to have low values = dark, high values = bright
846 void PixelReadConvert::ConvertFixGreyLevels()
851 uint32_t i; // to please M$VC6
856 if ( BitsAllocated == 8 )
858 uint8_t *deb = (uint8_t *)Raw;
859 for (i=0; i<RawSize; i++)
867 if ( BitsAllocated == 16 )
870 for (j=0; j<BitsStored-1; j++)
872 mask = (mask << 1) +1; // will be fff when BitsStored=12
875 uint16_t *deb = (uint16_t *)Raw;
876 for (i=0; i<RawSize/2; i++)
886 if ( BitsAllocated == 8 )
888 uint8_t smask8 = 255;
889 uint8_t *deb = (uint8_t *)Raw;
890 for (i=0; i<RawSize; i++)
892 *deb = smask8 - *deb;
897 if ( BitsAllocated == 16 )
899 uint16_t smask16 = 65535;
900 uint16_t *deb = (uint16_t *)Raw;
901 for (i=0; i<RawSize/2; i++)
903 *deb = smask16 - *deb;
912 * \brief Re-arrange the bits within the bytes.
913 * @return Boolean always true
915 bool PixelReadConvert::ConvertReArrangeBits() throw ( FormatError )
918 if ( BitsStored != BitsAllocated )
920 int l = (int)( RawSize / ( BitsAllocated / 8 ) );
921 if ( BitsAllocated == 16 )
923 // pmask : to mask the 'unused bits' (may contain overlays)
924 uint16_t pmask = 0xffff;
925 pmask = pmask >> ( BitsAllocated - BitsStored );
927 uint16_t *deb = (uint16_t*)Raw;
929 if ( !PixelSign ) // Pixels are unsigned
931 for(int i = 0; i<l; i++)
933 *deb = (*deb >> (BitsStored - HighBitPosition - 1)) & pmask;
937 else // Pixels are signed
939 // smask : to check the 'sign' when BitsStored != BitsAllocated
940 uint16_t smask = 0x0001;
941 smask = smask << ( 16 - (BitsAllocated - BitsStored + 1) );
942 // nmask : to propagate sign bit on negative values
943 int16_t nmask = (int16_t)0x8000;
944 nmask = nmask >> ( BitsAllocated - BitsStored - 1 );
946 std::cout << "BitsStored " << BitsStored
947 << " BitsAllocated " << BitsAllocated
949 std::cout << std::hex << "pmask " << pmask
950 << " smask " << smask
951 << " nmask " << nmask
954 for(int i = 0; i<l; i++)
956 *deb = *deb >> (BitsStored - HighBitPosition - 1);
969 else if ( BitsAllocated == 32 )
971 // pmask : to mask the 'unused bits' (may contain overlays)
972 uint32_t pmask = 0xffffffff;
973 pmask = pmask >> ( BitsAllocated - BitsStored );
975 uint32_t *deb = (uint32_t*)Raw;
979 for(int i = 0; i<l; i++)
981 *deb = (*deb >> (BitsStored - HighBitPosition - 1)) & pmask;
987 // smask : to check the 'sign' when BitsStored != BitsAllocated
988 uint32_t smask = 0x00000001;
989 smask = smask >> ( 32 - (BitsAllocated - BitsStored +1 ));
990 // nmask : to propagate sign bit on negative values
991 int32_t nmask = 0x80000000;
992 nmask = nmask >> ( BitsAllocated - BitsStored -1 );
994 for(int i = 0; i<l; i++)
996 *deb = *deb >> (BitsStored - HighBitPosition - 1);
1000 *deb = *deb & pmask;
1007 gdcmWarningMacro("Weird image (BitsAllocated !=8, 12, 16, 32)");
1008 throw FormatError( "Weird image !?" );
1015 * \brief Convert (Red plane, Green plane, Blue plane) to RGB pixels
1016 * \warning Works on all the frames at a time
1018 void PixelReadConvert::ConvertRGBPlanesToRGBPixels()
1020 gdcmWarningMacro("--> ConvertRGBPlanesToRGBPixels");
1022 uint8_t *localRaw = Raw;
1023 uint8_t *copyRaw = new uint8_t[ RawSize ];
1024 memmove( copyRaw, localRaw, RawSize );
1026 int l = XSize * YSize * ZSize;
1028 uint8_t *a = copyRaw;
1029 uint8_t *b = copyRaw + l;
1030 uint8_t *c = copyRaw + l + l;
1032 for (int j = 0; j < l; j++)
1034 *(localRaw++) = *(a++);
1035 *(localRaw++) = *(b++);
1036 *(localRaw++) = *(c++);
1042 * \brief Convert (cY plane, cB plane, cR plane) to RGB pixels
1043 * \warning Works on all the frames at a time
1045 void PixelReadConvert::ConvertYcBcRPlanesToRGBPixels()
1047 // Remarks for YBR newbees :
1048 // YBR_FULL works very much like RGB, i.e. three samples per pixel,
1049 // just the color space is YCbCr instead of RGB. This is particularly useful
1050 // for doppler ultrasound where most of the image is grayscale
1051 // (i.e. only populates the Y components) and Cb and Cr are mostly zero,
1052 // except for the few patches of color on the image.
1053 // On such images, RLE achieves a compression ratio that is much better
1054 // than the compression ratio on an equivalent RGB image.
1056 gdcmWarningMacro("--> ConvertYcBcRPlanesToRGBPixels");
1058 uint8_t *localRaw = Raw;
1059 uint8_t *copyRaw = new uint8_t[ RawSize ];
1060 memmove( copyRaw, localRaw, RawSize );
1062 // to see the tricks about YBR_FULL, YBR_FULL_422,
1063 // YBR_PARTIAL_422, YBR_ICT, YBR_RCT have a look at :
1064 // ftp://medical.nema.org/medical/dicom/final/sup61_ft.pdf
1065 // and be *very* affraid
1067 int l = XSize * YSize;
1068 int nbFrames = ZSize;
1070 uint8_t *a = copyRaw + 0;
1071 uint8_t *b = copyRaw + l;
1072 uint8_t *c = copyRaw + l+ l;
1075 /// We replaced easy to understand but time consuming floating point
1076 /// computations by the 'well known' integer computation counterpart
1078 /// http://lestourtereaux.free.fr/papers/data/yuvrgb.pdf
1079 /// for code optimisation.
1081 for ( int i = 0; i < nbFrames; i++ )
1083 for ( int j = 0; j < l; j++ )
1085 R = 38142 *(*a-16) + 52298 *(*c -128);
1086 G = 38142 *(*a-16) - 26640 *(*c -128) - 12845 *(*b -128);
1087 B = 38142 *(*a-16) + 66093 *(*b -128);
1096 if (R > 255) R = 255;
1097 if (G > 255) G = 255;
1098 if (B > 255) B = 255;
1100 *(localRaw++) = (uint8_t)R;
1101 *(localRaw++) = (uint8_t)G;
1102 *(localRaw++) = (uint8_t)B;
1111 /// \brief Deals with the color decoding i.e. handle:
1112 /// - R, G, B planes (as opposed to RGB pixels)
1113 /// - YBR (various) encodings.
1114 /// - LUT[s] (or "PALETTE COLOR").
1116 void PixelReadConvert::ConvertHandleColor()
1118 //////////////////////////////////
1119 // Deal with the color decoding i.e. handle:
1120 // - R, G, B planes (as opposed to RGB pixels)
1121 // - YBR (various) encodings.
1122 // - LUT[s] (or "PALETTE COLOR").
1124 // The classification in the color decoding schema is based on the blending
1125 // of two Dicom tags values:
1126 // * "Photometric Interpretation" for which we have the cases:
1127 // - [Photo A] MONOCHROME[1|2] pictures,
1128 // - [Photo B] RGB or YBR_FULL_422 (which acts as RGB),
1129 // - [Photo C] YBR_* (with the above exception of YBR_FULL_422)
1130 // - [Photo D] "PALETTE COLOR" which indicates the presence of LUT[s].
1131 // * "Planar Configuration" for which we have the cases:
1132 // - [Planar 0] 0 then Pixels are already RGB
1133 // - [Planar 1] 1 then we have 3 planes : R, G, B,
1134 // - [Planar 2] 2 then we have 1 gray Plane and 3 LUTs
1136 // Now in theory, one could expect some coherence when blending the above
1137 // cases. For example we should not encounter files belonging at the
1138 // time to case [Planar 0] and case [Photo D].
1139 // Alas, this was only theory ! Because in practice some odd (read ill
1140 // formated Dicom) files (e.g. gdcmData/US-PAL-8-10x-echo.dcm) we encounter:
1141 // - "Planar Configuration" = 0,
1142 // - "Photometric Interpretation" = "PALETTE COLOR".
1143 // Hence gdcm will use the folowing "heuristic" in order to be tolerant
1144 // towards Dicom-non-conformant files:
1145 // << whatever the "Planar Configuration" value might be, a
1146 // "Photometric Interpretation" set to "PALETTE COLOR" forces
1147 // a LUT intervention >>
1149 // Now we are left with the following handling of the cases:
1150 // - [Planar 0] OR [Photo A] no color decoding (since respectively
1151 // Pixels are already RGB and monochrome pictures have no color :),
1152 // - [Planar 1] AND [Photo B] handled with ConvertRGBPlanesToRGBPixels()
1153 // - [Planar 1] AND [Photo C] handled with ConvertYcBcRPlanesToRGBPixels()
1154 // - [Planar 2] OR [Photo D] requires LUT intervention.
1156 gdcmWarningMacro("--> ConvertHandleColor"
1157 << "Planar Configuration " << PlanarConfiguration );
1161 // [Planar 2] OR [Photo D]: LUT intervention done outside
1162 gdcmWarningMacro("--> RawRGB : LUT intervention done outside");
1166 if ( PlanarConfiguration == 1 )
1170 // [Planar 1] AND [Photo C] (remember YBR_FULL_422 acts as RGB)
1171 gdcmWarningMacro("--> YBRFull");
1172 ConvertYcBcRPlanesToRGBPixels();
1176 // [Planar 1] AND [Photo C]
1177 gdcmWarningMacro("--> YBRFull");
1178 ConvertRGBPlanesToRGBPixels();
1183 // When planarConf is 0, and RLELossless (forbidden by Dicom norm)
1184 // pixels need to be RGB-fyied anyway
1188 gdcmWarningMacro("--> RLE Lossless");
1189 ConvertRGBPlanesToRGBPixels();
1192 // In *normal *case, when planarConf is 0, pixels are already in RGB
1195 /// Computes the Pixels Size
1196 void PixelReadConvert::ComputeRawAndRGBSizes()
1198 int bitsAllocated = BitsAllocated;
1199 // Number of "Bits Allocated" is fixed to 16 when it's 12, since
1200 // in this case we will expand the image to 16 bits (see
1201 // \ref ReadAndDecompress12BitsTo16Bits() )
1202 if ( BitsAllocated == 12 )
1207 RawSize = XSize * YSize * ZSize
1208 * ( bitsAllocated / 8 )
1212 RGBSize = 3 * RawSize; // works for 8 and 16 bits per Pixel
1220 /// Allocates room for RGB Pixels
1221 void PixelReadConvert::AllocateRGB()
1225 RGB = new uint8_t[RGBSize];
1228 /// Allocates room for RAW Pixels
1229 void PixelReadConvert::AllocateRaw()
1233 Raw = new uint8_t[RawSize];
1236 //-----------------------------------------------------------------------------
1239 * \brief Print self.
1240 * @param indent Indentation string to be prepended during printing.
1241 * @param os Stream to print to.
1243 void PixelReadConvert::Print( std::ostream &os, std::string const &indent )
1246 << "--- Pixel information -------------------------"
1249 << "Pixel Data: offset " << PixelOffset
1250 << " x(" << std::hex << PixelOffset << std::dec
1251 << ") length " << PixelDataLength
1252 << " x(" << std::hex << PixelDataLength << std::dec
1253 << ")" << std::endl;
1255 if ( IsRLELossless )
1259 RLEInfo->Print( os, indent );
1263 gdcmWarningMacro("Set as RLE file but NO RLEinfo present.");
1267 if ( IsJPEG2000 || IsJPEGLossless || IsJPEGLossy || IsJPEGLS )
1271 JPEGInfo->Print( os, indent );
1275 gdcmWarningMacro("Set as JPEG file but NO JPEGinfo present.");
1280 //-----------------------------------------------------------------------------
1281 } // end namespace gdcm
1283 // NOTES on File internal calls
1285 // ---> GetImageData
1286 // ---> GetImageDataIntoVector
1287 // |---> GetImageDataIntoVectorRaw
1288 // | lut intervention
1290 // ---> GetImageDataRaw
1291 // ---> GetImageDataIntoVectorRaw