std_kernel_info.c

JPEG2000实现的源码

/*****************************************************************************/
/* Copyright 1998, Hewlett-Packard Company                                   */
/* All rights reserved                                                       */
/* File: "std_kernel_info.c"                                                 */
/* Description: Decomposition kernel loading and manipulation tools          */
/* Author: David Taubman                                                     */
/* Affiliation: Hewlett-Packard and                                          */
/*              The University of New South Wales, Australia                 */
/* Version: VM8.6                                                            */
/* Last Revised: 1 December, 2000                                            */
/*****************************************************************************/

/*****************************************************************************/
/* Modified by David Taubman to remove convolution as a means of importing   */
/* coefficients and to restore full support for different kernel types in    */
/* different directions and resolution levels.  Convolution is now a matter  */
/* of implementation choice, but has nothing to do with representation of    */
/* kernels in the codestream or import format.                               */
/*****************************************************************************/

#include <local_services.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <assert.h>
#include <ifc.h>
#include "std_kernel_info_local.h"
#include "std_coeffs.h"

static double p1_5x3_taps[2] = {p1_5X3_0,p1_5X3_1};
static double u1_5x3_taps[2] = {u1_5X3_0,u1_5X3_1};
static double p1_9x7_taps[2] = {p1_9X7_0,p1_9X7_1};
static double u1_9x7_taps[2] = {u1_9X7_0,u1_9X7_1};
static double p2_9x7_taps[2] = {p2_9X7_0,p2_9X7_1};
static double u2_9x7_taps[2] = {u2_9X7_0,u2_9X7_1};


/* ========================================================================= */
/* -------------------------- Internal Functions --------------------------- */
/* ========================================================================= */

/*****************************************************************************/
/* STATIC                  normalize_convolution_taps                        */
/*****************************************************************************/

static float
  normalize_convolution_taps(std_tap_info_ptr taps)

 /* Normalizes the kernel taps to have a unit absolute DC or Nyquist gain,
    as appropriate and returns the amount by which the filter taps were
    scaled in order to achieve this normalization. */

{
  float scale;
  float dc, nyq, tmp;
  int n, support;

  support = taps->neg_support + taps->pos_support + 1;
  for (dc=nyq=0.0F, n=0; n < support; n++)
    {
      dc += taps->taps[n];
      nyq += (n&1)?(-taps->taps[n]):(taps->taps[n]);
    }
  tmp=nyq; /* Christos Chrysafis (HPL)*/
  nyq = (nyq<0.0F)?(-nyq):nyq;
  dc = (dc<0.0F)?(dc):dc;
  scale = (dc>nyq)?dc:nyq;
  scale = 1.0F / scale;

  /* Fix for LH/HL discrepancy between VM8.5 and FDIS
  if ((tmp<0.0F) && (nyq>dc))
    scale = -scale;
  End of fix */

  for (n=0; n < support; n++)
    taps->taps[n] *= scale;
  return(scale);
}

/*****************************************************************************/
/* STATIC                    get_synthesis_from_analysis                     */
/*****************************************************************************/

static void
  get_synthesis_from_analysis(std_tap_info_ptr low, std_tap_info_ptr high)

 /* Fills in the low- and high-pass synthesis taps structures, which reside in
    the second entries of the 2-element arrays, `low' and `high', based on
    the low- and high-pass analysis taps, which reside in the first entries
    of the same arrays. */

{
  std_tap_info_ptr a_low, a_high, s_low, s_high;
  int num_taps, i;
  float *sp, *dp, a_low_dc, a_high_dc, s_low_dc, s_high_dc, scale;

  a_low = low; a_high = high; s_low = low+1; s_high = high+1;
  assert((a_low->taps != NULL) && (a_high->taps != NULL) &&
         (s_low->taps == NULL) && (s_high->taps == NULL));

  /* 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(INFO_MEM_KEY,sizeof(float)*(size_t) num_taps);
  sp = a_low->taps + a_low->neg_support;
  dp = s_high->taps + s_high->neg_support;
  a_low_dc = s_high_dc = 0.0F;
  for (i=-s_high->neg_support; i <= s_high->pos_support; i++)
    {
      dp[i] = sp[-i] * ((i&1)?-1.0F:1.0F);
      a_low_dc += sp[-i];
      s_high_dc += dp[i];
    }

  /* 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(INFO_MEM_KEY,sizeof(float)*(size_t) num_taps);
  sp = a_high->taps + a_high->neg_support;
  dp = s_low->taps + s_low->neg_support;
  a_high_dc = s_low_dc = 0.0F;
  for (i=-s_low->neg_support; i <= s_low->pos_support; i++)
    {
      dp[i] = sp[-i] * ((i&1)?-1.0F:1.0F);
      a_high_dc += sp[-i];
      s_low_dc += dp[i];
    }

  /* Now compute and apply the appropriate scaling factors for the
     synthesis filters. */

  scale = 2.0F / (a_low_dc*s_low_dc + a_high_dc*s_high_dc);
  num_taps = s_low->neg_support + s_low->pos_support + 1;
  for (i=0; i < num_taps; i++)
    s_low->taps[i] *= scale;
  num_taps = s_high->neg_support + s_high->pos_support + 1;
  for (i=0; i < num_taps; i++)
    s_high->taps[i] *= scale;

  /* Complete other fields in the synthesis taps structures. */

  s_low->int_taps = s_high->int_taps = NULL;
  s_low->int_rdx = s_high->int_rdx = -1;
}

/*****************************************************************************/
/* STATIC                   get_convolution_from_lifting                     */
/*****************************************************************************/

static void
  get_convolution_from_lifting(std_lifting_info_ptr lifting,
                               std_tap_info_ptr low, std_tap_info_ptr high)

 /* This function computes the low- and high-pass analysis taps from the
    lifting steps given in `lifting'.  Note that the lifting step taps and
    the analysis convolution taps all appear in inner-product order, rather
    than impulse response order.   The algorithm is very simple: we simply
    apply the lifting steps to a unit impulse at 0 and a unit impulse at 1,
    which yields the even and odd sub-sequences of the relevant impulse
    responses. */

{
  float even_buf_0[128], odd_buf_0[128], even_buf_1[128], odd_buf_1[128];
  std_tap_info_ptr step_taps;
  float *taps, sum;
  int n, k, step_idx, neg, pos;
  
  for (n=0; n < 128; n++)
    even_buf_0[n] = odd_buf_0[n] = even_buf_1[n] = odd_buf_1[n] = 0.0F;

  /* Apply lifting to an impulse at 0 (index 128 in deinterleaved buffers) */

  even_buf_0[64] = 1.0F;
  for (step_idx=0; step_idx < lifting->num_steps; step_idx++)
    {
      step_taps = lifting->steps + step_idx;
      neg = step_taps->neg_support;
      pos = step_taps->pos_support;
      taps = step_taps->taps + neg;
      if (step_idx & 1)
        { /* Update step. */
          for (n=neg; n < (128-pos); n++)
            {
              for (sum=0.0F, k=-neg; k <= pos; k++)
                sum += odd_buf_0[n+k] * taps[k];
              even_buf_0[n] += sum;
            }
        }
      else
        { /* Lifting step. */
          for (n=neg; n < (128-pos); n++)
            {
              for (sum=0.0F, k=-neg; k <= pos; k++)
                sum += even_buf_0[n+k] * taps[k];
              odd_buf_0[n] += sum;
            }
        }
    }
  for (n=0; n < 128; n++)
    {
      even_buf_0[n] *= lifting->lp_scale;
      odd_buf_0[n] *= lifting->hp_scale;
    }

  /* Apply lifting to an impulse at 1 (index 129 in deinterleaved buffers) */

  odd_buf_1[64] = 1.0F;
  for (step_idx=0; step_idx < lifting->num_steps; step_idx++)
    {
      step_taps = lifting->steps + step_idx;
      neg = step_taps->neg_support;
      pos = step_taps->pos_support;
      taps = step_taps->taps + neg;
      if (step_idx & 1)
        { /* Update step. */
          for (n=neg; n < (128-pos); n++)
            {
              for (sum=0.0F, k=-neg; k <= pos; k++)
                sum += odd_buf_1[n+k] * taps[k];
              even_buf_1[n] += sum;
            }
        }
      else
        { /* Lifting step. */
          for (n=neg; n < (128-pos); n++)
            {
              for (sum=0.0F, k=-neg; k <= pos; k++)
                sum += even_buf_1[n+k] * taps[k];
              odd_buf_1[n] += sum;
            }
        }
    }
  for (n=0; n < 128; n++)
    {
      even_buf_1[n] *= lifting->lp_scale;
      odd_buf_1[n] *= lifting->hp_scale;
    }

  /* Extract the low-pass convolution kernel. */

  neg = pos = 0; /* These will become the signed impulse response support. */
  for (n=0; n < 128; n++)
    {
      if (even_buf_0[n] != 0.0F)
        {
          k = n-64; k = k+k;
          if (k < neg)
            neg = k;
          else if (k > pos)
            pos = k;
        }
      if (even_buf_1[n] != 0.0F)
        {
          k = n-64; k = k+k-1;
          if (k < neg)
            neg = k;
          else if (k > pos)
            pos = k;
        }
    }
  neg = -neg;
  low->neg_support = pos;
  low->pos_support = neg;
  low->taps = (float *)
    local_malloc(INFO_MEM_KEY,sizeof(float)*(neg+pos+1));
  taps = low->taps + low->neg_support;
  for (n=-neg; n <= pos; n++)
    {
      if (n & 1)
        taps[-n] = even_buf_1[64+((n+1)>>1)];
      else
        taps[-n] = even_buf_0[64+(n>>1)];
    }
  low->int_rdx = -1;
  low->int_taps = NULL;

  /* Extract the high-pass convolution kernel. */

  neg = pos = 0; /* These will become the signed impulse response support. */
  for (n=0; n < 128; n++)
    {
      if (odd_buf_0[n] != 0.0F)
        {
          k = n-64; k = k+k+1;
          if (k < neg)
            neg = k;
          else if (k > pos)
            pos = k;
        }
      if (odd_buf_1[n] != 0.0F)
        {
          k = n-64; k = k+k;
          if (k < neg)
            neg = k;
          else if (k > pos)