]> Creatis software - gdcm.git/blob - src/gdcmPixelReadConvert.cxx
Comments
[gdcm.git] / src / gdcmPixelReadConvert.cxx
1 /*=========================================================================
2
3   Program:   gdcm
4   Module:    $RCSfile: gdcmPixelReadConvert.cxx,v $
5   Language:  C++
6   Date:      $Date: 2005/06/20 17:12:03 $
7   Version:   $Revision: 1.68 $
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 LUT 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 LUT 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 LUT 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 dans 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    // Build RGB Pixels
326    AllocateRGB();
327    
328    int j;
329    if( BitsAllocated <= 8)
330    {
331       uint8_t *localRGB = RGB;
332       for (size_t i = 0; i < RawSize; ++i )
333       {
334          j  = Raw[i] * 4;
335          *localRGB++ = LutRGBA[j];
336          *localRGB++ = LutRGBA[j+1];
337          *localRGB++ = LutRGBA[j+2];
338       }
339     }
340  
341     else  // deal with 16 bits pixels and 16 bits Palette color
342     {
343       uint16_t *localRGB = (uint16_t *)RGB;
344       for (size_t i = 0; i < RawSize/2; ++i )
345       {
346          j  = ((uint16_t *)Raw)[i] * 4;
347          *localRGB++ = ((uint16_t *)LutRGBA)[j];
348          *localRGB++ = ((uint16_t *)LutRGBA)[j+1];
349          *localRGB++ = ((uint16_t *)LutRGBA)[j+2];
350       } 
351     }
352  
353    return true;
354 }
355
356 //-----------------------------------------------------------------------------
357 // Protected
358
359 //-----------------------------------------------------------------------------
360 // Private
361 /**
362  * \brief Read from file a 12 bits per pixel image and decompress it
363  *        into a 16 bits per pixel image.
364  */
365 void PixelReadConvert::ReadAndDecompress12BitsTo16Bits( std::ifstream *fp )
366                throw ( FormatError )
367 {
368    int nbPixels = XSize * YSize;
369    uint16_t *localDecompres = (uint16_t*)Raw;
370
371    for( int p = 0; p < nbPixels; p += 2 )
372    {
373       uint8_t b0, b1, b2;
374
375       fp->read( (char*)&b0, 1);
376       if ( fp->fail() || fp->eof() )
377       {
378          throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
379                                 "Unfound first block" );
380       }
381
382       fp->read( (char*)&b1, 1 );
383       if ( fp->fail() || fp->eof())
384       {
385          throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
386                                 "Unfound second block" );
387       }
388
389       fp->read( (char*)&b2, 1 );
390       if ( fp->fail() || fp->eof())
391       {
392          throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
393                                 "Unfound second block" );
394       }
395
396       // Two steps are necessary to please VC++
397       //
398       // 2 pixels 12bit =     [0xABCDEF]
399       // 2 pixels 16bit = [0x0ABD] + [0x0FCE]
400       //                        A                     B                 D
401       *localDecompres++ =  ((b0 >> 4) << 8) + ((b0 & 0x0f) << 4) + (b1 & 0x0f);
402       //                        F                     C                 E
403       *localDecompres++ =  ((b2 & 0x0f) << 8) + ((b1 >> 4) << 4) + (b2 >> 4);
404
405       /// \todo JPR Troubles expected on Big-Endian processors ?
406    }
407 }
408
409 /**
410  * \brief     Reads from disk the Pixel Data of JPEG Dicom encapsulated
411  *            file and decompress it.
412  * @param     fp File Pointer
413  * @return    Boolean
414  */
415 bool PixelReadConvert::ReadAndDecompressJPEGFile( std::ifstream *fp )
416 {
417    if ( IsJPEG2000 )
418    {
419      // make sure this is the right JPEG compression
420      assert( !IsJPEGLossless || !IsJPEGLossy || !IsJPEGLS );
421      // FIXME this is really ugly but it seems I have to load the complete
422      // jpeg2000 stream to use jasper:
423       // I don't think we'll ever be able to deal with multiple fragments properly
424
425       unsigned long inputlength = 0;
426       JPEGFragment *jpegfrag = JPEGInfo->GetFirstFragment();
427       while( jpegfrag )
428       {
429          inputlength += jpegfrag->GetLength();
430          jpegfrag = JPEGInfo->GetNextFragment();
431       }
432       gdcmAssertMacro( inputlength != 0);
433       uint8_t *inputdata = new uint8_t[inputlength];
434       char *pinputdata = (char*)inputdata;
435       jpegfrag = JPEGInfo->GetFirstFragment();
436       while( jpegfrag )
437       {
438          fp->seekg( jpegfrag->GetOffset(), std::ios::beg);
439          fp->read(pinputdata, jpegfrag->GetLength());
440          pinputdata += jpegfrag->GetLength();
441          jpegfrag = JPEGInfo->GetNextFragment();
442       }
443       // Warning the inputdata buffer is delete in the function
444       if ( ! gdcm_read_JPEG2000_file( Raw, 
445           (char*)inputdata, inputlength ) )
446       {
447          return true;
448       }
449       // wow what happen, must be an error
450       return false;
451    }
452    else if ( IsJPEGLS )
453    {
454      // make sure this is the right JPEG compression
455      assert( !IsJPEGLossless || !IsJPEGLossy || !IsJPEG2000 );
456    // WARNING : JPEG-LS is NOT the 'classical' Jpeg Lossless : 
457    // [JPEG-LS is the basis for new lossless/near-lossless compression
458    // standard for continuous-tone images intended for JPEG2000. The standard
459    // is based on the LOCO-I algorithm (LOw COmplexity LOssless COmpression
460    // for Images) developed at Hewlett-Packard Laboratories]
461    //
462    // see http://datacompression.info/JPEGLS.shtml
463    //
464 #if 0
465    std::cerr << "count:" << JPEGInfo->GetFragmentCount() << std::endl;
466       unsigned long inputlength = 0;
467       JPEGFragment *jpegfrag = JPEGInfo->GetFirstFragment();
468       while( jpegfrag )
469       {
470          inputlength += jpegfrag->GetLength();
471          jpegfrag = JPEGInfo->GetNextFragment();
472       }
473       gdcmAssertMacro( inputlength != 0);
474       uint8_t *inputdata = new uint8_t[inputlength];
475       char *pinputdata = (char*)inputdata;
476       jpegfrag = JPEGInfo->GetFirstFragment();
477       while( jpegfrag )
478       {
479          fp->seekg( jpegfrag->GetOffset(), std::ios::beg);
480          fp->read(pinputdata, jpegfrag->GetLength());
481          pinputdata += jpegfrag->GetLength();
482          jpegfrag = JPEGInfo->GetNextFragment();
483       }  
484       
485   //fp->read((char*)Raw, PixelDataLength);
486
487   std::ofstream out("/tmp/jpegls.jpg");
488   out.write((char*)inputdata, inputlength);
489   out.close();
490   delete[] inputdata;
491 #endif
492
493       gdcmWarningMacro( "Sorry, JPEG-LS not yet taken into account" );
494       fp->seekg( JPEGInfo->GetFirstFragment()->GetOffset(), std::ios::beg);
495 //    if ( ! gdcm_read_JPEGLS_file( fp,Raw ) )
496          return false;
497    }
498    else
499    {
500      // make sure this is the right JPEG compression
501      assert( !IsJPEGLS || !IsJPEG2000 );
502      // Precompute the offset localRaw will be shifted with
503      int length = XSize * YSize * SamplesPerPixel;
504      int numberBytes = BitsAllocated / 8;
505
506      JPEGInfo->DecompressFromFile(fp, Raw, BitsStored, numberBytes, length );
507      return true;
508    }
509 }
510
511 /**
512  * \brief Build Red/Green/Blue/Alpha LUT from File
513  *         when (0028,0004),Photometric Interpretation = [PALETTE COLOR ]
514  *          and (0028,1101),(0028,1102),(0028,1102)
515  *            - xxx Palette Color Lookup Table Descriptor - are found
516  *          and (0028,1201),(0028,1202),(0028,1202)
517  *            - xxx Palette Color Lookup Table Data - are found
518  * \warning does NOT deal with :
519  *   0028 1100 Gray Lookup Table Descriptor (Retired)
520  *   0028 1221 Segmented Red Palette Color Lookup Table Data
521  *   0028 1222 Segmented Green Palette Color Lookup Table Data
522  *   0028 1223 Segmented Blue Palette Color Lookup Table Data
523  *   no known Dicom reader deals with them :-(
524  * @return a RGBA Lookup Table
525  */
526 void PixelReadConvert::BuildLUTRGBA()
527 {
528    if ( LutRGBA )
529    {
530       return;
531    }
532    // Not so easy : see
533    // http://www.barre.nom.fr/medical/dicom2/limitations.html#Color%20Lookup%20Tables
534                                                                                 
535    if ( ! IsPaletteColor )
536    {
537       return;
538    }
539                                                                                 
540    if (   LutRedDescriptor   == GDCM_UNFOUND
541        || LutGreenDescriptor == GDCM_UNFOUND
542        || LutBlueDescriptor  == GDCM_UNFOUND )
543    {
544       gdcmWarningMacro( "(At least) a LUT Descriptor is missing" );
545       return;
546    }
547
548    ////////////////////////////////////////////
549    // Extract the info from the LUT descriptors
550    int lengthR;   // Red LUT length in Bytes
551    int debR;      // Subscript of the first Lut Value
552    int nbitsR;    // Lut item size (in Bits)
553    int nbRead;    // nb of items in LUT descriptor (must be = 3)
554
555    nbRead = sscanf( LutRedDescriptor.c_str(),
556                         "%d\\%d\\%d",
557                         &lengthR, &debR, &nbitsR );
558    if( nbRead != 3 )
559    {
560       gdcmWarningMacro( "Wrong Red LUT descriptor" );
561    }                                                                                
562    int lengthG;  // Green LUT length in Bytes
563    int debG;     // Subscript of the first Lut Value
564    int nbitsG;   // Lut item size (in Bits)
565
566    nbRead = sscanf( LutGreenDescriptor.c_str(),
567                     "%d\\%d\\%d",
568                     &lengthG, &debG, &nbitsG );  
569    if( nbRead != 3 )
570    {
571       gdcmWarningMacro( "Wrong Green LUT descriptor" );
572    }
573                                                                                 
574    int lengthB;  // Blue LUT length in Bytes
575    int debB;     // Subscript of the first Lut Value
576    int nbitsB;   // Lut item size (in Bits)
577    nbRead = sscanf( LutRedDescriptor.c_str(),
578                     "%d\\%d\\%d",
579                     &lengthB, &debB, &nbitsB );
580    if( nbRead != 3 )
581    {
582       gdcmWarningMacro( "Wrong Blue LUT descriptor" );
583    }
584  
585    gdcmWarningMacro(" lengthR " << lengthR << " debR " 
586                  << debR << " nbitsR " << nbitsR);
587    gdcmWarningMacro(" lengthG " << lengthG << " debG " 
588                  << debG << " nbitsG " << nbitsG);
589    gdcmWarningMacro(" lengthB " << lengthB << " debB " 
590                  << debB << " nbitsB " << nbitsB);
591
592    if ( !lengthR ) // if = 2^16, this shall be 0 see : CP-143
593       lengthR=65536;
594    if( !lengthG ) // if = 2^16, this shall be 0
595       lengthG=65536;
596    if ( !lengthB ) // if = 2^16, this shall be 0
597       lengthB=65536; 
598                                                                                 
599    ////////////////////////////////////////////////////////
600
601    if ( ( ! LutRedData ) || ( ! LutGreenData ) || ( ! LutBlueData ) )
602    {
603       gdcmWarningMacro( "(At least) a LUT is missing" );
604       return;
605    }
606
607    // -------------------------------------------------------------
608    
609    if ( BitsAllocated <= 8)
610    {
611       // forge the 4 * 8 Bits Red/Green/Blue/Alpha LUT
612       LutRGBA = new uint8_t[ 1024 ]; // 256 * 4 (R, G, B, Alpha)
613       if ( !LutRGBA )
614          return;
615       LutItemNumber = 256;
616       LutItemSize   = 8;
617       memset( LutRGBA, 0, 1024 );
618                                                                                 
619       int mult;
620       if ( ( nbitsR == 16 ) && ( BitsAllocated == 8 ) )
621       {
622          // when LUT item size is different than pixel size
623          mult = 2; // high byte must be = low byte
624       }
625       else
626       {
627          // See PS 3.3-2003 C.11.1.1.2 p 619
628          mult = 1;
629       }
630                                                                                 
631       // if we get a black image, let's just remove the '+1'
632       // from 'i*mult+1' and check again
633       // if it works, we shall have to check the 3 Palettes
634       // to see which byte is ==0 (first one, or second one)
635       // and fix the code
636       // We give up the checking to avoid some (useless ?) overhead
637       // (optimistic asumption)
638       int i;
639       uint8_t *a;
640
641       //take "Subscript of the first Lut Value" (debR,debG,debB) into account!
642
643       a = LutRGBA + 0 + debR;
644       for( i=0; i < lengthR; ++i )
645       {
646          *a = LutRedData[i*mult+1];
647          a += 4;
648       }
649                                                                                 
650       a = LutRGBA + 1 + debG;
651       for( i=0; i < lengthG; ++i)
652       {
653          *a = LutGreenData[i*mult+1];
654          a += 4;
655       }
656                                                                                 
657       a = LutRGBA + 2 + debB;
658       for(i=0; i < lengthB; ++i)
659       {
660          *a = LutBlueData[i*mult+1];
661          a += 4;
662       }
663                                                                                 
664       a = LutRGBA + 3 ;
665       for(i=0; i < 256; ++i)
666       {
667          *a = 1; // Alpha component
668          a += 4;
669       }
670    }
671    else
672    {
673       // Probabely the same stuff is to be done for 16 Bits Pixels
674       // with 65536 entries LUT ?!?
675       // Still looking for accurate info on the web :-(
676
677       gdcmWarningMacro( "Sorry Palette Color Lookup Tables not yet dealt with"
678                          << " for 16 Bits Per Pixel images" );
679
680       // forge the 4 * 16 Bits Red/Green/Blue/Alpha LUT
681
682       LutRGBA = (uint8_t *)new uint16_t[ 65536*4 ]; // 2^16 * 4 (R, G, B, Alpha)
683       if ( !LutRGBA )
684          return;
685       memset( LutRGBA, 0, 65536*4*2 );  // 16 bits = 2 bytes ;-)
686
687       LutItemNumber = 65536;
688       LutItemSize   = 16;
689
690       int i;
691       uint16_t *a16;
692
693       //take "Subscript of the first Lut Value" (debR,debG,debB) into account!
694
695       a16 = (uint16_t*)LutRGBA + 0 + debR;
696       for( i=0; i < lengthR; ++i )
697       {
698          *a16 = ((uint16_t*)LutRedData)[i];
699          a16 += 4;
700       }
701                                                                               
702       a16 = (uint16_t*)LutRGBA + 1 + debG;
703       for( i=0; i < lengthG; ++i)
704       {
705          *a16 = ((uint16_t*)LutGreenData)[i];
706          a16 += 4;
707       }
708                                                                                 
709       a16 = (uint16_t*)LutRGBA + 2 + debB;
710       for(i=0; i < lengthB; ++i)
711       {
712          *a16 = ((uint16_t*)LutBlueData)[i];
713          a16 += 4;
714       }
715                                                                              
716       a16 = (uint16_t*)LutRGBA + 3 ;
717       for(i=0; i < 65536; ++i)
718       {
719          *a16 = 1; // Alpha component
720          a16 += 4;
721       }
722 /* Just to 'see' the LUT, at debug time
723
724       a16=(uint16_t*)LutRGBA;
725       for (int j=0;j<65536;j++)
726       {
727          std::cout << *a16     << " " << *(a16+1) << " "
728                    << *(a16+2) << " " << *(a16+3) << std::endl;
729          a16+=4;
730       }
731 */
732    }
733 }
734
735 /**
736  * \brief Swap the bytes, according to \ref SwapCode.
737  */
738 void PixelReadConvert::ConvertSwapZone()
739 {
740    unsigned int i;
741
742    if( BitsAllocated == 16 )
743    {
744       uint16_t *im16 = (uint16_t*)Raw;
745       switch( SwapCode )
746       {
747          case 1234:
748             break;
749          case 3412:
750          case 2143:
751          case 4321:
752             for( i = 0; i < RawSize / 2; i++ )
753             {
754                im16[i]= (im16[i] >> 8) | (im16[i] << 8 );
755             }
756             break;
757          default:
758             gdcmWarningMacro("SwapCode value (16 bits) not allowed.");
759       }
760    }
761    else if( BitsAllocated == 32 )
762    {
763       uint32_t s32;
764       uint16_t high;
765       uint16_t low;
766       uint32_t *im32 = (uint32_t*)Raw;
767       switch ( SwapCode )
768       {
769          case 1234:
770             break;
771          case 4321:
772             for( i = 0; i < RawSize / 4; i++ )
773             {
774                low     = im32[i] & 0x0000ffff;  // 4321
775                high    = im32[i] >> 16;
776                high    = ( high >> 8 ) | ( high << 8 );
777                low     = ( low  >> 8 ) | ( low  << 8 );
778                s32     = low;
779                im32[i] = ( s32 << 16 ) | high;
780             }
781             break;
782          case 2143:
783             for( i = 0; i < RawSize / 4; i++ )
784             {
785                low     = im32[i] & 0x0000ffff;   // 2143
786                high    = im32[i] >> 16;
787                high    = ( high >> 8 ) | ( high << 8 );
788                low     = ( low  >> 8 ) | ( low  << 8 );
789                s32     = high;
790                im32[i] = ( s32 << 16 ) | low;
791             }
792             break;
793          case 3412:
794             for( i = 0; i < RawSize / 4; i++ )
795             {
796                low     = im32[i] & 0x0000ffff; // 3412
797                high    = im32[i] >> 16;
798                s32     = low;
799                im32[i] = ( s32 << 16 ) | high;
800             }
801             break;
802          default:
803             gdcmWarningMacro("SwapCode value (32 bits) not allowed." );
804       }
805    }
806 }
807
808 /**
809  * \brief Deal with endianness i.e. re-arange bytes inside the integer
810  */
811 void PixelReadConvert::ConvertReorderEndianity()
812 {
813    if ( BitsAllocated != 8 )
814    {
815       ConvertSwapZone();
816    }
817
818    // Special kludge in order to deal with xmedcon broken images:
819    if ( BitsAllocated == 16
820      && BitsStored < BitsAllocated
821      && !PixelSign )
822    {
823       int l = (int)( RawSize / ( BitsAllocated / 8 ) );
824       uint16_t *deb = (uint16_t *)Raw;
825       for(int i = 0; i<l; i++)
826       {
827          if( *deb == 0xffff )
828          {
829            *deb = 0;
830          }
831          deb++;
832       }
833    }
834 }
835
836 /**
837  * \brief Deal with Grey levels i.e. re-arange them
838  *        to have low values = dark, high values = bright
839  */
840 void PixelReadConvert::ConvertFixGreyLevels()
841 {
842    if (!IsMonochrome1)
843       return;
844
845    uint32_t i; // to please M$VC6
846    int16_t j;
847
848    if (!PixelSign)
849    {
850       if ( BitsAllocated == 8 )
851       {
852          uint8_t *deb = (uint8_t *)Raw;
853          for (i=0; i<RawSize; i++)      
854          {
855             *deb = 255 - *deb;
856             deb++;
857          }
858          return;
859       }
860
861       if ( BitsAllocated == 16 )
862       {
863          uint16_t mask =1;
864          for (j=0; j<BitsStored-1; j++)
865          {
866             mask = (mask << 1) +1; // will be fff when BitsStored=12
867          }
868
869          uint16_t *deb = (uint16_t *)Raw;
870          for (i=0; i<RawSize/2; i++)      
871          {
872             *deb = mask - *deb;
873             deb++;
874          }
875          return;
876        }
877    }
878    else
879    {
880       if ( BitsAllocated == 8 )
881       {
882          uint8_t smask8 = 255;
883          uint8_t *deb = (uint8_t *)Raw;
884          for (i=0; i<RawSize; i++)      
885          {
886             *deb = smask8 - *deb;
887             deb++;
888          }
889          return;
890       }
891       if ( BitsAllocated == 16 )
892       {
893          uint16_t smask16 = 65535;
894          uint16_t *deb = (uint16_t *)Raw;
895          for (i=0; i<RawSize/2; i++)      
896          {
897             *deb = smask16 - *deb;
898             deb++;
899          }
900          return;
901       }
902    }
903 }
904
905 /**
906  * \brief  Re-arrange the bits within the bytes.
907  * @return Boolean always true
908  */
909 bool PixelReadConvert::ConvertReArrangeBits() throw ( FormatError )
910 {
911    if ( BitsStored != BitsAllocated )
912    {
913       int l = (int)( RawSize / ( BitsAllocated / 8 ) );
914       if ( BitsAllocated == 16 )
915       {
916          uint16_t mask = 0xffff;
917          mask = mask >> ( BitsAllocated - BitsStored );
918          uint16_t *deb = (uint16_t*)Raw;
919          for(int i = 0; i<l; i++)
920          {
921             *deb = (*deb >> (BitsStored - HighBitPosition - 1)) & mask;
922             deb++;
923          }
924       }
925       else if ( BitsAllocated == 32 )
926       {
927          uint32_t mask = 0xffffffff;
928          mask = mask >> ( BitsAllocated - BitsStored );
929          uint32_t *deb = (uint32_t*)Raw;
930          for(int i = 0; i<l; i++)
931          {
932             *deb = (*deb >> (BitsStored - HighBitPosition - 1)) & mask;
933             deb++;
934          }
935       }
936       else
937       {
938          gdcmWarningMacro("Weird image");
939          throw FormatError( "Weird image !?" );
940       }
941    }
942    return true;
943 }
944
945 /**
946  * \brief   Convert (Red plane, Green plane, Blue plane) to RGB pixels
947  * \warning Works on all the frames at a time
948  */
949 void PixelReadConvert::ConvertRGBPlanesToRGBPixels()
950 {
951    gdcmWarningMacro("ConvertRGBPlanesToRGBPixels");
952
953    uint8_t *localRaw = Raw;
954    uint8_t *copyRaw = new uint8_t[ RawSize ];
955    memmove( copyRaw, localRaw, RawSize );
956
957    int l = XSize * YSize * ZSize;
958
959    uint8_t *a = copyRaw;
960    uint8_t *b = copyRaw + l;
961    uint8_t *c = copyRaw + l + l;
962
963    for (int j = 0; j < l; j++)
964    {
965       *(localRaw++) = *(a++);
966       *(localRaw++) = *(b++);
967       *(localRaw++) = *(c++);
968    }
969    delete[] copyRaw;
970 }
971
972 /**
973  * \brief   Convert (cY plane, cB plane, cR plane) to RGB pixels
974  * \warning Works on all the frames at a time
975  */
976 void PixelReadConvert::ConvertYcBcRPlanesToRGBPixels()
977 {
978   // Remarks for YBR newbees :
979   // YBR_FULL works very much like RGB, i.e. three samples per pixel, 
980   // just the color space is YCbCr instead of RGB. This is particularly useful
981   // for doppler ultrasound where most of the image is grayscale 
982   // (i.e. only populates the Y components) and Cb and Cr are mostly zero,
983   // except for the few patches of color on the image.
984   // On such images, RLE achieves a compression ratio that is much better 
985   // than the compression ratio on an equivalent RGB image. 
986     
987    uint8_t *localRaw = Raw;
988    uint8_t *copyRaw = new uint8_t[ RawSize ];
989    memmove( copyRaw, localRaw, RawSize );
990
991    // to see the tricks about YBR_FULL, YBR_FULL_422,
992    // YBR_PARTIAL_422, YBR_ICT, YBR_RCT have a look at :
993    // ftp://medical.nema.org/medical/dicom/final/sup61_ft.pdf
994    // and be *very* affraid
995    //
996    int l        = XSize * YSize;
997    int nbFrames = ZSize;
998
999    uint8_t *a = copyRaw + 0;
1000    uint8_t *b = copyRaw + l;
1001    uint8_t *c = copyRaw + l+ l;
1002    int32_t R, G, B;
1003
1004    ///  We replaced easy to understand but time consuming floating point
1005    ///  computations by the 'well known' integer computation counterpart
1006    ///  Refer to :
1007    ///            http://lestourtereaux.free.fr/papers/data/yuvrgb.pdf
1008    ///  for code optimisation.
1009
1010    gdcmWarningMacro("ConvertYcBcRPlanesToRGBPixels");
1011
1012    for ( int i = 0; i < nbFrames; i++ )
1013    {
1014       for ( int j = 0; j < l; j++ )
1015       {
1016          R = 38142 *(*a-16) + 52298 *(*c -128);
1017          G = 38142 *(*a-16) - 26640 *(*c -128) - 12845 *(*b -128);
1018          B = 38142 *(*a-16) + 66093 *(*b -128);
1019
1020          R = (R+16384)>>15;
1021          G = (G+16384)>>15;
1022          B = (B+16384)>>15;
1023
1024          if (R < 0)   R = 0;
1025          if (G < 0)   G = 0;
1026          if (B < 0)   B = 0;
1027          if (R > 255) R = 255;
1028          if (G > 255) G = 255;
1029          if (B > 255) B = 255;
1030
1031          *(localRaw++) = (uint8_t)R;
1032          *(localRaw++) = (uint8_t)G;
1033          *(localRaw++) = (uint8_t)B;
1034          a++;
1035          b++;
1036          c++;
1037       }
1038    }
1039    delete[] copyRaw;
1040 }
1041
1042 /// \brief Deals with the color decoding i.e. handle:
1043 ///   - R, G, B planes (as opposed to RGB pixels)
1044 ///   - YBR (various) encodings.
1045 ///   - LUT[s] (or "PALETTE COLOR").
1046
1047 void PixelReadConvert::ConvertHandleColor()
1048 {
1049    //////////////////////////////////
1050    // Deal with the color decoding i.e. handle:
1051    //   - R, G, B planes (as opposed to RGB pixels)
1052    //   - YBR (various) encodings.
1053    //   - LUT[s] (or "PALETTE COLOR").
1054    //
1055    // The classification in the color decoding schema is based on the blending
1056    // of two Dicom tags values:
1057    // * "Photometric Interpretation" for which we have the cases:
1058    //  - [Photo A] MONOCHROME[1|2] pictures,
1059    //  - [Photo B] RGB or YBR_FULL_422 (which acts as RGB),
1060    //  - [Photo C] YBR_* (with the above exception of YBR_FULL_422)
1061    //  - [Photo D] "PALETTE COLOR" which indicates the presence of LUT[s].
1062    // * "Planar Configuration" for which we have the cases:
1063    //  - [Planar 0] 0 then Pixels are already RGB
1064    //  - [Planar 1] 1 then we have 3 planes : R, G, B,
1065    //  - [Planar 2] 2 then we have 1 gray Plane and 3 LUTs
1066    //
1067    // Now in theory, one could expect some coherence when blending the above
1068    // cases. For example we should not encounter files belonging at the
1069    // time to case [Planar 0] and case [Photo D].
1070    // Alas, this was only theory ! Because in practice some odd (read ill
1071    // formated Dicom) files (e.g. gdcmData/US-PAL-8-10x-echo.dcm) we encounter:
1072    //     - "Planar Configuration" = 0,
1073    //     - "Photometric Interpretation" = "PALETTE COLOR".
1074    // Hence gdcm will use the folowing "heuristic" in order to be tolerant
1075    // towards Dicom-non-conformant files:
1076    //   << whatever the "Planar Configuration" value might be, a
1077    //      "Photometric Interpretation" set to "PALETTE COLOR" forces
1078    //      a LUT intervention >>
1079    //
1080    // Now we are left with the following handling of the cases:
1081    // - [Planar 0] OR  [Photo A] no color decoding (since respectively
1082    //       Pixels are already RGB and monochrome pictures have no color :),
1083    // - [Planar 1] AND [Photo B] handled with ConvertRGBPlanesToRGBPixels()
1084    // - [Planar 1] AND [Photo C] handled with ConvertYcBcRPlanesToRGBPixels()
1085    // - [Planar 2] OR  [Photo D] requires LUT intervention.
1086
1087    gdcmWarningMacro("ConvertHandleColor");
1088
1089    if ( ! IsRawRGB() )
1090    {
1091       // [Planar 2] OR  [Photo D]: LUT intervention done outside
1092       return;
1093    }
1094                                                                                 
1095    if ( PlanarConfiguration == 1 )
1096    {
1097       if ( IsYBRFull )
1098       {
1099          // [Planar 1] AND [Photo C] (remember YBR_FULL_422 acts as RGB)
1100          ConvertYcBcRPlanesToRGBPixels();
1101       }
1102       else
1103       {
1104          // [Planar 1] AND [Photo C]
1105          ConvertRGBPlanesToRGBPixels();
1106       }
1107       return;
1108    }
1109                                                                                 
1110    // When planarConf is 0, and RLELossless (forbidden by Dicom norm)
1111    // pixels need to be RGB-fyied anyway
1112
1113    if (IsRLELossless)
1114    { 
1115      ConvertRGBPlanesToRGBPixels();
1116    }
1117
1118    // In *normal *case, when planarConf is 0, pixels are already in RGB
1119 }
1120
1121 /// Computes the Pixels Size
1122 void PixelReadConvert::ComputeRawAndRGBSizes()
1123 {
1124    int bitsAllocated = BitsAllocated;
1125    // Number of "Bits Allocated" is fixed to 16 when it's 12, since
1126    // in this case we will expand the image to 16 bits (see
1127    //    \ref ReadAndDecompress12BitsTo16Bits() )
1128    if (  BitsAllocated == 12 )
1129    {
1130       bitsAllocated = 16;
1131    }
1132                                                                                 
1133    RawSize =  XSize * YSize * ZSize
1134                      * ( bitsAllocated / 8 )
1135                      * SamplesPerPixel;
1136    if ( HasLUT )
1137    {
1138       RGBSize = 3 * RawSize; // works for 8 and 16 bits per Pixel
1139    }
1140    else
1141    {
1142       RGBSize = RawSize;
1143    }
1144 }
1145
1146 /// Allocates room for RGB Pixels
1147 void PixelReadConvert::AllocateRGB()
1148 {
1149   if ( RGB )
1150      delete [] RGB;
1151   RGB = new uint8_t[RGBSize];
1152 }
1153
1154 /// Allocates room for RAW Pixels
1155 void PixelReadConvert::AllocateRaw()
1156 {
1157   if ( Raw )
1158      delete [] Raw;
1159   Raw = new uint8_t[RawSize];
1160 }
1161
1162 //-----------------------------------------------------------------------------
1163 // Print
1164 /**
1165  * \brief        Print self.
1166  * @param indent Indentation string to be prepended during printing.
1167  * @param os     Stream to print to.
1168  */
1169 void PixelReadConvert::Print( std::ostream &os, std::string const &indent )
1170 {
1171    os << indent
1172       << "--- Pixel information -------------------------"
1173       << std::endl;
1174    os << indent
1175       << "Pixel Data: offset " << PixelOffset
1176       << " x(" << std::hex << PixelOffset << std::dec
1177       << ")   length " << PixelDataLength
1178       << " x(" << std::hex << PixelDataLength << std::dec
1179       << ")" << std::endl;
1180
1181    if ( IsRLELossless )
1182    {
1183       if ( RLEInfo )
1184       {
1185          RLEInfo->Print( os, indent );
1186       }
1187       else
1188       {
1189          gdcmWarningMacro("Set as RLE file but NO RLEinfo present.");
1190       }
1191    }
1192
1193    if ( IsJPEG2000 || IsJPEGLossless || IsJPEGLossy || IsJPEGLS )
1194    {
1195       if ( JPEGInfo )
1196       {
1197          JPEGInfo->Print( os, indent );
1198       }
1199       else
1200       {
1201          gdcmWarningMacro("Set as JPEG file but NO JPEGinfo present.");
1202       }
1203    }
1204 }
1205
1206 //-----------------------------------------------------------------------------
1207 } // end namespace gdcm
1208
1209 // NOTES on File internal calls
1210 // User
1211 // ---> GetImageData
1212 //     ---> GetImageDataIntoVector
1213 //        |---> GetImageDataIntoVectorRaw
1214 //        | lut intervention
1215 // User
1216 // ---> GetImageDataRaw
1217 //     ---> GetImageDataIntoVectorRaw
1218