spiht.c

该文件是小波信源编码SPIHT算法的C语言代码

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

/*
 * ----------------------------------------------------------------------------
 * 
 * Public License for the SPIHT Algorithm
 * Version 1.2, March 8, 2004
 * 
 * ----------------------------------------------------------------------------
 * 
 * The Set Partitioning In Hierarchical Trees (SPIHT) algorithm is protected
 * by US Patent #5,764,807 (issued June 9, 1998), US Patent #6,674,911 (issued
 * January 6, 2004), and other international patents and patents pending. An
 * implementation of the SPIHT algorithm is included herein with the gracious
 * permission of Dr. William A. Pearlman, President of PrimaComp, Inc.,
 * exclusive holder of patent rights. PrimaComp, Inc., has 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 Dr. Pearlman.
 * 
 * 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, Inc.,  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, Inc., 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.
 *   851 Maxwell Drive
 *   Schenectady, NY 12309
 *   U.S.A.
 *   [phone]
 *   email: [e-mail]
 *   http://www.cipr.rpi.edu/research/SPIHT
 * 
 * ----------------------------------------------------------------------------
*/


#include "libQccPack.h"


#define QCCSPIHT_TYPE_A 0
#define QCCSPIHT_TYPE_B 1
#define QCCSPIHT_TYPE_DELETE 2
#define QCCSPIHT_ENCODE 0
#define QCCSPIHT_DECODE 1

#define QCCSPIHT_MAXBITPLANES 128

typedef struct
{
  int row;
  int col;
  unsigned char type;
  unsigned char state;
} QccSPIHTCoefficientBlock;

typedef struct
{
  QccWAVSubbandPyramid original_coefficients;
  QccWAVSubbandPyramid reconstructed_coefficients;
  QccWAVSubbandPyramid squared_errors;
  double distortion;
  int num_pixels;
  int start_position;
  FILE *file;
} QccSPIHTDistortionTrace;


#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,
    /* Refinement Context  */
    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 void QccSPIHTRefineCoefficient(double *coefficient,
                                      int bit,
                                      double threshold)
{
  *coefficient +=
    (bit) ? (threshold / 2) : (-threshold / 2);
}


static void QccSPIHTDistortionTraceInitialize(QccSPIHTDistortionTrace
                                              *distortion_trace,
                                              QccWAVSubbandPyramid
                                              *original_coefficients,
                                              QccWAVSubbandPyramid
                                              *mask,
                                              QccBitBuffer *output_buffer,
                                              FILE *outfile)
{
  int row, col;

  QccWAVSubbandPyramidInitialize(&distortion_trace->original_coefficients);
  QccWAVSubbandPyramidInitialize(&distortion_trace->reconstructed_coefficients);
  QccWAVSubbandPyramidInitialize(&distortion_trace->squared_errors);

  QccWAVSubbandPyramidCopy(&distortion_trace->original_coefficients,
                           original_coefficients);
  QccWAVSubbandPyramidCopy(&distortion_trace->reconstructed_coefficients,
                           original_coefficients);
  QccWAVSubbandPyramidCopy(&distortion_trace->squared_errors,
                           original_coefficients);

  QccMatrixZero(distortion_trace->reconstructed_coefficients.matrix,
                distortion_trace->reconstructed_coefficients.num_rows,
                distortion_trace->reconstructed_coefficients.num_cols);

  distortion_trace->distortion = 0;
  distortion_trace->num_pixels = 0;
  distortion_trace->start_position = output_buffer->bit_cnt;

  if (mask == NULL)
    {
      for (row = 0; row < distortion_trace->squared_errors.num_rows; row++)
        for (col = 0; col < distortion_trace->squared_errors.num_cols; col++)
          {
            distortion_trace->squared_errors.matrix[row][col] *=
              distortion_trace->squared_errors.matrix[row][col];
            distortion_trace->distortion +=
              distortion_trace->squared_errors.matrix[row][col];
            distortion_trace->num_pixels++;
          }
    }
  else
    for (row = 0; row < distortion_trace->squared_errors.num_rows; row++)
      for (col = 0; col < distortion_trace->squared_errors.num_cols; col++)
        if (!QccAlphaTransparent(mask->matrix[row][col]))
          {
            distortion_trace->squared_errors.matrix[row][col] *=
              distortion_trace->squared_errors.matrix[row][col];
            distortion_trace->distortion +=
              distortion_trace->squared_errors.matrix[row][col];
            distortion_trace->num_pixels++;
          }
  
  distortion_trace->file = outfile;

  fprintf(outfile,
          "      Rate (bpp)\t      MSE\n");
  fprintf(outfile,
          "--------------------------------------------\n");
}


static void QccSPIHTDistortionTraceUpdate(QccSPIHTDistortionTrace
                                          *distortion_trace,
                                          int row,
                                          int col,
                                          QccBitBuffer *buffer)
{
  double diff;

  QccMatrix original_coefficients =
    distortion_trace->original_coefficients.matrix;
  QccMatrix reconstructed_coefficients =
    distortion_trace->reconstructed_coefficients.matrix;
  QccMatrix squared_errors =
    distortion_trace->squared_errors.matrix;

  distortion_trace->distortion -= squared_errors[row][col];
  diff = original_coefficients[row][col] - 
    reconstructed_coefficients[row][col];
  squared_errors[row][col] = diff * diff;
  distortion_trace->distortion += squared_errors[row][col];

  fprintf(distortion_trace->file,
          "%16.6f\t%16.6f\n",
          (double)(buffer->bit_cnt - distortion_trace->start_position) /
          distortion_trace->num_pixels,
          distortion_trace->distortion / distortion_trace->num_pixels);
}


static void QccSPIHTListClean(QccList *list)
{
  QccListNode *current_list_node;
  QccListNode *next_list_node;
  QccSPIHTCoefficientBlock *current_coefficient_block;

  current_list_node = list->start;
  while (current_list_node != NULL)
    {
      next_list_node = current_list_node->next;

      current_coefficient_block =
        QccSPIHTGetCoefficientBlockFromNode(current_list_node);

      if (current_coefficient_block->type == QCCSPIHT_TYPE_DELETE)
        QccListDeleteNode(list, current_list_node);

      current_list_node = next_list_node;
    }
}


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);
            }
          if (target_bit_cnt != QCCENT_ANYNUMBITS)
            if (buffer->bit_cnt >= target_bit_cnt)
              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 Error;
    }
  
  if (QccBitBufferPutInt(buffer, num_rows))
    {
      QccErrorAddMessage("(QccSPIHTEncodeHeader): Error calling QccBitBufferPutInt()");
      goto Error;
    }
  
  if (QccBitBufferPutInt(buffer, num_cols))
    {
      QccErrorAddMessage("(QccSPIHTEncodeHeader): Error calling QccBitBufferPutInt()");
      goto Error;
    }
  
  if (QccBitBufferPutDouble(buffer, image_mean))
    {
      QccErrorAddMessage("(QccSPIHTEncodeHeader): Error calling QccBitBufferPutDouble()");
      goto Error;
    }
  
  if (QccBitBufferPutInt(buffer, max_coefficient_bits))
    {
      QccErrorAddMessage("(QccSPIHTEncodeHeader): Error calling QccBitBufferPutInt()");
      goto Error;
    }
  
  if (QccBitBufferPutBit(buffer, arithmetic_coded))
    {
      QccErrorAddMessage("(QccSPIHTEncodeHeader): Error calling QccBitBufferPutBit()");
      goto Error;
    }
  
  return_value = 0;
  goto Return;
 Error:
  return_value = 1;
 Return:
  return(return_value);
}


static int QccSPIHTEncodeExtractSigns(QccWAVSubbandPyramid
                                      *image_subband_pyramid,
                                      QccWAVSubbandPyramid
                                      *mask_subband_pyramid,
                                      int **sign_array,
                                      int *max_coefficient_bits)
{
  double coefficient_magnitude;
  double max_coefficient = -MAXDOUBLE;
  int row, col;
  
  if (mask_subband_pyramid == NULL)
    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;