]> Creatis software - gdcm.git/blob - src/gdcmFile.cxx
* CLEANUP_ROUND for gdcmPixelConvert:
[gdcm.git] / src / gdcmFile.cxx
1   /*=========================================================================
2                                                                                 
3   Program:   gdcm
4   Module:    $RCSfile: gdcmFile.cxx,v $
5   Language:  C++
6   Date:      $Date: 2004/09/29 17:33:17 $
7   Version:   $Revision: 1.132 $
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 "gdcmFile.h"
20 #include "gdcmDebug.h"
21 #include "jpeg/ljpg/jpegless.h"
22
23 typedef std::pair<TagDocEntryHT::iterator,TagDocEntryHT::iterator> IterHT;
24
25 //-------------------------------------------------------------------------
26 // Constructor / Destructor
27 /**
28  * \brief Constructor dedicated to deal with the *pixels* area of a ACR/DICOMV3
29  *        file (gdcmHeader only deals with the ... header)
30  *        Opens (in read only and when possible) an existing file and checks
31  *        for DICOM compliance. Returns NULL on failure.
32  *        It will be up to the user to load the pixels into memory
33  *        (see GetImageData, GetImageDataRaw)
34  * \note  the in-memory representation of all available tags found in
35  *        the DICOM header is post-poned to first header information access.
36  *        This avoid a double parsing of public part of the header when
37  *        user sets an a posteriori shadow dictionary (efficiency can be
38  *        seen as a side effect).   
39  * @param header already built gdcmHeader
40  */
41 gdcmFile::gdcmFile(gdcmHeader *header)
42 {
43    Header     = header;
44    SelfHeader = false;
45    Initialise();
46 }
47
48 /**
49  * \brief Constructor dedicated to deal with the *pixels* area of a ACR/DICOMV3
50  *        file (gdcmHeader only deals with the ... header)
51  *        Opens (in read only and when possible) an existing file and checks
52  *        for DICOM compliance. Returns NULL on failure.
53  *        It will be up to the user to load the pixels into memory
54  *        (see GetImageData, GetImageDataRaw)
55  * \note  the in-memory representation of all available tags found in
56  *        the DICOM header is post-poned to first header information access.
57  *        This avoid a double parsing of public part of the header when
58  *        one sets an a posteriori shadow dictionary (efficiency can be
59  *        seen as a side effect).   
60  * @param filename file to be opened for parsing
61  */
62 gdcmFile::gdcmFile(std::string const & filename )
63 {
64    Header = new gdcmHeader( filename );
65    SelfHeader = true;
66    Initialise();
67 }
68
69 /**
70  * \brief Factorization for various forms of constructors.
71  */
72 void gdcmFile::Initialise()
73 {
74    if ( Header->IsReadable() )
75    {
76       ImageDataSizeRaw = ComputeDecompressedPixelDataSizeFromHeader();
77       if ( Header->HasLUT() )
78       {
79          ImageDataSize = 3 * ImageDataSizeRaw;
80       }
81       else
82       {
83          ImageDataSize = ImageDataSizeRaw;
84       }
85    }
86    SaveInitialValues();
87 }
88
89 /**
90  * \brief canonical destructor
91  * \note  If the gdcmHeader was created by the gdcmFile constructor,
92  *        it is destroyed by the gdcmFile
93  */
94 gdcmFile::~gdcmFile()
95
96    if( SelfHeader )
97    {
98       delete Header;
99    }
100    Header = 0;
101
102    DeleteInitialValues();
103 }
104
105 /**
106  * \brief Sets some initial values for the Constructor
107  * \warning not end user intended
108  */
109 void gdcmFile::SaveInitialValues()
110
111
112    PixelRead  = -1; // no ImageData read yet.
113    LastAllocatedPixelDataLength = 0;
114    Pixel_Data = 0;
115
116    InitialSpp = "";     
117    InitialPhotInt = "";
118    InitialPlanConfig = "";
119    InitialBitsAllocated = "";
120    InitialHighBit = "";
121   
122    InitialRedLUTDescr   = 0;
123    InitialGreenLUTDescr = 0;
124    InitialBlueLUTDescr  = 0;
125    InitialRedLUTData    = 0;
126    InitialGreenLUTData  = 0;
127    InitialBlueLUTData   = 0; 
128                 
129    if ( Header->IsReadable() )
130    {
131       // the following values *may* be modified 
132       // by gdcmFile::GetImageDataIntoVectorRaw
133       // we save their initial value.
134       InitialSpp           = Header->GetEntryByNumber(0x0028,0x0002);
135       InitialPhotInt       = Header->GetEntryByNumber(0x0028,0x0004);
136       InitialPlanConfig    = Header->GetEntryByNumber(0x0028,0x0006);
137       
138       InitialBitsAllocated = Header->GetEntryByNumber(0x0028,0x0100);
139       InitialHighBit       = Header->GetEntryByNumber(0x0028,0x0102);
140
141       // the following entries *may* be removed from the H table
142       // (NOT deleted ...) by gdcmFile::GetImageDataIntoVectorRaw  
143       // we keep a pointer on them.
144       InitialRedLUTDescr   = Header->GetDocEntryByNumber(0x0028,0x1101);
145       InitialGreenLUTDescr = Header->GetDocEntryByNumber(0x0028,0x1102);
146       InitialBlueLUTDescr  = Header->GetDocEntryByNumber(0x0028,0x1103);
147
148       InitialRedLUTData    = Header->GetDocEntryByNumber(0x0028,0x1201);
149       InitialGreenLUTData  = Header->GetDocEntryByNumber(0x0028,0x1202);
150       InitialBlueLUTData   = Header->GetDocEntryByNumber(0x0028,0x1203); 
151    }
152 }
153
154 /**
155  * \brief restores some initial values
156  * \warning not end user intended
157  */
158 void gdcmFile::RestoreInitialValues()
159 {   
160    if ( Header->IsReadable() )
161    {      
162       // the following values *may* have been modified 
163       // by gdcmFile::GetImageDataIntoVectorRaw
164       // we restore their initial value.
165       if ( InitialSpp != "")
166          Header->SetEntryByNumber(InitialSpp,0x0028,0x0002);
167       if ( InitialPhotInt != "")
168          Header->SetEntryByNumber(InitialPhotInt,0x0028,0x0004);
169       if ( InitialPlanConfig != "")
170
171          Header->SetEntryByNumber(InitialPlanConfig,0x0028,0x0006);
172       if ( InitialBitsAllocated != "")
173           Header->SetEntryByNumber(InitialBitsAllocated,0x0028,0x0100);
174       if ( InitialHighBit != "")
175           Header->SetEntryByNumber(InitialHighBit,0x0028,0x0102);
176                
177       // the following entries *may* be have been removed from the H table
178       // (NOT deleted ...) by gdcmFile::GetImageDataIntoVectorRaw  
179       // we restore them.
180
181       if (InitialRedLUTDescr)
182          Header->AddEntry(InitialRedLUTDescr);
183       if (InitialGreenLUTDescr)
184          Header->AddEntry(InitialGreenLUTDescr);
185       if (InitialBlueLUTDescr)
186          Header->AddEntry(InitialBlueLUTDescr);
187
188       if (InitialRedLUTData)
189          Header->AddEntry(InitialBlueLUTDescr);
190       if (InitialGreenLUTData)
191          Header->AddEntry(InitialGreenLUTData);
192       if (InitialBlueLUTData)
193          Header->AddEntry(InitialBlueLUTData);
194    }
195 }
196
197 /**
198  * \brief delete initial values (il they were saved)
199  *        of InitialLutDescriptors and InitialLutData
200  */
201 void gdcmFile::DeleteInitialValues()
202
203
204 // InitialLutDescriptors and InitialLutData
205 // will have to be deleted if the don't belong any longer
206 // to the Header H table when the header is deleted...
207
208    if ( InitialRedLUTDescr )           
209       delete InitialRedLUTDescr;
210   
211    if ( InitialGreenLUTDescr )
212       delete InitialGreenLUTDescr;
213       
214    if ( InitialBlueLUTDescr )      
215       delete InitialBlueLUTDescr; 
216        
217    if ( InitialRedLUTData )      
218       delete InitialRedLUTData;
219    
220    if ( InitialGreenLUTData != NULL)
221       delete InitialGreenLUTData;
222       
223    if ( InitialBlueLUTData != NULL)      
224       delete InitialBlueLUTData;      
225 }
226
227 //-----------------------------------------------------------------------------
228 // Print
229
230 //-----------------------------------------------------------------------------
231 // Public
232
233 /**
234  * \brief     computes the length (in bytes) we must ALLOCATE to receive the
235  *            image(s) pixels (multiframes taken into account) 
236  * \warning : it is NOT the group 7FE0 length
237  *          (no interest for compressed images).
238  */
239 int gdcmFile::ComputeDecompressedPixelDataSizeFromHeader()
240 {
241    // see PS 3.3-2003 : C.7.6.3.2.1  
242    // 
243    //   MONOCHROME1
244    //   MONOCHROME2
245    //   PALETTE COLOR
246    //   RGB
247    //   HSV  (Retired)
248    //   ARGB (Retired)
249    //   CMYK (Retired)
250    //   YBR_FULL
251    //   YBR_FULL_422 (no LUT, no Palette)
252    //   YBR_PARTIAL_422
253    //   YBR_ICT
254    //   YBR_RCT
255
256    // LUT's
257    // ex : gdcm-US-ALOKA-16.dcm
258    // 0028|1221 [OW]   [Segmented Red Palette Color Lookup Table Data]
259    // 0028|1222 [OW]   [Segmented Green Palette Color Lookup Table Data]  
260    // 0028|1223 [OW]   [Segmented Blue Palette Color Lookup Table Data]
261
262    // ex  : OT-PAL-8-face.dcm
263    // 0028|1201 [US]   [Red Palette Color Lookup Table Data]
264    // 0028|1202 [US]   [Green Palette Color Lookup Table Data]
265    // 0028|1203 [US]   [Blue Palette Color Lookup Table Data]
266
267    // Number of "Bits Allocated"
268    int numberBitsAllocated = Header->GetBitsAllocated();
269    if ( ( numberBitsAllocated == 0 ) || ( numberBitsAllocated == 12 ) )
270    {
271       numberBitsAllocated = 16;
272    } 
273
274    int DecompressedSize = Header->GetXSize()
275                         * Header->GetYSize() 
276                         * Header->GetZSize()
277                         * ( numberBitsAllocated / 8 )
278                         * Header->GetSamplesPerPixel();
279    
280    return DecompressedSize;
281 }
282
283 /**
284  * \brief   - Allocates necessary memory, 
285  *          - Reads the pixels from disk (uncompress if necessary),
286  *          - Transforms YBR pixels, if any, into RGB pixels
287  *          - Transforms 3 planes R, G, B, if any, into a single RGB Plane
288  *          - Transforms single Grey plane + 3 Palettes into a RGB Plane
289  *          - Copies the pixel data (image[s]/volume[s]) to newly allocated zone.
290  * @return  Pointer to newly allocated pixel data.
291  *          NULL if alloc fails 
292  */
293 uint8_t* gdcmFile::GetImageData()
294 {
295    // FIXME (Mathieu)
296    // I need to deallocate Pixel_Data before doing any allocation:
297    
298    if ( Pixel_Data )
299      if ( LastAllocatedPixelDataLength != ImageDataSize ) 
300         free(Pixel_Data);
301    if ( !Pixel_Data )
302       Pixel_Data = new uint8_t[ImageDataSize];
303     
304    if ( Pixel_Data )
305    {
306       LastAllocatedPixelDataLength = ImageDataSize;
307
308       // we load the pixels (and transform grey level + LUT into RGB)
309       GetImageDataIntoVector(Pixel_Data, ImageDataSize);
310
311       // We say the value *is* loaded.
312       GetHeader()->SetEntryByNumber( GDCM_BINLOADED,
313          GetHeader()->GetGrPixel(), GetHeader()->GetNumPixel());
314
315       // Will be 7fe0, 0010 in standard case
316       GetHeader()->SetEntryBinAreaByNumber( Pixel_Data, 
317          GetHeader()->GetGrPixel(), GetHeader()->GetNumPixel()); 
318    }      
319    PixelRead = 0; // no PixelRaw
320
321    return Pixel_Data;
322 }
323
324 /**
325  * \brief
326  *          Read the pixels from disk (uncompress if necessary),
327  *          Transforms YBR pixels, if any, into RGB pixels
328  *          Transforms 3 planes R, G, B, if any, into a single RGB Plane
329  *          Transforms single Grey plane + 3 Palettes into a RGB Plane   
330  *          Copies at most MaxSize bytes of pixel data to caller allocated
331  *          memory space.
332  * \warning This function allows people that want to build a volume
333  *          from an image stack *not to* have, first to get the image pixels, 
334  *          and then move them to the volume area.
335  *          It's absolutely useless for any VTK user since vtk chooses 
336  *          to invert the lines of an image, that is the last line comes first
337  *          (for some axis related reasons?). Hence he will have 
338  *          to load the image line by line, starting from the end.
339  *          VTK users have to call GetImageData
340  *     
341  * @param   destination Address (in caller's memory space) at which the
342  *          pixel data should be copied
343  * @param   maxSize Maximum number of bytes to be copied. When MaxSize
344  *          is not sufficient to hold the pixel data the copy is not
345  *          executed (i.e. no partial copy).
346  * @return  On success, the number of bytes actually copied. Zero on
347  *          failure e.g. MaxSize is lower than necessary.
348  */
349 size_t gdcmFile::GetImageDataIntoVector (void* destination, size_t maxSize)
350 {
351    GetImageDataIntoVectorRaw (destination, maxSize);
352    PixelRead = 0 ; // =0 : no ImageDataRaw 
353    if ( !Header->HasLUT() )
354    {
355       return ImageDataSize;
356    }
357                             
358    // from Lut R + Lut G + Lut B
359    uint8_t *newDest = new uint8_t[ImageDataSize];
360    uint8_t *a       = (uint8_t *)destination;
361    uint8_t *lutRGBA = Header->GetLUTRGBA();
362
363    if ( lutRGBA )
364    {
365       int j;
366       // move Gray pixels to temp area  
367       memmove(newDest, destination, ImageDataSizeRaw);
368       for (size_t i=0; i<ImageDataSizeRaw; ++i) 
369       {
370          // Build RGB Pixels
371          j    = newDest[i]*4;
372          *a++ = lutRGBA[j];
373          *a++ = lutRGBA[j+1];
374          *a++ = lutRGBA[j+2];
375       }
376       delete[] newDest;
377     
378    // now, it's an RGB image
379    // Lets's write it in the Header
380
381    // FIXME : Better use CreateOrReplaceIfExist ?
382
383    std::string spp = "3";        // Samples Per Pixel
384    Header->SetEntryByNumber(spp,0x0028,0x0002);
385    std::string rgb = "RGB ";     // Photometric Interpretation
386    Header->SetEntryByNumber(rgb,0x0028,0x0004);
387    std::string planConfig = "0"; // Planar Configuration
388    Header->SetEntryByNumber(planConfig,0x0028,0x0006);
389
390    }
391    else  // GetLUTRGBA() failed
392    { 
393       // (gdcm-US-ALOKA-16.dcm), contains Segmented xxx Palette Color 
394       // that are *more* than 65535 long ?!? 
395       // No idea how to manage such an image !
396       // Need to make RGB Pixels (?) from grey Pixels (?!) and Gray Lut  (!?!)
397       // It seems that *no Dicom Viewer* has any idea :-(
398         
399       std::string photomInterp = "MONOCHROME1 ";  // Photometric Interpretation
400       Header->SetEntryByNumber(photomInterp,0x0028,0x0004);
401    } 
402
403    /// \todo Drop Palette Color out of the Header?
404    return ImageDataSize; 
405 }
406
407 /**
408  * \brief   Allocates necessary memory, 
409  *          Transforms YBR pixels (if any) into RGB pixels
410  *          Transforms 3 planes R, G, B  (if any) into a single RGB Plane
411  *          Copies the pixel data (image[s]/volume[s]) to newly allocated zone. 
412  *          DOES NOT transform Grey plane + 3 Palettes into a RGB Plane
413  * @return  Pointer to newly allocated pixel data.
414  * \        NULL if alloc fails 
415  */
416 uint8_t* gdcmFile::GetImageDataRaw ()
417 {
418    size_t imgDataSize;
419    if ( Header->HasLUT() )
420       /// \todo Let gdcmHeader user a chance to get the right value
421       imgDataSize = ImageDataSizeRaw;
422    else 
423       imgDataSize = ImageDataSize;
424     
425    // FIXME (Mathieu)
426    // I need to deallocate Pixel_Data before doing any allocation:
427    
428    if ( Pixel_Data )
429       if ( LastAllocatedPixelDataLength != imgDataSize )
430          free(Pixel_Data);
431    if ( !Pixel_Data ) 
432       Pixel_Data = new uint8_t[imgDataSize];
433
434    if ( Pixel_Data )
435    {
436       LastAllocatedPixelDataLength = imgDataSize;
437       
438       // we load the pixels ( grey level or RGB, but NO transformation)
439        GetImageDataIntoVectorRaw(Pixel_Data, imgDataSize);
440
441       // We say the value *is* loaded.
442       GetHeader()->SetEntryByNumber( GDCM_BINLOADED,
443          GetHeader()->GetGrPixel(), GetHeader()->GetNumPixel());
444  
445       // will be 7fe0, 0010 in standard cases
446       GetHeader()->SetEntryBinAreaByNumber(Pixel_Data, 
447          GetHeader()->GetGrPixel(), GetHeader()->GetNumPixel());
448    } 
449    PixelRead = 1; // PixelRaw
450
451    return Pixel_Data;
452 }
453
454 /**
455  * \brief   Copies at most MaxSize bytes of pixel data to caller's
456  *          memory space.
457  * \warning This function was designed to avoid people that want to build
458  *          a volume from an image stack to need first to get the image pixels 
459  *          and then move them to the volume area.
460  *          It's absolutely useless for any VTK user since vtk chooses 
461  *          to invert the lines of an image, that is the last line comes first
462  *          (for some axis related reasons?). Hence he will have 
463  *          to load the image line by line, starting from the end.
464  *          VTK users hace to call GetImageData
465  * \warning DOES NOT transform the Grey Plane + Palette Color (if any) 
466  *                   into a single RGB Pixels Plane
467  *          the (VTK) user will manage the palettes
468  *     
469  * @param   destination Address (in caller's memory space) at which the
470  *          pixel data should be copied
471  * @param   maxSize Maximum number of bytes to be copied. When MaxSize
472  *          is not sufficient to hold the pixel data the copy is not
473  *          executed (i.e. no partial copy).
474  * @return  On success, the number of bytes actually copied. Zero on
475  *          failure e.g. MaxSize is lower than necessary.
476  */
477 size_t gdcmFile::GetImageDataIntoVectorRaw (void* destination, size_t maxSize)
478 {
479   // we save the initial values of the following
480   // in order to be able to restore the header in a disk-consistent state
481   // (if user asks twice to get the pixels from disk)
482
483    if ( PixelRead != -1 ) // File was "read" before
484    {  
485       RestoreInitialValues(); 
486    }
487    
488    PixelRead = 1 ; // PixelRaw
489     
490    if ( ImageDataSize > maxSize )
491    {
492       dbg.Verbose(0, "gdcmFile::GetImageDataIntoVector: pixel data bigger"
493                      "than caller's expected MaxSize");
494       return (size_t)0;
495    }
496
497    ReadPixelData( destination );
498
499    // Number of Bits Allocated for storing a Pixel
500    int numberBitsAllocated = Header->GetBitsAllocated();
501    if ( numberBitsAllocated == 0 )
502    {
503       numberBitsAllocated = 16;
504    }
505
506    // Number of Bits actually used
507    int numberBitsStored = Header->GetBitsStored();
508    if ( numberBitsStored == 0 )
509    {
510       numberBitsStored = numberBitsAllocated;
511    }
512
513    // High Bit Position
514    int highBitPosition = Header->GetHighBitPosition();
515    if ( highBitPosition == 0 )
516    {
517       highBitPosition = numberBitsAllocated - 1;
518    }
519
520    bool signedPixel = Header->IsSignedPixelData();
521
522 // SPLIT ME
523 ////////////////////////////////////////////////////////
524 // ENDIANITY SECTION: re-arange bytes inside the integer
525    if ( numberBitsAllocated != 8 )
526    {
527       SwapZone( destination, Header->GetSwapCode(), ImageDataSize,
528                 numberBitsAllocated );
529    }
530      
531    // to avoid pb with some xmedcon breakers images 
532    if (  ( numberBitsAllocated == 16 )
533       && ( numberBitsStored < numberBitsAllocated )
534       && ( ! signedPixel ) )
535    {
536       int l = (int)(ImageDataSize / (numberBitsAllocated/8));
537       uint16_t *deb = (uint16_t *)destination;
538       for(int i = 0; i<l; i++)
539       {
540          if( *deb == 0xffff )
541          {
542            *deb = 0;
543          }
544          deb++;   
545       }
546    }
547
548 // SPLIT ME
549 //////////////////////////////////
550 // re arange bits inside the bytes
551    if ( numberBitsStored != numberBitsAllocated )
552    {
553       int l = (int)(ImageDataSize / (numberBitsAllocated/8));
554       if ( numberBitsAllocated == 16 )
555       {
556          uint16_t mask = 0xffff;
557          mask = mask >> ( numberBitsAllocated - numberBitsStored );
558          uint16_t *deb = (uint16_t *)destination;
559          for(int i = 0; i<l; i++)
560          {
561             *deb = (*deb >> (numberBitsStored - highBitPosition - 1)) & mask;
562             deb++;
563          }
564       }
565       else if ( numberBitsAllocated == 32 )
566       {
567          uint32_t mask = 0xffffffff;
568          mask         = mask >> ( numberBitsAllocated - numberBitsStored );
569          uint32_t *deb = (uint32_t *)destination;
570          for(int i = 0; i<l; i++)
571          {
572             *deb = (*deb >> (numberBitsStored - highBitPosition - 1)) & mask;
573             deb++;
574          }
575       }
576       else
577       {
578          dbg.Verbose(0, "gdcmFile::GetImageDataIntoVector: weird image");
579          return 0;
580       }
581    }
582
583 #ifdef GDCM_DEBUG
584    FILE*  DebugFile;
585    DebugFile = fopen( "SpuriousFile.RAW", "wb" );
586    fwrite( PixelConvertor.GetUncompressed(),
587            PixelConvertor.GetUncompressedsSize(),
588            1, DebugFile );
589    fclose( DebugFile );
590 #endif //GDCM_DEBUG
591
592 // SPLIT ME
593 //////////////////////////////////
594 // Deal with the color
595    
596    // Monochrome pictures don't require color intervention
597    if ( Header->IsMonochrome() )
598    {
599       return ImageDataSize; 
600    }
601       
602    // Planar configuration = 0 : Pixels are already RGB
603    // Planar configuration = 1 : 3 planes : R, G, B
604    // Planar configuration = 2 : 1 gray Plane + 3 LUT
605
606    // Well ... supposed to be !
607    // See US-PAL-8-10x-echo.dcm: PlanarConfiguration=0,
608    //                            PhotometricInterpretation=PALETTE COLOR
609    // and heuristic has to be found :-( 
610
611    int planConf = Header->GetPlanarConfiguration();  // 0028,0006
612
613    // Whatever Planar Configuration is, 
614    // "PALETTE COLOR " implies that we deal with the palette. 
615    if ( Header->IsPaletteColor() )
616    {
617       planConf = 2;
618    }
619
620    switch ( planConf )
621    {
622       case 0:
623          // Pixels are already RGB
624          break;
625       case 1:
626          if ( Header->IsYBRFull() )
627          {
628             // Warning : YBR_FULL_422 acts as RGB
629             //         : we need to make RGB Pixels from Planes Y,cB,cR
630
631             // to see the tricks about YBR_FULL, YBR_FULL_422, 
632             // YBR_PARTIAL_422, YBR_ICT, YBR_RCT have a look at :
633             // ftp://medical.nema.org/medical/dicom/final/sup61_ft.pdf
634             // and be *very* affraid
635             //
636             int l        = Header->GetXSize() * Header->GetYSize();
637             int nbFrames = Header->GetZSize();
638
639             uint8_t* newDest = new uint8_t[ImageDataSize];
640             uint8_t* x       = newDest;
641             uint8_t* a       = (uint8_t*)destination;
642             uint8_t* b       = a + l;
643             uint8_t* c       = b + l;
644             double R,G,B;
645
646             /// \todo : Replace by the 'well known' integer computation
647             /// counterpart
648             /// see http://lestourtereaux.free.fr/papers/data/yuvrgb.pdf
649             /// for code optimisation
650     
651             for (int i = 0; i < nbFrames; i++)
652             {
653                for (int j = 0; j < l; j++)
654                {
655                   R = 1.164 *(*a-16) + 1.596 *(*c -128) + 0.5;
656                   G = 1.164 *(*a-16) - 0.813 *(*c -128) - 0.392 *(*b -128) + 0.5;
657                   B = 1.164 *(*a-16) + 2.017 *(*b -128) + 0.5;
658
659                   if (R < 0.0)   R = 0.0;
660                   if (G < 0.0)   G = 0.0;
661                   if (B < 0.0)   B = 0.0;
662                   if (R > 255.0) R = 255.0;
663                   if (G > 255.0) G = 255.0;
664                   if (B > 255.0) B = 255.0;
665
666                   *(x++) = (uint8_t)R;
667                   *(x++) = (uint8_t)G;
668                   *(x++) = (uint8_t)B;
669                   a++; b++; c++;  
670                }
671             }
672             memmove(destination, newDest, ImageDataSize);
673             delete[] newDest;
674          }
675          else
676          {
677             // need to make RGB Pixels from R,G,B Planes
678             // (all the Frames at a time)
679
680             int l = Header->GetXSize() * Header->GetYSize() * Header->GetZSize();
681
682             uint8_t *newDest = new uint8_t[ImageDataSize];
683             uint8_t *x       = newDest;
684             uint8_t *a       = (uint8_t *)destination;
685             uint8_t *b       = a + l;
686             uint8_t *c       = b + l;
687
688             for (int j = 0; j < l; j++)
689             {
690                *(x++) = *(a++);
691                *(x++) = *(b++);
692                *(x++) = *(c++);
693             }
694             memmove(destination, newDest, ImageDataSize);
695             delete[] newDest;
696          }
697          break;
698       case 2:                      
699          // Palettes were found
700          // Let the user deal with them !
701          return ImageDataSize;
702    }
703    // now, it's an RGB image
704    // Lets's write it in the Header
705  
706    // Droping Palette Color out of the Header
707    // has been moved to the Write process.
708
709    // TODO : move 'values' modification to the write process
710    //      : save also (in order to be able to restore)
711    //      : 'high bit' -when not equal to 'bits stored' + 1
712    //      : 'bits allocated', when it's equal to 12 ?!
713
714    std::string spp = "3";            // Samples Per Pixel
715    std::string photInt = "RGB ";     // Photometric Interpretation
716    std::string planConfig = "0";     // Planar Configuration
717      
718    Header->SetEntryByNumber(spp,0x0028,0x0002);
719    Header->SetEntryByNumber(photInt,0x0028,0x0004);
720    Header->SetEntryByNumber(planConfig,0x0028,0x0006);
721  
722    return ImageDataSize; 
723 }
724
725 /**
726  * \brief   Points the internal Pixel_Data pointer to the callers inData
727  *          image representation, BUT WITHOUT COPYING THE DATA.
728  *          'image' Pixels are presented as C-like 2D arrays : line per line.
729  *          'volume'Pixels are presented as C-like 3D arrays : plane per plane 
730  * \warning Since the pixels are not copied, it is the caller's responsability
731  *          not to deallocate it's data before gdcm uses them (e.g. with
732  *          the Write() method.
733  * @param inData user supplied pixel area
734  * @param expectedSize total image size, in Bytes
735  *
736  * @return boolean
737  */
738 bool gdcmFile::SetImageData(uint8_t* inData, size_t expectedSize)
739 {
740    Header->SetImageDataSize( expectedSize );
741 // FIXME : if already allocated, memory leak !
742    Pixel_Data     = inData;
743    ImageDataSize = ImageDataSizeRaw = expectedSize;
744    PixelRead     = 1;
745 // FIXME : 7fe0, 0010 IS NOT set ...
746    return true;
747 }
748
749 /**
750  * \brief Writes on disk A SINGLE Dicom file
751  *        NO test is performed on  processor "Endiannity".
752  *        It's up to the user to call his Reader properly
753  * @param fileName name of the file to be created
754  *                 (any already existing file is over written)
755  * @return false if write fails
756  */
757
758 bool gdcmFile::WriteRawData(std::string const & fileName)
759 {
760    FILE* fp1 = fopen(fileName.c_str(), "wb");
761    if (fp1 == NULL)
762    {
763       printf("Fail to open (write) file [%s] \n", fileName.c_str());
764       return false;
765    }
766    fwrite (Pixel_Data, ImageDataSize, 1, fp1);
767    fclose (fp1);
768
769    return true;
770 }
771
772 /**
773  * \brief Writes on disk A SINGLE Dicom file, 
774  *        using the Implicit Value Representation convention
775  *        NO test is performed on  processor "Endiannity".
776  * @param fileName name of the file to be created
777  *                 (any already existing file is overwritten)
778  * @return false if write fails
779  */
780
781 bool gdcmFile::WriteDcmImplVR (std::string const & fileName)
782 {
783    return WriteBase(fileName, gdcmImplicitVR);
784 }
785
786 /**
787 * \brief Writes on disk A SINGLE Dicom file, 
788  *        using the Explicit Value Representation convention
789  *        NO test is performed on  processor "Endiannity". * @param fileName name of the file to be created
790  *                 (any already existing file is overwritten)
791  * @return false if write fails
792  */
793
794 bool gdcmFile::WriteDcmExplVR (std::string const & fileName)
795 {
796    return WriteBase(fileName, gdcmExplicitVR);
797 }
798
799 /**
800  * \brief Writes on disk A SINGLE Dicom file, 
801  *        using the ACR-NEMA convention
802  *        NO test is performed on  processor "Endiannity".
803  *        (a l'attention des logiciels cliniques 
804  *        qui ne prennent en entrée QUE des images ACR ...
805  * \warning if a DICOM_V3 header is supplied,
806  *         groups < 0x0008 and shadow groups are ignored
807  * \warning NO TEST is performed on processor "Endiannity".
808  * @param fileName name of the file to be created
809  *                 (any already existing file is overwritten)
810  * @return false if write fails
811  */
812
813 bool gdcmFile::WriteAcr (std::string const & fileName)
814 {
815    return WriteBase(fileName, gdcmACR);
816 }
817
818 //-----------------------------------------------------------------------------
819 // Protected
820 /**
821  * \brief NOT a end user inteded function
822  *        (used by WriteDcmExplVR, WriteDcmImplVR, WriteAcr, etc)
823  * @param fileName name of the file to be created
824  *                 (any already existing file is overwritten)
825  * @param  type file type (ExplicitVR, ImplicitVR, ...)
826  * @return false if write fails
827  */
828 bool gdcmFile::WriteBase (std::string const & fileName, FileType type)
829 {
830    if ( PixelRead == -1 && type != gdcmExplicitVR)
831    {
832       return false;
833    }
834
835    FILE* fp1 = fopen(fileName.c_str(), "wb");
836    if (fp1 == NULL)
837    {
838       printf("Failed to open (write) File [%s] \n", fileName.c_str());
839       return false;
840    }
841
842    if ( type == gdcmImplicitVR || type == gdcmExplicitVR )
843    {
844       // writing Dicom File Preamble
845       uint8_t* filePreamble = new uint8_t[128];
846       memset(filePreamble, 0, 128);
847       fwrite(filePreamble, 128, 1, fp1);
848       fwrite("DICM", 4, 1, fp1);
849
850       delete[] filePreamble;
851    }
852
853    // --------------------------------------------------------------
854    // Special Patch to allow gdcm to re-write ACR-LibIDO formated images
855    //
856    // if recognition code tells us we dealt with a LibIDO image
857    // we reproduce on disk the switch between lineNumber and columnNumber
858    // just before writting ...
859    
860    /// \todo the best trick would be *change* the recognition code
861    ///       but pb expected if user deals with, e.g. COMPLEX images
862
863    std::string rows, columns; 
864    if ( Header->GetFileType() == gdcmACR_LIBIDO)
865    {
866       rows    = Header->GetEntryByNumber(0x0028, 0x0010);
867       columns = Header->GetEntryByNumber(0x0028, 0x0011);
868
869       Header->SetEntryByNumber(columns,  0x0028, 0x0010);
870       Header->SetEntryByNumber(rows   ,  0x0028, 0x0011);
871    }
872    // ----------------- End of Special Patch ----------------
873       
874    uint16_t grPixel  = Header->GetGrPixel();
875    uint16_t numPixel = Header->GetNumPixel();;
876           
877    gdcmDocEntry* PixelElement = 
878       GetHeader()->GetDocEntryByNumber(grPixel, numPixel);  
879  
880    if ( PixelRead == 1 )
881    {
882       // we read pixel 'as is' (no tranformation LUT -> RGB)
883       PixelElement->SetLength( ImageDataSizeRaw );
884    }
885    else if ( PixelRead == 0 )
886    {
887       // we tranformed GrayLevel pixels + LUT into RGB Pixel
888       PixelElement->SetLength( ImageDataSize );
889    }
890  
891    Header->Write(fp1, type);
892
893    // --------------------------------------------------------------
894    // Special Patch to allow gdcm to re-write ACR-LibIDO formated images
895    // 
896    // ...and we restore the Header to be Dicom Compliant again 
897    // just after writting
898
899    if ( Header->GetFileType() == gdcmACR_LIBIDO )
900    {
901       Header->SetEntryByNumber(rows   , 0x0028, 0x0010);
902       Header->SetEntryByNumber(columns, 0x0028, 0x0011);
903    }
904    // ----------------- End of Special Patch ----------------
905    
906    // fwrite(Pixel_Data, ImageDataSize, 1, fp1);  // should be useless, now
907    fclose (fp1);
908
909    return true;
910 }
911
912 //-----------------------------------------------------------------------------
913 // Private
914 /**
915  * \brief   Swap the bytes, according to swap code.
916  * \warning not end user intended
917  * @param   im area to deal with
918  * @param   swap swap code
919  * @param   lgr Area Length
920  * @param   nb Pixels Bit number 
921  */
922 void gdcmFile::SwapZone(void* im, int swap, int lgr, int nb)
923 {
924    int i;
925
926    if( nb == 16 )
927    {
928       uint16_t* im16 = (uint16_t*)im;
929       switch( swap )
930       {
931          case 0:
932          case 12:
933          case 1234:
934             break;
935          case 21:
936          case 3412:
937          case 2143:
938          case 4321:
939             for(i=0; i < lgr/2; i++)
940             {
941                im16[i]= (im16[i] >> 8) | (im16[i] << 8 );
942             }
943             break;
944          default:
945             std::cout << "SWAP value (16 bits) not allowed :i" << swap << 
946             std::endl;
947       }
948    }
949    else if( nb == 32 )
950    {
951       uint32_t s32;
952       uint16_t fort, faible;
953       uint32_t* im32 = (uint32_t*)im;
954       switch ( swap )
955       {
956          case 0:
957          case 1234:
958             break;
959          case 4321:
960             for(i = 0; i < lgr/4; i++)
961             {
962                faible  = im32[i] & 0x0000ffff;  // 4321
963                fort    = im32[i] >> 16;
964                fort    = ( fort >> 8   ) | ( fort << 8 );
965                faible  = ( faible >> 8 ) | ( faible << 8);
966                s32     = faible;
967                im32[i] = ( s32 << 16 ) | fort;
968             }
969             break;
970          case 2143:
971             for(i = 0; i < lgr/4; i++)
972             {
973                faible  = im32[i] & 0x0000ffff;   // 2143
974                fort    = im32[i] >> 16;
975                fort    = ( fort >> 8 ) | ( fort << 8 );
976                faible  = ( faible >> 8) | ( faible << 8);
977                s32     = fort; 
978                im32[i] = ( s32 << 16 ) | faible;
979             }
980             break;
981          case 3412:
982             for(i = 0; i < lgr/4; i++)
983             {
984                faible  = im32[i] & 0x0000ffff; // 3412
985                fort    = im32[i] >> 16;
986                s32     = faible;
987                im32[i] = ( s32 << 16 ) | fort;
988             }
989             break;
990          default:
991             std::cout << "SWAP value (32 bits) not allowed : " << swap << 
992             std::endl;
993       }
994    }
995 }
996
997 /**
998  * \brief   Read pixel data from disk (optionaly decompressing) into the
999  *          caller specified memory location.
1000  * @param   destination where the pixel data should be stored.
1001  *
1002  */
1003 bool gdcmFile::ReadPixelData(void* destination) 
1004 {
1005    FILE* fp = Header->OpenFile();
1006
1007    if ( !fp )
1008    {
1009       return false;
1010    }
1011    if ( fseek(fp, Header->GetPixelOffset(), SEEK_SET) == -1 )
1012    {
1013       Header->CloseFile();
1014       return false;
1015    }
1016
1017    // ----------------------  Compacted File (12 Bits Per Pixel)
1018    // unpack 12 Bits pixels into 16 Bits pixels
1019    // 2 pixels 12bit =     [0xABCDEF]
1020    // 2 pixels 16bit = [0x0ABD] + [0x0FCE]
1021    
1022    if ( Header->GetBitsAllocated() == 12 )
1023    {
1024       int nbPixels = Header->GetXSize() * Header->GetYSize();
1025       uint8_t b0, b1, b2;
1026       
1027       uint16_t* pdestination = (uint16_t*)destination;    
1028       for(int p = 0; p < nbPixels; p += 2 )
1029       {
1030          fread(&b0,1,1,fp);
1031          fread(&b1,1,1,fp);
1032          fread(&b2,1,1,fp);      
1033
1034          //Two steps is necessary to please VC++
1035          *pdestination++ =  ((b0 >> 4) << 8) + ((b0 & 0x0f) << 4) + (b1 & 0x0f);
1036          //                     A                     B                 D
1037          *pdestination++ =  ((b2 & 0x0f) << 8) + ((b1 >> 4) << 4) + (b2 >> 4);
1038          //                     F                     C                 E
1039   
1040          // Troubles expected on Big-Endian processors ?
1041       }
1042
1043       Header->CloseFile();
1044       return true;
1045    }
1046
1047    // ----------------------  Uncompressed File
1048    if ( !Header->IsDicomV3()                             ||
1049         Header->IsImplicitVRLittleEndianTransferSyntax() ||
1050         Header->IsExplicitVRLittleEndianTransferSyntax() ||
1051         Header->IsExplicitVRBigEndianTransferSyntax()    ||
1052         Header->IsDeflatedExplicitVRLittleEndianTransferSyntax() )
1053    {
1054       size_t ItemRead = fread(destination, Header->GetPixelAreaLength(), 1, fp);
1055       Header->CloseFile();
1056       if ( ItemRead != 1 )
1057       {
1058          return false;
1059       }
1060       else
1061       {
1062          return true;
1063       }
1064    }
1065
1066    // ---------------------- Run Length Encoding
1067    if ( Header->IsRLELossLessTransferSyntax() )
1068    {
1069       bool res = gdcm_read_RLE_file (fp,destination);
1070       Header->CloseFile();
1071       return res; 
1072    }  
1073     
1074    // --------------- SingleFrame/Multiframe JPEG Lossless/Lossy/2000 
1075    int nb;
1076    std::string str_nb = Header->GetEntryByNumber(0x0028,0x0100);
1077    if ( str_nb == GDCM_UNFOUND )
1078    {
1079       nb = 16;
1080    }
1081    else
1082    {
1083       nb = atoi( str_nb.c_str() );
1084       if ( nb == 12 )
1085       {
1086          nb = 16;  // ?? 12 should be ACR-NEMA only
1087       }
1088    }
1089
1090    int nBytes= nb/8;
1091    int taille = Header->GetXSize() * Header->GetYSize()  
1092                 * Header->GetSamplesPerPixel();
1093    long fragmentBegining; // for ftell, fseek
1094
1095    bool jpg2000     = Header->IsJPEG2000();
1096    bool jpgLossless = Header->IsJPEGLossless();
1097
1098    bool res = true;
1099    uint16_t ItemTagGr, ItemTagEl;
1100    int ln;  
1101    
1102    //  Position on begining of Jpeg Pixels
1103    
1104    fread(&ItemTagGr,2,1,fp);  // Reading (fffe) : Item Tag Gr
1105    fread(&ItemTagEl,2,1,fp);  // Reading (e000) : Item Tag El
1106    if(Header->GetSwapCode())
1107    {
1108       ItemTagGr = Header->SwapShort(ItemTagGr);
1109       ItemTagEl = Header->SwapShort(ItemTagEl);
1110    }
1111
1112    fread(&ln,4,1,fp);
1113    if( Header->GetSwapCode() )
1114    {
1115       ln = Header->SwapLong( ln );    // Basic Offset Table Item length
1116    }
1117
1118    if ( ln != 0 )
1119    {
1120       // What is it used for ?!?
1121       uint8_t* BasicOffsetTableItemValue = new uint8_t[ln+1];
1122       fread(BasicOffsetTableItemValue,ln,1,fp);
1123       //delete[] BasicOffsetTableItemValue;
1124    }
1125
1126    // first Fragment initialisation
1127    fread(&ItemTagGr,2,1,fp);  // Reading (fffe) : Item Tag Gr
1128    fread(&ItemTagEl,2,1,fp);  // Reading (e000) : Item Tag El
1129    if( Header->GetSwapCode() )
1130    {
1131       ItemTagGr = Header->SwapShort( ItemTagGr );
1132       ItemTagEl = Header->SwapShort( ItemTagEl );
1133    }
1134
1135    // parsing fragments until Sequence Delim. Tag found
1136    while ( ItemTagGr == 0xfffe && ItemTagEl != 0xe0dd )
1137    {
1138       // --- for each Fragment
1139       fread(&ln,4,1,fp);
1140       if( Header->GetSwapCode() )
1141       {
1142          ln = Header->SwapLong(ln);    // Fragment Item length
1143       }
1144       fragmentBegining = ftell( fp );
1145
1146       if ( jpg2000 )
1147       {
1148       // JPEG 2000 :    call to ???
1149       res = gdcm_read_JPEG2000_file (fp,destination);  // Not Yet written 
1150       // ------------------------------------- endif (JPEG2000)
1151       }
1152       else if (jpgLossless)
1153       {
1154          // JPEG LossLess : call to xmedcom Lossless JPEG
1155          // Reading Fragment pixels
1156          JPEGLosslessDecodeImage (fp, (uint16_t*)destination,
1157                Header->GetPixelSize() * 8 * Header->GetSamplesPerPixel(), ln);
1158          res = 1; // in order not to break the loop
1159   
1160       } // ------------------------------------- endif (JPEGLossless)
1161       else
1162       {
1163          // JPEG Lossy : call to IJG 6b
1164          if ( Header->GetBitsStored() == 8)
1165          {
1166             // Reading Fragment pixels
1167             res = gdcm_read_JPEG_file (fp,destination);
1168          }
1169          else
1170          {
1171             // Reading Fragment pixels
1172             res = gdcm_read_JPEG_file12 (fp,destination);
1173          }
1174          // ------------------------------------- endif (JPEGLossy)
1175       }
1176
1177       if ( !res )
1178       {
1179          break;
1180       }
1181                
1182       // location in user's memory 
1183       // for next fragment (if any) 
1184       destination = (uint8_t*)destination + taille * nBytes;
1185
1186       fseek(fp,fragmentBegining, SEEK_SET); // To be sure we start 
1187       fseek(fp,ln,SEEK_CUR);                // at the begining of next fragment
1188       
1189       ItemTagGr = ItemTagEl = 0;
1190       fread(&ItemTagGr,2,1,fp);  // Reading (fffe) : Item Tag Gr
1191       fread(&ItemTagEl,2,1,fp);  // Reading (e000) : Item Tag El
1192       if( Header->GetSwapCode() )
1193       {
1194          ItemTagGr = Header->SwapShort( ItemTagGr );
1195          ItemTagEl = Header->SwapShort( ItemTagEl );
1196       }
1197    }
1198    // endWhile parsing fragments until Sequence Delim. Tag found    
1199
1200    Header->CloseFile();
1201    return res;
1202 }
1203 //-----------------------------------------------------------------------------