hplx_info.c

JPEG2000 EBCOT算法源码

/*****************************************************************************/
/* Copyright 1998, Hewlett-Packard Company                                   */
/* All rights reserved                                                       */
/* Author: David Taubman                                                     */
/* Version: V2.0                                                             */
/* Last Revised: 9/22/98                                                     */
/*****************************************************************************/
#include "local_heap.h"
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <assert.h>
#include "line_block_ifc.h"
#include "hplx_coeffs.h"


#define OVFL_ALLOWANCE 1.2

/* ========================================================================= */
/* ------------------- Implementation of my_filter_info -------------------- */
/* ========================================================================= */

typedef
  struct tap_info {
    int neg_support;
    int pos_support;
    float *taps;
  } tap_info, *tap_info_ptr;

typedef
  struct band_info {
    float step_size;
    float step_wmse;
    int nearest_exponent;
    int extra_lsbs;
  } band_info, *band_info_ptr;

typedef
  struct level_info {
    tap_info hor_low[2];    /* First entry in each array is for analysis;  */
    tap_info hor_high[2];   /* second entry in array is for synthesis. */
    tap_info vert_low[2];
    tap_info vert_high[2];
    int min_band, max_band; /* Min and max allowable band indices. */
    band_info bands[4];
  } level_info, *level_info_ptr;

typedef
  struct component_info {
    int num_levels;
    level_info_ptr levels;
    char *id_copy;
    int ideal_quant;
    int radix_quant;
  } component_info, *component_info_ptr;

typedef
  struct my_filter_info {
    filter_info_obj base;
    int num_components;
    int num_levels;
    component_info_ptr components;
  } my_filter_info_obj, *my_filter_info_ref;
  /* The `levels' array contains one entry for each resolution level, plus
     one additional level (at the end of the array) for holding information
     about the image itself. */

/*****************************************************************************/
/* STATIC                        copy_tap_info                               */
/*****************************************************************************/

static void
  copy_tap_info(tap_info_ptr src, tap_info_ptr dest)
{
  int num_taps, i;

  dest->neg_support = src->neg_support;
  dest->pos_support = src->pos_support;
  num_taps = src->neg_support + src->pos_support + 1;
  dest->taps = (float *)
    local_malloc(sizeof(float) * (size_t) num_taps);
  for (i=0; i < num_taps; i++)
    dest->taps[i] = src->taps[i];
}

/*****************************************************************************/
/* STATIC                      check_built_in_gain                           */
/*****************************************************************************/

static void
  check_built_in_gain(double sum, char *id)
{
  sum = fabs(sum) - 1.0;
  if ((sum < -0.001) || (sum > 0.001))
    {
      fprintf(stderr,"Build in filter set, \"%s\", is incorrectly "
              "defined!\n  Built in filters must have unit DC and "
              "Nyquist gain!\n",id);
      exit(-1);
    }
}

/*****************************************************************************/
/* STATIC                       set_built_in_taps                            */
/*****************************************************************************/

static int
  set_built_in_taps(tap_info_ptr hor_low, tap_info_ptr hor_high,
                    tap_info_ptr vert_low, tap_info_ptr vert_high,
                    char *id)

 /* Returns 1 if the `id' string was able to be interpreted as a built-in
    filter specification.  Returns 0 otherwise.  "0" is interpreted as
    the Daubechies 9/7 biorthogonal filters for the moment.  The function
    fills out the low- and high-pass analysis taps for horizontal and
    vertical filters.  Normalization and generation of synthesis taps are
    left to the `complete_1d_filter_bank' function.  Note that the
    analysis taps are organized for inner products (reverse order to
    filtering), as explained in the interface definition header file. */

{
  static double d97_low[9] =
    {ha4, ha3, ha2, ha1, ha0, ha1, ha2, ha3, ha4};
  static double d97_high[7] =
    {ga3, ga2, ga1, ga0, ga1, ga2, ga3};
  int i;
  double sum;

  if (strcmp(id,"0") == 0)
    {
      hor_low->neg_support = hor_low->pos_support = 4; /* 9-tap. */
      hor_low->taps = (float *)
        local_malloc(sizeof(float)*9);
      for (i=0, sum=0.0; i < 9; i++)
        sum += (hor_low->taps[i] = (float) d97_low[i]);
      check_built_in_gain(sum,id);
      hor_high->neg_support = hor_high->pos_support = 3; /* 7-tap. */
      hor_high->taps = (float *)
        local_malloc(sizeof(float)*7);
      for (i=0, sum=0.0; i < 7; i++)
        sum += (hor_high->taps[i] = (float) d97_high[i]) * ((i&1)?1.0:-1.0);
      check_built_in_gain(sum,id);
      copy_tap_info(hor_low,vert_low);
      copy_tap_info(hor_high,vert_high);
      return(1);
    }
  return(0);
}

/*****************************************************************************/
/* STATIC                        read_tap_info                              */
/*****************************************************************************/

static int
  read_tap_info(tap_info_ptr low, tap_info_ptr high, FILE *fp)

 /* Reads support sizes and tap values for the low- and high-pass wavelet
    filters from the supplied file.  Returns 0 if there are no more filters
    in the file.  Otherwise, returns 1.  Two file formats are supported:
    1) VM format
      The first line should contain four integers of the form
      <SL> <SH> <SH> <SL>
      where SL denotes the length of the low-pass analysis filter and
      SH denotes the length of the high-pass analysis filter.  These should
      be followed by the low-pass analysis tap values (one per line) and then
      the high-pass analysis tap values (one per line).  In this case, the
      remainder of the file is discarded since synthesis filters can always
      be constructed from the analysis filters.
    2) Local format
      This file format is introduced by the appearance of an initial line
      containing the string, "<analysis low>".
      Support sizes (negative and positive support distances) for the low-pass
      filter should appear on a single line, followed by the tap values with
      one line for each tap.  This should be followed by a line containing
      the string, "<analysis high>", followed by support sizes and tap values
      for the high-pass filter.  Filter taps and support sizes are
      quoted for inner product (waveform interpretation), rather than
      filter impulse responses.  For other aspects of the interpretation
      consult the comments in "line_block_ifc.h" appearing with the
      `filter_info__get_taps' function. */

{
  tap_info_ptr taps;
  char line[257];
  int neg, pos, num, n;
  int format, sl1, sh1, sl2, sh2;

  format = 0;
  do {
    if (fgets(line,256,fp) == NULL)
      return(0);
    if (strncmp(line,"<analysis low>",strlen("<analysis low>")) == 0)
      format = 2;
    else if ((sscanf(line,"%d %d %d %d",&sl1,&sh1,&sh2,&sl2) == 4) &&
             (sl1 == sl2) && (sh1 == sh2))
      format = 1;
  } while (!format);
  if (format == 1)
    { /* VM format. */
      taps = low;
      neg = (sl1 - 1) / 2;
      pos = sl1 / 2;
      taps->neg_support = neg;
      taps->pos_support = pos;
      num = neg + pos + 1;
      assert(num == sl1);
      taps->taps = (float *)
        local_malloc(sizeof(float) * (size_t) num);
      for (n=0; n < num; n++)
        if ((fgets(line,256,fp) == NULL) ||
            (sscanf(line,"%f",taps->taps+n) != 1))
          {
            fprintf(stderr,"Syntax error in wavelet filter file!!\n"
                    "  Insufficient tap values found for filter.\n");
            exit(-1);
          }
      taps = high;
      neg = sh1 / 2;
      pos = (sh1 - 1) / 2;
      taps->neg_support = neg;
      taps->pos_support = pos;
      num = neg + pos + 1;
      assert(num == sh1);
      taps->taps = (float *)
        local_malloc(sizeof(float) * (size_t) num);
      for (n=0; n < num; n++)
        if ((fgets(line,256,fp) == NULL) ||
            (sscanf(line,"%f",taps->taps+n) != 1))
          {
            fprintf(stderr,"Syntax error in wavelet filter file!!\n"
                    "  Insufficient tap values found for filter.\n");
            exit(-1);
          }
      while (fgets(line,256,fp) != NULL);
    }
  else
    { /* Non-VM format. */
      taps = low;
      if ((fgets(line,256,fp) == NULL) ||
          (sscanf(line,"%d %d",&neg,&pos) != 2) ||
          (neg < 0) || (pos < 0))
        {
          fprintf(stderr,"Syntax error in wavelet filter file!!\n"
                  "  Each set of filter taps must be preceded by two +ve "
                  "integers,\n  representing the negative and positive "
                  "support sizes.\n");
          exit(-1);
        }
      taps->neg_support = neg;
      taps->pos_support = pos;
      num = neg+pos+1;
      taps->taps = (float *)
        local_malloc(sizeof(float) * (size_t) num);
      for (n=0; n < num; n++)
        if ((fgets(line,256,fp) == NULL) ||
            (sscanf(line,"%f",taps->taps+n) != 1))
          {
            fprintf(stderr,"Syntax error in wavelet filter file!!\n"
                    "  Insufficient tap values found for filter.\n");
            exit(-1);
          }

      taps = high;
      if ((fgets(line,256,fp) == NULL) ||
          (strncmp(line,"<analysis high>",strlen("<analysis high>")) != 0))
        {
          fprintf(stderr,"Syntax error in wavelet filter file!!\n"
                  "  Expected to find a string of the form,\n"
                  "  \"<analysis high>\" before the high-pass tap values!\n");
          exit(-1);
        }
      if ((fgets(line,256,fp) == NULL) ||
          (sscanf(line,"%d %d",&neg,&pos) != 2) ||
          (neg < 0) || (pos < 0))
        {
          fprintf(stderr,"Syntax error in wavelet filter file!!\n"
                  "  Each set of filter taps must be preceded by two +ve "
                  "integers,\n  representing the negative and positive "
                  "support sizes.\n");
          exit(-1);
        }
      taps->neg_support = neg;
      taps->pos_support = pos;
      num = neg+pos+1;
      taps->taps = (float *)
        local_malloc(sizeof(float) * (size_t) num);
      for (n=0; n < num; n++)
        if ((fgets(line,256,fp) == NULL) ||
            (sscanf(line,"%f",taps->taps+n) != 1))
          {
            fprintf(stderr,"Syntax error in wavelet filter file!!\n"
                    "  Insufficient tap values found for filter.\n");
            exit(-1);
          }
    }
  return(1);
}

/*****************************************************************************/
/* STATIC                     read_taps_from_file                            */
/*****************************************************************************/

static void
  read_taps_from_file(tap_info_ptr hor_low, tap_info_ptr hor_high,