src/gdcmFile.cxx

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