lra.c

JPEG2000实现的源码

        /* The actual work */
        lra_band->neg_dprime = lambda / ( lra_band->energy_weight * 
                                          lra_band->variance);
        /* Clamp neg_dprime from above at 1000 */
        if( lra_band->neg_dprime > 1000.0)
          lra_band->neg_dprime = 1000.0;
      }
      /* Done with bands for this level */
    }
    /* Done with levels for this component */
  }
  /* Done with components for this image */
}

/*****************************************************************************/
/* STATIC                compute_constraint                                  */
/*****************************************************************************/
/* 
 * Given the neg_dprime values (which should have been computed in a
 * previous call to min_subb_functionals()) and size weights, this
 * function determines the theoretical total rate, or what some folks
 * might call the expected rate.  
 */
static float 
  compute_constraint (lra_stats_ref lra_stats, lra_model_data *model_data)
{
  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                */
  float rate = 0;                   /* The total rate for the image */

  /*
   * Loop over all components, levels, and bands,
   * 
   */
  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;

        /* Estimate the rate for this subband */
        lra_band->rate = 
          R_of_dprime(lra_band->codebook - 1, /* "-1": Zero based arrays */
                      lra_band->neg_dprime, model_data);

        /* Threshold rates which are too small to consider (why?) */
        if(lra_band->rate < LRA_MINALLOWABLERATE)
          lra_band->rate = 0.0;

        /* Accumulate the expected rate for the entire image */
        rate += lra_band->rate * lra_band->size_weight;
      }
      /* Done with bands for this level */
    }
    /* Done with levels for this component */
  }
  /* Done with components for this image */

  /* Return the rate in bpp */
  return(rate);
}

/*****************************************************************************/
/* STATIC                do_deltas                                           */
/*****************************************************************************/
/* 
 * Determine the step sizes based on the -dprime values (which depend
 * upon the lambda value which should yield the desired rate).
 */
static void
  do_deltas( lra_stats_ref lra_stats, lra_model_data *model_data)
{
  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                */
  float sigma;                      /* For debug */
  float delta;                      /* for debug */
  float min_step, max_step;

  max_step = exponent_mantissa_to_step2(-31, 0);
  min_step = exponent_mantissa_to_step2(0, 0);

  /*
   * Loop over all components, levels, and bands,
   * 
   */
  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;

        delta = Delta_of_dprime(lra_band->codebook - 1, /* "-1": Zero based arrays */
                                lra_band->neg_dprime, model_data);
        sigma = (float) sqrt( lra_band->variance );
        lra_band->step_size = delta * sigma;
        lra_band->step_size *= 2.0F;
        if ((lra_band->step_size > max_step) || (lra_band->step_size == 0.0F)) 
          lra_band->step_size = max_step * 2.0F;
        if (lra_band->step_size < min_step) lra_band->step_size = min_step;
      }
      /* Done with bands for this level */
    }
    /* Done with levels for this component */
  }
  /* Done with components for this image */
}

/*****************************************************************************/
/* STATIC                allocate_rate                                       */
/*****************************************************************************/
/* 
 * Given the file rate required, this function estimates the step sizes.
 */
static void
  allocate_rate(lra_stats_ref lra_stats, float desired_rate, lra_model_data *model_data)
{
  double lambda1;                  /* Current upper bound on lambda    */
  double lambda2;                  /* Midpoint of lambda1 and lambda3  */
  double lambda3;                  /* Current lower bound on lambda    */
  double rate1;                    /* Rate corresponding to lambda1    */
  double rate2;                    /* Rate corresponding to lambda2    */
  double rate3;                    /* Rate corresponding to lambda3    */
  double tolerance  = 0.00000001F;  /* Rate tolerance                   */
  double lambda_tol = 0.00000001F;  /* lambda tolerance                 */
  long   lambda_found = 0;         /* Have we found lambda? Flag       */
  long   count=0;                  /* Bisection search iteration count */

  /* 
   * Compute upper bound on lamba 
   */
  lambda1 = 10.0;
  min_subb_functionals(lra_stats, (float) lambda1);
  rate1 = compute_constraint(lra_stats, model_data);
  while (rate1 > desired_rate) {
    lambda1 *= 2;
    min_subb_functionals(lra_stats, (float) lambda1);
    rate1 = compute_constraint(lra_stats, model_data);
    if(lambda1 > 1.0e+30) break;
  }
  
  /* 
   * Compute lower bound on lambda 
   */
  lambda3 = 0.1;
  min_subb_functionals(lra_stats, (float) lambda3);
  rate3 = compute_constraint(lra_stats, model_data);
  while (rate3 <= desired_rate) {
    lambda3 /= 2;
    min_subb_functionals(lra_stats, (float) lambda3);
    rate3 = compute_constraint(lra_stats, model_data);
    if(lambda3 < 1.0e-10) break;
  }

  /*
   * Bisection search to find the lambda which yields required rate 
   */
  lambda2 = (lambda1 + lambda3) / 2.0F;
  min_subb_functionals(lra_stats, (float) lambda2);
  rate2 = compute_constraint(lra_stats, model_data);

  /* This isn't usually part of a bisection search, but since we've
   * already done the work to compute the rate, we might as well check
   * to see if we're done...  */
  if ((fabs (rate2 - desired_rate) <= tolerance) && (desired_rate >= rate2)) {
    /* Compute the step sizes */
    do_deltas(lra_stats, model_data);
    return;
  }
  do {
    /* Pick the upper half interval as our next test interval */
    if (desired_rate < rate2) {
      lambda3 = lambda2;
      lambda2 = (lambda1 + lambda3) / 2.0F;
      min_subb_functionals(lra_stats, (float) lambda2);
      rate2 = compute_constraint(lra_stats, model_data);
    }
    /* Pick the lower half interval as our next test interval */
    else if (desired_rate > rate2) {
      lambda1 = lambda2;
      lambda2 = (lambda1 + lambda3) / 2.0F;
      min_subb_functionals(lra_stats, (float) lambda2);
      rate2 = compute_constraint(lra_stats, model_data);
    }

    /* Are we done? */
    if ((fabs (desired_rate - rate2) <= tolerance) && (desired_rate >= rate2))
      lambda_found = 1;
    else if ((lambda1 - lambda3) < lambda_tol) {
      lambda_found = 1;
      /* Estimate lambda on high side ==> rate a little on the low side */
      lambda2 = lambda1;
      min_subb_functionals(lra_stats, (float) lambda2);
      rate2 = compute_constraint(lra_stats, model_data);
    }
    count++;

    /* To avoid infinite loops, increase tolerances */
    if (count > 1000) {  
      lambda_tol *= 10.0;
      tolerance *= 10.0;
      count = 0;
    }
  } while (!lambda_found);

  /* Compute the step sizes corresponding to the lambda determined */
  do_deltas(lra_stats, model_data);
}

/*****************************************************************************/
/* EXTERN                   find_step_sizes                                  */
/*****************************************************************************/
/* 
 * Use Lagrangian Rate Allocation to compute the step sizes, given the
 * (estimated) requested bit rate and the subband statistics 
 */
void  
find_step_sizes(lra_stats_ref lra_stats,
                float rate_target_adj,   /* Overall file rate desired - header rate */
                float rate_achieved_adj, /* Rate attained just now - header rate    */
                int   iterations,        /* Current iteration                       */
                int quantizer)           /* LRA_TCQ or LRA_SQ                       */
{
  float rate_new;              /* Estimated rate which, if fed into LRA, will   */
                               /* hopefully produce step sizes which yield a    */
                               /* rate of rate_target_adj                       */
  static float rate_new_prev = 0.0;
  static float stepSizeCkSum[200];
  lra_model_data *model_data;
  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                */
  float stepSizeCkSum2;
  static float rate_new_diff;
  static float rate_new_diff_prev;
  static float rate_new_min = 0.0;
  int exponent, mantissa;
  int stepSizeNotFound;
  int i;

  /* 
   * Setup LRA object
   */
  /* Setup the parameters (for all codebooks) for Scalar Quantization*/
  model_data = set_models(quantizer);
  /* Determine the "codebook" to be used for each subband, based on kurtosis */
  select_codebooks(lra_stats);

  /* 
   * Do the actual step size computations and save to file
   */
  /* Estimate required bit rate by affine correction from previous attempts */
  rate_new = affine_rate_estimate(rate_new_prev, rate_achieved_adj, rate_target_adj, iterations);
  if (rate_new_min == 0.0) rate_new_min = rate_new;
  if ((rate_new > 0.0) && (rate_new < rate_new_min)) rate_new_min = rate_new;
  if (rate_new < 0.0) {
    rate_new = rate_new_min/2.0F;
    rate_new_min = rate_new;
  }
  rate_new_diff_prev = rate_new_diff;
  rate_new_diff = rate_new - rate_new_prev;
  if (rate_new_diff == 0.0) rate_new_diff = rate_new_diff_prev;

  /* Compute the step sizes and store in lra_quantizer_obj */
  stepSizeNotFound = 1;
  while (stepSizeNotFound) {
    allocate_rate(lra_stats, rate_new, model_data);
    stepSizeCkSum2 = 0.0;

    /* Compute stepsize check sum as result of allocate_rate() call */
    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;
          step_to_exponent_mantissa2(lra_band->step_size, &exponent, &mantissa);
          exponent += (IMPLEMENTATION_PRECISION - 2);
          stepSizeCkSum2 += exponent_mantissa_to_step2(exponent, mantissa);
        }
        /* Done with bands for this level */
      }
      /* Done with levels for this component */
    }

    /* See if stepsizes have repeated */
    stepSizeNotFound = 0;
    for (i=0; i<iterations; i++) {
      if ((stepSizeCkSum[i] == stepSizeCkSum2) && (!stepSizeNotFound)) {
        rate_new += rate_new_diff;
        if (rate_new < 0.0) {
          rate_new = rate_new_min/2.0F;
          rate_new_min = rate_new;
        }
        stepSizeNotFound = 1;
      }
    }
  }

  rate_new_prev = rate_new;
  stepSizeCkSum[iterations] = stepSizeCkSum2;
  local_free(model_data);

}