spiht.c

该源码是SPIHT即分层树集合划分算法的c代码,希望对搞压缩编码的朋友们有用.

/*
 *
 * QccPack: Quantization, compression, and coding utilities
 * Copyright (C) 1997-2001  James E. Fowler
 * 
 */

/*
 * ----------------------------------------------------------------------------
 * 
 * Public License for the SPIHT Algorithm
 * Version 1.0, September 22, 1999
 * 
 * ----------------------------------------------------------------------------
 * 
 * The Set Partitioning In Hierarchical Trees (SPIHT) algorithm is protected
 * by US Patent #5,764,807 (issued June 9, 1998) and other international
 * patents and patents pending. An implementation of the SPIHT algorithm is
 * included herein with the permission of PrimaComp, Inc., exclusive holder
 * of patent rights. PrimaComp has graciously granted the following license
 * governing the terms and conditions for use, copying, distribution, and
 * modification of the SPIHT algorithm implementation contained herein
 * (hereafter referred to as "the SPIHT source code").
 * 
 * 0. Use of the SPIHT source code, including any executable-program or
 *    linkable-library form resulting from its compilation, is restricted to
 *    solely academic or non-commercial research activities.
 * 
 * 1. Any other use, including, but not limited to, use in the development
 *    of a commercial product, use in a commercial application, or commercial
 *    distribution, is prohibited by this license. Such acts require a separate
 *    license directly from PrimaComp, Inc.
 * 
 * 2. For academic and non-commercial purposes, this license does not restrict
 *    use; copying, distribution, and modification are permitted under the
 *    terms of the GNU General Public License as published by the Free Software
 *    Foundation, with the further restriction that the terms of the present
 *    license shall also apply to all subsequent copies, distributions,
 *    or modifications of the SPIHT source code.
 * 
 * NO WARRANTY
 * 
 * 3. PrimaComp dislaims all warranties, expressed or implied, including
 *    without limitation any warranty whatsoever as to the fitness for a
 *    particular use or the merchantability of the SPIHT source code.
 * 
 * 4. In no event shall PrimaComp be liable for any loss of profits, loss of 
 *    business, loss of use or loss of data, nor for indirect, special,
 *    incidental or consequential damages of any kind related to use of the
 *    SPIHT source code.
 * 
 * 
 * END OF TERMS AND CONDITIONS
 * 
 * 
 * Persons desiring to license the SPIHT algorithm for commercial purposes or
 * for uses otherwise prohibited by this license may wish to contact
 * PrimaComp, Inc., regarding the possibility of negotiating such licenses:
 * 
 *   PrimaComp, Inc.
 *   1223 Peoples Ave.
 *   Troy, New York  12180
 *   USA
 *   [phone]
 *   [e-mail]
 *   http://www.primacomp.com
 * 
 * ----------------------------------------------------------------------------
 */


#include "libQccPack.h"


void QccSPIHTHeader(void)
{
  QccString header;

  QccStringSprintf(header,
		   "%s\n\nQccPackSPIHT Module Version %d.%d %s,\n%s\n\nfor:",
		   QCCSPIHT_LICENSE,
		   QCCSPIHT_MAJORVERSION,
		   QCCSPIHT_MINORVERSION,
		   QCCSPIHT_DATE,
		   QCCSPIHT_COPYRIGHT);

  QccSetUserHeader(header);
}


#define QCCSPIHT_LASTBIT 4

static const int QccSPIHTArithmeticContexts[] =
{
  /*  LIP Context  */
  16, 8, 8, 4, 8, 4, 4, 2, 8, 4, 4, 2, 4, 2, 2,
  /*  LIS Type A Context  */
  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
  /*  LIS Type A Offspring Context  */
  16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16,
  /*  LIS Type B Context  */
  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
  /* Sign Context  */
  2,
  2
};
#define QCCSPIHT_LIP_CONTEXT 0
#define QCCSPIHT_LIS_TYPEA_CONTEXT (QCCSPIHT_LIP_CONTEXT + 15)
#define QCCSPIHT_LIS_TYPEA_OFFSPRING_CONTEXT (QCCSPIHT_LIS_TYPEA_CONTEXT + 16)
#define QCCSPIHT_LIS_TYPEB_CONTEXT (QCCSPIHT_LIS_TYPEA_OFFSPRING_CONTEXT + 16)
#define QCCSPIHT_SIGN_CONTEXT (QCCSPIHT_LIS_TYPEB_CONTEXT + 16)
#define QCCSPIHT_REFINEMENT_CONTEXT (QCCSPIHT_SIGN_CONTEXT + 1)
#define QCCSPIHT_NUM_CONTEXTS (QCCSPIHT_REFINEMENT_CONTEXT + 1)


static QccSPIHTCoefficientBlock
*QccSPIHTGetCoefficientBlockFromNode(QccListNode *node)
{
  return((QccSPIHTCoefficientBlock *)(node->value));
}


static int QccSPIHTInputOutput(QccBitBuffer *buffer,
			       int *symbol,
			       int method,
			       QccENTArithmeticModel *model,
			       int target_bit_cnt)
{
  
  if (method == QCCSPIHT_ENCODE)
    {
      if (model == NULL)
	{
	  if (QccBitBufferPutBit(buffer, *symbol))
	    {
	      QccErrorAddMessage("(QccSPIHTInputOutput): Error calling QccBitBufferPutBit()");
	      return(1);
	    }
	}
      else
	{
	  if (QccENTArithmeticEncode(symbol, 1,
				     model, buffer))
	    {
	      QccErrorAddMessage("(QccSPIHTInputOutput): Error calling QccENTArithmeticEncode()");
	      return(1);
	    }
	}
      if (buffer->bit_cnt >= target_bit_cnt)
	{
	  QccBitBufferFlush(buffer);
	  return(1);
	}
    }
  else
    if (model == NULL)
      {
	if (QccBitBufferGetBit(buffer, symbol))
	  {
	    QccErrorAddMessage("(QccSPIHTInputOutput): Error calling QccBitBufferGetBit()");
	    return(1);
	  }
      }
    else
      if (QccENTArithmeticDecode(buffer,
				 model,
				 symbol, 1))
	{
	  QccErrorAddMessage("(QccSPIHTInputOutput): Error calling QccENTArithmeticDecode()");
	  return(1);
	}

  return(0);
}


int QccSPIHTEncodeHeader(QccBitBuffer *buffer, 
			 int num_levels, 
			 int num_rows, int num_cols,
			 double image_mean,
			 int max_coefficient_bits,
			 int arithmetic_coded)
{
  int return_value;
  
  if (QccBitBufferPutChar(buffer, (unsigned char)num_levels))
    {
      QccErrorAddMessage("(QccSPIHTEncodeHeader): Error calling QccBitBufferPuChar()");
      goto QccError;
    }
  
  if (QccBitBufferPutInt(buffer, num_rows))
    {
      QccErrorAddMessage("(QccSPIHTEncodeHeader): Error calling QccBitBufferPutInt()");
      goto QccError;
    }
  
  if (QccBitBufferPutInt(buffer, num_cols))
    {
      QccErrorAddMessage("(QccSPIHTEncodeHeader): Error calling QccBitBufferPutInt()");
      goto QccError;
    }
  
  if (QccBitBufferPutDouble(buffer, image_mean))
    {
      QccErrorAddMessage("(QccSPIHTEncodeHeader): Error calling QccBitBufferPutDouble()");
      goto QccError;
    }
  
  if (QccBitBufferPutChar(buffer, (unsigned char)max_coefficient_bits))
    {
      QccErrorAddMessage("(QccSPIHTEncodeHeader): Error calling QccBitBufferPutChar()");
      goto QccError;
    }
  
  if (QccBitBufferPutBit(buffer, arithmetic_coded))
    {
      QccErrorAddMessage("(QccSPIHTEncodeHeader): Error calling QccBitBufferPutBit()");
      goto QccError;
    }
  
  return_value = 0;
  goto QccReturn;
 QccError:
  return_value = 1;
 QccReturn:
  return(return_value);
}


static int QccSPIHTEncodeDWT(QccWAVSubbandPyramid *image_subband_pyramid,
			     int **sign_array,
			     const QccIMGImageComponent *image,
			     int num_levels,
			     double *image_mean,
			     int *max_coefficient_bits,
			     QccWAVSubbandPyramid *mask_subband_pyramid,
			     const QccIMGImageComponent *mask,
			     const QccWAVWavelet *wavelet,
			     const QccWAVPerceptualWeights
			     *perceptual_weights)
{
  double coefficient_magnitude;
  double max_coefficient = -MAXFLOAT;
  int row, col;
  
  if (QccMatrixCopy(image_subband_pyramid->matrix, image->image,
		    image->num_rows, image->num_cols))
    {
      QccErrorAddMessage("(QccSPIHTEncodeDWT): Error calling QccMatrixCopy()");
      return(1);
    }
  
  if (mask != NULL)
    {
      if (QccMatrixCopy(mask_subband_pyramid->matrix, mask->image,
			mask->num_rows, mask->num_cols))
	{
	  QccErrorAddMessage("(QccSPIHTEncodeDWT): Error calling QccMatrixCopy()");
	  return(1);
	}

      for (row = 0; row < image_subband_pyramid->num_rows; row++)
	for (col = 0; col < image_subband_pyramid->num_cols; col++)
	  if (QccAlphaTransparent(mask_subband_pyramid->matrix[row][col]))
	    image_subband_pyramid->matrix[row][col] = 0;
    }

  if (QccWAVSubbandPyramidSubtractMean(image_subband_pyramid,
				       image_mean,
				       NULL))
    {
      QccErrorAddMessage("(QccSPIHTEncodeDWT): Error calling QccWAVSubbandPyramidSubtractMean()");
      return(1);
    }
  
  if (mask != NULL)
    {
      if (QccWAVSubbandPyramidShapeAdaptiveDWT(image_subband_pyramid,
					       mask_subband_pyramid,
					       num_levels,
					       wavelet))
	{
	  QccErrorAddMessage("(QccSPIHTEncodeDWT): Error calling QccWAVSubbandPyramidShapeAdaptiveDWT()");
	  return(1);
	}
    }
  else
    if (QccWAVSubbandPyramidDWT(image_subband_pyramid,
				num_levels,
				wavelet))
      {
	QccErrorAddMessage("(QccSPIHTEncodeDWT): Error calling QccWAVSubbandPyramidDWT()");
	return(1);
      }
  
  if (perceptual_weights != NULL)
    if (QccWAVPerceptualWeightsApply(image_subband_pyramid,
				     perceptual_weights))
      {
	QccErrorAddMessage("(QccSPIHTEncodeDWT): Error calling QccWAVPerceptualWeightsApply()");
	return(1);
      }
  
  for (row = 0; row < image_subband_pyramid->num_rows; row++)
    for (col = 0; col < image_subband_pyramid->num_cols; col++)
      {
	coefficient_magnitude = fabs(image_subband_pyramid->matrix[row][col]);
	if (image_subband_pyramid->matrix[row][col] != coefficient_magnitude)
	  {
	    image_subband_pyramid->matrix[row][col] = coefficient_magnitude;
	    sign_array[row][col] = 1;
	  }
	else
	  sign_array[row][col] = 0;
	if (coefficient_magnitude > max_coefficient)
	  max_coefficient = coefficient_magnitude;
      }
  
  *max_coefficient_bits = (int)floor(QccMathLog2(max_coefficient));
  
  return(0);
}


static int QccSPIHTDecodeInverseDWT(QccWAVSubbandPyramid
				    *image_subband_pyramid,
				    QccWAVSubbandPyramid
				    *mask_subband_pyramid,
				    int **sign_array,
				    QccIMGImageComponent *image,
				    double image_mean,
				    const QccWAVWavelet *wavelet,
				    const QccWAVPerceptualWeights
				    *perceptual_weights)
{
  int row, col;
  
  for (row = 0; row < image_subband_pyramid->num_rows; row++)
    for (col = 0; col < image_subband_pyramid->num_cols; col++)
      if (sign_array[row][col])
	image_subband_pyramid->matrix[row][col] *= -1;
  
  if (perceptual_weights != NULL)
    if (QccWAVPerceptualWeightsRemove(image_subband_pyramid,
				      perceptual_weights))
      {
	QccErrorAddMessage("(QccSPIHTDecodeInverseDWT): Error calling QccWAVPerceptualWeightsRemove()");
	return(1);
      }
  
  if (mask_subband_pyramid != NULL)
    {
      if (QccWAVSubbandPyramidInverseShapeAdaptiveDWT(image_subband_pyramid,
						      mask_subband_pyramid,
						      wavelet))
	{
	  QccErrorAddMessage("(QccSPIHTDecodeInverseDWT): Error calling QccWAVSubbandPyramidInverseShapeAdaptiveDWT()");
	  return(1);
	}
    }
  else
    if (QccWAVSubbandPyramidInverseDWT(image_subband_pyramid,
				       wavelet))
      {
	QccErrorAddMessage("(QccSPIHTDecodeInverseDWT): Error calling QccWAVSubbandPyramidInverseDWT()");
	return(1);
      }
  
  if (QccWAVSubbandPyramidAddMean(image_subband_pyramid,
				  image_mean))
    {
      QccErrorAddMessage("(QccSPIHTDecodeInverseDWT): Error calling QccWAVSubbandPyramidAddMean()");
      return(1);
    }
  
  if (mask_subband_pyramid != NULL)
    for (row = 0; row < image_subband_pyramid->num_cols; row++)
      for (col = 0; col < image_subband_pyramid->num_cols; col++)
	if (QccAlphaTransparent(mask_subband_pyramid->matrix[row][col]))
	  image_subband_pyramid->matrix[row][col] = 0;

  if (QccMatrixCopy(image->image, image_subband_pyramid->matrix,
		    image->num_rows, image->num_cols))
    {
      QccErrorAddMessage("(QccSPIHTDecodeInverseDWT): Error calling QccMatrixCopy()");
      return(1);
    }
  
  return(0);
}


static int QccSPIHTAlgorithmInitialize(QccWAVSubbandPyramid
				       *subband_pyramid,
				       QccList *LIP,
				       QccList *LIS,
				       int block_size)
{
  int subband_num_rows, subband_num_cols;
  int row, col;
  QccSPIHTCoefficientBlock coefficient_block;
  QccListNode *new_node;
  
  if (QccWAVSubbandPyramidSubbandResolution(subband_pyramid,
					    0,
					    &subband_num_rows,
					    &subband_num_cols))
    {
      QccErrorAddMessage("(QccSPIHTAlgorithmInitialize): Error calling QccWAVSubbandPyramidSubbandResolution()");
      return(1);
    }
  
  for (row = 0; row < subband_num_rows; row += block_size)
    for (col = 0; col < subband_num_cols; col += block_size)
      {
	coefficient_block.row = row;
	coefficient_block.col = col;
	coefficient_block.type = QCCSPIHT_TYPE_A;
	coefficient_block.state = 0;
	if ((new_node =
	     QccListCreateNode(sizeof(QccSPIHTCoefficientBlock),
			       (void *)(&coefficient_block))) == NULL)
	  {
	    QccErrorAddMessage("(QccSPIHTAlgorithmInitialize): Error calling QccListCreateNode()");
	    return(1);
	  }
	if (QccListAppendNode(LIP, new_node))
	  {
	    QccErrorAddMessage("(QccSPIHTAlgorithmInitialize): Error calling QccListAppendNode()");
	    return(1);
	  }
	
	if ((new_node =
	     QccListCreateNode(sizeof(QccSPIHTCoefficientBlock),
			       (void *)(&coefficient_block))) == NULL)
	  {
	    QccErrorAddMessage("(QccSPIHTAlgorithmInitialize): Error calling QccListCreateNode()");
	    return(1);
	  }

	if (QccListAppendNode(LIS, new_node))
	  {
	    QccErrorAddMessage("(QccSPIHTAlgorithmInitialize): Error calling QccListAppendNode()");
	    return(1);
	  }
      }
  
  return(0);
}


static int QccSPIHTGetOffspring(QccWAVSubbandPyramid
				*subband_pyramid,
				int row, int col,
				int *offspring_row,
				int *offspring_col,
				int block_size)
{
  int baseband_num_rows;
  int baseband_num_cols;
  int num_offspring;
  
  if ((2*row >= subband_pyramid->num_rows) ||
      (2*col >= subband_pyramid->num_cols))
    return(0);
  
  QccWAVSubbandPyramidSubbandResolution(subband_pyramid,
					0,
					&baseband_num_rows,
					&baseband_num_cols);
  
  if ((row < baseband_num_rows) &&
      (col < baseband_num_cols))
    {
      num_offspring = 3;
      offspring_row[0] = row;
      offspring_col[0] = col + baseband_num_cols;
      offspring_row[1] = row + baseband_num_rows;
      offspring_col[1] = col;
      offspring_row[2] = row + baseband_num_rows;
      offspring_col[2] = col + baseband_num_cols;
    }
  else
    {
      num_offspring = 4;
      offspring_row[0] = row * 2;
      offspring_col[0] = col * 2;
      offspring_row[1] = row * 2 + block_size;
      offspring_col[1] = col * 2;
      offspring_row[2] = row * 2;
      offspring_col[2] = col * 2 + block_size;
      offspring_row[3] = row * 2 + block_size;
      offspring_col[3] = col * 2 + block_size;
    }
  
  return(num_offspring);
}


static double QccSPIHTMaximumCoefficient(QccWAVSubbandPyramid
					 *subband_pyramid,
					 int base_row, int base_col,
					 int block_size)
{
  int row, col;
  double maximum_coefficient = -MAXFLOAT;
  
  for (row = 0; row < block_size; row++)
    for (col = 0; col < block_size; col++)
      if (subband_pyramid->matrix[base_row + row][base_col + col] >
	  maximum_coefficient)
	maximum_coefficient =
	  subband_pyramid->matrix[base_row + row][base_col + col];
  
  return(maximum_coefficient);
}


static int QccSPIHTDescendantsSignificance(QccWAVSubbandPyramid
					   *subband_pyramid,
					   int row, int col,
					   double threshold,
					   int block_size)
{
  int offspring_row[4];
  int offspring_col[4];
  int num_offspring;
  int offspring;
  
  num_offspring = 
    QccSPIHTGetOffspring(subband_pyramid,
			 row, col,
			 offspring_row, offspring_col,
			 block_size);
  
  if (!num_offspring)
    return(-1);
  
  
  for (offspring = 0; offspring < num_offspring; offspring++)
    if (QccSPIHTMaximumCoefficient(subband_pyramid,
				   offspring_row[offspring],
				   offspring_col[offspring],
				   block_size) >=
	threshold)
      return(1);
  
  for (offspring = 0; offspring < num_offspring; offspring++)
    if (QccSPIHTDescendantsSignificance(subband_pyramid,
					offspring_row[offspring],
					offspring_col[offspring],
					threshold,
					block_size) == 1)
      return(1);
  
  return(0);
}


static int QccSPIHTLSignificance(QccWAVSubbandPyramid
				 *subband_pyramid,
				 int row, int col,
				 double threshold,
				 int block_size)
{