]> Creatis software - gdcm.git/blob - src/gdcmPixelReadConvert.cxx
to avoid warnings
[gdcm.git] / src / gdcmPixelReadConvert.cxx
1 /*=========================================================================
2
3   Program:   gdcm
4   Module:    $RCSfile: gdcmPixelReadConvert.cxx,v $
5   Language:  C++
6   Date:      $Date: 2005/06/23 09:18:37 $
7   Version:   $Revision: 1.72 $
8                                                                                 
9   Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
10   l'Image). All rights reserved. See Doc/License.txt or
11   http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
12                                                                                 
13      This software is distributed WITHOUT ANY WARRANTY; without even
14      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
15      PURPOSE.  See the above copyright notices for more information.
16                                                                                 
17 =========================================================================*/
18
19 #include "gdcmPixelReadConvert.h"
20 #include "gdcmDebug.h"
21 #include "gdcmFile.h"
22 #include "gdcmGlobal.h"
23 #include "gdcmTS.h"
24 #include "gdcmDocEntry.h"
25 #include "gdcmRLEFramesInfo.h"
26 #include "gdcmJPEGFragmentsInfo.h"
27
28 #include <fstream>
29 #include <stdio.h> //for sscanf
30
31 namespace gdcm
32 {
33
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))
39
40 //-----------------------------------------------------------------------------
41 // Constructor / Destructor
42 /// Constructor
43 PixelReadConvert::PixelReadConvert() 
44 {
45    RGB          = 0;
46    RGBSize      = 0;
47    Raw          = 0;
48    RawSize      = 0;
49    LutRGBA      = 0;
50    LutRedData   = 0;
51    LutGreenData = 0;
52    LutBlueData  = 0;
53 }
54
55 /// Canonical Destructor
56 PixelReadConvert::~PixelReadConvert() 
57 {
58    Squeeze();
59 }
60
61 //-----------------------------------------------------------------------------
62 // Public
63 /**
64  * \brief Predicate to know whether the image[s] (once Raw) is RGB.
65  * \note See comments of \ref ConvertHandleColor
66  */
67 bool PixelReadConvert::IsRawRGB()
68 {
69    if (   IsMonochrome
70        || PlanarConfiguration == 2
71        || IsPaletteColor )
72    {
73       return false;
74    }
75    return true;
76 }
77 /**
78  * \brief Gets various usefull informations from the file header
79  * @param file gdcm::File pointer
80  */
81 void PixelReadConvert::GrabInformationsFromFile( File *file )
82 {
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 )
87    {
88       BitsAllocated = 16;
89    }
90
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 )
95    {
96       BitsStored = BitsAllocated;
97    }
98
99    // High Bit Position, defaulted to "Bits Allocated" - 1
100    HighBitPosition = file->GetHighBitPosition();
101    if ( HighBitPosition == 0 )
102    {
103       HighBitPosition = BitsAllocated - 1;
104    }
105
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();
114    IsRaw =
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;
121
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);
128
129    PixelOffset     = file->GetPixelOffset();
130    PixelDataLength = file->GetPixelAreaLength();
131    RLEInfo         = file->GetRLEInfo();
132    JPEGInfo        = file->GetJPEGInfo();
133
134    IsMonochrome    = file->IsMonochrome();
135    IsMonochrome1   = file->IsMonochrome1();
136    IsPaletteColor  = file->IsPaletteColor();
137    IsYBRFull       = file->IsYBRFull();
138
139    PlanarConfiguration = file->GetPlanarConfiguration();
140
141    /////////////////////////////////////////////////////////////////
142    // LUT section:
143    HasLUT = file->HasLUT();
144    if ( HasLUT )
145    {
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 );
150    
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.
164    
165       // //// Red round
166       file->LoadEntryBinArea(0x0028, 0x1201);
167       LutRedData = (uint8_t*)file->GetEntryBinArea( 0x0028, 0x1201 );
168       if ( ! LutRedData )
169       {
170          gdcmWarningMacro( "Unable to read Red Palette Color Lookup Table data" );
171       }
172
173       // //// Green round:
174       file->LoadEntryBinArea(0x0028, 0x1202);
175       LutGreenData = (uint8_t*)file->GetEntryBinArea(0x0028, 0x1202 );
176       if ( ! LutGreenData)
177       {
178          gdcmWarningMacro( "Unable to read Green Palette Color Lookup Table data" );
179       }
180
181       // //// Blue round:
182       file->LoadEntryBinArea(0x0028, 0x1203);
183       LutBlueData = (uint8_t*)file->GetEntryBinArea( 0x0028, 0x1203 );
184       if ( ! LutBlueData )
185       {
186          gdcmWarningMacro( "Unable to read Blue Palette Color Lookup Table data" );
187       }
188    }
189
190    ComputeRawAndRGBSizes();
191 }
192
193 /// \brief Reads from disk and decompresses Pixels
194 bool PixelReadConvert::ReadAndDecompressPixelData( std::ifstream *fp )
195 {
196    // ComputeRawAndRGBSizes is already made by 
197    // ::GrabInformationsFromfile. So, the structure sizes are
198    // correct
199    Squeeze();
200
201    //////////////////////////////////////////////////
202    //// First stage: get our hands on the Pixel Data.
203    if ( !fp )
204    {
205       gdcmWarningMacro( "Unavailable file pointer." );
206       return false;
207    }
208
209    fp->seekg( PixelOffset, std::ios::beg );
210    if ( fp->fail() || fp->eof() )
211    {
212       gdcmWarningMacro( "Unable to find PixelOffset in file." );
213       return false;
214    }
215
216    AllocateRaw();
217
218    //////////////////////////////////////////////////
219    //// Second stage: read from disk and decompress.
220    if ( BitsAllocated == 12 )
221    {
222       ReadAndDecompress12BitsTo16Bits( fp);
223    }
224    else if ( IsRaw )
225    {
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 )
231       {
232          gdcmWarningMacro( "Mismatch between PixelReadConvert : "
233                             << PixelDataLength << " and RawSize : " << RawSize );
234       }
235       if ( PixelDataLength > RawSize )
236       {
237          fp->read( (char*)Raw, RawSize);
238       }
239       else
240       {
241          fp->read( (char*)Raw, PixelDataLength);
242       }
243
244       if ( fp->fail() || fp->eof())
245       {
246          gdcmWarningMacro( "Reading of Raw pixel data failed." );
247          return false;
248       }
249    } 
250    else if ( IsRLELossless )
251    {
252       if ( ! RLEInfo->DecompressRLEFile( fp, Raw, XSize, YSize, ZSize, BitsAllocated ) )
253       {
254          gdcmWarningMacro( "RLE decompressor failed." );
255          return false;
256       }
257    }
258    else if ( IsMPEG )
259    {
260       //gdcmWarningMacro( "Sorry, MPEG not yet taken into account" );
261       //return false;
262       //ReadMPEGFile(fp, Raw, PixelDataLength); // fp has already been seek to start of mpeg
263       return true;
264    }
265    else
266    {
267       // Default case concerns JPEG family
268       if ( ! ReadAndDecompressJPEGFile( fp ) )
269       {
270          gdcmWarningMacro( "JPEG decompressor failed." );
271          return false;
272       }
273    }
274
275    ////////////////////////////////////////////
276    //// Third stage: twigle the bytes and bits.
277    ConvertReorderEndianity();
278    ConvertReArrangeBits();
279    ConvertFixGreyLevels();
280    ConvertHandleColor();
281
282    return true;
283 }
284
285 /// Deletes Pixels Area
286 void PixelReadConvert::Squeeze() 
287 {
288    if ( RGB )
289       delete [] RGB;
290    RGB = 0;
291
292    if ( Raw )
293       delete [] Raw;
294    Raw = 0;
295
296    if ( LutRGBA )
297       delete [] LutRGBA;
298    LutRGBA = 0;
299 }
300
301 /**
302  * \brief Build the RGB image from the Raw image and the LUTs.
303  */
304 bool PixelReadConvert::BuildRGBImage()
305 {
306    if ( RGB )
307    {
308       // The job is already done.
309       return true;
310    }
311
312    if ( ! Raw )
313    {
314       // The job can't be done
315       return false;
316    }
317
318    BuildLUTRGBA();
319    if ( ! LutRGBA )
320    {
321       // The job can't be done
322       return false;
323    }
324
325    gdcmWarningMacro( "--> BuildRGBImage" );
326                                                                                 
327    // Build RGB Pixels
328    AllocateRGB();
329    
330    int j;
331    if ( BitsAllocated <= 8 )
332    {
333       uint8_t *localRGB = RGB;
334       for (size_t i = 0; i < RawSize; ++i )
335       {
336          j  = Raw[i] * 4;
337          *localRGB++ = LutRGBA[j];
338          *localRGB++ = LutRGBA[j+1];
339          *localRGB++ = LutRGBA[j+2];
340       }
341     }
342  
343     else  // deal with 16 bits pixels and 16 bits Palette color
344     {
345       uint16_t *localRGB = (uint16_t *)RGB;
346       for (size_t i = 0; i < RawSize/2; ++i )
347       {
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];
352       } 
353     }
354  
355    return true;
356 }
357
358 //-----------------------------------------------------------------------------
359 // Protected
360
361 //-----------------------------------------------------------------------------
362 // Private
363 /**
364  * \brief Read from file a 12 bits per pixel image and decompress it
365  *        into a 16 bits per pixel image.
366  */
367 void PixelReadConvert::ReadAndDecompress12BitsTo16Bits( std::ifstream *fp )
368                throw ( FormatError )
369 {
370    int nbPixels = XSize * YSize;
371    uint16_t *localDecompres = (uint16_t*)Raw;
372
373    for( int p = 0; p < nbPixels; p += 2 )
374    {
375       uint8_t b0, b1, b2;
376
377       fp->read( (char*)&b0, 1);
378       if ( fp->fail() || fp->eof() )
379       {
380          throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
381                                 "Unfound first block" );
382       }
383
384       fp->read( (char*)&b1, 1 );
385       if ( fp->fail() || fp->eof())
386       {
387          throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
388                                 "Unfound second block" );
389       }
390
391       fp->read( (char*)&b2, 1 );
392       if ( fp->fail() || fp->eof())
393       {
394          throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
395                                 "Unfound second block" );
396       }
397
398       // Two steps are necessary to please VC++
399       //
400       // 2 pixels 12bit =     [0xABCDEF]
401       // 2 pixels 16bit = [0x0ABD] + [0x0FCE]
402       //                        A                     B                 D
403       *localDecompres++ =  ((b0 >> 4) << 8) + ((b0 & 0x0f) << 4) + (b1 & 0x0f);
404       //                        F                     C                 E
405       *localDecompres++ =  ((b2 & 0x0f) << 8) + ((b1 >> 4) << 4) + (b2 >> 4);
406
407       /// \todo JPR Troubles expected on Big-Endian processors ?
408    }
409 }
410
411 /**
412  * \brief     Reads from disk the Pixel Data of JPEG Dicom encapsulated
413  *            file and decompress it.
414  * @param     fp File Pointer
415  * @return    Boolean
416  */
417 bool PixelReadConvert::ReadAndDecompressJPEGFile( std::ifstream *fp )
418 {
419    if ( IsJPEG2000 )
420    {
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
426
427       unsigned long inputlength = 0;
428       JPEGFragment *jpegfrag = JPEGInfo->GetFirstFragment();
429       while( jpegfrag )
430       {
431          inputlength += jpegfrag->GetLength();
432          jpegfrag = JPEGInfo->GetNextFragment();
433       }
434       gdcmAssertMacro( inputlength != 0);
435       uint8_t *inputdata = new uint8_t[inputlength];
436       char *pinputdata = (char*)inputdata;
437       jpegfrag = JPEGInfo->GetFirstFragment();
438       while( jpegfrag )
439       {
440          fp->seekg( jpegfrag->GetOffset(), std::ios::beg);
441          fp->read(pinputdata, jpegfrag->GetLength());
442          pinputdata += jpegfrag->GetLength();
443          jpegfrag = JPEGInfo->GetNextFragment();
444       }
445       // Warning the inputdata buffer is delete in the function
446       if ( ! gdcm_read_JPEG2000_file( Raw, 
447           (char*)inputdata, inputlength ) )
448       {
449          return true;
450       }
451       // wow what happen, must be an error
452       return false;
453    }
454    else if ( IsJPEGLS )
455    {
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]
463    //
464    // see http://datacompression.info/JPEGLS.shtml
465    //
466 #if 0
467    std::cerr << "count:" << JPEGInfo->GetFragmentCount() << std::endl;
468       unsigned long inputlength = 0;
469       JPEGFragment *jpegfrag = JPEGInfo->GetFirstFragment();
470       while( jpegfrag )
471       {
472          inputlength += jpegfrag->GetLength();
473          jpegfrag = JPEGInfo->GetNextFragment();
474       }
475       gdcmAssertMacro( inputlength != 0);
476       uint8_t *inputdata = new uint8_t[inputlength];
477       char *pinputdata = (char*)inputdata;
478       jpegfrag = JPEGInfo->GetFirstFragment();
479       while( jpegfrag )
480       {
481          fp->seekg( jpegfrag->GetOffset(), std::ios::beg);
482          fp->read(pinputdata, jpegfrag->GetLength());
483          pinputdata += jpegfrag->GetLength();
484          jpegfrag = JPEGInfo->GetNextFragment();
485       }  
486       
487   //fp->read((char*)Raw, PixelDataLength);
488
489   std::ofstream out("/tmp/jpegls.jpg");
490   out.write((char*)inputdata, inputlength);
491   out.close();
492   delete[] inputdata;
493 #endif
494
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 ) )
498          return false;
499    }
500    else
501    {
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;
507
508      JPEGInfo->DecompressFromFile(fp, Raw, BitsStored, numberBytes, length );
509      return true;
510    }
511 }
512
513 /**
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
527  */
528 void PixelReadConvert::BuildLUTRGBA()
529 {
530    if ( LutRGBA )
531    {
532       return;
533    }
534    // Not so easy : see
535    // http://www.barre.nom.fr/medical/dicom2/limitations.html#Color%20Lookup%20Tables
536                                                                                 
537    if ( ! IsPaletteColor )
538    {
539       return;
540    }
541                                                                                 
542    if (   LutRedDescriptor   == GDCM_UNFOUND
543        || LutGreenDescriptor == GDCM_UNFOUND
544        || LutBlueDescriptor  == GDCM_UNFOUND )
545    {
546       gdcmWarningMacro( "(At least) a LUT Descriptor is missing" );
547       return;
548    }
549
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)
556
557    nbRead = sscanf( LutRedDescriptor.c_str(),
558                         "%d\\%d\\%d",
559                         &lengthR, &debR, &nbitsR );
560    if ( nbRead != 3 )
561    {
562       gdcmWarningMacro( "Wrong Red LUT descriptor" );
563    }                                                                                
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)
567
568    nbRead = sscanf( LutGreenDescriptor.c_str(),
569                     "%d\\%d\\%d",
570                     &lengthG, &debG, &nbitsG );  
571    if ( nbRead != 3 )
572    {
573       gdcmWarningMacro( "Wrong Green LUT descriptor" );
574    }
575                                                                                 
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(),
580                     "%d\\%d\\%d",
581                     &lengthB, &debB, &nbitsB );
582    if ( nbRead != 3 )
583    {
584       gdcmWarningMacro( "Wrong Blue LUT descriptor" );
585    }
586  
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);
593
594    if ( !lengthR ) // if = 2^16, this shall be 0 see : CP-143
595       lengthR=65536;
596    if ( !lengthG ) // if = 2^16, this shall be 0
597       lengthG=65536;
598    if ( !lengthB ) // if = 2^16, this shall be 0
599       lengthB=65536; 
600                                                                                 
601    ////////////////////////////////////////////////////////
602
603    if ( ( ! LutRedData ) || ( ! LutGreenData ) || ( ! LutBlueData ) )
604    {
605       gdcmWarningMacro( "(At least) a LUT is missing" );
606       return;
607    }
608
609    // -------------------------------------------------------------
610    
611    if ( BitsAllocated <= 8 )
612    {
613       // forge the 4 * 8 Bits Red/Green/Blue/Alpha LUT
614       LutRGBA = new uint8_t[ 1024 ]; // 256 * 4 (R, G, B, Alpha)
615       if ( !LutRGBA )
616          return;
617       LutItemNumber = 256;
618       LutItemSize   = 8;
619       memset( LutRGBA, 0, 1024 );
620                                                                                 
621       int mult;
622       if ( ( nbitsR == 16 ) && ( BitsAllocated == 8 ) )
623       {
624          // when LUT item size is different than pixel size
625          mult = 2; // high byte must be = low byte
626       }
627       else
628       {
629          // See PS 3.3-2003 C.11.1.1.2 p 619
630          mult = 1;
631       }
632                                                                                 
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)
637       // and fix the code
638       // We give up the checking to avoid some (useless ?) overhead
639       // (optimistic asumption)
640       int i;
641       uint8_t *a;
642
643       //take "Subscript of the first Lut Value" (debR,debG,debB) into account!
644
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?
648
649       a = LutRGBA + 0 + debR;
650       for( i=0; i < lengthR; ++i )
651       {
652          *a = LutRedData[i*mult+1]; 
653          a += 4;
654       }
655                                                                                 
656       a = LutRGBA + 1 + debG;
657       for( i=0; i < lengthG; ++i)
658       {
659          *a = LutGreenData[i*mult+1];
660          a += 4;
661       }
662                                                                                 
663       a = LutRGBA + 2 + debB;
664       for(i=0; i < lengthB; ++i)
665       {
666          *a = LutBlueData[i*mult+1];
667          a += 4;
668       }
669                                     
670       a = LutRGBA + 3 ;
671       for(i=0; i < 256; ++i)
672       {
673          *a = 1; // Alpha component
674          a += 4;
675       }
676    }
677    else
678    {
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 :-(
682
683       gdcmWarningMacro( "Sorry Palette Color Lookup Tables not yet dealt with"
684                          << " for 16 Bits Per Pixel images" );
685
686       // forge the 4 * 16 Bits Red/Green/Blue/Alpha LUT
687
688       LutRGBA = (uint8_t *)new uint16_t[ 65536*4 ]; // 2^16 * 4 (R, G, B, Alpha)
689       if ( !LutRGBA )
690          return;
691       memset( LutRGBA, 0, 65536*4*2 );  // 16 bits = 2 bytes ;-)
692
693       LutItemNumber = 65536;
694       LutItemSize   = 16;
695
696       int i;
697       uint16_t *a16;
698
699       //take "Subscript of the first Lut Value" (debR,debG,debB) into account!
700
701       a16 = (uint16_t*)LutRGBA + 0 + debR;
702       for( i=0; i < lengthR; ++i )
703       {
704          *a16 = ((uint16_t*)LutRedData)[i];
705          a16 += 4;
706       }
707                                                                               
708       a16 = (uint16_t*)LutRGBA + 1 + debG;
709       for( i=0; i < lengthG; ++i)
710       {
711          *a16 = ((uint16_t*)LutGreenData)[i];
712          a16 += 4;
713       }
714                                                                                 
715       a16 = (uint16_t*)LutRGBA + 2 + debB;
716       for(i=0; i < lengthB; ++i)
717       {
718          *a16 = ((uint16_t*)LutBlueData)[i];
719          a16 += 4;
720       }
721                                                                              
722       a16 = (uint16_t*)LutRGBA + 3 ;
723       for(i=0; i < 65536; ++i)
724       {
725          *a16 = 1; // Alpha component
726          a16 += 4;
727       }
728 /* Just to 'see' the LUT, at debug time
729
730       a16=(uint16_t*)LutRGBA;
731       for (int j=0;j<65536;j++)
732       {
733          std::cout << *a16     << " " << *(a16+1) << " "
734                    << *(a16+2) << " " << *(a16+3) << std::endl;
735          a16+=4;
736       }
737 */
738    }
739 }
740
741 /**
742  * \brief Swap the bytes, according to \ref SwapCode.
743  */
744 void PixelReadConvert::ConvertSwapZone()
745 {
746    unsigned int i;
747
748    if ( BitsAllocated == 16 )
749    {
750       uint16_t *im16 = (uint16_t*)Raw;
751       switch( SwapCode )
752       {
753          case 1234:
754             break;
755          case 3412:
756          case 2143:
757          case 4321:
758             for( i = 0; i < RawSize / 2; i++ )
759             {
760                im16[i]= (im16[i] >> 8) | (im16[i] << 8 );
761             }
762             break;
763          default:
764             gdcmWarningMacro("SwapCode value (16 bits) not allowed.");
765       }
766    }
767    else if ( BitsAllocated == 32 )
768    {
769       uint32_t s32;
770       uint16_t high;
771       uint16_t low;
772       uint32_t *im32 = (uint32_t*)Raw;
773       switch ( SwapCode )
774       {
775          case 1234:
776             break;
777          case 4321:
778             for( i = 0; i < RawSize / 4; i++ )
779             {
780                low     = im32[i] & 0x0000ffff;  // 4321
781                high    = im32[i] >> 16;
782                high    = ( high >> 8 ) | ( high << 8 );
783                low     = ( low  >> 8 ) | ( low  << 8 );
784                s32     = low;
785                im32[i] = ( s32 << 16 ) | high;
786             }
787             break;
788          case 2143:
789             for( i = 0; i < RawSize / 4; i++ )
790             {
791                low     = im32[i] & 0x0000ffff;   // 2143
792                high    = im32[i] >> 16;
793                high    = ( high >> 8 ) | ( high << 8 );
794                low     = ( low  >> 8 ) | ( low  << 8 );
795                s32     = high;
796                im32[i] = ( s32 << 16 ) | low;
797             }
798             break;
799          case 3412:
800             for( i = 0; i < RawSize / 4; i++ )
801             {
802                low     = im32[i] & 0x0000ffff; // 3412
803                high    = im32[i] >> 16;
804                s32     = low;
805                im32[i] = ( s32 << 16 ) | high;
806             }
807             break;
808          default:
809             gdcmWarningMacro("SwapCode value (32 bits) not allowed." );
810       }
811    }
812 }
813
814 /**
815  * \brief Deal with endianness i.e. re-arange bytes inside the integer
816  */
817 void PixelReadConvert::ConvertReorderEndianity()
818 {
819    if ( BitsAllocated != 8 )
820    {
821       ConvertSwapZone();
822    }
823
824    // Special kludge in order to deal with xmedcon broken images:
825    if ( BitsAllocated == 16
826      && BitsStored < BitsAllocated
827      && !PixelSign )
828    {
829       int l = (int)( RawSize / ( BitsAllocated / 8 ) );
830       uint16_t *deb = (uint16_t *)Raw;
831       for(int i = 0; i<l; i++)
832       {
833          if ( *deb == 0xffff )
834          {
835            *deb = 0;
836          }
837          deb++;
838       }
839    }
840 }
841
842 /**
843  * \brief Deal with Grey levels i.e. re-arange them
844  *        to have low values = dark, high values = bright
845  */
846 void PixelReadConvert::ConvertFixGreyLevels()
847 {
848    if (!IsMonochrome1)
849       return;
850
851    uint32_t i; // to please M$VC6
852    int16_t j;
853
854    if (!PixelSign)
855    {
856       if ( BitsAllocated == 8 )
857       {
858          uint8_t *deb = (uint8_t *)Raw;
859          for (i=0; i<RawSize; i++)      
860          {
861             *deb = 255 - *deb;
862             deb++;
863          }
864          return;
865       }
866
867       if ( BitsAllocated == 16 )
868       {
869          uint16_t mask =1;
870          for (j=0; j<BitsStored-1; j++)
871          {
872             mask = (mask << 1) +1; // will be fff when BitsStored=12
873          }
874
875          uint16_t *deb = (uint16_t *)Raw;
876          for (i=0; i<RawSize/2; i++)      
877          {
878             *deb = mask - *deb;
879             deb++;
880          }
881          return;
882        }
883    }
884    else
885    {
886       if ( BitsAllocated == 8 )
887       {
888          uint8_t smask8 = 255;
889          uint8_t *deb = (uint8_t *)Raw;
890          for (i=0; i<RawSize; i++)      
891          {
892             *deb = smask8 - *deb;
893             deb++;
894          }
895          return;
896       }
897       if ( BitsAllocated == 16 )
898       {
899          uint16_t smask16 = 65535;
900          uint16_t *deb = (uint16_t *)Raw;
901          for (i=0; i<RawSize/2; i++)      
902          {
903             *deb = smask16 - *deb;
904             deb++;
905          }
906          return;
907       }
908    }
909 }
910
911 /**
912  * \brief  Re-arrange the bits within the bytes.
913  * @return Boolean always true
914  */
915 bool PixelReadConvert::ConvertReArrangeBits() throw ( FormatError )
916 {
917
918    if ( BitsStored != BitsAllocated )
919    {
920       int l = (int)( RawSize / ( BitsAllocated / 8 ) );
921       if ( BitsAllocated == 16 )
922       {
923          // pmask : to mask the 'unused bits' (may contain overlays)
924          uint16_t pmask = 0xffff;
925          pmask = pmask >> ( BitsAllocated - BitsStored );
926
927          uint16_t *deb = (uint16_t*)Raw;
928
929          if ( !PixelSign )  // Pixels are unsigned
930          {
931             for(int i = 0; i<l; i++)
932             {   
933                *deb = (*deb >> (BitsStored - HighBitPosition - 1)) & pmask;
934                deb++;
935             }
936          }
937          else // Pixels are signed
938          {
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 = 0x8000;  
944             nmask = nmask >> ( BitsAllocated - BitsStored - 1 );
945 /*
946 std::cout  << "BitsStored "     << BitsStored 
947            << " BitsAllocated " << BitsAllocated 
948            << std::endl;
949 std::cout  << std::hex  << "pmask " << pmask
950            << " smask " << smask
951            << " nmask " << nmask
952            << std::endl;
953 */ 
954             for(int i = 0; i<l; i++)
955             {
956                *deb = *deb >> (BitsStored - HighBitPosition - 1);
957                if ( *deb & smask )
958                {
959                   *deb = *deb | nmask;
960                }
961                else
962                {
963                   *deb = *deb & pmask;
964                }
965                deb++;
966             }
967          }
968       }
969       else if ( BitsAllocated == 32 )
970       {
971          // pmask : to mask the 'unused bits' (may contain overlays)
972          uint32_t pmask = 0xffffffff;
973          pmask = pmask >> ( BitsAllocated - BitsStored );
974
975          uint32_t *deb = (uint32_t*)Raw;
976
977          if ( !PixelSign )
978          {
979             for(int i = 0; i<l; i++)
980             {             
981                *deb = (*deb >> (BitsStored - HighBitPosition - 1)) & pmask;
982                deb++;
983             }
984          }
985          else
986          {
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 );
993
994             for(int i = 0; i<l; i++)
995             {
996                *deb = *deb >> (BitsStored - HighBitPosition - 1);
997                if ( *deb & smask )
998                   *deb = *deb | nmask;
999                else
1000                   *deb = *deb & pmask;
1001                deb++;
1002             }
1003          }
1004       }
1005       else
1006       {
1007          gdcmWarningMacro("Weird image (BitsAllocated !=8, 12, 16, 32)");
1008          throw FormatError( "Weird image !?" );
1009       }
1010    }
1011    return true;
1012 }
1013
1014 /**
1015  * \brief   Convert (Red plane, Green plane, Blue plane) to RGB pixels
1016  * \warning Works on all the frames at a time
1017  */
1018 void PixelReadConvert::ConvertRGBPlanesToRGBPixels()
1019 {
1020    gdcmWarningMacro("--> ConvertRGBPlanesToRGBPixels");
1021
1022    uint8_t *localRaw = Raw;
1023    uint8_t *copyRaw = new uint8_t[ RawSize ];
1024    memmove( copyRaw, localRaw, RawSize );
1025
1026    int l = XSize * YSize * ZSize;
1027
1028    uint8_t *a = copyRaw;
1029    uint8_t *b = copyRaw + l;
1030    uint8_t *c = copyRaw + l + l;
1031
1032    for (int j = 0; j < l; j++)
1033    {
1034       *(localRaw++) = *(a++);
1035       *(localRaw++) = *(b++);
1036       *(localRaw++) = *(c++);
1037    }
1038    delete[] copyRaw;
1039 }
1040
1041 /**
1042  * \brief   Convert (cY plane, cB plane, cR plane) to RGB pixels
1043  * \warning Works on all the frames at a time
1044  */
1045 void PixelReadConvert::ConvertYcBcRPlanesToRGBPixels()
1046 {
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. 
1055  
1056    gdcmWarningMacro("--> ConvertYcBcRPlanesToRGBPixels");
1057    
1058    uint8_t *localRaw = Raw;
1059    uint8_t *copyRaw = new uint8_t[ RawSize ];
1060    memmove( copyRaw, localRaw, RawSize );
1061
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
1066    //
1067    int l        = XSize * YSize;
1068    int nbFrames = ZSize;
1069
1070    uint8_t *a = copyRaw + 0;
1071    uint8_t *b = copyRaw + l;
1072    uint8_t *c = copyRaw + l+ l;
1073    int32_t R, G, B;
1074
1075    ///  We replaced easy to understand but time consuming floating point
1076    ///  computations by the 'well known' integer computation counterpart
1077    ///  Refer to :
1078    ///            http://lestourtereaux.free.fr/papers/data/yuvrgb.pdf
1079    ///  for code optimisation.
1080
1081    for ( int i = 0; i < nbFrames; i++ )
1082    {
1083       for ( int j = 0; j < l; j++ )
1084       {
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);
1088
1089          R = (R+16384)>>15;
1090          G = (G+16384)>>15;
1091          B = (B+16384)>>15;
1092
1093          if (R < 0)   R = 0;
1094          if (G < 0)   G = 0;
1095          if (B < 0)   B = 0;
1096          if (R > 255) R = 255;
1097          if (G > 255) G = 255;
1098          if (B > 255) B = 255;
1099
1100          *(localRaw++) = (uint8_t)R;
1101          *(localRaw++) = (uint8_t)G;
1102          *(localRaw++) = (uint8_t)B;
1103          a++;
1104          b++;
1105          c++;
1106       }
1107    }
1108    delete[] copyRaw;
1109 }
1110
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").
1115
1116 void PixelReadConvert::ConvertHandleColor()
1117 {
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").
1123    //
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
1135    //
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 >>
1148    //
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.
1155
1156    gdcmWarningMacro("--> ConvertHandleColor"
1157                     << "Planar Configuration " << PlanarConfiguration );
1158
1159    if ( ! IsRawRGB() )
1160    {
1161       // [Planar 2] OR  [Photo D]: LUT intervention done outside
1162       gdcmWarningMacro("--> RawRGB : LUT intervention done outside");
1163       return;
1164    }
1165                                                                                 
1166    if ( PlanarConfiguration == 1 )
1167    {
1168       if ( IsYBRFull )
1169       {
1170          // [Planar 1] AND [Photo C] (remember YBR_FULL_422 acts as RGB)
1171          gdcmWarningMacro("--> YBRFull");
1172          ConvertYcBcRPlanesToRGBPixels();
1173       }
1174       else
1175       {
1176          // [Planar 1] AND [Photo C]
1177          gdcmWarningMacro("--> YBRFull");
1178          ConvertRGBPlanesToRGBPixels();
1179       }
1180       return;
1181    }
1182                                                                                 
1183    // When planarConf is 0, and RLELossless (forbidden by Dicom norm)
1184    // pixels need to be RGB-fyied anyway
1185
1186    if (IsRLELossless)
1187    { 
1188      gdcmWarningMacro("--> RLE Lossless");
1189      ConvertRGBPlanesToRGBPixels();
1190    }
1191
1192    // In *normal *case, when planarConf is 0, pixels are already in RGB
1193 }
1194
1195 /// Computes the Pixels Size
1196 void PixelReadConvert::ComputeRawAndRGBSizes()
1197 {
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 )
1203    {
1204       bitsAllocated = 16;
1205    }
1206                                                                                 
1207    RawSize =  XSize * YSize * ZSize
1208                      * ( bitsAllocated / 8 )
1209                      * SamplesPerPixel;
1210    if ( HasLUT )
1211    {
1212       RGBSize = 3 * RawSize; // works for 8 and 16 bits per Pixel
1213    }
1214    else
1215    {
1216       RGBSize = RawSize;
1217    }
1218 }
1219
1220 /// Allocates room for RGB Pixels
1221 void PixelReadConvert::AllocateRGB()
1222 {
1223   if ( RGB )
1224      delete [] RGB;
1225   RGB = new uint8_t[RGBSize];
1226 }
1227
1228 /// Allocates room for RAW Pixels
1229 void PixelReadConvert::AllocateRaw()
1230 {
1231   if ( Raw )
1232      delete [] Raw;
1233   Raw = new uint8_t[RawSize];
1234 }
1235
1236 //-----------------------------------------------------------------------------
1237 // Print
1238 /**
1239  * \brief        Print self.
1240  * @param indent Indentation string to be prepended during printing.
1241  * @param os     Stream to print to.
1242  */
1243 void PixelReadConvert::Print( std::ostream &os, std::string const &indent )
1244 {
1245    os << indent
1246       << "--- Pixel information -------------------------"
1247       << std::endl;
1248    os << indent
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;
1254
1255    if ( IsRLELossless )
1256    {
1257       if ( RLEInfo )
1258       {
1259          RLEInfo->Print( os, indent );
1260       }
1261       else
1262       {
1263          gdcmWarningMacro("Set as RLE file but NO RLEinfo present.");
1264       }
1265    }
1266
1267    if ( IsJPEG2000 || IsJPEGLossless || IsJPEGLossy || IsJPEGLS )
1268    {
1269       if ( JPEGInfo )
1270       {
1271          JPEGInfo->Print( os, indent );
1272       }
1273       else
1274       {
1275          gdcmWarningMacro("Set as JPEG file but NO JPEGinfo present.");
1276       }
1277    }
1278 }
1279
1280 //-----------------------------------------------------------------------------
1281 } // end namespace gdcm
1282
1283 // NOTES on File internal calls
1284 // User
1285 // ---> GetImageData
1286 //     ---> GetImageDataIntoVector
1287 //        |---> GetImageDataIntoVectorRaw
1288 //        | lut intervention
1289 // User
1290 // ---> GetImageDataRaw
1291 //     ---> GetImageDataIntoVectorRaw
1292