]> Creatis software - gdcm.git/blob - gdcmPixelReadConvert.cxx
3e4a1e346c8b0db4cc26618b18ec9bc0e4759363
[gdcm.git] / gdcmPixelReadConvert.cxx
1 /*=========================================================================
2  
3   Program:   gdcm
4   Module:    $RCSfile: gdcmPixelReadConvert.cxx,v $
5   Language:  C++
6   Date:      $Date: 2005/09/07 14:12:23 $
7   Version:   $Revision: 1.77 $
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    RLEInfo      = 0;
54    JPEGInfo     = 0;
55    UserFunction = 0;
56    FileInternal = 0;
57 }
58
59 /// Canonical Destructor
60 PixelReadConvert::~PixelReadConvert() 
61 {
62    Squeeze();
63 }
64
65 //-----------------------------------------------------------------------------
66 // Public
67 /**
68  * \brief Predicate to know whether the image[s] (once Raw) is RGB.
69  * \note See comments of \ref ConvertHandleColor
70  */
71 bool PixelReadConvert::IsRawRGB()
72 {
73    if (   IsMonochrome
74        || PlanarConfiguration == 2
75        || IsPaletteColor )
76    {
77       return false;
78    }
79    return true;
80 }
81 /**
82  * \brief Gets various usefull informations from the file header
83  * @param file gdcm::File pointer
84  */
85 void PixelReadConvert::GrabInformationsFromFile( File *file )
86 {
87    // Number of Bits Allocated for storing a Pixel is defaulted to 16
88    // when absent from the file.
89    BitsAllocated = file->GetBitsAllocated();
90    if ( BitsAllocated == 0 )
91    {
92       BitsAllocated = 16;
93    }
94
95    // Number of "Bits Stored", defaulted to number of "Bits Allocated"
96    // when absent from the file.
97    BitsStored = file->GetBitsStored();
98    if ( BitsStored == 0 )
99    {
100       BitsStored = BitsAllocated;
101    }
102
103    // High Bit Position, defaulted to "Bits Allocated" - 1
104    HighBitPosition = file->GetHighBitPosition();
105    if ( HighBitPosition == 0 )
106    {
107       HighBitPosition = BitsAllocated - 1;
108    }
109
110    XSize           = file->GetXSize();
111    YSize           = file->GetYSize();
112    ZSize           = file->GetZSize();
113    SamplesPerPixel = file->GetSamplesPerPixel();
114    //PixelSize       = file->GetPixelSize();  Useless
115    PixelSign       = file->IsSignedPixelData();
116    SwapCode        = file->GetSwapCode();
117    std::string ts  = file->GetTransferSyntax();
118    IsRaw =
119         ( ! file->IsDicomV3() )
120      || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ImplicitVRLittleEndian
121      || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ImplicitVRLittleEndianDLXGE
122      || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ExplicitVRLittleEndian
123      || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ExplicitVRBigEndian
124      || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::DeflatedExplicitVRLittleEndian;
125
126    IsMPEG          = Global::GetTS()->IsMPEG(ts);
127    IsJPEG2000      = Global::GetTS()->IsJPEG2000(ts);
128    IsJPEGLS        = Global::GetTS()->IsJPEGLS(ts);
129    IsJPEGLossy     = Global::GetTS()->IsJPEGLossy(ts);
130    IsJPEGLossless  = Global::GetTS()->IsJPEGLossless(ts);
131    IsRLELossless   = Global::GetTS()->IsRLELossless(ts);
132
133    PixelOffset     = file->GetPixelOffset();
134    PixelDataLength = file->GetPixelAreaLength();
135    RLEInfo         = file->GetRLEInfo();
136    JPEGInfo        = file->GetJPEGInfo();
137
138    IsMonochrome    = file->IsMonochrome();
139    IsMonochrome1   = file->IsMonochrome1();
140    IsPaletteColor  = file->IsPaletteColor();
141    IsYBRFull       = file->IsYBRFull();
142
143    PlanarConfiguration = file->GetPlanarConfiguration();
144
145    /////////////////////////////////////////////////////////////////
146    // LUT section:
147    HasLUT = file->HasLUT();
148    if ( HasLUT )
149    {
150       // Just in case some access to a File element requires disk access.
151       LutRedDescriptor   = file->GetEntryValue( 0x0028, 0x1101 );
152       LutGreenDescriptor = file->GetEntryValue( 0x0028, 0x1102 );
153       LutBlueDescriptor  = file->GetEntryValue( 0x0028, 0x1103 );
154    
155       // Depending on the value of Document::MAX_SIZE_LOAD_ELEMENT_VALUE
156       // [ refer to invocation of Document::SetMaxSizeLoadEntry() in
157       // Document::Document() ], the loading of the value (content) of a
158       // [Bin|Val]Entry occurence migth have been hindered (read simply NOT
159       // loaded). Hence, we first try to obtain the LUTs data from the file
160       // and when this fails we read the LUTs data directly from disk.
161       // \TODO Reading a [Bin|Val]Entry directly from disk is a kludge.
162       //       We should NOT bypass the [Bin|Val]Entry class. Instead
163       //       an access to an UNLOADED content of a [Bin|Val]Entry occurence
164       //       (e.g. BinEntry::GetBinArea()) should force disk access from
165       //       within the [Bin|Val]Entry class itself. The only problem
166       //       is that the [Bin|Val]Entry is unaware of the FILE* is was
167       //       parsed from. Fix that. FIXME.
168    
169       // //// Red round
170       file->LoadEntryBinArea(0x0028, 0x1201);
171       LutRedData = (uint8_t*)file->GetEntryBinArea( 0x0028, 0x1201 );
172       if ( ! LutRedData )
173       {
174          gdcmWarningMacro( "Unable to read Red Palette Color Lookup Table data" );
175       }
176
177       // //// Green round:
178       file->LoadEntryBinArea(0x0028, 0x1202);
179       LutGreenData = (uint8_t*)file->GetEntryBinArea(0x0028, 0x1202 );
180       if ( ! LutGreenData)
181       {
182          gdcmWarningMacro( "Unable to read Green Palette Color Lookup Table data" );
183       }
184
185       // //// Blue round:
186       file->LoadEntryBinArea(0x0028, 0x1203);
187       LutBlueData = (uint8_t*)file->GetEntryBinArea( 0x0028, 0x1203 );
188       if ( ! LutBlueData )
189       {
190          gdcmWarningMacro( "Unable to read Blue Palette Color Lookup Table data" );
191       }
192    }
193    FileInternal = file;   
194
195    ComputeRawAndRGBSizes();
196 }
197
198 /// \brief Reads from disk and decompresses Pixels
199 bool PixelReadConvert::ReadAndDecompressPixelData( std::ifstream *fp )
200 {
201    // ComputeRawAndRGBSizes is already made by 
202    // ::GrabInformationsFromfile. So, the structure sizes are
203    // correct
204    Squeeze();
205
206    //////////////////////////////////////////////////
207    //// First stage: get our hands on the Pixel Data.
208    if ( !fp )
209    {
210       gdcmWarningMacro( "Unavailable file pointer." );
211       return false;
212    }
213
214    fp->seekg( PixelOffset, std::ios::beg );
215    if ( fp->fail() || fp->eof() )
216    {
217       gdcmWarningMacro( "Unable to find PixelOffset in file." );
218       return false;
219    }
220
221    AllocateRaw();
222
223    //////////////////////////////////////////////////
224    //// Second stage: read from disk and decompress.
225    if ( BitsAllocated == 12 )
226    {
227       ReadAndDecompress12BitsTo16Bits( fp);
228    }
229    else if ( IsRaw )
230    {
231       // This problem can be found when some obvious informations are found
232       // after the field containing the image data. In this case, these
233       // bad data are added to the size of the image (in the PixelDataLength
234       // variable). But RawSize is the right size of the image !
235       if ( PixelDataLength != RawSize )
236       {
237          gdcmWarningMacro( "Mismatch between PixelReadConvert : "
238                             << PixelDataLength << " and RawSize : " << RawSize );
239       }
240       if ( PixelDataLength > RawSize )
241       {
242          fp->read( (char*)Raw, RawSize);
243       }
244       else
245       {
246          fp->read( (char*)Raw, PixelDataLength);
247       }
248
249       if ( fp->fail() || fp->eof())
250       {
251          gdcmWarningMacro( "Reading of Raw pixel data failed." );
252          return false;
253       }
254    } 
255    else if ( IsRLELossless )
256    {
257       if ( ! RLEInfo->DecompressRLEFile( fp, Raw, XSize, YSize, ZSize, BitsAllocated ) )
258       {
259          gdcmWarningMacro( "RLE decompressor failed." );
260          return false;
261       }
262    }
263    else if ( IsMPEG )
264    {
265       //gdcmWarningMacro( "Sorry, MPEG not yet taken into account" );
266       //return false;
267       // fp has already been seek to start of mpeg
268       //ReadMPEGFile(fp, Raw, PixelDataLength); 
269       return true;
270    }
271    else
272    {
273       // Default case concerns JPEG family
274       if ( ! ReadAndDecompressJPEGFile( fp ) )
275       {
276          gdcmWarningMacro( "JPEG decompressor failed." );
277          return false;
278       }
279    }
280
281    ////////////////////////////////////////////
282    //// Third stage: twigle the bytes and bits.
283    ConvertReorderEndianity();
284    ConvertReArrangeBits();
285    ConvertFixGreyLevels();
286    if (UserFunction) // user is allowed to Mirror, TopDown, Rotate,...the image
287       UserFunction( Raw, FileInternal);
288    ConvertHandleColor();
289
290    return true;
291 }
292
293 /// Deletes Pixels Area
294 void PixelReadConvert::Squeeze() 
295 {
296    if ( RGB )
297       delete [] RGB;
298    RGB = 0;
299
300    if ( Raw )
301       delete [] Raw;
302    Raw = 0;
303
304    if ( LutRGBA )
305       delete [] LutRGBA;
306    LutRGBA = 0;
307 }
308
309 /**
310  * \brief Build the RGB image from the Raw image and the LUTs.
311  */
312 bool PixelReadConvert::BuildRGBImage()
313 {
314    if ( RGB )
315    {
316       // The job is already done.
317       return true;
318    }
319
320    if ( ! Raw )
321    {
322       // The job can't be done
323       return false;
324    }
325
326    BuildLUTRGBA();
327    if ( ! LutRGBA )
328    {
329       // The job can't be done
330       return false;
331    }
332
333    gdcmWarningMacro( "--> BuildRGBImage" );
334                                                                                 
335    // Build RGB Pixels
336    AllocateRGB();
337    
338    int j;
339    if ( BitsAllocated <= 8 )
340    {
341       uint8_t *localRGB = RGB;
342       for (size_t i = 0; i < RawSize; ++i )
343       {
344          j  = Raw[i] * 4;
345          *localRGB++ = LutRGBA[j];
346          *localRGB++ = LutRGBA[j+1];
347          *localRGB++ = LutRGBA[j+2];
348       }
349     }
350  
351     else  // deal with 16 bits pixels and 16 bits Palette color
352     {
353       uint16_t *localRGB = (uint16_t *)RGB;
354       for (size_t i = 0; i < RawSize/2; ++i )
355       {
356          j  = ((uint16_t *)Raw)[i] * 4;
357          *localRGB++ = ((uint16_t *)LutRGBA)[j];
358          *localRGB++ = ((uint16_t *)LutRGBA)[j+1];
359          *localRGB++ = ((uint16_t *)LutRGBA)[j+2];
360       } 
361     }
362  
363    return true;
364 }
365
366 //-----------------------------------------------------------------------------
367 // Protected
368
369 //-----------------------------------------------------------------------------
370 // Private
371 /**
372  * \brief Read from file a 12 bits per pixel image and decompress it
373  *        into a 16 bits per pixel image.
374  */
375 void PixelReadConvert::ReadAndDecompress12BitsTo16Bits( std::ifstream *fp )
376                throw ( FormatError )
377 {
378    int nbPixels = XSize * YSize;
379    uint16_t *localDecompres = (uint16_t*)Raw;
380
381    for( int p = 0; p < nbPixels; p += 2 )
382    {
383       uint8_t b0, b1, b2;
384
385       fp->read( (char*)&b0, 1);
386       if ( fp->fail() || fp->eof() )
387       {
388          throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
389                                 "Unfound first block" );
390       }
391
392       fp->read( (char*)&b1, 1 );
393       if ( fp->fail() || fp->eof())
394       {
395          throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
396                                 "Unfound second block" );
397       }
398
399       fp->read( (char*)&b2, 1 );
400       if ( fp->fail() || fp->eof())
401       {
402          throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
403                                 "Unfound second block" );
404       }
405
406       // Two steps are necessary to please VC++
407       //
408       // 2 pixels 12bit =     [0xABCDEF]
409       // 2 pixels 16bit = [0x0ABD] + [0x0FCE]
410       //                        A                     B                 D
411       *localDecompres++ =  ((b0 >> 4) << 8) + ((b0 & 0x0f) << 4) + (b1 & 0x0f);
412       //                        F                     C                 E
413       *localDecompres++ =  ((b2 & 0x0f) << 8) + ((b1 >> 4) << 4) + (b2 >> 4);
414
415       /// \todo JPR Troubles expected on Big-Endian processors ?
416    }
417 }
418
419 /**
420  * \brief     Reads from disk the Pixel Data of JPEG Dicom encapsulated
421  *            file and decompress it.
422  * @param     fp File Pointer
423  * @return    Boolean
424  */
425 bool PixelReadConvert::ReadAndDecompressJPEGFile( std::ifstream *fp )
426 {
427    if ( IsJPEG2000 )
428    {
429      // make sure this is the right JPEG compression
430      assert( !IsJPEGLossless || !IsJPEGLossy || !IsJPEGLS );
431      // FIXME this is really ugly but it seems I have to load the complete
432      // jpeg2000 stream to use jasper:
433      // I don't think we'll ever be able to deal with multiple fragments properly
434
435       unsigned long inputlength = 0;
436       JPEGFragment *jpegfrag = JPEGInfo->GetFirstFragment();
437       while( jpegfrag )
438       {
439          inputlength += jpegfrag->GetLength();
440          jpegfrag = JPEGInfo->GetNextFragment();
441       }
442       gdcmAssertMacro( inputlength != 0);
443       uint8_t *inputdata = new uint8_t[inputlength];
444       char *pinputdata = (char*)inputdata;
445       jpegfrag = JPEGInfo->GetFirstFragment();
446       while( jpegfrag )
447       {
448          fp->seekg( jpegfrag->GetOffset(), std::ios::beg);
449          fp->read(pinputdata, jpegfrag->GetLength());
450          pinputdata += jpegfrag->GetLength();
451          jpegfrag = JPEGInfo->GetNextFragment();
452       }
453       // Warning the inputdata buffer is delete in the function
454       if ( ! gdcm_read_JPEG2000_file( Raw, 
455           (char*)inputdata, inputlength ) )
456       {
457          return true;
458       }
459       // wow what happen, must be an error
460       return false;
461    }
462    else if ( IsJPEGLS )
463    {
464      // make sure this is the right JPEG compression
465      assert( !IsJPEGLossless || !IsJPEGLossy || !IsJPEG2000 );
466    // WARNING : JPEG-LS is NOT the 'classical' Jpeg Lossless : 
467    // [JPEG-LS is the basis for new lossless/near-lossless compression
468    // standard for continuous-tone images intended for JPEG2000. The standard
469    // is based on the LOCO-I algorithm (LOw COmplexity LOssless COmpression
470    // for Images) developed at Hewlett-Packard Laboratories]
471    //
472    // see http://datacompression.info/JPEGLS.shtml
473    //
474 #if 0
475    std::cerr << "count:" << JPEGInfo->GetFragmentCount() << std::endl;
476       unsigned long inputlength = 0;
477       JPEGFragment *jpegfrag = JPEGInfo->GetFirstFragment();
478       while( jpegfrag )
479       {
480          inputlength += jpegfrag->GetLength();
481          jpegfrag = JPEGInfo->GetNextFragment();
482       }
483       gdcmAssertMacro( inputlength != 0);
484       uint8_t *inputdata = new uint8_t[inputlength];
485       char *pinputdata = (char*)inputdata;
486       jpegfrag = JPEGInfo->GetFirstFragment();
487       while( jpegfrag )
488       {
489          fp->seekg( jpegfrag->GetOffset(), std::ios::beg);
490          fp->read(pinputdata, jpegfrag->GetLength());
491          pinputdata += jpegfrag->GetLength();
492          jpegfrag = JPEGInfo->GetNextFragment();
493       }  
494       
495   //fp->read((char*)Raw, PixelDataLength);
496
497   std::ofstream out("/tmp/jpegls.jpg");
498   out.write((char*)inputdata, inputlength);
499   out.close();
500   delete[] inputdata;
501 #endif
502
503       gdcmWarningMacro( "Sorry, JPEG-LS not yet taken into account" );
504       fp->seekg( JPEGInfo->GetFirstFragment()->GetOffset(), std::ios::beg);
505 //    if ( ! gdcm_read_JPEGLS_file( fp,Raw ) )
506          return false;
507    }
508    else
509    {
510      // make sure this is the right JPEG compression
511      assert( !IsJPEGLS || !IsJPEG2000 );
512      // Precompute the offset localRaw will be shifted with
513      int length = XSize * YSize * SamplesPerPixel;
514      int numberBytes = BitsAllocated / 8;
515
516      JPEGInfo->DecompressFromFile(fp, Raw, BitsStored, numberBytes, length );
517      return true;
518    }
519 }
520
521 /**
522  * \brief Build Red/Green/Blue/Alpha LUT from File
523  *         when (0028,0004),Photometric Interpretation = [PALETTE COLOR ]
524  *          and (0028,1101),(0028,1102),(0028,1102)
525  *            - xxx Palette Color Lookup Table Descriptor - are found
526  *          and (0028,1201),(0028,1202),(0028,1202)
527  *            - xxx Palette Color Lookup Table Data - are found
528  * \warning does NOT deal with :
529  *   0028 1100 Gray Lookup Table Descriptor (Retired)
530  *   0028 1221 Segmented Red Palette Color Lookup Table Data
531  *   0028 1222 Segmented Green Palette Color Lookup Table Data
532  *   0028 1223 Segmented Blue Palette Color Lookup Table Data
533  *   no known Dicom reader deals with them :-(
534  * @return a RGBA Lookup Table
535  */
536 void PixelReadConvert::BuildLUTRGBA()
537 {
538    if ( LutRGBA )
539    {
540       return;
541    }
542    // Not so easy : see
543    // http://www.barre.nom.fr/medical/dicom2/limitations.html#Color%20Lookup%20Tables
544                                                                                 
545    if ( ! IsPaletteColor )
546    {
547       return;
548    }
549                                                                                 
550    if (   LutRedDescriptor   == GDCM_UNFOUND
551        || LutGreenDescriptor == GDCM_UNFOUND
552        || LutBlueDescriptor  == GDCM_UNFOUND )
553    {
554       gdcmWarningMacro( "(At least) a LUT Descriptor is missing" );
555       return;
556    }
557
558    ////////////////////////////////////////////
559    // Extract the info from the LUT descriptors
560    int lengthR;   // Red LUT length in Bytes
561    int debR;      // Subscript of the first Lut Value
562    int nbitsR;    // Lut item size (in Bits)
563    int nbRead;    // nb of items in LUT descriptor (must be = 3)
564
565    nbRead = sscanf( LutRedDescriptor.c_str(),
566                         "%d\\%d\\%d",
567                         &lengthR, &debR, &nbitsR );
568    if ( nbRead != 3 )
569    {
570       gdcmWarningMacro( "Wrong Red LUT descriptor" );
571    }                                                                                
572    int lengthG;  // Green LUT length in Bytes
573    int debG;     // Subscript of the first Lut Value
574    int nbitsG;   // Lut item size (in Bits)
575
576    nbRead = sscanf( LutGreenDescriptor.c_str(),
577                     "%d\\%d\\%d",
578                     &lengthG, &debG, &nbitsG );  
579    if ( nbRead != 3 )
580    {
581       gdcmWarningMacro( "Wrong Green LUT descriptor" );
582    }
583                                                                                 
584    int lengthB;  // Blue LUT length in Bytes
585    int debB;     // Subscript of the first Lut Value
586    int nbitsB;   // Lut item size (in Bits)
587    nbRead = sscanf( LutRedDescriptor.c_str(),
588                     "%d\\%d\\%d",
589                     &lengthB, &debB, &nbitsB );
590    if ( nbRead != 3 )
591    {
592       gdcmWarningMacro( "Wrong Blue LUT descriptor" );
593    }
594  
595    gdcmWarningMacro(" lengthR " << lengthR << " debR " 
596                  << debR << " nbitsR " << nbitsR);
597    gdcmWarningMacro(" lengthG " << lengthG << " debG " 
598                  << debG << " nbitsG " << nbitsG);
599    gdcmWarningMacro(" lengthB " << lengthB << " debB " 
600                  << debB << " nbitsB " << nbitsB);
601
602    if ( !lengthR ) // if = 2^16, this shall be 0 see : CP-143
603       lengthR=65536;
604    if ( !lengthG ) // if = 2^16, this shall be 0
605       lengthG=65536;
606    if ( !lengthB ) // if = 2^16, this shall be 0
607       lengthB=65536; 
608                                                                                 
609    ////////////////////////////////////////////////////////
610
611    if ( ( ! LutRedData ) || ( ! LutGreenData ) || ( ! LutBlueData ) )
612    {
613       gdcmWarningMacro( "(At least) a LUT is missing" );
614       return;
615    }
616
617    // -------------------------------------------------------------
618    
619    if ( BitsAllocated <= 8 )
620    {
621       // forge the 4 * 8 Bits Red/Green/Blue/Alpha LUT
622       LutRGBA = new uint8_t[ 1024 ]; // 256 * 4 (R, G, B, Alpha)
623       if ( !LutRGBA )
624          return;
625       LutItemNumber = 256;
626       LutItemSize   = 8;
627       memset( LutRGBA, 0, 1024 );
628                                                                                 
629       int mult;
630       if ( ( nbitsR == 16 ) && ( BitsAllocated == 8 ) )
631       {
632          // when LUT item size is different than pixel size
633          mult = 2; // high byte must be = low byte
634       }
635       else
636       {
637          // See PS 3.3-2003 C.11.1.1.2 p 619
638          mult = 1;
639       }
640                                                                                 
641       // if we get a black image, let's just remove the '+1'
642       // from 'i*mult+1' and check again
643       // if it works, we shall have to check the 3 Palettes
644       // to see which byte is ==0 (first one, or second one)
645       // and fix the code
646       // We give up the checking to avoid some (useless ?) overhead
647       // (optimistic asumption)
648       int i;
649       uint8_t *a;
650
651       //take "Subscript of the first Lut Value" (debR,debG,debB) into account!
652
653       //FIXME :  +1 : to get 'low value' byte
654       //         Trouble expected on Big Endian Processors ?
655       //         16 BIts Per Pixel Palette Color to be swapped?
656
657       a = LutRGBA + 0 + debR;
658       for( i=0; i < lengthR; ++i )
659       {
660          *a = LutRedData[i*mult+1]; 
661          a += 4;
662       }
663                                                                                 
664       a = LutRGBA + 1 + debG;
665       for( i=0; i < lengthG; ++i)
666       {
667          *a = LutGreenData[i*mult+1];
668          a += 4;
669       }
670                                                                                 
671       a = LutRGBA + 2 + debB;
672       for(i=0; i < lengthB; ++i)
673       {
674          *a = LutBlueData[i*mult+1];
675          a += 4;
676       }
677                                     
678       a = LutRGBA + 3 ;
679       for(i=0; i < 256; ++i)
680       {
681          *a = 1; // Alpha component
682          a += 4;
683       }
684    }
685    else
686    {
687       // Probabely the same stuff is to be done for 16 Bits Pixels
688       // with 65536 entries LUT ?!?
689       // Still looking for accurate info on the web :-(
690
691       gdcmWarningMacro( "Sorry Palette Color Lookup Tables not yet dealt with"
692                          << " for 16 Bits Per Pixel images" );
693
694       // forge the 4 * 16 Bits Red/Green/Blue/Alpha LUT
695
696       LutRGBA = (uint8_t *)new uint16_t[ 65536*4 ]; // 2^16 * 4 (R, G, B, Alpha)
697       if ( !LutRGBA )
698          return;
699       memset( LutRGBA, 0, 65536*4*2 );  // 16 bits = 2 bytes ;-)
700
701       LutItemNumber = 65536;
702       LutItemSize   = 16;
703
704       int i;
705       uint16_t *a16;
706
707       //take "Subscript of the first Lut Value" (debR,debG,debB) into account!
708
709       a16 = (uint16_t*)LutRGBA + 0 + debR;
710       for( i=0; i < lengthR; ++i )
711       {
712          *a16 = ((uint16_t*)LutRedData)[i];
713          a16 += 4;
714       }
715                                                                               
716       a16 = (uint16_t*)LutRGBA + 1 + debG;
717       for( i=0; i < lengthG; ++i)
718       {
719          *a16 = ((uint16_t*)LutGreenData)[i];
720          a16 += 4;
721       }
722                                                                                 
723       a16 = (uint16_t*)LutRGBA + 2 + debB;
724       for(i=0; i < lengthB; ++i)
725       {
726          *a16 = ((uint16_t*)LutBlueData)[i];
727          a16 += 4;
728       }
729                                                                              
730       a16 = (uint16_t*)LutRGBA + 3 ;
731       for(i=0; i < 65536; ++i)
732       {
733          *a16 = 1; // Alpha component
734          a16 += 4;
735       }
736 /* Just to 'see' the LUT, at debug time
737
738       a16=(uint16_t*)LutRGBA;
739       for (int j=0;j<65536;j++)
740       {
741          std::cout << *a16     << " " << *(a16+1) << " "
742                    << *(a16+2) << " " << *(a16+3) << std::endl;
743          a16+=4;
744       }
745 */
746    }
747 }
748
749 /**
750  * \brief Swap the bytes, according to \ref SwapCode.
751  */
752 void PixelReadConvert::ConvertSwapZone()
753 {
754    unsigned int i;
755
756    if ( BitsAllocated == 16 )
757    {
758       uint16_t *im16 = (uint16_t*)Raw;
759       switch( SwapCode )
760       {
761          case 1234:
762             break;
763          case 3412:
764          case 2143:
765          case 4321:
766             for( i = 0; i < RawSize / 2; i++ )
767             {
768                im16[i]= (im16[i] >> 8) | (im16[i] << 8 );
769             }
770             break;
771          default:
772             gdcmWarningMacro("SwapCode value (16 bits) not allowed.");
773       }
774    }
775    else if ( BitsAllocated == 32 )
776    {
777       uint32_t s32;
778       uint16_t high;
779       uint16_t low;
780       uint32_t *im32 = (uint32_t*)Raw;
781       switch ( SwapCode )
782       {
783          case 1234:
784             break;
785          case 4321:
786             for( i = 0; i < RawSize / 4; i++ )
787             {
788                low     = im32[i] & 0x0000ffff;  // 4321
789                high    = im32[i] >> 16;
790                high    = ( high >> 8 ) | ( high << 8 );
791                low     = ( low  >> 8 ) | ( low  << 8 );
792                s32     = low;
793                im32[i] = ( s32 << 16 ) | high;
794             }
795             break;
796          case 2143:
797             for( i = 0; i < RawSize / 4; i++ )
798             {
799                low     = im32[i] & 0x0000ffff;   // 2143
800                high    = im32[i] >> 16;
801                high    = ( high >> 8 ) | ( high << 8 );
802                low     = ( low  >> 8 ) | ( low  << 8 );
803                s32     = high;
804                im32[i] = ( s32 << 16 ) | low;
805             }
806             break;
807          case 3412:
808             for( i = 0; i < RawSize / 4; i++ )
809             {
810                low     = im32[i] & 0x0000ffff; // 3412
811                high    = im32[i] >> 16;
812                s32     = low;
813                im32[i] = ( s32 << 16 ) | high;
814             }
815             break;
816          default:
817             gdcmWarningMacro("SwapCode value (32 bits) not allowed." );
818       }
819    }
820 }
821
822 /**
823  * \brief Deal with endianness i.e. re-arange bytes inside the integer
824  */
825 void PixelReadConvert::ConvertReorderEndianity()
826 {
827    if ( BitsAllocated != 8 )
828    {
829       ConvertSwapZone();
830    }
831
832    // Special kludge in order to deal with xmedcon broken images:
833    if ( BitsAllocated == 16
834      && BitsStored < BitsAllocated
835      && !PixelSign )
836    {
837       int l = (int)( RawSize / ( BitsAllocated / 8 ) );
838       uint16_t *deb = (uint16_t *)Raw;
839       for(int i = 0; i<l; i++)
840       {
841          if ( *deb == 0xffff )
842          {
843            *deb = 0;
844          }
845          deb++;
846       }
847    }
848 }
849
850 /**
851  * \brief Deal with Grey levels i.e. re-arange them
852  *        to have low values = dark, high values = bright
853  */
854 void PixelReadConvert::ConvertFixGreyLevels()
855 {
856    if (!IsMonochrome1)
857       return;
858
859    uint32_t i; // to please M$VC6
860    int16_t j;
861
862    if (!PixelSign)
863    {
864       if ( BitsAllocated == 8 )
865       {
866          uint8_t *deb = (uint8_t *)Raw;
867          for (i=0; i<RawSize; i++)      
868          {
869             *deb = 255 - *deb;
870             deb++;
871          }
872          return;
873       }
874
875       if ( BitsAllocated == 16 )
876       {
877          uint16_t mask =1;
878          for (j=0; j<BitsStored-1; j++)
879          {
880             mask = (mask << 1) +1; // will be fff when BitsStored=12
881          }
882
883          uint16_t *deb = (uint16_t *)Raw;
884          for (i=0; i<RawSize/2; i++)      
885          {
886             *deb = mask - *deb;
887             deb++;
888          }
889          return;
890        }
891    }
892    else
893    {
894       if ( BitsAllocated == 8 )
895       {
896          uint8_t smask8 = 255;
897          uint8_t *deb = (uint8_t *)Raw;
898          for (i=0; i<RawSize; i++)      
899          {
900             *deb = smask8 - *deb;
901             deb++;
902          }
903          return;
904       }
905       if ( BitsAllocated == 16 )
906       {
907          uint16_t smask16 = 65535;
908          uint16_t *deb = (uint16_t *)Raw;
909          for (i=0; i<RawSize/2; i++)      
910          {
911             *deb = smask16 - *deb;
912             deb++;
913          }
914          return;
915       }
916    }
917 }
918
919 /**
920  * \brief  Re-arrange the bits within the bytes.
921  * @return Boolean always true
922  */
923 bool PixelReadConvert::ConvertReArrangeBits() throw ( FormatError )
924 {
925
926    if ( BitsStored != BitsAllocated )
927    {
928       int l = (int)( RawSize / ( BitsAllocated / 8 ) );
929       if ( BitsAllocated == 16 )
930       {
931          // pmask : to mask the 'unused bits' (may contain overlays)
932          uint16_t pmask = 0xffff;
933          pmask = pmask >> ( BitsAllocated - BitsStored );
934
935          uint16_t *deb = (uint16_t*)Raw;
936
937          if ( !PixelSign )  // Pixels are unsigned
938          {
939             for(int i = 0; i<l; i++)
940             {   
941                *deb = (*deb >> (BitsStored - HighBitPosition - 1)) & pmask;
942                deb++;
943             }
944          }
945          else // Pixels are signed
946          {
947             // smask : to check the 'sign' when BitsStored != BitsAllocated
948             uint16_t smask = 0x0001;
949             smask = smask << ( 16 - (BitsAllocated - BitsStored + 1) );
950             // nmask : to propagate sign bit on negative values
951             int16_t nmask = (int16_t)0x8000;  
952             nmask = nmask >> ( BitsAllocated - BitsStored - 1 );
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