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