gdcmpixelreadconvert.cxx

DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.

 
   return true;
}

//-----------------------------------------------------------------------------
// Protected

//-----------------------------------------------------------------------------
// Private
/**
 * \brief Read from file a 12 bits per pixel image and decompress it
 *        into a 16 bits per pixel image.
 */
void PixelReadConvert::ReadAndDecompress12BitsTo16Bits( std::ifstream *fp )
               throw ( FormatError )
{
   int nbPixels = XSize * YSize;
   uint16_t *localDecompres = (uint16_t*)Raw;

   for( int p = 0; p < nbPixels; p += 2 )
   {
      uint8_t b0, b1, b2;

      fp->read( (char*)&b0, 1);
      if ( fp->fail() || fp->eof() )
      {
         throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
                                "Unfound first block" );
      }

      fp->read( (char*)&b1, 1 );
      if ( fp->fail() || fp->eof())
      {
         throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
                                "Unfound second block" );
      }

      fp->read( (char*)&b2, 1 );
      if ( fp->fail() || fp->eof())
      {
         throw FormatError( "PixelReadConvert::ReadAndDecompress12BitsTo16Bits()",
                                "Unfound second block" );
      }

      // Two steps are necessary to please VC++
      //
      // 2 pixels 12bit =     [0xABCDEF]
      // 2 pixels 16bit = [0x0ABD] + [0x0FCE]
      //                        A                     B                 D
      *localDecompres++ =  ((b0 >> 4) << 8) + ((b0 & 0x0f) << 4) + (b1 & 0x0f);
      //                        F                     C                 E
      *localDecompres++ =  ((b2 & 0x0f) << 8) + ((b1 >> 4) << 4) + (b2 >> 4);

      /// \todo JPR Troubles expected on Big-Endian processors ?
   }
}

/**
 * \brief     Reads from disk the Pixel Data of JPEG Dicom encapsulated
 *            file and decompress it.
 * @param     fp File Pointer
 * @return    Boolean
 */
bool PixelReadConvert::ReadAndDecompressJPEGFile( std::ifstream *fp )
{
   if ( IsJPEG2000 )
   {
     // make sure this is the right JPEG compression
     assert( !IsJPEGLossless || !IsJPEGLossy || !IsJPEGLS );
     // FIXME this is really ugly but it seems I have to load the complete
     // jpeg2000 stream to use jasper:
     // I don't think we'll ever be able to deal with multiple fragments properly

     if( ZSize == 1 )
       {
       unsigned long inputlength = 0;
       JPEGFragment *jpegfrag = JPEGInfo->GetFirstFragment();
       while( jpegfrag )
         {
         inputlength += jpegfrag->GetLength();
         jpegfrag = JPEGInfo->GetNextFragment();
         }
       gdcmAssertMacro( inputlength != 0);
       uint8_t *inputdata = new uint8_t[inputlength];
       char *pinputdata = (char*)inputdata;
       jpegfrag = JPEGInfo->GetFirstFragment();
       while( jpegfrag )
         {
         fp->seekg( jpegfrag->GetOffset(), std::ios::beg);
         fp->read(pinputdata, jpegfrag->GetLength());
         pinputdata += jpegfrag->GetLength();
         jpegfrag = JPEGInfo->GetNextFragment();
         }
       // Warning the inputdata buffer is delete in the function
       if ( gdcm_read_JPEG2000_file( Raw,
           (char*)inputdata, inputlength ) )
         {
         return true;
         }
       // wow what happen, must be an error
       gdcmWarningMacro( "gdcm_read_JPEG2000_file() failed "); 
       return false;
       }
     else
       {
       if( (unsigned int)ZSize != JPEGInfo->GetFragmentCount() )
         {
         gdcmErrorMacro( "Sorry GDCM does not handle this type of fragments" );
         return false;
         }
       // Hopefully every dicom fragment is *exactly* the j2k stream
       JPEGFragment *jpegfrag = JPEGInfo->GetFirstFragment();
       char *praw = (char*)Raw;
       while( jpegfrag )
         {
         unsigned long inputlength = jpegfrag->GetLength();
         char *inputdata = new char[inputlength];
         fp->seekg( jpegfrag->GetOffset(), std::ios::beg);
         fp->read(inputdata, jpegfrag->GetLength());
         // Warning the inputdata buffer is delete in the function
         gdcm_read_JPEG2000_file( praw, 
           inputdata, inputlength) ;
         praw += XSize*YSize*SamplesPerPixel*(BitsAllocated/8);
         jpegfrag = JPEGInfo->GetNextFragment();
         }
       return true;
       }
   }
   else if ( IsJPEGLS )
   {
     // make sure this is the right JPEG compression
     assert( !IsJPEGLossless || !IsJPEGLossy || !IsJPEG2000 );
   // WARNING : JPEG-LS is NOT the 'classical' Jpeg Lossless : 
   // [JPEG-LS is the basis for new lossless/near-lossless compression
   // standard for continuous-tone images intended for JPEG2000. The standard
   // is based on the LOCO-I algorithm (LOw COmplexity LOssless COmpression
   // for Images) developed at Hewlett-Packard Laboratories]
   //
   // see http://datacompression.info/JPEGLS.shtml
   //
#if 0
   std::cerr << "count:" << JPEGInfo->GetFragmentCount() << std::endl;
      unsigned long inputlength = 0;
      JPEGFragment *jpegfrag = JPEGInfo->GetFirstFragment();
      while( jpegfrag )
      {
         inputlength += jpegfrag->GetLength();
         jpegfrag = JPEGInfo->GetNextFragment();
      }
      gdcmAssertMacro( inputlength != 0);
      uint8_t *inputdata = new uint8_t[inputlength];
      char *pinputdata = (char*)inputdata;
      jpegfrag = JPEGInfo->GetFirstFragment();
      while( jpegfrag )
      {
         fp->seekg( jpegfrag->GetOffset(), std::ios::beg);
         fp->read(pinputdata, jpegfrag->GetLength());
         pinputdata += jpegfrag->GetLength();
         jpegfrag = JPEGInfo->GetNextFragment();
      }  
      
  //fp->read((char*)Raw, PixelDataLength);

  std::ofstream out("/tmp/jpegls.jpg");
  out.write((char*)inputdata, inputlength);
  out.close();
  delete[] inputdata;
#endif

      gdcmWarningMacro( "Sorry, JPEG-LS not yet taken into account" );
      fp->seekg( JPEGInfo->GetFirstFragment()->GetOffset(), std::ios::beg);
//    if ( ! gdcm_read_JPEGLS_file( fp,Raw ) )
         return false;
   }
   else if( JPEGInfo )
   {
     // make sure this is the right JPEG compression
     assert( !IsJPEGLS || !IsJPEG2000 );
     // Precompute the offset localRaw will be shifted with
     int length = XSize * YSize * SamplesPerPixel;
     int numberBytes = BitsAllocated / 8;

     JPEGInfo->DecompressFromFile(fp, Raw, BitsStored, numberBytes, length );
     return true;
   }
   //else (not sure how get there...), must be one of those crazy DICOM file
   return false;
}

/**
 * \brief Build Red/Green/Blue/Alpha LUT from File when :
 *         - (0028,0004) : Photometric Interpretation == [PALETTE COLOR ]
 *         and
 *         - (0028,1101),(0028,1102),(0028,1102)
 *            xxx Palette Color Lookup Table Descriptor are found
 *          and
 *         - (0028,1201),(0028,1202),(0028,1202)
 *           xxx Palette Color Lookup Table Data - are found
 * \warning does NOT deal with :
 *   - 0028 1100 Gray Lookup Table Descriptor (Retired)
 *   - 0028 1221 Segmented Red Palette Color Lookup Table Data
 *   - 0028 1222 Segmented Green Palette Color Lookup Table Data
 *   - 0028 1223 Segmented Blue Palette Color Lookup Table Data
 *   no known Dicom reader deals with them :-(
 * @return a RGBA Lookup Table
 */
void PixelReadConvert::BuildLUTRGBA()
{

   // Note to code reviewers :
   // The problem is *much more* complicated, since a lot of manufacturers
   // Don't follow the norm :
   // have a look at David Clunie's remark at the end of this .cxx file.
   if ( LutRGBA )
   
   {
      return;
   }
   // Not so easy : see
   // http://www.barre.nom.fr/medical/dicom2/limitations.html#Color%20Lookup%20Tables
                                                                                
   if ( ! IsPaletteColor )
   {
      return;
   }
                                                                                
   if (   LutRedDescriptor   == GDCM_UNFOUND
       || LutGreenDescriptor == GDCM_UNFOUND
       || LutBlueDescriptor  == GDCM_UNFOUND )
   {
      gdcmWarningMacro( "(At least) a LUT Descriptor is missing" );
      return;
   }

   ////////////////////////////////////////////
   // Extract the info from the LUT descriptors
   int lengthR;   // Red LUT length in Bytes
   int debR;      // Subscript of the first Lut Value
   int nbitsR;    // Lut item size (in Bits)
   int nbRead;    // nb of items in LUT descriptor (must be = 3)

   nbRead = sscanf( LutRedDescriptor.c_str(),
                        "%d\\%d\\%d",
                        &lengthR, &debR, &nbitsR );
   if ( nbRead != 3 )
   {
      gdcmWarningMacro( "Wrong Red LUT descriptor" );
   }                                                                                
   int lengthG;  // Green LUT length in Bytes
   int debG;     // Subscript of the first Lut Value
   int nbitsG;   // Lut item size (in Bits)

   nbRead = sscanf( LutGreenDescriptor.c_str(),
                    "%d\\%d\\%d",
                    &lengthG, &debG, &nbitsG );  
   if ( nbRead != 3 )
   {
      gdcmWarningMacro( "Wrong Green LUT descriptor" );
   }
                                                                                
   int lengthB;  // Blue LUT length in Bytes
   int debB;     // Subscript of the first Lut Value
   int nbitsB;   // Lut item size (in Bits)
   nbRead = sscanf( LutRedDescriptor.c_str(),
                    "%d\\%d\\%d",
                    &lengthB, &debB, &nbitsB );
   if ( nbRead != 3 )
   {
      gdcmWarningMacro( "Wrong Blue LUT descriptor" );
   }
 
   gdcmDebugMacro(" lengthR " << lengthR << " debR " 
                << debR << " nbitsR " << nbitsR);
   gdcmDebugMacro(" lengthG " << lengthG << " debG " 
                << debG << " nbitsG " << nbitsG);
   gdcmDebugMacro(" lengthB " << lengthB << " debB " 
                << debB << " nbitsB " << nbitsB);

   if ( !lengthR ) // if = 2^16, this shall be 0 see : CP-143
      lengthR=65536;
   if ( !lengthG ) // if = 2^16, this shall be 0
      lengthG=65536;
   if ( !lengthB ) // if = 2^16, this shall be 0
      lengthB=65536; 
                                                                                
   ////////////////////////////////////////////////////////

   if ( ( ! LutRedData ) || ( ! LutGreenData ) || ( ! LutBlueData ) )
   {
      gdcmWarningMacro( "(At least) a LUT is missing" );
      return;
   }

   // -------------------------------------------------------------
   
   if ( BitsAllocated <= 8 )
   {
      // forge the 4 * 8 Bits Red/Green/Blue/Alpha LUT
      LutRGBA = new uint8_t[ 1024 ]; // 256 * 4 (R, G, B, Alpha)
      if ( !LutRGBA )
         return;
      LutItemNumber = 256;
      LutItemSize   = 8;
      memset( LutRGBA, 0, 1024 );
                                                                                
      int mult;
      if ( ( nbitsR == 16 ) && ( BitsAllocated == 8 ) )
      {
         // when LUT item size is different than pixel size
         mult = 2; // high byte must be = low byte
      }
      else
      {
         // See PS 3.3-2003 C.11.1.1.2 p 619
         mult = 1;
      }
                                                                                
      // if we get a black image, let's just remove the '+1'
      // from 'i*mult+1' and check again
      // if it works, we shall have to check the 3 Palettes
      // to see which byte is ==0 (first one, or second one)
      // and fix the code
      // We give up the checking to avoid some (useless ?) overhead
      // (optimistic asumption)
      int i;
      uint8_t *a;

      //take "Subscript of the first Lut Value" (debR,debG,debB) into account!

      //FIXME :  +1 : to get 'low value' byte
      //         Trouble expected on Big Endian Processors ?
      //         16 BIts Per Pixel Palette Color to be swapped?

      a = LutRGBA + 0 + debR;
      for( i=0; i < lengthR; ++i )
      {
         *a = LutRedData[i*mult+1]; 
         a += 4;
      }
                                                                                
      a = LutRGBA + 1 + debG;
      for( i=0; i < lengthG; ++i)
      {
         *a = LutGreenData[i*mult+1];
         a += 4;
      }
                                                                                
      a = LutRGBA + 2 + debB;
      for(i=0; i < lengthB; ++i)
      {
         *a = LutBlueData[i*mult+1];
         a += 4;
      }
                                    
      a = LutRGBA + 3 ;
      for(i=0; i < 256; ++i)
      {
         *a = 1; // Alpha component
         a += 4;
      }
   }
   else
   {
      // Probabely the same stuff is to be done for 16 Bits Pixels
      // with 65536 entries LUT ?!?
      // Still looking for accurate info on the web :-(

      gdcmWarningMacro( "Sorry Palette Color Lookup Tables not yet dealt with"
                         << " for 16 Bits Per Pixel images" );

      // forge the 4 * 16 Bits Red/Green/Blue/Alpha LUT

      LutRGBA = (uint8_t *)new uint16_t[ 65536*4 ]; // 2^16 * 4 (R, G, B, Alpha)
      if ( !LutRGBA )
         return;
      memset( LutRGBA, 0, 65536*4*2 );  // 16 bits = 2 bytes ;-)

      LutItemNumber = 65536;
      LutItemSize   = 16;

      int i;
      uint16_t *a16;

      //take "Subscript of the first Lut Value" (debR,debG,debB) into account!