]> Creatis software - gdcm.git/blob - src/gdcmPixelReadConvert.cxx
Some constructors (at least one ...) use BitsAllocated =! BitsStored for
[gdcm.git] / src / gdcmPixelReadConvert.cxx
1 /*=========================================================================
2
3   Program:   gdcm
4   Module:    $RCSfile: gdcmPixelReadConvert.cxx,v $
5   Language:  C++
6   Date:      $Date: 2005/06/22 07:52:32 $
7   Version:   $Revision: 1.69 $
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    if ( BitsStored != BitsAllocated )
918    {
919       int l = (int)( RawSize / ( BitsAllocated / 8 ) );
920       if ( BitsAllocated == 16 
921       {
922          // pmask : to mask the 'unused bits' (may contain overlays)
923          uint16_t pmask = 0xffff;
924          pmask = pmask >> ( BitsAllocated - BitsStored );
925
926          uint16_t *deb = (uint16_t*)Raw;
927
928          if ( !PixelSign )
929          {
930             for(int i = 0; i<l; i++)
931             {             
932                *deb = (*deb >> (BitsStored - HighBitPosition - 1)) & pmask;
933                deb++;
934             }
935          }
936          else
937          {
938             // smask : to check the 'sign' when BitsStored != BitsAllocated
939             uint16_t smask = 0x8000;
940             smask = smask >> ( BitsAllocated - BitsStored );
941             // nmask : to propagate sign bit on negative values
942             int16_t nmask = 0x8000;  
943             nmask = nmask >> ( BitsAllocated - BitsStored );
944
945             for(int i = 0; i<l; i++)
946             {
947                *deb = *deb >> (BitsStored - HighBitPosition - 1);
948                if ( *deb & smask )
949                   *deb = *deb | nmask;
950                else
951                   *deb = *deb & pmask;
952                deb++;
953             }
954          }
955       }
956       else if ( BitsAllocated == 32 )
957       {
958          // pmask : to mask the 'unused bits' (may contain overlays)
959          uint32_t pmask = 0xffffffff;
960          pmask = pmask >> ( BitsAllocated - BitsStored );
961
962          uint32_t *deb = (uint32_t*)Raw;
963
964          if ( !PixelSign )
965          {
966             for(int i = 0; i<l; i++)
967             {             
968                *deb = (*deb >> (BitsStored - HighBitPosition - 1)) & pmask;
969                deb++;
970             }
971          }
972          else
973          {
974             // smask : to check the 'sign' when BitsStored != BitsAllocated
975             uint32_t smask = 0x80000000;
976             smask = smask >> ( BitsAllocated - BitsStored );
977             // nmask : to propagate sign bit on negative values
978             int32_t nmask = 0x80000000;   
979             nmask = nmask >> ( BitsAllocated - BitsStored )
980
981             for(int i = 0; i<l; i++)
982             {
983                *deb = *deb >> (BitsStored - HighBitPosition - 1);
984                if ( *deb & smask )
985                   *deb = *deb | nmask;
986                else
987                   *deb = *deb & pmask;
988                deb++;
989             }
990          }
991       }
992       else
993       {
994          gdcmWarningMacro("Weird image (BitsAllocated !=8, 12, 16, 32)");
995          throw FormatError( "Weird image !?" );
996       }
997    }
998    return true;
999 }
1000
1001 /**
1002  * \brief   Convert (Red plane, Green plane, Blue plane) to RGB pixels
1003  * \warning Works on all the frames at a time
1004  */
1005 void PixelReadConvert::ConvertRGBPlanesToRGBPixels()
1006 {
1007    gdcmWarningMacro("--> ConvertRGBPlanesToRGBPixels");
1008
1009    uint8_t *localRaw = Raw;
1010    uint8_t *copyRaw = new uint8_t[ RawSize ];
1011    memmove( copyRaw, localRaw, RawSize );
1012
1013    int l = XSize * YSize * ZSize;
1014
1015    uint8_t *a = copyRaw;
1016    uint8_t *b = copyRaw + l;
1017    uint8_t *c = copyRaw + l + l;
1018
1019    for (int j = 0; j < l; j++)
1020    {
1021       *(localRaw++) = *(a++);
1022       *(localRaw++) = *(b++);
1023       *(localRaw++) = *(c++);
1024    }
1025    delete[] copyRaw;
1026 }
1027
1028 /**
1029  * \brief   Convert (cY plane, cB plane, cR plane) to RGB pixels
1030  * \warning Works on all the frames at a time
1031  */
1032 void PixelReadConvert::ConvertYcBcRPlanesToRGBPixels()
1033 {
1034   // Remarks for YBR newbees :
1035   // YBR_FULL works very much like RGB, i.e. three samples per pixel, 
1036   // just the color space is YCbCr instead of RGB. This is particularly useful
1037   // for doppler ultrasound where most of the image is grayscale 
1038   // (i.e. only populates the Y components) and Cb and Cr are mostly zero,
1039   // except for the few patches of color on the image.
1040   // On such images, RLE achieves a compression ratio that is much better 
1041   // than the compression ratio on an equivalent RGB image. 
1042  
1043    gdcmWarningMacro("--> ConvertYcBcRPlanesToRGBPixels");
1044    
1045    uint8_t *localRaw = Raw;
1046    uint8_t *copyRaw = new uint8_t[ RawSize ];
1047    memmove( copyRaw, localRaw, RawSize );
1048
1049    // to see the tricks about YBR_FULL, YBR_FULL_422,
1050    // YBR_PARTIAL_422, YBR_ICT, YBR_RCT have a look at :
1051    // ftp://medical.nema.org/medical/dicom/final/sup61_ft.pdf
1052    // and be *very* affraid
1053    //
1054    int l        = XSize * YSize;
1055    int nbFrames = ZSize;
1056
1057    uint8_t *a = copyRaw + 0;
1058    uint8_t *b = copyRaw + l;
1059    uint8_t *c = copyRaw + l+ l;
1060    int32_t R, G, B;
1061
1062    ///  We replaced easy to understand but time consuming floating point
1063    ///  computations by the 'well known' integer computation counterpart
1064    ///  Refer to :
1065    ///            http://lestourtereaux.free.fr/papers/data/yuvrgb.pdf
1066    ///  for code optimisation.
1067
1068    for ( int i = 0; i < nbFrames; i++ )
1069    {
1070       for ( int j = 0; j < l; j++ )
1071       {
1072          R = 38142 *(*a-16) + 52298 *(*c -128);
1073          G = 38142 *(*a-16) - 26640 *(*c -128) - 12845 *(*b -128);
1074          B = 38142 *(*a-16) + 66093 *(*b -128);
1075
1076          R = (R+16384)>>15;
1077          G = (G+16384)>>15;
1078          B = (B+16384)>>15;
1079
1080          if (R < 0)   R = 0;
1081          if (G < 0)   G = 0;
1082          if (B < 0)   B = 0;
1083          if (R > 255) R = 255;
1084          if (G > 255) G = 255;
1085          if (B > 255) B = 255;
1086
1087          *(localRaw++) = (uint8_t)R;
1088          *(localRaw++) = (uint8_t)G;
1089          *(localRaw++) = (uint8_t)B;
1090          a++;
1091          b++;
1092          c++;
1093       }
1094    }
1095    delete[] copyRaw;
1096 }
1097
1098 /// \brief Deals with the color decoding i.e. handle:
1099 ///   - R, G, B planes (as opposed to RGB pixels)
1100 ///   - YBR (various) encodings.
1101 ///   - LUT[s] (or "PALETTE COLOR").
1102
1103 void PixelReadConvert::ConvertHandleColor()
1104 {
1105    //////////////////////////////////
1106    // Deal with the color decoding i.e. handle:
1107    //   - R, G, B planes (as opposed to RGB pixels)
1108    //   - YBR (various) encodings.
1109    //   - LUT[s] (or "PALETTE COLOR").
1110    //
1111    // The classification in the color decoding schema is based on the blending
1112    // of two Dicom tags values:
1113    // * "Photometric Interpretation" for which we have the cases:
1114    //  - [Photo A] MONOCHROME[1|2] pictures,
1115    //  - [Photo B] RGB or YBR_FULL_422 (which acts as RGB),
1116    //  - [Photo C] YBR_* (with the above exception of YBR_FULL_422)
1117    //  - [Photo D] "PALETTE COLOR" which indicates the presence of LUT[s].
1118    // * "Planar Configuration" for which we have the cases:
1119    //  - [Planar 0] 0 then Pixels are already RGB
1120    //  - [Planar 1] 1 then we have 3 planes : R, G, B,
1121    //  - [Planar 2] 2 then we have 1 gray Plane and 3 LUTs
1122    //
1123    // Now in theory, one could expect some coherence when blending the above
1124    // cases. For example we should not encounter files belonging at the
1125    // time to case [Planar 0] and case [Photo D].
1126    // Alas, this was only theory ! Because in practice some odd (read ill
1127    // formated Dicom) files (e.g. gdcmData/US-PAL-8-10x-echo.dcm) we encounter:
1128    //     - "Planar Configuration" = 0,
1129    //     - "Photometric Interpretation" = "PALETTE COLOR".
1130    // Hence gdcm will use the folowing "heuristic" in order to be tolerant
1131    // towards Dicom-non-conformant files:
1132    //   << whatever the "Planar Configuration" value might be, a
1133    //      "Photometric Interpretation" set to "PALETTE COLOR" forces
1134    //      a LUT intervention >>
1135    //
1136    // Now we are left with the following handling of the cases:
1137    // - [Planar 0] OR  [Photo A] no color decoding (since respectively
1138    //       Pixels are already RGB and monochrome pictures have no color :),
1139    // - [Planar 1] AND [Photo B] handled with ConvertRGBPlanesToRGBPixels()
1140    // - [Planar 1] AND [Photo C] handled with ConvertYcBcRPlanesToRGBPixels()
1141    // - [Planar 2] OR  [Photo D] requires LUT intervention.
1142
1143    gdcmWarningMacro("--> ConvertHandleColor"
1144                     << "Planar Configuration " << PlanarConfiguration );
1145
1146    if ( ! IsRawRGB() )
1147    {
1148       // [Planar 2] OR  [Photo D]: LUT intervention done outside
1149       gdcmWarningMacro("--> RawRGB : LUT intervention done outside");
1150       return;
1151    }
1152                                                                                 
1153    if ( PlanarConfiguration == 1 )
1154    {
1155       if ( IsYBRFull )
1156       {
1157          // [Planar 1] AND [Photo C] (remember YBR_FULL_422 acts as RGB)
1158          gdcmWarningMacro("--> YBRFull");
1159          ConvertYcBcRPlanesToRGBPixels();
1160       }
1161       else
1162       {
1163          // [Planar 1] AND [Photo C]
1164          gdcmWarningMacro("--> YBRFull");
1165          ConvertRGBPlanesToRGBPixels();
1166       }
1167       return;
1168    }
1169                                                                                 
1170    // When planarConf is 0, and RLELossless (forbidden by Dicom norm)
1171    // pixels need to be RGB-fyied anyway
1172
1173    if (IsRLELossless)
1174    { 
1175      gdcmWarningMacro("--> RLE Lossless");
1176      ConvertRGBPlanesToRGBPixels();
1177    }
1178
1179    // In *normal *case, when planarConf is 0, pixels are already in RGB
1180 }
1181
1182 /// Computes the Pixels Size
1183 void PixelReadConvert::ComputeRawAndRGBSizes()
1184 {
1185    int bitsAllocated = BitsAllocated;
1186    // Number of "Bits Allocated" is fixed to 16 when it's 12, since
1187    // in this case we will expand the image to 16 bits (see
1188    //    \ref ReadAndDecompress12BitsTo16Bits() )
1189    if (  BitsAllocated == 12 )
1190    {
1191       bitsAllocated = 16;
1192    }
1193                                                                                 
1194    RawSize =  XSize * YSize * ZSize
1195                      * ( bitsAllocated / 8 )
1196                      * SamplesPerPixel;
1197    if ( HasLUT )
1198    {
1199       RGBSize = 3 * RawSize; // works for 8 and 16 bits per Pixel
1200    }
1201    else
1202    {
1203       RGBSize = RawSize;
1204    }
1205 }
1206
1207 /// Allocates room for RGB Pixels
1208 void PixelReadConvert::AllocateRGB()
1209 {
1210   if ( RGB )
1211      delete [] RGB;
1212   RGB = new uint8_t[RGBSize];
1213 }
1214
1215 /// Allocates room for RAW Pixels
1216 void PixelReadConvert::AllocateRaw()
1217 {
1218   if ( Raw )
1219      delete [] Raw;
1220   Raw = new uint8_t[RawSize];
1221 }
1222
1223 //-----------------------------------------------------------------------------
1224 // Print
1225 /**
1226  * \brief        Print self.
1227  * @param indent Indentation string to be prepended during printing.
1228  * @param os     Stream to print to.
1229  */
1230 void PixelReadConvert::Print( std::ostream &os, std::string const &indent )
1231 {
1232    os << indent
1233       << "--- Pixel information -------------------------"
1234       << std::endl;
1235    os << indent
1236       << "Pixel Data: offset " << PixelOffset
1237       << " x(" << std::hex << PixelOffset << std::dec
1238       << ")   length " << PixelDataLength
1239       << " x(" << std::hex << PixelDataLength << std::dec
1240       << ")" << std::endl;
1241
1242    if ( IsRLELossless )
1243    {
1244       if ( RLEInfo )
1245       {
1246          RLEInfo->Print( os, indent );
1247       }
1248       else
1249       {
1250          gdcmWarningMacro("Set as RLE file but NO RLEinfo present.");
1251       }
1252    }
1253
1254    if ( IsJPEG2000 || IsJPEGLossless || IsJPEGLossy || IsJPEGLS )
1255    {
1256       if ( JPEGInfo )
1257       {
1258          JPEGInfo->Print( os, indent );
1259       }
1260       else
1261       {
1262          gdcmWarningMacro("Set as JPEG file but NO JPEGinfo present.");
1263       }
1264    }
1265 }
1266
1267 //-----------------------------------------------------------------------------
1268 } // end namespace gdcm
1269
1270 // NOTES on File internal calls
1271 // User
1272 // ---> GetImageData
1273 //     ---> GetImageDataIntoVector
1274 //        |---> GetImageDataIntoVectorRaw
1275 //        | lut intervention
1276 // User
1277 // ---> GetImageDataRaw
1278 //     ---> GetImageDataIntoVectorRaw
1279