src/gdcmPixelConvert.cxx

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: gdcmPixelConvert.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2004/10/13 04:05:04 $
   7   Version:   $Revision: 1.10 $
   8
   9   Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
  10   l'Image). All rights reserved. See Doc/License.txt or
  11   http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
  12
  13      This software is distributed WITHOUT ANY WARRANTY; without even
  14      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  15      PURPOSE.  See the above copyright notices for more information.
  16
  17 =========================================================================*/
  18
  19 //////////////////   TEMPORARY NOTE
  20 // look for "fixMem" and convert that to a member of this class
  21 // Removing the prefix fixMem and dealing with allocations should do the trick
  22 //
  23 // grep PIXELCONVERT everywhere and clean up !
  24
  25 #include "gdcmDebug.h"
  26 #include "gdcmPixelConvert.h"
  27
  28 // External JPEG decompression
  29
  30 // for JPEGLosslessDecodeImage
  31 #include "jpeg/ljpg/jpegless.h"
  32
  33 namespace gdcm
  34 {
  35
  36 #define str2num(str, typeNum) *((typeNum *)(str))
  37
  38 // For JPEG 2000, body in file gdcmJpeg2000.cxx
  39 bool gdcm_read_JPEG2000_file (FILE* fp, void* image_buffer);
  40
  41 // For JPEG 8 Bits, body in file gdcmJpeg8.cxx
  42 bool gdcm_read_JPEG_file     (FILE* fp, void* image_buffer);
  43
  44 // For JPEG 12 Bits, body in file gdcmJpeg12.cxx
  45 bool gdcm_read_JPEG_file12   (FILE* fp, void* image_buffer);
  46
  47
  48
  49 //-----------------------------------------------------------------------------
  50 // Constructor / Destructor
  51 PixelConvert::PixelConvert()
  52 {
  53    RGB = 0;
  54    RGBSize = 0;
  55    Decompressed = 0;
  56    DecompressedSize = 0;
  57 }
  58
  59 void PixelConvert::Squeeze()
  60 {
  61    if ( RGB ) {
  62       delete [] RGB;
  63    }
  64    if ( Decompressed ) {
  65       delete [] Decompressed;
  66    }
  67 }
  68
  69 PixelConvert::~PixelConvert()
  70 {
  71    Squeeze();
  72 }
  73
  74 void PixelConvert::AllocateRGB()
  75 {
  76   if ( RGB ) {
  77      delete [] RGB;
  78   }
  79   RGB = new uint8_t[RGBSize];
  80 }
  81
  82 void PixelConvert::AllocateDecompressed()
  83 {
  84   if ( Decompressed ) {
  85      delete [] Decompressed;
  86   }
  87   Decompressed = new uint8_t[ DecompressedSize ];
  88 }
  89
  90 /**
  91  * \brief Read from file a 12 bits per pixel image and uncompress it
  92  *        into a 16 bits per pixel image.
  93  */
  94 void PixelConvert::Decompress12BitsTo16Bits(
  95                   uint8_t* pixelZone,
  96                   FILE* filePtr)
  97                throw ( FormatError )
  98 {
  99    int nbPixels = XSize * YSize;
 100    uint16_t* destination = (uint16_t*)pixelZone;
 101
 102    for( int p = 0; p < nbPixels; p += 2 )
 103    {
 104       uint8_t b0, b1, b2;
 105       size_t ItemRead;
 106
 107       ItemRead = fread( &b0, 1, 1, filePtr);
 108       if ( ItemRead != 1 )
 109       {
 110          throw FormatError( "File::Decompress12BitsTo16Bits()",
 111                                 "Unfound first block" );
 112       }
 113
 114       ItemRead = fread( &b1, 1, 1, filePtr);
 115       if ( ItemRead != 1 )
 116       {
 117          throw FormatError( "File::Decompress12BitsTo16Bits()",
 118                                 "Unfound second block" );
 119       }
 120
 121       ItemRead = fread( &b2, 1, 1, filePtr);
 122       if ( ItemRead != 1 )
 123       {
 124          throw FormatError( "File::Decompress12BitsTo16Bits()",
 125                                 "Unfound second block" );
 126       }
 127
 128       // Two steps are necessary to please VC++
 129       //
 130       // 2 pixels 12bit =     [0xABCDEF]
 131       // 2 pixels 16bit = [0x0ABD] + [0x0FCE]
 132       //                     A                     B                 D
 133       *destination++ =  ((b0 >> 4) << 8) + ((b0 & 0x0f) << 4) + (b1 & 0x0f);
 134       //                     F                     C                 E
 135       *destination++ =  ((b2 & 0x0f) << 8) + ((b1 >> 4) << 4) + (b2 >> 4);
 136
 137       /// \todo JPR Troubles expected on Big-Endian processors ?
 138    }
 139 }
 140
 141 /**
 142  * \brief     Try to deal with RLE 16 Bits.
 143  *            We assume the RLE has allready been parsed and loaded in
 144  *            Uncompressed (through \ref ReadAndDecompressJPEGFile ).
 145  *            We here need to make 16 Bits Pixels from Low Byte and
 146  *            High Byte 'Planes'...(for what it may mean)
 147  * @return    Boolean
 148  */
 149 bool PixelConvert::UncompressRLE16BitsFromRLE8Bits(
 150                        int NumberOfFrames,
 151                        uint8_t* fixMemUncompressed )
 152 {
 153    size_t PixelNumber = XSize * YSize;
 154    size_t fixMemUncompressedSize = XSize * YSize * NumberOfFrames;
 155
 156    // We assumed Uncompressed contains the decoded RLE pixels but as
 157    // 8 bits per pixel. In order to convert those pixels to 16 bits
 158    // per pixel we cannot work in place within Uncompressed and hence
 159    // we copy Uncompressed in a safe place, say OldUncompressed.
 160
 161    uint8_t* OldUncompressed = new uint8_t[ fixMemUncompressedSize * 2 ];
 162    memmove( OldUncompressed, fixMemUncompressed, fixMemUncompressedSize * 2);
 163
 164    uint8_t* x = fixMemUncompressed;
 165    uint8_t* a = OldUncompressed;
 166    uint8_t* b = a + PixelNumber;
 167
 168    for ( int i = 0; i < NumberOfFrames; i++ )
 169    {
 170       for ( unsigned int j = 0; j < PixelNumber; j++ )
 171       {
 172          *(x++) = *(a++);
 173          *(x++) = *(b++);
 174       }
 175    }
 176
 177    delete[] OldUncompressed;
 178
 179    /// \todo check that operator new []didn't fail, and sometimes return false
 180    return true;
 181 }
 182
 183 /**
 184  * \brief Implementation of the RLE decoding algorithm for uncompressing
 185  *        a RLE fragment. [refer to PS 3.5-2003, section G.3.2 p 86]
 186  */
 187 bool PixelConvert::ReadAndUncompressRLEFragment( uint8_t* decodedZone,
 188                                                long fragmentSize,
 189                                                long uncompressedSegmentSize,
 190                                                FILE* fp )
 191 {
 192    int8_t count;
 193    long numberOfOutputBytes = 0;
 194    long numberOfReadBytes = 0;
 195
 196    while( numberOfOutputBytes < uncompressedSegmentSize )
 197    {
 198       fread( &count, 1, 1, fp );
 199       numberOfReadBytes += 1;
 200       if ( count >= 0 )
 201       // Note: count <= 127 comparison is always true due to limited range
 202       //       of data type int8_t [since the maximum of an exact width
 203       //       signed integer of width N is 2^(N-1) - 1, which for int8_t
 204       //       is 127].
 205       {
 206          fread( decodedZone, count + 1, 1, fp);
 207          numberOfReadBytes += count + 1;
 208          decodedZone         += count + 1;
 209          numberOfOutputBytes += count + 1;
 210       }
 211       else
 212       {
 213          if ( ( count <= -1 ) && ( count >= -127 ) )
 214          {
 215             int8_t newByte;
 216             fread( &newByte, 1, 1, fp);
 217             numberOfReadBytes += 1;
 218             for( int i = 0; i < -count + 1; i++ )
 219             {
 220                decodedZone[i] = newByte;
 221             }
 222             decodedZone         += -count + 1;
 223             numberOfOutputBytes += -count + 1;
 224          }
 225       }
 226       // if count = 128 output nothing
 227
 228       if ( numberOfReadBytes > fragmentSize )
 229       {
 230          dbg.Verbose(0, "File::gdcm_read_RLE_fragment: we read more "
 231                         "bytes than the segment size.");
 232          return false;
 233       }
 234    }
 235    return true;
 236 }
 237
 238 /**
 239  * \brief     Reads from disk the Pixel Data of 'Run Length Encoded'
 240  *            Dicom encapsulated file and uncompress it.
 241  * @param     fp already open File Pointer
 242  * @param     image_buffer destination Address (in caller's memory space)
 243  *            at which the pixel data should be copied
 244  * @return    Boolean
 245  */
 246 bool PixelConvert::ReadAndDecompressRLEFile(
 247                           void* image_buffer,
 248                           FILE* fp )
 249 {
 250    uint8_t* im = (uint8_t*)image_buffer;
 251    long uncompressedSegmentSize = XSize * YSize;
 252
 253
 254    // Loop on the frame[s]
 255    for( RLEFramesInfo::RLEFrameList::iterator
 256         it  = RLEInfo->Frames.begin();
 257         it != RLEInfo->Frames.end();
 258       ++it )
 259    {
 260       // Loop on the fragments
 261       for( int k = 1; k <= (*it)->NumberFragments; k++ )
 262       {
 263          fseek( fp, (*it)->Offset[k] ,SEEK_SET );
 264          (void)PixelConvert::ReadAndUncompressRLEFragment(
 265                                  (uint8_t*) im, (*it)->Length[k],
 266                                  uncompressedSegmentSize, fp );
 267          im += uncompressedSegmentSize;
 268       }
 269    }
 270
 271    if ( BitsAllocated == 16 )
 272    {
 273       // Try to deal with RLE 16 Bits
 274       (void)UncompressRLE16BitsFromRLE8Bits( ZSize,
 275                                              (uint8_t*) image_buffer);
 276    }
 277
 278    return true;
 279 }
 280
 281 /**
 282  * \brief   Swap the bytes, according to swap code.
 283  * \warning not end user intended
 284  * @param   im area to deal with
 285  */
 286 void PixelConvert::SwapZone( uint8_t* im )
 287 {
 288    unsigned int i;
 289
 290    if( BitsAllocated == 16 )
 291    {
 292       uint16_t* im16 = (uint16_t*)im;
 293       switch( SwapCode )
 294       {
 295          case 0:
 296          case 12:
 297          case 1234:
 298             break;
 299          case 21:
 300          case 3412:
 301          case 2143:
 302          case 4321:
 303             for( i = 0; i < DecompressedSize / 2; i++ )
 304             {
 305                im16[i]= (im16[i] >> 8) | (im16[i] << 8 );
 306             }
 307             break;
 308          default:
 309             dbg.Verbose( 0, "PixelConvert::SwapZone: SwapCode value "
 310                             "(16 bits) not allowed." );
 311       }
 312    }
 313    else if( BitsAllocated == 32 )
 314    {
 315       uint32_t s32;
 316       uint16_t fort, faible;
 317       uint32_t* im32 = (uint32_t*)im;
 318       switch ( SwapCode )
 319       {
 320          case 0:
 321          case 1234:
 322             break;
 323          case 4321:
 324             for( i = 0; i < DecompressedSize / 4; i++ )
 325             {
 326                faible  = im32[i] & 0x0000ffff;  // 4321
 327                fort    = im32[i] >> 16;
 328                fort    = ( fort >> 8   ) | ( fort << 8 );
 329                faible  = ( faible >> 8 ) | ( faible << 8);
 330                s32     = faible;
 331                im32[i] = ( s32 << 16 ) | fort;
 332             }
 333             break;
 334          case 2143:
 335             for( i = 0; i < DecompressedSize / 4; i++ )
 336             {
 337                faible  = im32[i] & 0x0000ffff;   // 2143
 338                fort    = im32[i] >> 16;
 339                fort    = ( fort >> 8 ) | ( fort << 8 );
 340                faible  = ( faible >> 8) | ( faible << 8);
 341                s32     = fort;
 342                im32[i] = ( s32 << 16 ) | faible;
 343             }
 344             break;
 345          case 3412:
 346             for( i = 0; i < DecompressedSize / 4; i++ )
 347             {
 348                faible  = im32[i] & 0x0000ffff; // 3412
 349                fort    = im32[i] >> 16;
 350                s32     = faible;
 351                im32[i] = ( s32 << 16 ) | fort;
 352             }
 353             break;
 354          default:
 355             dbg.Verbose( 0, "PixelConvert::SwapZone: SwapCode value "
 356                             "(32 bits) not allowed." );
 357       }
 358    }
 359 }
 360
 361 /**
 362  * \brief Deal with endianity i.e. re-arange bytes inside the integer
 363  */
 364 void PixelConvert::ReorderEndianity( uint8_t* pixelZone )
 365 {
 366    if ( BitsAllocated != 8 )
 367    {
 368       SwapZone( pixelZone );
 369    }
 370
 371    // Special kludge in order to deal with xmedcon broken images:
 372    if (  ( BitsAllocated == 16 )
 373        && ( BitsStored < BitsAllocated )
 374        && ( ! PixelSign ) )
 375    {
 376       int l = (int)( DecompressedSize / ( BitsAllocated / 8 ) );
 377       uint16_t *deb = (uint16_t *)pixelZone;
 378       for(int i = 0; i<l; i++)
 379       {
 380          if( *deb == 0xffff )
 381          {
 382            *deb = 0;
 383          }
 384          deb++;
 385       }
 386    }
 387 }
 388
 389 /**
 390  * \brief     Reads from disk the Pixel Data of JPEG Dicom encapsulated
 391  &            file and uncompress it.
 392  * @param     fp already open File Pointer
 393  * @param     destination Where decompressed fragments should end up
 394  * @return    Boolean
 395  */
 396 bool PixelConvert::ReadAndDecompressJPEGFile(
 397                           uint8_t* destination,
 398                           FILE* fp )
 399 {
 400    // Loop on the fragment[s]
 401    for( JPEGFragmentsInfo::JPEGFragmentsList::iterator
 402         it  = JPEGInfo->Fragments.begin();
 403         it != JPEGInfo->Fragments.end();
 404       ++it )
 405    {
 406       fseek( fp, (*it)->Offset, SEEK_SET );
 407
 408       if ( IsJPEG2000 )
 409       {
 410          if ( ! gdcm_read_JPEG2000_file( fp, destination ) )
 411          {
 412             return false;
 413          }
 414       }
 415       else if ( IsJPEGLossless )
 416       {
 417          // JPEG LossLess : call to xmedcom Lossless JPEG
 418          JPEGLosslessDecodeImage( fp,
 419                                   (uint16_t*)destination,
 420                                   PixelSize * 8 * SamplesPerPixel,
 421                                   (*it)->Length );
 422       }
 423       else if ( BitsStored == 8)
 424       {
 425          // JPEG Lossy : call to IJG 6b
 426          if ( ! gdcm_read_JPEG_file ( fp, destination ) )
 427          {
 428             return false;
 429          }
 430       }
 431       else if ( BitsStored == 12)
 432       {
 433          // Reading Fragment pixels
 434          if ( ! gdcm_read_JPEG_file12 ( fp, destination ) )
 435          {
 436             return false;
 437          }
 438       }
 439       else
 440       {
 441          // other JPEG lossy not supported
 442          dbg.Error(" File::ReadAndDecompressJPEGFile: unknown jpeg lossy "
 443                    " compression ");
 444          return false;
 445       }
 446
 447       // Advance to next free location in destination
 448       // for next fragment decompression (if any)
 449       int length = XSize * YSize * SamplesPerPixel;
 450       int numberBytes = BitsAllocated / 8;
 451
 452       destination = (uint8_t*)destination + length * numberBytes;
 453
 454    }
 455    return true;
 456 }
 457
 458 /**
 459  * \brief  Re-arrange the bits within the bytes.
 460  * @param  fp already open File Pointer
 461  * @param  destination Where decompressed fragments should end up
 462  * @return Boolean
 463  */
 464 bool PixelConvert::ReArrangeBits( uint8_t* pixelZone )
 465      throw ( FormatError )
 466 {
 467    if ( BitsStored != BitsAllocated )
 468    {
 469       int l = (int)( DecompressedSize / ( BitsAllocated / 8 ) );
 470       if ( BitsAllocated == 16 )
 471       {
 472          uint16_t mask = 0xffff;
 473          mask = mask >> ( BitsAllocated - BitsStored );
 474          uint16_t* deb = (uint16_t*)pixelZone;
 475          for(int i = 0; i<l; i++)
 476          {
 477             *deb = (*deb >> (BitsStored - HighBitPosition - 1)) & mask;
 478             deb++;
 479          }
 480       }
 481       else if ( BitsAllocated == 32 )
 482       {
 483          uint32_t mask = 0xffffffff;
 484          mask = mask >> ( BitsAllocated - BitsStored );
 485          uint32_t* deb = (uint32_t*)pixelZone;
 486          for(int i = 0; i<l; i++)
 487          {
 488             *deb = (*deb >> (BitsStored - HighBitPosition - 1)) & mask;
 489             deb++;
 490          }
 491       }
 492       else
 493       {
 494          dbg.Verbose(0, "PixelConvert::ReArrangeBits: weird image");
 495          throw FormatError( "File::ReArrangeBits()",
 496                                 "weird image !?" );
 497       }
 498    }
 499    return true; //???
 500 }
 501
 502 /**
 503  * \brief   Convert (Y plane, cB plane, cR plane) to RGB pixels
 504  * \warning Works on all the frames at a time
 505  */
 506 void PixelConvert::ConvertYcBcRPlanesToRGBPixels(
 507                            uint8_t* destination,
 508                            size_t imageDataSize )
 509 {
 510    uint8_t* oldPixelZone = new uint8_t[ imageDataSize ];
 511    memmove( oldPixelZone, destination, imageDataSize );
 512
 513    // to see the tricks about YBR_FULL, YBR_FULL_422,
 514    // YBR_PARTIAL_422, YBR_ICT, YBR_RCT have a look at :
 515    // ftp://medical.nema.org/medical/dicom/final/sup61_ft.pdf
 516    // and be *very* affraid
 517    //
 518    int l        = XSize * YSize;
 519    int nbFrames = ZSize;
 520
 521    uint8_t* a = oldPixelZone;
 522    uint8_t* b = oldPixelZone + l;
 523    uint8_t* c = oldPixelZone + l + l;
 524    double R, G, B;
 525
 526    /// \todo : Replace by the 'well known' integer computation
 527    ///         counterpart. Refer to
 528    ///            http://lestourtereaux.free.fr/papers/data/yuvrgb.pdf
 529    ///         for code optimisation.
 530
 531    for ( int i = 0; i < nbFrames; i++ )
 532    {
 533       for ( int j = 0; j < l; j++ )
 534       {
 535          R = 1.164 *(*a-16) + 1.596 *(*c -128) + 0.5;
 536          G = 1.164 *(*a-16) - 0.813 *(*c -128) - 0.392 *(*b -128) + 0.5;
 537          B = 1.164 *(*a-16) + 2.017 *(*b -128) + 0.5;
 538
 539          if (R < 0.0)   R = 0.0;
 540          if (G < 0.0)   G = 0.0;
 541          if (B < 0.0)   B = 0.0;
 542          if (R > 255.0) R = 255.0;
 543          if (G > 255.0) G = 255.0;
 544          if (B > 255.0) B = 255.0;
 545
 546          *(destination++) = (uint8_t)R;
 547          *(destination++) = (uint8_t)G;
 548          *(destination++) = (uint8_t)B;
 549          a++;
 550          b++;
 551          c++;
 552       }
 553    }
 554    delete[] oldPixelZone;
 555 }
 556
 557 /**
 558  * \brief   Convert (Red plane, Green plane, Blue plane) to RGB pixels
 559  * \warning Works on all the frames at a time
 560  */
 561 void PixelConvert::ConvertRGBPlanesToRGBPixels(
 562                           uint8_t* destination,
 563                           size_t imageDataSize )
 564 {
 565    uint8_t* oldPixelZone = new uint8_t[ imageDataSize ];
 566    memmove( oldPixelZone, destination, imageDataSize );
 567
 568    int l = XSize * YSize * ZSize;
 569
 570    uint8_t* a = oldPixelZone;
 571    uint8_t* b = oldPixelZone + l;
 572    uint8_t* c = oldPixelZone + l + l;
 573
 574    for (int j = 0; j < l; j++)
 575    {
 576       *(destination++) = *(a++);
 577       *(destination++) = *(b++);
 578       *(destination++) = *(c++);
 579    }
 580    delete[] oldPixelZone;
 581 }
 582
 583 bool PixelConvert::ReadAndDecompressPixelData( void* destination, FILE* fp )
 584 {
 585    if ( !fp )
 586    {
 587       return false;
 588    }
 589
 590    if ( fseek(fp, PixelOffset, SEEK_SET) == -1 )
 591    {
 592       return false;
 593    }
 594
 595    if ( BitsAllocated == 12 )
 596    {
 597       Decompress12BitsTo16Bits( (uint8_t*)destination, fp);
 598       return true;
 599    }
 600
 601    //////////// Decompressed File
 602    if ( IsUncompressed )
 603    {
 604       size_t ItemRead = fread( destination, PixelDataLength, 1, fp);
 605       if ( ItemRead != 1 )
 606       {
 607          return false;
 608       }
 609       else
 610       {
 611          return true;
 612       }
 613    }
 614
 615    ///////////// Run Length Encoding
 616    if ( IsRLELossless )
 617    {
 618       return ReadAndDecompressRLEFile( destination, fp );
 619    }
 620
 621    ///////////// SingleFrame/Multiframe JPEG Lossless/Lossy/2000
 622    return ReadAndDecompressJPEGFile( (uint8_t*)destination, fp );
 623 }
 624
 625 void PixelConvert::ComputeDecompressedImageDataSize()
 626 {
 627    int bitsAllocated;
 628    // Number of "Bits Allocated" is fixed to 16 when it's 12, since
 629    // in this case we will expand the image to 16 bits (see
 630    //    \ref Decompress12BitsTo16Bits() )
 631    if (  BitsAllocated == 12 )
 632    {
 633       bitsAllocated = 16;
 634    }
 635
 636    DecompressedSize = XSize * YSize * ZSize
 637                     * ( bitsAllocated / 8 )
 638                     * SamplesPerPixel;
 639 }
 640
 641 } // end namespace gdcm