lra.c

JPEG2000实现的源码

/*****************************************************************************/
/* Copyright 1999, Fujifilm Software California Inc.                         */
/* Copyright 1999 Science Applications International Corporation (SAIC).     */
/* Copyright 1995 University of Arizona, Arizona Board of Regents.           */
/* All rights reserved                                                       */
/* File: "lra.c"                                                             */
/* Description: Lagrangian Rate Allocation implementation                    */
/* Author: Troy Chinen                                                       */
/* Affiliation: Fujifilm Software California Inc.                            */
/* Version: VM8.0                                                            */
/* Last Revised: 20 July, 2000                                               */
/*****************************************************************************/
/* 
 * Lagrangian Rate Allocation (LRA) is currently implemented using the
 * modules lra_stats and lra.  The lra_stats routines work within the
 * quantizer object (SQ,TCQ,or Mask) to compute and save subband
 * statistics to files.  The lra routines (called at the rate control
 * level) compute quantization step sizes based on these statistics
 * files.  
 */

/*****************************************************************************/
/* Modified for general wavelet decompositions.                              */
/* Copyright 2000 Science Applications International Corporation (SAIC).     */
/* Copyright 1995 University of Arizona, Arizona Board of Regents.           */
/* All Rights Reserved for modified parts.                                   */
/*****************************************************************************/

#include <local_services.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <assert.h>
#include <ifc.h>
#include <stdio.h>
#include <float.h>
#include "lra_local.h"

#define STEP_MANTISSA_BITS 11 /* `forward_info' and `reverse_info' objects. */


/*****************************************************************************/
/* STATIC                   exponent_mantissa_to_step2                        */
/*****************************************************************************/

static float
  exponent_mantissa_to_step2(int exponent, int mantissa)

 /* Computes a relative or absolute step size as 2{-E} * (1+2^{-Mbits)*M),
    where E is the `exponent' value, M is the `mantissa' value and
    Mbits is the STEP_MANTISSA_BITS constant. */

{
  float step;

  assert((mantissa >= 0) && (mantissa < (1<<STEP_MANTISSA_BITS)));
  step = 1.0F + ((float) mantissa) / ((float)(1<<STEP_MANTISSA_BITS));
  while (exponent > 0)
    {
      step *= 0.5F;
      exponent--;
    }
  while (exponent < 0)
    {
      step *= 2.0F;
      exponent++;
    }
  return(step);
}

/*****************************************************************************/
/* STATIC                   step_to_exponent_mantissa2                        */
/*****************************************************************************/

static void
  step_to_exponent_mantissa2(float step, int *exponent, int *mantissa)

 /* Computes legal `exponent' and `mantissa' values which evaluate as closely
    as possible to `step' when supplied to `exponent_mantissa_to_step2',
    as defined above. */

{
  int exp_val, mant_val;

  exp_val = 0;
  while (step >= 2.0F)
    {
      step *= 0.5F;
      exp_val--;
    }
  while (step < 1.0F)
    {
      step *= 2.0F;
      exp_val++;
    }
  mant_val = (int) floor((step-1.0)*((float)(1<<STEP_MANTISSA_BITS)) + 0.5);
  if (mant_val == (1<<STEP_MANTISSA_BITS))
    {
      mant_val = 0;
      exp_val--;
    }
  assert((mant_val >= 0) && (mant_val < (1<<STEP_MANTISSA_BITS)));
  *exponent = exp_val;
  *mantissa = mant_val;
}

/*****************************************************************************/
/* STATIC               affine_rate_estimate                                 */
/*****************************************************************************/
/* 
 * Estimates a rate for LRA via linear correction based on previous attempts.
 * (Maintains state data from past attempts)
 *
 * This function works with adjusted rates which means that they all
 * rates have the header rate subtracted already (typically a very
 * very small quantity) Now that I think about it, this is probably
 * not such a great idea since the header size is constant...  
 *
 * The rationale for the x and y notation is that I'm viewing LRA as a 
 * function: you put in a "rate requested" (which often, it seems, has
 * nothing whatever to do with a rate) and out pops step sizes which, on
 * the next iteration, will lead to a "rate achieved".  Hence,
 * LRA("rate requested") = "rate achieved".  The x's refer to rates
 * requested and the y's refer to rates achieved.
 *
 * Here's what this function does:
 * 
 * Iteration 0 (first time through)
 *  output = rate_target_adj.  This will be exactly correct if the mapping
 *           LRA(rate requested) = rate achieved is the identity function.
 *           This is probably the best we can do since the rate achieved
 *           at this point corresponds to fake step sizes.
 *
 * Iteration 1
 *  output = the fraction (rate_target_adj / rate_achieved_adj)
 *           times the rate_achieved_adj from the previous iteration.  This will
 *           be exactly correct if the mapping LRA(rate requested) = rate achieved
 *           is linear (no offset).
 *
 * Iteration 2 and higher
 *  output = an affine estimate based on the past two attempts.  This
 *           will be exactly correct if the mapping LRA(rate
 *           requested) = rate achieved is affine.  
 */
static float
  affine_rate_estimate(float rate_new_prev, float rate_achieved_adj, float rate_target_adj, int iterations)
{
  static float y_0 = 0;        /* Rate achieved the previous iteration         */
  static float y_1 = 0;        /* Rate achieved this iteration                 */
  static float x_0 = 0;        /* Rate requested the previous iteration        */
  static float x_1 = 0;        /* Rate requested this iteration                */
  static float x_2 = 0;        /* Rate to use for current step size generation */
  float     y_star;            /* Just so the formulas look nicer              */

  /* Aging.  That is, set z_{n+1} = z_n */
  x_0 = x_1;
  x_1 = x_2;
x_1 = rate_new_prev;
  y_0 = y_1;
  y_1 = rate_achieved_adj;
  y_star = rate_target_adj;

  /* Compute requested rate, handling bootstrap cases */
  if      ( 0==iterations ) {

    /* Try guessing the target bit rate */
    x_2 = y_star;
  }
  else if ( 1==iterations ) {

    /* Because previous iteration had fake steps */
    y_0 = 0;

    /* Affine correction (reduces to linear) */
    x_2 = ((y_star - y_0) / (y_1 - y_0)) * (x_1 - x_0) + x_0;
  }
  else if (y_0 == y_1) {
    x_2 += (x_1 - x_0);
  }
  else {

    /* Affine correction */
    x_2 = ((y_star - y_0) / (y_1 - y_0)) * (x_1 - x_0) + x_0;
  }

  return(x_2);
}

/*****************************************************************************/
/* STATIC               set_models                                           */
/*****************************************************************************/
/* 
 * Initializes the parameter arrays based on whether the quantizer will be
 * Trellis Coded Quantization (TCQ) or Scalar Quantization (SQ)
 * 
 * Lagrangian Rate Allocation uses curve fit parameters which depend
 * on the pdf of the pixels in the subband; our goal is to estimate
 * this pdf.  In this implementation, we assume a that a GGD (which
 * depends on a single parameter alpha) generates our pixels
 * intensities.  Our task thus reduces to estimating this alpha which
 * we estimate based on the subband's kurtosis.  As a practical
 * matter, we partition the range of alpha into five intervals which of
 * course results in 5 GGD types, which in turn results in 5 sets of
 * curve fit parameters.  The interval/class/GGDtype is called a 
 * "codebook" value because it determines which set of parameters are
 * used for a particular subband.
 */
static lra_model_data *set_models(int quantizer)
{

  lra_model_data *model_data;

  model_data = (lra_model_data *) local_malloc(LRA_MEM_KEY, 
                                               sizeof(lra_model_data));
  if(quantizer == LRA_TCQ) {
    model_data->rmh[0] = (float) -1.660964;
    model_data->rmh[1] = (float) -1.660964;
    model_data->rmh[2] = (float) -1.660964;
    model_data->rmh[3] = (float) -1.660964;
    model_data->rmh[4] = (float) -1.660964;
    
    model_data->rbh[0] = (float) -0.2984774;
    model_data->rbh[1] = (float) 0.0764613;
    model_data->rbh[2] = (float) 0.2143643;
    model_data->rbh[3] = (float) 0.3023283;
    model_data->rbh[4] = (float) 0.3186138;
    
    model_data->rmlin[0] = (float) 0.0563047;
    model_data->rmlin[1] = (float)  0.0;
    model_data->rmlin[2] = (float) -0.195004;
    model_data->rmlin[3] = (float) -0.3344047;
    model_data->rmlin[4] = (float) -1.491667;
    
    model_data->rblin[0] = (float) 0.1480394;
    model_data->rblin[1] = (float) 0.0;
    model_data->rblin[2] = (float) -0.124035;
    model_data->rblin[3] = (float) -0.1526097;
    model_data->rblin[4] = (float) -0.3311460;
    
    model_data->rpw[0] = (float) 2.832414e+02;
    model_data->rpw[1] = (float) 14.772349;
    model_data->rpw[2] = (float) 598.09130;
    model_data->rpw[3] = (float) 32.754752;
    model_data->rpw[4] = (float) 70.80316;
    
    model_data->rbs[0] = (float) 72.078073;
    model_data->rbs[1] = (float) 2.254300;
    model_data->rbs[2] = (float) 70.18852;
    model_data->rbs[3] = (float) 1.2153120;
    model_data->rbs[4] = (float) 1.326688;
    
    model_data->rzps[0] = (float) -0.0937507;
    model_data->rzps[1] = (float)  0.0460343;
    model_data->rzps[2] = (float) 0.04866951;
    model_data->rzps[3] = (float) 0.075033;
    model_data->rzps[4] = (float) -0.003975;
    
    model_data->rlimit[0] = (float) 5e-03;
    model_data->rlimit[1] = (float) 1000.0;
    model_data->rlimit[2] = (float) 0.15;
    model_data->rlimit[3] = (float)  0.3;
    model_data->rlimit[4] = (float) 0.6;
    
    model_data->rhyplim[0] = (float) 5e-03;
    model_data->rhyplim[1] = (float) 6.1906222e-03;
    model_data->rhyplim[2] = (float) 0.15;
    model_data->rhyplim[3] = (float) 0.3;
    model_data->rhyplim[4] = (float) 0.6;
    
    model_data->dmh[0] = (float) 0.5;
    model_data->dmh[1] = (float) 0.5;
    model_data->dmh[2] = (float) 0.5;
    model_data->dmh[3] = (float) 0.5;
    model_data->dmh[4] = (float) 0.5;
    
    model_data->dbh[0] = (float) 0.224935;
    model_data->dbh[1] = (float) 0.224935;
    model_data->dbh[2] = (float) 0.224935;
    model_data->dbh[3] = (float) 0.224935;
    model_data->dbh[4] = (float) 0.224935;
    
    model_data->dmlin[0] = (float) 0.027569;
    model_data->dmlin[1] = (float) 0.02369;
    model_data->dmlin[2] = (float) 0.031087;
    model_data->dmlin[3] = (float) 0.0212835;
    model_data->dmlin[4] = (float) 0.04727;
    
    model_data->dblin[0] = (float) 0.10964;
    model_data->dblin[1] = (float) 0.082817;
    model_data->dblin[2] = (float) 0.092541;
    model_data->dblin[3] = (float) 0.062685;
    model_data->dblin[4] = (float) 0.10813;
    
    model_data->dpw[0] = (float) 8.5537e+02;
    model_data->dpw[1] = (float) 1.02504e+08;
    model_data->dpw[2] = (float) 1.5239e+03;
    model_data->dpw[3] = (float) 3.2113e+08;
    model_data->dpw[4] = (float) 6.8090e+03;
    
    model_data->dbs[0] = (float) 2.9323e+02;
    model_data->dbs[1] = (float) 3.26066e+07;
    model_data->dbs[2] = (float) 3.9963e+02;
    model_data->dbs[3] = (float) 8.1289e+07;
    model_data->dbs[4] = (float) 1.8067e+03;
    
    model_data->dzps[0] = (float) -1.5067;
    model_data->dzps[1] = (float) -1.13293;
    model_data->dzps[2] = (float) -0.8759;
    model_data->dzps[3] = (float) -0.5922;
    model_data->dzps[4] = (float) 0.58176;
    
    model_data->dlimit[0] = (float) 1e-04;
    model_data->dlimit[1] = (float) 4e-04;
    model_data->dlimit[2] = (float) 1e-03;
    model_data->dlimit[3] = (float) 1e-03;
    model_data->dlimit[4] = (float) 4e-03;
    
    model_data->dhyplim[0] = (float) 1e-04;
    model_data->dhyplim[1] = (float) 4e-04;
    model_data->dhyplim[2] = (float) 1e-03;
    model_data->dhyplim[3] = (float) 1e-03;
    model_data->dhyplim[4] = (float) 4e-03;