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