src/gdcmFile.cxx

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