]> Creatis software - gdcm.git/blob - src/gdcmPixelReadConvert.cxx
Should avoid compile error or dash8.kitware/GDCM-g++-64bits
[gdcm.git] / src / gdcmPixelReadConvert.cxx
1 /*=========================================================================
2  
3   Program:   gdcm
4   Module:    $RCSfile: gdcmPixelReadConvert.cxx,v $
5   Language:  C++
6   Date:      $Date: 2005/10/26 13:16:45 $
7   Version:   $Revision: 1.91 $
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    
118    if (! file->IsDicomV3() )  // Should be ACR-NEMA file
119    {
120       IsRaw = true;
121       // Don't loose time checking unexistant Transfer Syntax !
122       IsPrivateGETransferSyntax = IsMPEG
123              = IsJPEG2000 = IsJPEGLS = IsJPEGLossy  
124              = IsJPEGLossless = IsRLELossless 
125              = false;
126    }
127    else
128    {
129       std::string ts = file->GetTransferSyntax();
130
131    
132 //   if ( ts == GDCM_UNKNOWN )
133 //   {
134 //     gdcmErrorMacro( "Could someone tell me how in the world could this happen !" );
135 // -->
136 //--> on ALL acr-nema images ! JPRx
137 //-->
138 //     abort(); // DO NOT REMOVE.  WE SHOULD NEVER READ SUCH IMAGE EVER (only gdcm can write such broekn dicom file)
139 //   }
140
141       IsRaw =
142            ( ! file->IsDicomV3() )  // Should be ACR-NEMA file
143         || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ImplicitVRLittleEndian
144         || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ImplicitVRBigEndianPrivateGE
145         || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ExplicitVRLittleEndian
146         || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ExplicitVRBigEndian
147         || Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::DeflatedExplicitVRLittleEndian;
148      
149       IsPrivateGETransferSyntax = 
150                 ( Global::GetTS()->GetSpecialTransferSyntax(ts) == TS::ImplicitVRBigEndianPrivateGE );
151
152       IsMPEG          = Global::GetTS()->IsMPEG(ts);
153       IsJPEG2000      = Global::GetTS()->IsJPEG2000(ts);
154       IsJPEGLS        = Global::GetTS()->IsJPEGLS(ts);
155       IsJPEGLossy     = Global::GetTS()->IsJPEGLossy(ts);
156       IsJPEGLossless  = Global::GetTS()->IsJPEGLossless(ts);
157       IsRLELossless   = Global::GetTS()->IsRLELossless(ts);
158    }
159
160    PixelOffset     = file->GetPixelOffset();
161    PixelDataLength = file->GetPixelAreaLength();
162    RLEInfo         = file->GetRLEInfo();
163    JPEGInfo        = file->GetJPEGInfo();
164
165    IsMonochrome    = file->IsMonochrome();
166    IsMonochrome1   = file->IsMonochrome1();
167    IsPaletteColor  = file->IsPaletteColor();
168    IsYBRFull       = file->IsYBRFull();
169
170    PlanarConfiguration = file->GetPlanarConfiguration();
171
172    /////////////////////////////////////////////////////////////////
173    // LUT section:
174    HasLUT = file->HasLUT();
175    if ( HasLUT )
176    {
177       // Just in case some access to a File element requires disk access.
178       LutRedDescriptor   = file->GetEntryString( 0x0028, 0x1101 );
179       LutGreenDescriptor = file->GetEntryString( 0x0028, 0x1102 );
180       LutBlueDescriptor  = file->GetEntryString( 0x0028, 0x1103 );
181    
182       // The following comment is probabely meaningless, since LUT are *always*
183       // loaded at parsing time, whatever their length is.
184          
185       // Depending on the value of Document::MAX_SIZE_LOAD_ELEMENT_VALUE
186       // [ refer to invocation of Document::SetMaxSizeLoadEntry() in
187       // Document::Document() ], the loading of the value (content) of a
188       // [Bin|Val]Entry occurence migth have been hindered (read simply NOT
189       // loaded). Hence, we first try to obtain the LUTs data from the file
190       // and when this fails we read the LUTs data directly from disk.
191       // \TODO Reading a [Bin|Val]Entry directly from disk is a kludge.
192       //       We should NOT bypass the [Bin|Val]Entry class. Instead
193       //       an access to an UNLOADED content of a [Bin|Val]Entry occurence
194       //       (e.g. DataEntry::GetBinArea()) should force disk access from
195       //       within the [Bin|Val]Entry class itself. The only problem
196       //       is that the [Bin|Val]Entry is unaware of the FILE* is was
197       //       parsed from. Fix that. FIXME.
198    
199       // //// Red round
200       file->LoadEntryBinArea(0x0028, 0x1201);
201       LutRedData = (uint8_t*)file->GetEntryBinArea( 0x0028, 0x1201 );
202       if ( ! LutRedData )
203       {
204          gdcmWarningMacro("Unable to read Red Palette Color Lookup Table data");
205       }
206
207       // //// Green round:
208       file->LoadEntryBinArea(0x0028, 0x1202);
209       LutGreenData = (uint8_t*)file->GetEntryBinArea(0x0028, 0x1202 );
210       if ( ! LutGreenData)
211       {
212          gdcmWarningMacro("Unable to read Green Palette Color Lookup Table data");
213       }
214
215       // //// Blue round:
216       file->LoadEntryBinArea(0x0028, 0x1203);
217       LutBlueData = (uint8_t*)file->GetEntryBinArea( 0x0028, 0x1203 );
218       if ( ! LutBlueData )
219       {
220          gdcmWarningMacro("Unable to read Blue Palette Color Lookup Table data");
221       }
222    }
223    FileInternal = file;   
224
225    ComputeRawAndRGBSizes();
226 }
227
228 /// \brief Reads from disk and decompresses Pixels
229 bool PixelReadConvert::ReadAndDecompressPixelData( std::ifstream *fp )
230 {
231    // ComputeRawAndRGBSizes is already made by 
232    // ::GrabInformationsFromfile. So, the structure sizes are
233    // correct
234    Squeeze();
235
236    //////////////////////////////////////////////////
237    //// First stage: get our hands on the Pixel Data.
238    if ( !fp )
239    {
240       gdcmWarningMacro( "Unavailable file pointer." );
241       return false;
242    }
243
244    fp->seekg( PixelOffset, std::ios::beg );
245    if ( fp->fail() || fp->eof() )
246    {
247       gdcmWarningMacro( "Unable to find PixelOffset in file." );
248       return false;
249    }
250
251    AllocateRaw();
252
253    //////////////////////////////////////////////////
254    //// Second stage: read from disk and decompress.
255    if ( BitsAllocated == 12 )
256    {
257       ReadAndDecompress12BitsTo16Bits( fp);
258    }
259    else if ( IsRaw )
260    {
261       // This problem can be found when some obvious informations are found
262       // after the field containing the image data. In this case, these
263       // bad data are added to the size of the image (in the PixelDataLength
264       // variable). But RawSize is the right size of the image !
265       if ( PixelDataLength != RawSize )
266       {
267          gdcmWarningMacro( "Mismatch between PixelReadConvert : "
268                           << PixelDataLength << " and RawSize : " << RawSize );
269       }
270       if ( PixelDataLength > RawSize )
271       {
272          fp->read( (char*)Raw, RawSize);
273       }
274       else
275       {
276          fp->read( (char*)Raw, PixelDataLength);
277       }
278
279       if ( fp->fail() || fp->eof())
280       {
281          gdcmWarningMacro( "Reading of Raw pixel data failed." );
282          return false;
283       }
284    } 
285    else if ( IsRLELossless )
286    {
287       if ( ! RLEInfo->DecompressRLEFile
288                                ( fp, Raw, XSize, YSize, ZSize, BitsAllocated ) )
289       {
290          gdcmWarningMacro( "RLE decompressor failed." );
291          return false;
292       }
293    }
294    else if ( IsMPEG )
295    {
296       //gdcmWarningMacro( "Sorry, MPEG not yet taken into account" );
297       //return false;
298       // fp has already been seek to start of mpeg
299       //ReadMPEGFile(fp, Raw, PixelDataLength); 
300       return true;
301    }
302    else
303    {
304       // Default case concerns JPEG family
305       if ( ! ReadAndDecompressJPEGFile( fp ) )
306       {
307          gdcmWarningMacro( "JPEG decompressor failed." );
308          return false;
309       }
310    }
311
312    ////////////////////////////////////////////
313    //// Third stage: twigle the bytes and bits.
314    ConvertReorderEndianity();
315    ConvertReArrangeBits();
316    ConvertFixGreyLevels();
317    if (UserFunction) // user is allowed to Mirror, TopDown, Rotate,...the image
318       UserFunction( Raw, FileInternal);
319    ConvertHandleColor();
320
321    return true;
322 }
323
324 /// Deletes Pixels Area
325 void PixelReadConvert::Squeeze() 
326 {
327    if ( RGB )
328       delete [] RGB;
329    RGB = 0;
330
331    if ( Raw )
332       delete [] Raw;
333    Raw = 0;
334
335    if ( LutRGBA )
336       delete [] LutRGBA;
337    LutRGBA = 0;
338 }
339
340 /**
341  * \brief Build the RGB image from the Raw image and the LUTs.
342  */
343 bool PixelReadConvert::BuildRGBImage()
344 {
345    if ( RGB )
346    {
347       // The job is already done.
348       return true;
349    }
350
351    if ( ! Raw )
352    {
353       // The job can't be done
354       return false;
355    }
356
357    BuildLUTRGBA();
358    if ( ! LutRGBA )
359    {
360       // The job can't be done
361       return false;
362    }
363
364    gdcmWarningMacro( "--> BuildRGBImage" );
365                                                                                 
366    // Build RGB Pixels
367    AllocateRGB();
368    
369    int j;
370    if ( BitsAllocated <= 8 )
371    {
372       uint8_t *localRGB = RGB;
373       for (size_t i = 0; i < RawSize; ++i )
374       {
375          j  = Raw[i] * 4;
376          *localRGB++ = LutRGBA[j];
377          *localRGB++ = LutRGBA[j+1];
378          *localRGB++ = LutRGBA[j+2];
379       }
380     }
381  
382     else  // deal with 16 bits pixels and 16 bits Palette color
383     {
384       uint16_t *localRGB = (uint16_t *)RGB;
385       for (size_t i = 0; i < RawSize/2; ++i )
386       {
387          j  = ((uint16_t *)Raw)[i] * 4;
388          *localRGB++ = ((uint16_t *)LutRGBA)[j];
389          *localRGB++ = ((uint16_t *)LutRGBA)[j+1];
390          *localRGB++ = ((uint16_t *)LutRGBA)[j+2];
391       } 
392     }
393  
394    return true;
395 }
396
397 //-----------------------------------------------------------------------------
398 // Protected
399
400 //-----------------------------------------------------------------------------
401 // Private
402 /**
403  * \brief Read from file a 12 bits per pixel image and decompress it
404  *        into a 16 bits per pixel image.
405  */
406 void PixelReadConvert::ReadAndDecompress12BitsTo16Bits( std::ifstream *fp )
407                throw ( FormatError )
408 {
409    int nbPixels = XSize * YSize;
410    uint16_t *localDecompres = (uint16_t*)Raw;
411
412    for( int p = 0; p < nbPixels; p += 2 )
413    {
414       uint8_t b0, b1, b2;
415
416       fp->read( (char*)&b0, 1);
417       if ( fp->fail() || fp->eof() )
418       {
419          throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
420                                 "Unfound first block" );
421       }
422
423       fp->read( (char*)&b1, 1 );
424       if ( fp->fail() || fp->eof())
425       {
426          throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
427                                 "Unfound second block" );
428       }
429
430       fp->read( (char*)&b2, 1 );
431       if ( fp->fail() || fp->eof())
432       {
433          throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
434                                 "Unfound second block" );
435       }
436
437       // Two steps are necessary to please VC++
438       //
439       // 2 pixels 12bit =     [0xABCDEF]
440       // 2 pixels 16bit = [0x0ABD] + [0x0FCE]
441       //                        A                     B                 D
442       *localDecompres++ =  ((b0 >> 4) << 8) + ((b0 & 0x0f) << 4) + (b1 & 0x0f);
443       //                        F                     C                 E
444       *localDecompres++ =  ((b2 & 0x0f) << 8) + ((b1 >> 4) << 4) + (b2 >> 4);
445
446       /// \todo JPR Troubles expected on Big-Endian processors ?
447    }
448 }
449
450 /**
451  * \brief     Reads from disk the Pixel Data of JPEG Dicom encapsulated
452  *            file and decompress it.
453  * @param     fp File Pointer
454  * @return    Boolean
455  */
456 bool PixelReadConvert::ReadAndDecompressJPEGFile( std::ifstream *fp )
457 {
458    if ( IsJPEG2000 )
459    {
460      // make sure this is the right JPEG compression
461      assert( !IsJPEGLossless || !IsJPEGLossy || !IsJPEGLS );
462      // FIXME this is really ugly but it seems I have to load the complete
463      // jpeg2000 stream to use jasper:
464      // I don't think we'll ever be able to deal with multiple fragments properly
465
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       // Warning the inputdata buffer is delete in the function
485       if ( ! gdcm_read_JPEG2000_file( Raw, 
486           (char*)inputdata, inputlength ) )
487       {
488          return true;
489       }
490       // wow what happen, must be an error
491       return false;
492    }
493    else if ( IsJPEGLS )
494    {
495      // make sure this is the right JPEG compression
496      assert( !IsJPEGLossless || !IsJPEGLossy || !IsJPEG2000 );
497    // WARNING : JPEG-LS is NOT the 'classical' Jpeg Lossless : 
498    // [JPEG-LS is the basis for new lossless/near-lossless compression
499    // standard for continuous-tone images intended for JPEG2000. The standard
500    // is based on the LOCO-I algorithm (LOw COmplexity LOssless COmpression
501    // for Images) developed at Hewlett-Packard Laboratories]
502    //
503    // see http://datacompression.info/JPEGLS.shtml
504    //
505 #if 0
506    std::cerr << "count:" << JPEGInfo->GetFragmentCount() << std::endl;
507       unsigned long inputlength = 0;
508       JPEGFragment *jpegfrag = JPEGInfo->GetFirstFragment();
509       while( jpegfrag )
510       {
511          inputlength += jpegfrag->GetLength();
512          jpegfrag = JPEGInfo->GetNextFragment();
513       }
514       gdcmAssertMacro( inputlength != 0);
515       uint8_t *inputdata = new uint8_t[inputlength];
516       char *pinputdata = (char*)inputdata;
517       jpegfrag = JPEGInfo->GetFirstFragment();
518       while( jpegfrag )
519       {
520          fp->seekg( jpegfrag->GetOffset(), std::ios::beg);
521          fp->read(pinputdata, jpegfrag->GetLength());
522          pinputdata += jpegfrag->GetLength();
523          jpegfrag = JPEGInfo->GetNextFragment();
524       }  
525       
526   //fp->read((char*)Raw, PixelDataLength);
527
528   std::ofstream out("/tmp/jpegls.jpg");
529   out.write((char*)inputdata, inputlength);
530   out.close();
531   delete[] inputdata;
532 #endif
533
534       gdcmWarningMacro( "Sorry, JPEG-LS not yet taken into account" );
535       fp->seekg( JPEGInfo->GetFirstFragment()->GetOffset(), std::ios::beg);
536 //    if ( ! gdcm_read_JPEGLS_file( fp,Raw ) )
537          return false;
538    }
539    else
540    {
541      // make sure this is the right JPEG compression
542      assert( !IsJPEGLS || !IsJPEG2000 );
543      // Precompute the offset localRaw will be shifted with
544      int length = XSize * YSize * SamplesPerPixel;
545      int numberBytes = BitsAllocated / 8;
546
547      JPEGInfo->DecompressFromFile(fp, Raw, BitsStored, numberBytes, length );
548      return true;
549    }
550 }
551
552 /**
553  * \brief Build Red/Green/Blue/Alpha LUT from File when :
554  *         - (0028,0004) : Photometric Interpretation == [PALETTE COLOR ]
555  *         and
556  *         - (0028,1101),(0028,1102),(0028,1102)
557  *            xxx Palette Color Lookup Table Descriptor are found
558  *          and 
559  *         - (0028,1201),(0028,1202),(0028,1202)
560  *           xxx Palette Color Lookup Table Data - are found
561  * \warning does NOT deal with :
562  *   - 0028 1100 Gray Lookup Table Descriptor (Retired)
563  *   - 0028 1221 Segmented Red Palette Color Lookup Table Data
564  *   - 0028 1222 Segmented Green Palette Color Lookup Table Data
565  *   - 0028 1223 Segmented Blue Palette Color Lookup Table Data
566  *   no known Dicom reader deals with them :-(
567  * @return a RGBA Lookup Table
568  */
569 void PixelReadConvert::BuildLUTRGBA()
570 {
571
572    // Note to code reviewers :
573    // The problem is *much more* complicated, since a lot of manufacturers
574    // Don't follow the norm :
575    // have a look at David Clunie's remark at the end of this .cxx file.
576    if ( LutRGBA )
577    
578    {
579       return;
580    }
581    // Not so easy : see
582    // http://www.barre.nom.fr/medical/dicom2/limitations.html#Color%20Lookup%20Tables
583                                                                                 
584    if ( ! IsPaletteColor )
585    {
586       return;
587    }
588                                                                                 
589    if (   LutRedDescriptor   == GDCM_UNFOUND
590        || LutGreenDescriptor == GDCM_UNFOUND
591        || LutBlueDescriptor  == GDCM_UNFOUND )
592    {
593       gdcmWarningMacro( "(At least) a LUT Descriptor is missing" );
594       return;
595    }
596
597    ////////////////////////////////////////////
598    // Extract the info from the LUT descriptors
599    int lengthR;   // Red LUT length in Bytes
600    int debR;      // Subscript of the first Lut Value
601    int nbitsR;    // Lut item size (in Bits)
602    int nbRead;    // nb of items in LUT descriptor (must be = 3)
603
604    nbRead = sscanf( LutRedDescriptor.c_str(),
605                         "%d\\%d\\%d",
606                         &lengthR, &debR, &nbitsR );
607    if ( nbRead != 3 )
608    {
609       gdcmWarningMacro( "Wrong Red LUT descriptor" );
610    }                                                                                
611    int lengthG;  // Green LUT length in Bytes
612    int debG;     // Subscript of the first Lut Value
613    int nbitsG;   // Lut item size (in Bits)
614
615    nbRead = sscanf( LutGreenDescriptor.c_str(),
616                     "%d\\%d\\%d",
617                     &lengthG, &debG, &nbitsG );  
618    if ( nbRead != 3 )
619    {
620       gdcmWarningMacro( "Wrong Green LUT descriptor" );
621    }
622                                                                                 
623    int lengthB;  // Blue LUT length in Bytes
624    int debB;     // Subscript of the first Lut Value
625    int nbitsB;   // Lut item size (in Bits)
626    nbRead = sscanf( LutRedDescriptor.c_str(),
627                     "%d\\%d\\%d",
628                     &lengthB, &debB, &nbitsB );
629    if ( nbRead != 3 )
630    {
631       gdcmWarningMacro( "Wrong Blue LUT descriptor" );
632    }
633  
634    gdcmWarningMacro(" lengthR " << lengthR << " debR " 
635                  << debR << " nbitsR " << nbitsR);
636    gdcmWarningMacro(" lengthG " << lengthG << " debG " 
637                  << debG << " nbitsG " << nbitsG);
638    gdcmWarningMacro(" lengthB " << lengthB << " debB " 
639                  << debB << " nbitsB " << nbitsB);
640
641    if ( !lengthR ) // if = 2^16, this shall be 0 see : CP-143
642       lengthR=65536;
643    if ( !lengthG ) // if = 2^16, this shall be 0
644       lengthG=65536;
645    if ( !lengthB ) // if = 2^16, this shall be 0
646       lengthB=65536; 
647                                                                                 
648    ////////////////////////////////////////////////////////
649
650    if ( ( ! LutRedData ) || ( ! LutGreenData ) || ( ! LutBlueData ) )
651    {
652       gdcmWarningMacro( "(At least) a LUT is missing" );
653       return;
654    }
655
656    // -------------------------------------------------------------
657    
658    if ( BitsAllocated <= 8 )
659    {
660       // forge the 4 * 8 Bits Red/Green/Blue/Alpha LUT
661       LutRGBA = new uint8_t[ 1024 ]; // 256 * 4 (R, G, B, Alpha)
662       if ( !LutRGBA )
663          return;
664       LutItemNumber = 256;
665       LutItemSize   = 8;
666       memset( LutRGBA, 0, 1024 );
667                                                                                 
668       int mult;
669       if ( ( nbitsR == 16 ) && ( BitsAllocated == 8 ) )
670       {
671          // when LUT item size is different than pixel size
672          mult = 2; // high byte must be = low byte
673       }
674       else
675       {
676          // See PS 3.3-2003 C.11.1.1.2 p 619
677          mult = 1;
678       }
679                                                                                 
680       // if we get a black image, let's just remove the '+1'
681       // from 'i*mult+1' and check again
682       // if it works, we shall have to check the 3 Palettes
683       // to see which byte is ==0 (first one, or second one)
684       // and fix the code
685       // We give up the checking to avoid some (useless ?) overhead
686       // (optimistic asumption)
687       int i;
688       uint8_t *a;
689
690       //take "Subscript of the first Lut Value" (debR,debG,debB) into account!
691
692       //FIXME :  +1 : to get 'low value' byte
693       //         Trouble expected on Big Endian Processors ?
694       //         16 BIts Per Pixel Palette Color to be swapped?
695
696       a = LutRGBA + 0 + debR;
697       for( i=0; i < lengthR; ++i )
698       {
699          *a = LutRedData[i*mult+1]; 
700          a += 4;
701       }
702                                                                                 
703       a = LutRGBA + 1 + debG;
704       for( i=0; i < lengthG; ++i)
705       {
706          *a = LutGreenData[i*mult+1];
707          a += 4;
708       }
709                                                                                 
710       a = LutRGBA + 2 + debB;
711       for(i=0; i < lengthB; ++i)
712       {
713          *a = LutBlueData[i*mult+1];
714          a += 4;
715       }
716                                     
717       a = LutRGBA + 3 ;
718       for(i=0; i < 256; ++i)
719       {
720          *a = 1; // Alpha component
721          a += 4;
722       }
723    }
724    else
725    {
726       // Probabely the same stuff is to be done for 16 Bits Pixels
727       // with 65536 entries LUT ?!?
728       // Still looking for accurate info on the web :-(
729
730       gdcmWarningMacro( "Sorry Palette Color Lookup Tables not yet dealt with"
731                          << " for 16 Bits Per Pixel images" );
732
733       // forge the 4 * 16 Bits Red/Green/Blue/Alpha LUT
734
735       LutRGBA = (uint8_t *)new uint16_t[ 65536*4 ]; // 2^16 * 4 (R, G, B, Alpha)
736       if ( !LutRGBA )
737          return;
738       memset( LutRGBA, 0, 65536*4*2 );  // 16 bits = 2 bytes ;-)
739
740       LutItemNumber = 65536;
741       LutItemSize   = 16;
742
743       int i;
744       uint16_t *a16;
745
746       //take "Subscript of the first Lut Value" (debR,debG,debB) into account!
747
748       a16 = (uint16_t*)LutRGBA + 0 + debR;
749       for( i=0; i < lengthR; ++i )
750       {
751          *a16 = ((uint16_t*)LutRedData)[i];
752          a16 += 4;
753       }
754                                                                               
755       a16 = (uint16_t*)LutRGBA + 1 + debG;
756       for( i=0; i < lengthG; ++i)
757       {
758          *a16 = ((uint16_t*)LutGreenData)[i];
759          a16 += 4;
760       }
761                                                                                 
762       a16 = (uint16_t*)LutRGBA + 2 + debB;
763       for(i=0; i < lengthB; ++i)
764       {
765          *a16 = ((uint16_t*)LutBlueData)[i];
766          a16 += 4;
767       }
768                                                                              
769       a16 = (uint16_t*)LutRGBA + 3 ;
770       for(i=0; i < 65536; ++i)
771       {
772          *a16 = 1; // Alpha component
773          a16 += 4;
774       }
775 /* Just to 'see' the LUT, at debug time
776 // Don't remove this commented out code.
777
778       a16=(uint16_t*)LutRGBA;
779       for (int j=0;j<65536;j++)
780       {
781          std::cout << *a16     << " " << *(a16+1) << " "
782                    << *(a16+2) << " " << *(a16+3) << std::endl;
783          a16+=4;
784       }
785 */
786    }
787 }
788
789 /**
790  * \brief Swap the bytes, according to \ref SwapCode.
791  */
792 void PixelReadConvert::ConvertSwapZone()
793 {
794    unsigned int i;
795    uint16_t localSwapCode = SwapCode;
796    
797    // If this file is 'ImplicitVR BigEndian PrivateGE Transfer Syntax', 
798    // then the header is in little endian format and the pixel data is in 
799    // big endian format.  When reading the header, GDCM has already established
800    // a byte swapping code suitable for this machine to read the
801    // header. In TS::ImplicitVRBigEndianPrivateGE, this code will need
802    // to be switched in order to read the pixel data.  This must be
803    // done REGARDLESS of the processor endianess!
804    //
805    // Example:  Assume we are on a little endian machine.  When
806    // GDCM reads the header, the header will match the machine
807    // endianess and the swap code will be established as a no-op.
808    // When GDCM reaches the pixel data, it will need to switch the
809    // swap code to do big endian to little endian conversion.
810    //
811    // Now, assume we are on a big endian machine.  When GDCM reads the
812    // header, the header will be recognized as a different endianess
813    // than the machine endianess, and a swap code will be established
814    // to convert from little endian to big endian.  When GDCM readers
815    // the pixel data, the pixel data endianess will now match the
816    // machine endianess.  But we currently have a swap code that
817    // converts from little endian to big endian.  In this case, we
818    // need to switch the swap code to a no-op.
819    //
820    // Therefore, in either case, if the file is in
821    // 'ImplicitVR BigEndian PrivateGE Transfer Syntax', then GDCM needs to switch
822    // the byte swapping code when entering the pixel data.
823
824 /* //Let me check something.
825    //I wait for the Dashboard !   
826    if ( IsPrivateGETransferSyntax )
827    {
828       // PrivateGETransferSyntax only exists for 'true' Dicom images
829       // we assume there is no 'exotic' 32 bits endianess!
830       switch (localSwapCode)
831       {
832          case 1234:
833             localSwapCode = 4321;
834             break;
835          case 4321:
836             localSwapCode = 1234;
837             break;
838       }  
839    }
840 */   
841    if ( BitsAllocated == 16 )
842    {
843       uint16_t *im16 = (uint16_t*)Raw;
844       switch( localSwapCode )
845       {
846          case 1234:
847             break;
848          case 3412:
849          case 2143:
850          case 4321:
851             for( i = 0; i < RawSize / 2; i++ )
852             {
853                im16[i]= (im16[i] >> 8) | (im16[i] << 8 );
854             }
855             break;
856          default:
857             gdcmWarningMacro("SwapCode value (16 bits) not allowed.");
858       }
859    }
860    else if ( BitsAllocated == 32 )
861    {
862       uint32_t s32;
863       uint16_t high;
864       uint16_t low;
865       uint32_t *im32 = (uint32_t*)Raw;
866       switch ( localSwapCode )
867       {
868          case 1234:
869             break;
870          case 4321:
871             for( i = 0; i < RawSize / 4; i++ )
872             {
873                low     = im32[i] & 0x0000ffff;  // 4321
874                high    = im32[i] >> 16;
875                high    = ( high >> 8 ) | ( high << 8 );
876                low     = ( low  >> 8 ) | ( low  << 8 );
877                s32     = low;
878                im32[i] = ( s32 << 16 ) | high;
879             }
880             break;
881          case 2143:
882             for( i = 0; i < RawSize / 4; i++ )
883             {
884                low     = im32[i] & 0x0000ffff;   // 2143
885                high    = im32[i] >> 16;
886                high    = ( high >> 8 ) | ( high << 8 );
887                low     = ( low  >> 8 ) | ( low  << 8 );
888                s32     = high;
889                im32[i] = ( s32 << 16 ) | low;
890             }
891             break;
892          case 3412:
893             for( i = 0; i < RawSize / 4; i++ )
894             {
895                low     = im32[i] & 0x0000ffff; // 3412
896                high    = im32[i] >> 16;
897                s32     = low;
898                im32[i] = ( s32 << 16 ) | high;
899             }
900             break;
901          default:
902             gdcmWarningMacro("SwapCode value (32 bits) not allowed." );
903       }
904    }
905 }
906
907 /**
908  * \brief Deal with endianness i.e. re-arange bytes inside the integer
909  */
910 void PixelReadConvert::ConvertReorderEndianity()
911 {
912    if ( BitsAllocated != 8 )
913    {
914       ConvertSwapZone();
915    }
916
917    // Special kludge in order to deal with xmedcon broken images:
918    if ( BitsAllocated == 16
919      && BitsStored < BitsAllocated
920      && !PixelSign )
921    {
922       int l = (int)( RawSize / ( BitsAllocated / 8 ) );
923       uint16_t *deb = (uint16_t *)Raw;
924       for(int i = 0; i<l; i++)
925       {
926          if ( *deb == 0xffff )
927          {
928            *deb = 0;
929          }
930          deb++;
931       }
932    }
933 }
934
935 /**
936  * \brief Deal with Grey levels i.e. re-arange them
937  *        to have low values = dark, high values = bright
938  */
939 void PixelReadConvert::ConvertFixGreyLevels()
940 {
941    if (!IsMonochrome1)
942       return;
943
944    uint32_t i; // to please M$VC6
945    int16_t j;
946
947    if (!PixelSign)
948    {
949       if ( BitsAllocated == 8 )
950       {
951          uint8_t *deb = (uint8_t *)Raw;
952          for (i=0; i<RawSize; i++)      
953          {
954             *deb = 255 - *deb;
955             deb++;
956          }
957          return;
958       }
959
960       if ( BitsAllocated == 16 )
961       {
962          uint16_t mask =1;
963          for (j=0; j<BitsStored-1; j++)
964          {
965             mask = (mask << 1) +1; // will be fff when BitsStored=12
966          }
967
968          uint16_t *deb = (uint16_t *)Raw;
969          for (i=0; i<RawSize/2; i++)      
970          {
971             *deb = mask - *deb;
972             deb++;
973          }
974          return;
975        }
976    }
977    else
978    {
979       if ( BitsAllocated == 8 )
980       {
981          uint8_t smask8 = 255;
982          uint8_t *deb = (uint8_t *)Raw;
983          for (i=0; i<RawSize; i++)      
984          {
985             *deb = smask8 - *deb;
986             deb++;
987          }
988          return;
989       }
990       if ( BitsAllocated == 16 )
991       {
992          uint16_t smask16 = 65535;
993          uint16_t *deb = (uint16_t *)Raw;
994          for (i=0; i<RawSize/2; i++)      
995          {
996             *deb = smask16 - *deb;
997             deb++;
998          }
999          return;
1000       }
1001    }
1002 }
1003
1004 /**
1005  * \brief  Re-arrange the bits within the bytes.
1006  * @return Boolean always true
1007  */
1008 bool PixelReadConvert::ConvertReArrangeBits() throw ( FormatError )
1009 {
1010
1011    if ( BitsStored != BitsAllocated )
1012    {
1013       int l = (int)( RawSize / ( BitsAllocated / 8 ) );
1014       if ( BitsAllocated == 16 )
1015       {
1016          // pmask : to mask the 'unused bits' (may contain overlays)
1017          uint16_t pmask = 0xffff;
1018          pmask = pmask >> ( BitsAllocated - BitsStored );
1019
1020          uint16_t *deb = (uint16_t*)Raw;
1021
1022          if ( !PixelSign )  // Pixels are unsigned
1023          {
1024             for(int i = 0; i<l; i++)
1025             {   
1026                *deb = (*deb >> (BitsStored - HighBitPosition - 1)) & pmask;
1027                deb++;
1028             }
1029          }
1030          else // Pixels are signed
1031          {
1032             // smask : to check the 'sign' when BitsStored != BitsAllocated
1033             uint16_t smask = 0x0001;
1034             smask = smask << ( 16 - (BitsAllocated - BitsStored + 1) );
1035             // nmask : to propagate sign bit on negative values
1036             int16_t nmask = (int16_t)0x8000;  
1037             nmask = nmask >> ( BitsAllocated - BitsStored - 1 );
1038  
1039             for(int i = 0; i<l; i++)
1040             {
1041                *deb = *deb >> (BitsStored - HighBitPosition - 1);
1042                if ( *deb & smask )
1043                {
1044                   *deb = *deb | nmask;
1045                }
1046                else
1047                {
1048                   *deb = *deb & pmask;
1049                }
1050                deb++;
1051             }
1052          }
1053       }
1054       else if ( BitsAllocated == 32 )
1055       {
1056          // pmask : to mask the 'unused bits' (may contain overlays)
1057          uint32_t pmask = 0xffffffff;
1058          pmask = pmask >> ( BitsAllocated - BitsStored );
1059
1060          uint32_t *deb = (uint32_t*)Raw;
1061
1062          if ( !PixelSign )
1063          {
1064             for(int i = 0; i<l; i++)
1065             {             
1066                *deb = (*deb >> (BitsStored - HighBitPosition - 1)) & pmask;
1067                deb++;
1068             }
1069          }
1070          else
1071          {
1072             // smask : to check the 'sign' when BitsStored != BitsAllocated
1073             uint32_t smask = 0x00000001;
1074             smask = smask >> ( 32 - (BitsAllocated - BitsStored +1 ));
1075             // nmask : to propagate sign bit on negative values
1076             int32_t nmask = 0x80000000;   
1077             nmask = nmask >> ( BitsAllocated - BitsStored -1 );
1078
1079             for(int i = 0; i<l; i++)
1080             {
1081                *deb = *deb >> (BitsStored - HighBitPosition - 1);
1082                if ( *deb & smask )
1083                   *deb = *deb | nmask;
1084                else
1085                   *deb = *deb & pmask;
1086                deb++;
1087             }
1088          }
1089       }
1090       else
1091       {
1092          gdcmWarningMacro("Weird image (BitsAllocated !=8, 12, 16, 32)");
1093          throw FormatError( "Weird image !?" );
1094       }
1095    }
1096    return true;
1097 }
1098
1099 /**
1100  * \brief   Convert (Red plane, Green plane, Blue plane) to RGB pixels
1101  * \warning Works on all the frames at a time
1102  */
1103 void PixelReadConvert::ConvertRGBPlanesToRGBPixels()
1104 {
1105    gdcmWarningMacro("--> ConvertRGBPlanesToRGBPixels");
1106
1107    uint8_t *localRaw = Raw;
1108    uint8_t *copyRaw = new uint8_t[ RawSize ];
1109    memmove( copyRaw, localRaw, RawSize );
1110
1111    int l = XSize * YSize * ZSize;
1112
1113    uint8_t *a = copyRaw;
1114    uint8_t *b = copyRaw + l;
1115    uint8_t *c = copyRaw + l + l;
1116
1117    for (int j = 0; j < l; j++)
1118    {
1119       *(localRaw++) = *(a++);
1120       *(localRaw++) = *(b++);
1121       *(localRaw++) = *(c++);
1122    }
1123    delete[] copyRaw;
1124 }
1125
1126 /**
1127  * \brief   Convert (cY plane, cB plane, cR plane) to RGB pixels
1128  * \warning Works on all the frames at a time
1129  */
1130 void PixelReadConvert::ConvertYcBcRPlanesToRGBPixels()
1131 {
1132   // Remarks for YBR newbees :
1133   // YBR_FULL works very much like RGB, i.e. three samples per pixel, 
1134   // just the color space is YCbCr instead of RGB. This is particularly useful
1135   // for doppler ultrasound where most of the image is grayscale 
1136   // (i.e. only populates the Y components) and Cb and Cr are mostly zero,
1137   // except for the few patches of color on the image.
1138   // On such images, RLE achieves a compression ratio that is much better 
1139   // than the compression ratio on an equivalent RGB image. 
1140  
1141    gdcmWarningMacro("--> ConvertYcBcRPlanesToRGBPixels");
1142    
1143    uint8_t *localRaw = Raw;
1144    uint8_t *copyRaw = new uint8_t[ RawSize ];
1145    memmove( copyRaw, localRaw, RawSize );
1146
1147    // to see the tricks about YBR_FULL, YBR_FULL_422,
1148    // YBR_PARTIAL_422, YBR_ICT, YBR_RCT have a look at :
1149    // ftp://medical.nema.org/medical/dicom/final/sup61_ft.pdf
1150    // and be *very* affraid
1151    //
1152    int l        = XSize * YSize;
1153    int nbFrames = ZSize;
1154
1155    uint8_t *a = copyRaw + 0;
1156    uint8_t *b = copyRaw + l;
1157    uint8_t *c = copyRaw + l+ l;
1158    int32_t R, G, B;
1159
1160    ///  We replaced easy to understand but time consuming floating point
1161    ///  computations by the 'well known' integer computation counterpart
1162    ///  Refer to :
1163    ///            http://lestourtereaux.free.fr/papers/data/yuvrgb.pdf
1164    ///  for code optimisation.
1165
1166    for ( int i = 0; i < nbFrames; i++ )
1167    {
1168       for ( int j = 0; j < l; j++ )
1169       {
1170          R = 38142 *(*a-16) + 52298 *(*c -128);
1171          G = 38142 *(*a-16) - 26640 *(*c -128) - 12845 *(*b -128);
1172          B = 38142 *(*a-16) + 66093 *(*b -128);
1173
1174          R = (R+16384)>>15;
1175          G = (G+16384)>>15;
1176          B = (B+16384)>>15;
1177
1178          if (R < 0)   R = 0;
1179          if (G < 0)   G = 0;
1180          if (B < 0)   B = 0;
1181          if (R > 255) R = 255;
1182          if (G > 255) G = 255;
1183          if (B > 255) B = 255;
1184
1185          *(localRaw++) = (uint8_t)R;
1186          *(localRaw++) = (uint8_t)G;
1187          *(localRaw++) = (uint8_t)B;
1188          a++;
1189          b++;
1190          c++;
1191       }
1192    }
1193    delete[] copyRaw;
1194 }
1195
1196 /// \brief Deals with the color decoding i.e. handle:
1197 ///   - R, G, B planes (as opposed to RGB pixels)
1198 ///   - YBR (various) encodings.
1199 ///   - LUT[s] (or "PALETTE COLOR").
1200
1201 void PixelReadConvert::ConvertHandleColor()
1202 {
1203    //////////////////////////////////
1204    // Deal with the color decoding i.e. handle:
1205    //   - R, G, B planes (as opposed to RGB pixels)
1206    //   - YBR (various) encodings.
1207    //   - LUT[s] (or "PALETTE COLOR").
1208    //
1209    // The classification in the color decoding schema is based on the blending
1210    // of two Dicom tags values:
1211    // * "Photometric Interpretation" for which we have the cases:
1212    //  - [Photo A] MONOCHROME[1|2] pictures,
1213    //  - [Photo B] RGB or YBR_FULL_422 (which acts as RGB),
1214    //  - [Photo C] YBR_* (with the above exception of YBR_FULL_422)
1215    //  - [Photo D] "PALETTE COLOR" which indicates the presence of LUT[s].
1216    // * "Planar Configuration" for which we have the cases:
1217    //  - [Planar 0] 0 then Pixels are already RGB
1218    //  - [Planar 1] 1 then we have 3 planes : R, G, B,
1219    //  - [Planar 2] 2 then we have 1 gray Plane and 3 LUTs
1220    //
1221    // Now in theory, one could expect some coherence when blending the above
1222    // cases. For example we should not encounter files belonging at the
1223    // time to case [Planar 0] and case [Photo D].
1224    // Alas, this was only theory ! Because in practice some odd (read ill
1225    // formated Dicom) files (e.g. gdcmData/US-PAL-8-10x-echo.dcm) we encounter:
1226    //     - "Planar Configuration" = 0,
1227    //     - "Photometric Interpretation" = "PALETTE COLOR".
1228    // Hence gdcm will use the folowing "heuristic" in order to be tolerant
1229    // towards Dicom-non-conformant files:
1230    //   << whatever the "Planar Configuration" value might be, a
1231    //      "Photometric Interpretation" set to "PALETTE COLOR" forces
1232    //      a LUT intervention >>
1233    //
1234    // Now we are left with the following handling of the cases:
1235    // - [Planar 0] OR  [Photo A] no color decoding (since respectively
1236    //       Pixels are already RGB and monochrome pictures have no color :),
1237    // - [Planar 1] AND [Photo B] handled with ConvertRGBPlanesToRGBPixels()
1238    // - [Planar 1] AND [Photo C] handled with ConvertYcBcRPlanesToRGBPixels()
1239    // - [Planar 2] OR  [Photo D] requires LUT intervention.
1240
1241    gdcmWarningMacro("--> ConvertHandleColor"
1242                     << "Planar Configuration " << PlanarConfiguration );
1243
1244    if ( ! IsRawRGB() )
1245    {
1246       // [Planar 2] OR  [Photo D]: LUT intervention done outside
1247       gdcmWarningMacro("--> RawRGB : LUT intervention done outside");
1248       return;
1249    }
1250                                                                                 
1251    if ( PlanarConfiguration == 1 )
1252    {
1253       if ( IsYBRFull )
1254       {
1255          // [Planar 1] AND [Photo C] (remember YBR_FULL_422 acts as RGB)
1256          gdcmWarningMacro("--> YBRFull");
1257          ConvertYcBcRPlanesToRGBPixels();
1258       }
1259       else
1260       {
1261          // [Planar 1] AND [Photo C]
1262          gdcmWarningMacro("--> YBRFull");
1263          ConvertRGBPlanesToRGBPixels();
1264       }
1265       return;
1266    }
1267                                                                                 
1268    // When planarConf is 0, and RLELossless (forbidden by Dicom norm)
1269    // pixels need to be RGB-fyied anyway
1270
1271    if (IsRLELossless)
1272    { 
1273      gdcmWarningMacro("--> RLE Lossless");
1274      ConvertRGBPlanesToRGBPixels();
1275    }
1276
1277    // In *normal *case, when planarConf is 0, pixels are already in RGB
1278 }
1279
1280 /// Computes the Pixels Size
1281 void PixelReadConvert::ComputeRawAndRGBSizes()
1282 {
1283    int bitsAllocated = BitsAllocated;
1284    // Number of "Bits Allocated" is fixed to 16 when it's 12, since
1285    // in this case we will expand the image to 16 bits (see
1286    //    \ref ReadAndDecompress12BitsTo16Bits() )
1287    if (  BitsAllocated == 12 )
1288    {
1289       bitsAllocated = 16;
1290    }
1291                                                                                 
1292    RawSize =  XSize * YSize * ZSize
1293                      * ( bitsAllocated / 8 )
1294                      * SamplesPerPixel;
1295    if ( HasLUT )
1296    {
1297       RGBSize = 3 * RawSize; // works for 8 and 16 bits per Pixel
1298    }
1299    else
1300    {
1301       RGBSize = RawSize;
1302    }
1303 }
1304
1305 /// Allocates room for RGB Pixels
1306 void PixelReadConvert::AllocateRGB()
1307 {
1308   if ( RGB )
1309      delete [] RGB;
1310   RGB = new uint8_t[RGBSize];
1311 }
1312
1313 /// Allocates room for RAW Pixels
1314 void PixelReadConvert::AllocateRaw()
1315 {
1316   if ( Raw )
1317      delete [] Raw;
1318   Raw = new uint8_t[RawSize];
1319 }
1320
1321 //-----------------------------------------------------------------------------
1322 // Print
1323 /**
1324  * \brief        Print self.
1325  * @param indent Indentation string to be prepended during printing.
1326  * @param os     Stream to print to.
1327  */
1328 void PixelReadConvert::Print( std::ostream &os, std::string const &indent )
1329 {
1330    os << indent
1331       << "--- Pixel information -------------------------"
1332       << std::endl;
1333    os << indent
1334       << "Pixel Data: offset " << PixelOffset
1335       << " x(" << std::hex << PixelOffset << std::dec
1336       << ")   length " << PixelDataLength
1337       << " x(" << std::hex << PixelDataLength << std::dec
1338       << ")" << std::endl;
1339
1340    if ( IsRLELossless )
1341    {
1342       if ( RLEInfo )
1343       {
1344          RLEInfo->Print( os, indent );
1345       }
1346       else
1347       {
1348          gdcmWarningMacro("Set as RLE file but NO RLEinfo present.");
1349       }
1350    }
1351
1352    if ( IsJPEG2000 || IsJPEGLossless || IsJPEGLossy || IsJPEGLS )
1353    {
1354       if ( JPEGInfo )
1355       {
1356          JPEGInfo->Print( os, indent );
1357       }
1358       else
1359       {
1360          gdcmWarningMacro("Set as JPEG file but NO JPEGinfo present.");
1361       }
1362    }
1363 }
1364
1365 //-----------------------------------------------------------------------------
1366 } // end namespace gdcm
1367
1368 // Note to developpers :
1369 // Here is a very detailled post from David Clunie, on the troubles caused 
1370 // 'non standard' LUT and LUT description
1371 // We shall have to take it into accound in our code.
1372 // Some day ...
1373
1374
1375 /*
1376 Subject: Problem with VOI LUTs in Agfa and Fuji CR and GE DX images, was Re: VOI LUT issues
1377 Date: Sun, 06 Feb 2005 17:13:40 GMT
1378 From: David Clunie <dclunie@dclunie.com>
1379 Reply-To: dclunie@dclunie.com
1380 Newsgroups: comp.protocols.dicom
1381 References: <1107553502.040221.189550@o13g2000cwo.googlegroups.com>
1382
1383 > THE LUT that comes with [my] image claims to be 16-bit, but none of the
1384 > values goes higher than 4095.  That being said, though, none of my
1385 > original pixel values goes higher than that, either.  I have read
1386 > elsewhere on this group that when that happens you are supposed to
1387 > adjust the LUT.  Can someone be more specific?  There was a thread from
1388 > 2002 where Marco and David were mentioning doing precisely that.
1389 >
1390 > Thanks
1391 >
1392 > -carlos rodriguez
1393
1394
1395 You have encountered the well known "we know what the standard says but
1396 we are going to ignore it and do what we have been doing for almost
1397 a decade regardless" CR vendor bug. Agfa started this, but they are not
1398 the only vendor doing this now; GE and Fuji may have joined the club.
1399
1400 Sadly, one needs to look at the LUT Data, figure out what the maximum
1401 value actually encoded is, and find the next highest power of 2 (e.g.
1402 212 in this case), to figure out what the range of the data is
1403 supposed to be. I have assumed that if the maximum value in the LUT
1404 data is less than a power of 2 minus 1 (e.g. 0xebc) then the intent
1405 of the vendor was not to use the maximum available grayscale range
1406 of the display (e.g. the maximum is 0xfff in this case). An alternative
1407 would be to scale to the actual maximum rather than a power of two.
1408
1409 Very irritating, and in theory not totally reliable if one really
1410 intended the full 16 bits and only used, say 15, but that is extremely
1411 unlikely since everything would be too dark, and this heuristic
1412 seems to work OK.
1413
1414 There has never been anything in the standard that describes having
1415 to go through these convolutions. Since the only value in the
1416 standard that describes the bit depth of the LUT values is LUT
1417 Descriptor value 3 and that is (usually) always required to be
1418 either 8 or 16, it mystifies me how the creators' of these images
1419 imagine that the receiver is going to divine the range that is intended. Further, the standard is quite explicit that this 3rd
1420 value defines the range of LUT values, but as far as I am aware, all
1421 the vendors are ignoring the standard and indeed sending a third value
1422 of 16 in all cases.
1423
1424 This problem is not confined to CR, and is also seen with DX products.
1425
1426 Typically I have seen:
1427
1428 - Agfa CR, which usually (always ?) sends LUTs, values up to 0x0fff
1429 - Fuji CR, which occasionally send LUTs, values up to 0x03ff
1430 - GE DX, for presentation, which always have LUTs, up to 0x3fff
1431
1432 Swissray, Siemens, Philips, Canon and Kodak never seem to send VOI LUTs
1433 at this point (which is a whole other problem). Note that the presence
1434 or absence of a VOI LUT as opposed to window values may be configurable
1435 on the modality in some cases, and I have just looked at what I happen
1436 to have received from a myriad of sites over whose configuration I have
1437 no control. This may be why the majority of Fuji images have no VOI LUTs,
1438 but a few do (or it may be the Siemens system that these Fuji images went
1439 through that perhaps added it). I do have some test Hologic DX images that
1440 are not from a clinical site that do actually get this right (a value
1441 of 12 for the third value and a max of 0xfff).
1442
1443 Since almost every vendor that I have encountered that encodes LUTs
1444 makes this mistake, perhaps it is time to amend the standard to warn
1445 implementor's of receivers and/or sanction this bad behavior. We have
1446 talked about this in the past in WG 6 but so far everyone has been
1447 reluctant to write into the standard such a comment. Maybe it is time
1448 to try again, since if one is not aware of this problem, one cannot
1449 effectively implement display using VOI LUTs, and there is a vast
1450 installed base to contend with.
1451
1452 I did not check presentation states, in which VOI LUTs could also be
1453 encountered, for the prevalence of this mistake, nor did I look at the
1454 encoding of Modality LUT's, which are unusual. Nor did I check digital
1455 mammography images. I would be interested to hear from anyone who has.
1456
1457 David
1458
1459 PS. The following older thread in this newsgroup discusses this:
1460
1461 "http://groups-beta.google.com/group/comp.protocols.dicom/browse_frm/t hread/6a033444802a35fc/0f0a9a1e35c1468e?q=voi+lut&_done=%2Fgroup%2Fcom p.protocols.dicom%2Fsearch%3Fgroup%3Dcomp.protocols.dicom%26q%3Dvoi+lu t%26qt_g%3D1%26searchnow%3DSearch+this+group%26&_doneTitle=Back+to+Sea rch&&d#0f0a9a1e35c1468e"
1462
1463 PPS. From a historical perspective, the following may be of interest.
1464
1465 In the original standard in 1993, all that was said about this was a
1466 reference to the corresponding such where Modality LUTs are described
1467 that said:
1468
1469 "The third value specifies the number of bits for each entry in the
1470 LUT Data. It shall take the value 8 or 16. The LUT Data shall be stored
1471 in a format equivalent to 8 or 16 bits allocated and high bit equal
1472 1-bits allocated."
1473
1474 Since the high bit hint was not apparently explicit enough, a very
1475 early CP, CP 15 (submitted by Agfa as it happens), replaced this with:
1476
1477 "The third value conveys the range of LUT entry values. It shall take
1478 the value 8 or 16, corresponding with the LUT entry value range of
1479 256 or 65536.
1480
1481 Note:    The third value is not required for describing the
1482     LUT data and is only included for informational usage
1483     and for maintaining compatibility with ACRNEMA 2.0.
1484
1485 The LUT Data contains the LUT entry values."
1486
1487 That is how it read in the 1996, 1998 and 1999 editions.
1488
1489 By the 2000 edition, Supplement 33 that introduced presentation states
1490 extensively reworked this entire section and tried to explain this in
1491 different words:
1492
1493 "The output range is from 0 to 2^n-1 where n is the third value of LUT
1494 Descriptor. This range is always unsigned."
1495
1496 and also added a note to spell out what the output range meant in the
1497 VOI LUT section:
1498
1499 "9. The output of the Window Center/Width or VOI LUT transformation
1500 is either implicitly scaled to the full range of the display device
1501 if there is no succeeding transformation defined, or implicitly scaled
1502 to the full input range of the succeeding transformation step (such as
1503 the Presentation LUT), if present. See C.11.6.1."
1504
1505 It still reads this way in the 2004 edition.
1506
1507 Note that LUTs in other applications than the general VOI LUT allow for
1508 values other than 8 or 16 in the third value of LUT descriptor to permit
1509 ranges other than 0 to 255 or 65535.
1510
1511 In addition, the DX Image Module specializes the VOI LUT
1512 attributes as follows, in PS 3.3 section C.8.11.3.1.5 (added in Sup 32):
1513
1514 "The third value specifies the number of bits for each entry in the LUT
1515 Data (analogous to ìbits storedî). It shall be between 10-16. The LUT
1516 Data shall be stored in a format equivalent to 16 ìbits allocatedî and
1517 ìhigh bitî equal to ìbits storedî - 1. The third value conveys the range
1518 of LUT entry values. These unsigned LUT entry values shall range between
1519 0 and 2^n-1, where n is the third value of the LUT Descriptor."
1520
1521 So in the case of the GE DX for presentation images, the third value of
1522 LUT descriptor is allowed to be and probably should be 14 rather than 16.
1523
1524 */