src/gdcmFile.cxx

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