hplx_info.c

JPEG2000 EBCOT算法源码

                      tap_info_ptr vert_low, tap_info_ptr vert_high,
                      char *filename)
{
  FILE *fp;
  int success;

  fp = fopen(filename,"r");
  if (fp == NULL)
    {
      fprintf(stderr,"Unable to open filter tap information file, \"%s\"!\n",
              filename);
      exit(-1);
    }
  success = read_tap_info(hor_low,hor_high,fp);
  if (success && !read_tap_info(vert_low,vert_high,fp))
    {
      copy_tap_info(hor_low,vert_low);
      copy_tap_info(hor_high,vert_high);
    }
  if (!success)
    {
      fprintf(stderr,"Insufficent filters specified in file, \"%s\"!\n",
              filename);
      exit(-1);
    }
  fclose(fp);
}

/*****************************************************************************/
/* STATIC                 complete_1d_filter_bank                            */
/*****************************************************************************/

static void
  complete_1d_filter_bank(tap_info_ptr a_low, tap_info_ptr a_high,
                          tap_info_ptr s_low, tap_info_ptr s_high)

 /* On entry, the low and high pass analysis taps must be available.
    The function generates corresponding synthesis filters and
    normalizes all filters in an appropriate manner.  Note that the
    analysis taps appear in inner product order (opposite order to
    impulse responses of filters), as explained in the interface definition
    header file. */

{
  int num_taps, i;
  float *sp, *dp;
  double lp_factor, hp_factor, hp_sign;

  /* Compute high-pass synthesis from low-pass analysis. */

  s_high->pos_support = a_low->neg_support;
  s_high->neg_support = a_low->pos_support;
  num_taps = s_high->pos_support + s_high->neg_support + 1;
  s_high->taps = (float *)
    local_malloc(sizeof(float)*(size_t) num_taps);
  sp = a_low->taps + a_low->neg_support;
  dp = s_high->taps + s_high->neg_support;
  for (i=-s_high->neg_support; i <= s_high->pos_support; i++)
    dp[i] = sp[-i] * ((i&1)?-1.0F:1.0F);

  /* Compute low-pass synthesis from high-pass analysis. */

  s_low->pos_support = a_high->neg_support;
  s_low->neg_support = a_high->pos_support;
  num_taps = s_low->pos_support + s_low->neg_support + 1;
  s_low->taps = (float *)
    local_malloc(sizeof(float)*(size_t) num_taps);
  sp = a_high->taps + a_high->neg_support;
  dp = s_low->taps + s_low->neg_support;
  for (i=-s_low->neg_support; i <= s_low->pos_support; i++)
    dp[i] = sp[-i] * ((i&1)?-1.0F:1.0F);

  /* Now normalize the filters for unit DC and Nyquist gain. */

  num_taps = a_low->neg_support + a_low->pos_support + 1;
  for (sp=a_low->taps, lp_factor=0.0, i=num_taps; i > 0; i--)
    lp_factor += *(sp++);
  lp_factor = 1.0 / lp_factor;
  for (sp=a_low->taps, i=num_taps; i > 0; i--)
    *(sp++) *= (float) lp_factor;

  num_taps = s_low->neg_support + s_low->pos_support + 1;
  for (sp=s_low->taps, hp_factor=0.0, i=num_taps; i > 0; i--)
    hp_factor += *(sp++);  /* synthesis LP at DC is analysis HP at Nyquist */
  hp_sign = (hp_factor < 0.0)?-1.0:1.0;
  hp_factor = 1.0 / hp_factor;
  for (sp=s_low->taps, i=num_taps; i > 0; i--)
    *(sp++) *= (float)(2.0*hp_factor);

  num_taps = a_high->neg_support + a_high->pos_support + 1;
  for (sp=a_high->taps, i=num_taps; i > 0; i--)
    *(sp++) *= (float) (hp_sign * hp_factor);

  num_taps = s_high->neg_support + s_high->pos_support + 1;
  for (sp=s_high->taps, i=num_taps; i > 0; i--)
    *(sp++) *= (float)(2.0 * lp_factor * hp_sign);
}

/*****************************************************************************/
/* STATIC                     get_nearest_exponent                           */
/*****************************************************************************/

static int
  get_nearest_exponent(float step)

 /* Returns the largest integer, s, such that 2^s is no larger than
    `step' * OVFL_ALLOWANCE (OVFL_ALLOWANCE is typically a little larger
    than 1.0).  Note that the integer may be negative or positive. */

{
  int shift;
  float val;

  shift = 0;
  val = 1.0F / step;
  while (val < 1.0F)
    {
      shift++;
      val *= 2.0F;
    }
  while (val > (float) OVFL_ALLOWANCE)
    {
      shift--;
      val *= 0.5F;
    }
  return(shift);
}

/*****************************************************************************/
/* STATIC                         get_extra_lsbs                             */
/*****************************************************************************/

static int
  get_extra_lsbs(float step, float nom_range)

 /* Returns the largest integer, E, such that the ratio between
    `nom_range' and (`step' * 2^{-E}) does not exceed 2^R * OVFL_ALLOWANCE,
    where R is the STD_NOMINAL_RANGE_BITS value.  The value returned here
    may be zero or even negative, but the caller will make any necessary
    adjustments to the quantization step size in order to ensure that the
    number of extra LSB's will ultimately be greater than zero. */

{
  int extra;
  float val;

  extra = 0;
  val = nom_range / (step * (float)(1<<STD_NOMINAL_RANGE_BITS));
  while (val < 1.0F)
    {
      extra++;
      val *= 2.0F;
    }
  while (val > (float) OVFL_ALLOWANCE)
    {
      extra--;
      val *= 0.5F;
    }
  return(extra); 
}

/*****************************************************************************/
/* STATIC                     convolve_and_get_norm                          */
/*****************************************************************************/

static float
  convolve_and_get_norm(tap_info_ptr synth, tap_info_ptr in, tap_info_ptr out,
                        int interp_factor)

 /* This function is key to the determination of synthesis waveforms and their
    L2-norms.  Essentially, the function takes the synthesis waveform provided
    via the `synth' argument, interpolates by a `interp_factor', inserting
    `interp_factor'-1 zeros between the entries, and convolves with the
    filter represented by `in', writing the result into `out' and returning
    the L2-norm of the final set of waveform coefficients. */

{
  int out_neg, out_pos, in_taps, out_taps, synth_taps, i, k;
  float *sp, *op, *ip, *opp, val;
  double sum;

  out_neg = synth->neg_support*interp_factor + in->neg_support;
  out_pos = synth->pos_support*interp_factor + in->pos_support;
  out->neg_support = out_neg;
  out->pos_support = out_pos;
  out_taps = out_neg+out_pos+1;
  out->taps = (float *)
    local_realloc(out->taps,sizeof(float)*(size_t) out_taps);

  /* Perform interpolated convolution, as described above. */
  
  for (i=0; i < out_taps; i++)
    out->taps[i] = 0.0F; /* Ready to accumulate. */
  synth_taps = synth->neg_support + synth->pos_support + 1;
  in_taps = in->neg_support + in->pos_support + 1;
  op = NULL;
  for (sp=synth->taps, opp=out->taps,
       i=0; i < synth_taps; i++, opp+=interp_factor)
    {
      val = *(sp++);
      for (k=0, ip=in->taps, op=opp; k < in_taps; k++)
        *(op++) += val * *(ip++);
    }
  assert(op == (out->taps + out_taps));

  /* Compute L2-norm. */

  for (sum=0.0, i=0; i < out_taps; i++)
    {
      val = out->taps[i];
      sum += val*val;
    }
  return((float) sqrt(sum));
}

/*****************************************************************************/
/* STATIC                       __initialize                                 */
/*****************************************************************************/

static void
  __initialize(filter_info_ref base, int decomposition, int levels,
               int components, char *component_ids[],
               int component_flags[], float component_base_steps[],
               int *max_bitplanes, int argc, char *argv[])
{
  my_filter_info_ref self = (my_filter_info_ref) base;
  component_info_ptr comp_info;
  char *id;
  level_info_ptr lev;
  int in_bitplanes, out_bitplanes;
  float base_step, nom_range, mse, factor;
  float hor_low_norm, hor_high_norm, vert_low_norm, vert_high_norm, ideal_norm;
  float save_hor_low_norm, save_vert_low_norm;
  tap_info conv_h1, conv_h2, conv_v1, conv_v2, tmp_taps;
  int i, n, comp, min_extra_lsbs;

  assert((levels >= 0) && (components > 0));
  if (decomposition != DECOMPOSITION__MALLAT)
    {
      fprintf(stderr,"The `filter_info' object currently only supports "
              "Mallat style\n  decompositions!\n");
      exit(-1);
    }
  in_bitplanes = *max_bitplanes;
  out_bitplanes = IMPLEMENTATION_PRECISION-2;
  if (in_bitplanes == 0)
    in_bitplanes = out_bitplanes;
  *max_bitplanes = out_bitplanes;
  min_extra_lsbs = 1 - (in_bitplanes-out_bitplanes);
  self->num_levels = levels;
  self->num_components = components;
  self->components = (component_info_ptr)
    local_malloc(sizeof(component_info)*(size_t) components);
  memset(self->components,0,sizeof(component_info)*(size_t) components);
  for (comp=0; comp < components; comp++)
    {
      comp_info = self->components + comp;
      comp_info->num_levels = levels;
      if (component_flags[comp] & FILTER_INFO__IDEAL_QUANT)
        comp_info->ideal_quant = 1;
      if (component_flags[comp] & FILTER_INFO__RADIX_QUANT)
        comp_info->radix_quant = 1;
      base_step = component_base_steps[comp];
      base_step *= (float)(1<<STD_NOMINAL_RANGE_BITS);
      comp_info->levels = (level_info_ptr)
        local_malloc(sizeof(level_info)*(size_t)(levels+1));
      for (lev=comp_info->levels, n=0; n < levels; n++, lev++)
        {
          lev->min_band = (n==0)?0:1;
          lev->max_band = 3;
        }
      lev->min_band = lev->max_band = 0;

      /* First, initialize the image-level information in the last
        `level_info' record. */

      lev->bands[0].step_size = base_step;
      nom_range = (float)(1<<STD_NOMINAL_RANGE_BITS);
      mse = base_step / nom_range;
      mse *= mse;
      lev->bands[0].step_wmse = mse;
      lev->bands[0].nearest_exponent = get_nearest_exponent(base_step);
      lev->bands[0].extra_lsbs = get_extra_lsbs(base_step,nom_range);
      while (lev->bands[0].extra_lsbs < min_extra_lsbs)
        {
          lev->bands[0].extra_lsbs++;
          lev->bands[0].step_size *= 2.0F;
          lev->bands[0].nearest_exponent++;
          lev->bands[0].step_wmse *= 4.0F;
        }
      if (levels == 0)
        continue;

      /* Now fill out the filter tap information. */

      id = component_ids[comp];
      assert(id != NULL);
      lev = comp_info->levels;
      if (!set_built_in_taps(lev->hor_low,lev->hor_high,
                             lev->vert_low,lev->vert_high,id))
        read_taps_from_file(lev->hor_low,lev->hor_high,
                            lev->vert_low,lev->vert_high,id);
      complete_1d_filter_bank(lev->hor_low,lev->hor_high,
                              lev->hor_low+1,lev->hor_high+1);