lra.c

JPEG2000实现的源码

    model_data->dmhyp[0] = (float) 0.211731;
    model_data->dmhyp[1] = (float) 0.30278;
    model_data->dmhyp[2] = (float) 0.390313;
    model_data->dmhyp[3] = (float) 0.6517865;
    model_data->dmhyp[4] = (float) 15.37833;
  }
  else if (quantizer == LRA_SQ) { /* VM2.1 scalar quantization? */
    model_data->rmh[0] = (float) -1.6609640e+00;
    model_data->rmh[1] = (float) -1.6609640e+00;
    model_data->rmh[2] = (float) -1.6609640e+00;
    model_data->rmh[3] = (float) -1.6609640e+00;
    model_data->rmh[4] = (float) -1.6609640e+00;
    
    model_data->rbh[0] = (float) -1.0260913e-01;
    model_data->rbh[1] = (float) 2.7437747e-01;
    model_data->rbh[2] = (float) 4.2484655e-01;
    model_data->rbh[3] = (float) 5.0952118e-01;
    model_data->rbh[4] = (float)  5.2848715e-01;
    
    model_data->rmlin[0] = (float) 4.0600905e-01;
    model_data->rmlin[1] = (float) -9.4199861e-02;
    model_data->rmlin[2] = (float) -4.3148024e-01;
    model_data->rmlin[3] = (float) -7.0342040e-01;
    model_data->rmlin[4] = (float) -1.5321467e+00;
    
    model_data->rblin[0] = (float) 3.1856767e-01;
    model_data->rblin[1] = (float) -1.3058530e-01;
    model_data->rblin[2] = (float) -3.3616060e-01;
    model_data->rblin[3] = (float) -4.5592897e-01;
    model_data->rblin[4] = (float) -6.0792320e-01;
    
    model_data->rpw[0] = (float) 3.2224528e+00;
    model_data->rpw[1] = (float) 4.0915129e+00;
    model_data->rpw[2] = (float) 3.8673443e+00;
    model_data->rpw[3] = (float) 2.5888557e+00;
    model_data->rpw[4] = (float) 1.4162744e+00;
    
    model_data->rbs[0] = (float) 5.9115108e-01;
    model_data->rbs[1] = (float) 4.9344464e-01;
    model_data->rbs[2] = (float) 3.8587481e-01;
    model_data->rbs[3] = (float) 1.5083935e-01;
    model_data->rbs[4] = (float) 3.5027385e-02;
    
    model_data->rzps[0] = (float) 1.7208469e-01;
    model_data->rzps[1] = (float) 8.1925205e-02;
    model_data->rzps[2] = (float) 4.2383903e-02;
    model_data->rzps[3] = (float) 2.2666454e-02;
    model_data->rzps[4] = (float) -2.4877324e-02;
    
    model_data->rlimit[0] = (float) 1.6419807e-01;
    model_data->rlimit[1] = (float) 4.1090561e-02;
    model_data->rlimit[2] = (float) 1.6630799e-01;
    model_data->rlimit[3] = (float) 2.2482261e-01;
    model_data->rlimit[4] = (float) 4.0107102e-01;
    
    model_data->rhyplim[0] = (float) 3.8100757e-01;
    model_data->rhyplim[1] = (float) 3.8768640e-01;
    model_data->rhyplim[2] = (float)  4.5342879e-01;
    model_data->rhyplim[3] = (float) 7.4443556e-01;
    model_data->rhyplim[4] = (float) 8.7115472e-01;
    
    model_data->dmh[0] = (float) 5.0000000e-01;
    model_data->dmh[1] = (float) 5.0000000e-01;
    model_data->dmh[2] = (float) 5.0000000e-01;
    model_data->dmh[3] = (float) 5.0000000e-01;
    model_data->dmh[4] = (float) 5.0000000e-01;
    
    model_data->dbh[0] = (float) 4.6866289e-01;
    model_data->dbh[1] = (float) 4.6866289e-01;
    model_data->dbh[2] = (float) 4.6866289e-01;
    model_data->dbh[3] = (float) 4.6866289e-01;
    model_data->dbh[4] = (float) 4.6866289e-01;
    
    model_data->dmlin[0] = (float) 4.9816881e-02;
    model_data->dmlin[1] = (float) 8.3733463e-02;
    model_data->dmlin[2] = (float) 6.4246442e-02;
    model_data->dmlin[3] = (float) 4.3910965e-02;
    model_data->dmlin[4] = (float) 3.6370467e-02;
    
    model_data->dblin[0] = (float) 1.1955573e-01;
    model_data->dblin[1] = (float) 1.6224159e-01;
    model_data->dblin[2] = (float) 9.4898809e-02;
    model_data->dblin[3] = (float) 6.3973807e-02;
    model_data->dblin[4] = (float) 5.4645144e-02;
    
    model_data->dpw[0] = (float) 1.6961542e+01;
    model_data->dpw[1] = (float) 8.5298580e+01;
    model_data->dpw[2] = (float) 8.2061860e+01;
    model_data->dpw[3] = (float) 3.2483845e+01;
    model_data->dpw[4] = (float) 2.9983854e+01;
    
    model_data->dbs[0] = (float) 3.4745758e+00;
    model_data->dbs[1] = (float) 1.0733488e+01;
    model_data->dbs[2] = (float) 4.8907865e+00;
    model_data->dbs[3] = (float) 3.8051394e+00;
    model_data->dbs[4] = (float) 5.3635380e+00;
    
    model_data->dzps[0] = (float) -6.3578223e-01;
    model_data->dzps[1] = (float) -5.2664118e-01;
    model_data->dzps[2] = (float) -6.0012445e-01;
    model_data->dzps[3] = (float)  -5.1308439e-01;
    model_data->dzps[4] = (float) -2.7996784e-02;
    
    model_data->dlimit[0] = (float) 3.9819528e-03;
    model_data->dlimit[1] = (float) 1.1545279e-02;
    model_data->dlimit[2] = (float) 3.3334501e-02;
    model_data->dlimit[3] = (float) 3.4922214e-02;
    model_data->dlimit[4] = (float) 3.1444217e-02;
    
    model_data->dhyplim[0] = (float) 3.9819528e-03;
    model_data->dhyplim[1] = (float) 1.1545279e-02;
    model_data->dhyplim[2] = (float) 3.3334501e-02;
    model_data->dhyplim[3] = (float)  3.4922214e-02;
    model_data->dhyplim[4] = (float) 3.1444217e-02;
    
    model_data->dmhyp[0] = (float) 1.8509322e-01;
    model_data->dmhyp[1] = (float) 1.4690569e-01;
    model_data->dmhyp[2] = (float) 2.8372985e-01;
    model_data->dmhyp[3] = (float) 4.0509855e-01;
    model_data->dmhyp[4] = (float)  2.0850623e+00;
  }
  return(model_data);
}

/*****************************************************************************/
/* STATIC               select_codebook                                      */
/*****************************************************************************/
/* 
 * Estimates the Generalized Gaussian Density (GGD) class as a
 * function of the kurtosis via look up table.  This is done for each
 * subband.
 * 
 * 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 void
  select_codebooks(lra_stats_ref lra_stats)
{
  lra_component_info_ptr lra_comp;  /* Convenience pointer         */
  lra_level_info_ptr     lra_lev;   /* Convenience pointer         */
  lra_band_info_ptr      lra_band;  /* Convenience pointer         */
  int c;                            /* Component level counter     */
  int n;                            /* Level counter               */
  int b;                            /* Band counter                */
  int cb;                           /* Code book counter           */

  /* The following is the kurtosis discriminant function. It defines
   * the boundaries between the five quantization distributions on the
   * basis of kurtosis ).  NOTE THAT I HAVE MODIFIED THIS ARRAY: I
   * added the float versions of +infinity and -infinity to facilitate
   * clean looking interval checking--Troy
   */
  static float kdf[6] =            
  {
    (float) +1.e37,
    (float) 17.4250,
    (float) 4.8810, 
    (float) 7.8250, 
    (float) 3.3810,
    (float) -1.e37
  };

  /*
   * Loop over all components, levels, and bands, computing the
   * codebook via table lookup.
   */
  for (c=0; c < lra_stats->num_components; c++) {
    lra_comp = lra_stats->components + c;
    for (n=0; n <= lra_comp->num_levels; n++) {   
      lra_lev = lra_comp->levels + n;
      for (b=lra_lev->min_band; b <= lra_lev->max_band; b++) {
        lra_band = lra_lev->bands + b;        
        if (!lra_band->valid_band)
          continue;
        
        /* Test each codebook interval... */
        for( cb = 1; cb < 6; cb++) 
          /* If the kurtosis lies in the cb-th codebook interval */
          if ( kdf[cb] < lra_band->kurtosis  &&  lra_band->kurtosis <= kdf[cb-1] )
            /* Then we've found our codebook */
            lra_band->codebook = cb;
      }
      /* Done with bands for this level */
    }
    /* Done with levels for this component */
  }
  /* Done with components for this image */
}

/*****************************************************************************/
/* STATIC               R_of_dprime                                          */
/*****************************************************************************/
/* 
 * Uses the curve fit parameters (which depend upon TCQ vs SQ, e.g.)
 * to compute rate as a function of -dprime, the negative of the
 * derivative of the distortion (NOT dprime).
 *
 * FYI:
 * -dprime = lambda /
 * (energy_weight*color_weight*variance) for each subband.
 * 
 * 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 float 
  R_of_dprime(long alpha, float dprime, lra_model_data *model_data)
{
  double rh, rlin = 0.0, rhyper = 0.0, temp, temp2, m;

  temp = log10 (dprime);
  rh = model_data->rmh[alpha] * temp + model_data->rbh[alpha];
  if (dprime >= model_data->rlimit[alpha]){
    rlin = model_data->rmlin[alpha] * temp + model_data->rblin[alpha];
  }
  if (dprime >= model_data->rhyplim[alpha]) {
    m = -(model_data->rmh[alpha] + model_data->rmlin[alpha]);
    temp2 = (temp + model_data->rbs[alpha] - model_data->rzps[alpha]) / model_data->rbs[alpha];
    temp2 = 1 + pow (temp2, model_data->rpw[alpha]);
    rhyper = m * model_data->rbs[alpha] * ( pow (temp2, 1.0 / model_data->rpw[alpha]) - 1);
  }
  temp2 = rh + rlin + rhyper;
  return ((float) temp2);
}

/*****************************************************************************/
/* STATIC               Delta_of_dprime                                      */
/*****************************************************************************/
/* 
 * Uses the curve fit parameters (which depend on TCQ vs. SQ, e.g.)
 * to compute delta as a function of -dprime, the negative of the
 * derivative of the distortion (NOT dprime).
 *
 * FYI:
 * -dprime = lambda /
 * (energy_weight*color_weight*variance) for each subband.
 * 
 * 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 float 
  Delta_of_dprime (long alpha, float dprime, lra_model_data *model_data)
{
  double dh, dlin = 0.0, dhyper = 0.0, temp, temp2;
  double d_sum;
  int d_sum_max;

  temp = log10 (dprime);
  dh = model_data->dmh[alpha] * temp + model_data->dbh[alpha];
  if (dprime >= model_data->dlimit[alpha]){
    dlin = model_data->dmlin[alpha] * temp + model_data->dblin[alpha];
  }
  if (dprime >= model_data->dhyplim[alpha]) {
    temp2 = (temp + model_data->dbs[alpha] - model_data->dzps[alpha]) / model_data->dbs[alpha];
    temp2 = 1 + pow (temp2, model_data->dpw[alpha]);
    dhyper = model_data->dmhyp[alpha] * model_data->dbs[alpha] * (pow (temp2, 1.0 / model_data->dpw[alpha]) - 1);
  }
  d_sum_max = (int)(log10(FLT_MAX) - 1);
  d_sum = dh + dlin + dhyper;
  if (d_sum > d_sum_max) d_sum = d_sum_max; /* Avoid overflow */
  return ((float) pow (10.0, d_sum));
}

/*****************************************************************************/
/* STATIC                min_subb_functionals                                */
/*****************************************************************************/
/* 
 * Minimizes the Lagrangian functional in the Shoham & Gersho
 * algorithm for all subbands.  Ultimately, this boils down to the
 * following simple computation for each subband:
 *
 *     -dprime = lambda / (energy_weight*color_weight*variance) 
 *
 * These values are stored as a field in the lra_quantizer_obj itself. 
 */
static void 
  min_subb_functionals (lra_stats_ref lra_stats, float lambda)
{
  lra_component_info_ptr lra_comp;  /* Convenience pointer         */
  lra_level_info_ptr     lra_lev;   /* Convenience pointer         */
  lra_band_info_ptr      lra_band;  /* Convenience pointer         */
  int c;                            /* Component level counter     */
  int n;                            /* Level counter               */
  int b;                            /* Band counter                */

  /*
   * Loop over all components, levels, and bands, computing the
   * -dprime values.
   */
  for (c=0; c < lra_stats->num_components; c++) {
    lra_comp = lra_stats->components + c;
    for (n=0; n <= lra_comp->num_levels; n++) {   
      lra_lev = lra_comp->levels + n;
      for (b=lra_lev->min_band; b <= lra_lev->max_band; b++) {
        lra_band = lra_lev->bands + b;        
        if (!lra_band->valid_band)
          continue;