kernels.cpp

JPEG2000的C++实现代码

/*****************************************************************************/
// File: kernels.cpp [scope = CORESYS/DWT-KERNELS]
// Version: Kakadu, V2.2
// Author: David Taubman
// Last Revised: 20 June, 2001
/*****************************************************************************/
// Copyright 2001, David Taubman, The University of New South Wales (UNSW)
// The copyright owner is Unisearch Ltd, Australia (commercial arm of UNSW)
// Neither this copyright statement, nor the licensing details below
// may be removed from this file or dissociated from its contents.
/*****************************************************************************/
// Licensee: Book Owner
// License number: 99999
// The Licensee has been granted a NON-COMMERCIAL license to the contents of
// this source file, said Licensee being the owner of a copy of the book,
// "JPEG2000: Image Compression Fundamentals, Standards and Practice," by
// Taubman and Marcellin (Kluwer Academic Publishers, 2001).  A brief summary
// of the license appears below.  This summary is not to be relied upon in
// preference to the full text of the license agreement, which was accepted
// upon breaking the seal of the compact disc accompanying the above-mentioned
// book.
// 1. The Licensee has the right to Non-Commercial Use of the Kakadu software,
//    Version 2.2, including distribution of one or more Applications built
//    using the software, provided such distribution is not for financial
//    return.
// 2. The Licensee has the right to personal use of the Kakadu software,
//    Version 2.2.
// 3. The Licensee has the right to distribute Reusable Code (including
//    source code and dynamically or statically linked libraries) to a Third
//    Party, provided the Third Party possesses a license to use the Kakadu
//    software, Version 2.2, and provided such distribution is not for
//    financial return.
/******************************************************************************
Description:
   Implements the services defined by "kdu_kernels.h"
******************************************************************************/

#include <assert.h>
#include <math.h>
#include "kdu_elementary.h"
#include "kdu_messaging.h"
#include "kdu_params.h"
#include "kdu_kernels.h"


/* ========================================================================= */
/*                                kdu_kernels                                */
/* ========================================================================= */

/*****************************************************************************/
/*                          kdu_kernels::kdu_kernels                         */
/*****************************************************************************/

kdu_kernels::kdu_kernels(int kernel_id, bool reversible)
{
  this->kernel_id = kernel_id;
  this->downshifts = NULL;
  if (kernel_id == Ckernels_W5X3)
    {
      num_steps = 2;
      lifting_factors = new float[num_steps];
      lifting_factors[0] = -0.5F;
      lifting_factors[1] = 0.25F;
      if (reversible)
        {
          downshifts = new int[num_steps];
          downshifts[0] = 1;
          downshifts[1] = 2;
        }
    }
  else if (kernel_id == Ckernels_W9X7)
    {
      num_steps = 4;
      lifting_factors = new float[num_steps];
      if (reversible)
        { kdu_error e; e << "The W9X7 kernel may not be used for reversible "
          "compression!"; }
      lifting_factors[0] = (float) -1.586134342;
      lifting_factors[1] = (float) -0.052980118;
      lifting_factors[2] = (float) 0.882911075;
      lifting_factors[3] = (float)  0.443506852;
    }
  else
    { kdu_error e; e << "Illegal DWT kernel ID used to construct a "
      "`kdu_kernels' object."; }

  // Now let's derive all the remaining quantities.

  low_analysis_L   = num_steps;   // These lengths may be pessimistic
  high_analysis_L  = num_steps-1; // if one or more of the lifting factors
  low_synthesis_L  = num_steps-1; // is equal to 0.
  high_synthesis_L = num_steps;
  
  low_analysis_taps = (new float[2*low_analysis_L+1])+low_analysis_L;
  high_analysis_taps = (new float[2*high_analysis_L+1])+high_analysis_L;
  low_synthesis_taps = (new float[2*low_synthesis_L+1])+low_synthesis_L;
  high_synthesis_taps = (new float[2*high_synthesis_L+1])+high_synthesis_L;

  // Initialize the vector expansion buffers.

  max_expansion_levels = 4;
  work_L = num_steps+1; // Allow for placing the input impulse at n=1 or n=0.
  for (int d=1; d < max_expansion_levels; d++)
    work_L = work_L*2 + num_steps;
  work1 = (new float[2*work_L+1]) + work_L;
  work2 = (new float[2*work_L+1]) + work_L;
  bibo_step_gains = new double[num_steps];

  // Deduce synthesis impulse responses, without scaling factors.

  int n, k;

  for (n=0; n <= work_L; n++)
    work1[n] = work1[-n] = 0.0F;
  work1[0]=1.0F; // Simulate an impulse in the low-pass subband.
  for (k=num_steps-1; k >= 0; k--)
    {
      if (k&1)
        n = -(num_steps & (~1)); // Smallest even integer >= -num_steps.
      else
        n = ((-num_steps) & (~1)) + 1; // Smallest odd integer >= -num_steps.
      for (; n <= num_steps; n+=2)
        work1[n] -= lifting_factors[k]*(work1[n-1]+work1[n+1]);
    }
  for (n=0; n <= low_synthesis_L; n++)
    low_synthesis_taps[n] = low_synthesis_taps[-n] = work1[n];

  for (n=0; n <= work_L; n++)
    work1[n] = work1[-n] = 0.0F;
  work1[1]=1.0F; // Simulate an impulse in the high-pass subband.
  for (k=num_steps-1; k >= 0; k--)
    {
      if (k&1)
        n = -(num_steps & (~1)); // Smallest even integer >= -num_steps.
      else
        n = ((-num_steps) & (~1)) + 1; // Smallest odd integer >= -num_steps.
      for (; n <= (num_steps+1); n+=2)
        work1[n] -= lifting_factors[k]*(work1[n-1]+work1[n+1]);
    }
  for (n=0; n <= high_synthesis_L; n++)
    high_synthesis_taps[n] = high_synthesis_taps[-n] = work1[n+1];

  // Deduce analysis kernels from synthesis kernels.

  float sign_flip;
  for (sign_flip=1.0F, n=0; n <= low_analysis_L; n++, sign_flip=-sign_flip)
    low_analysis_taps[n] = low_analysis_taps[-n] =
      sign_flip*high_synthesis_taps[n];
  for (sign_flip=1.0F, n=0; n <= high_analysis_L; n++, sign_flip=-sign_flip)
    high_analysis_taps[n] = high_analysis_taps[-n] =
      sign_flip*low_synthesis_taps[n];

  // Deduce scaling factors and normalize filter taps.

  if (reversible)
    { low_scale = high_scale = 1.0F; return; }

  float gain;

  gain=low_analysis_taps[0];
  for (n=1; n <= low_analysis_L; n++)
    gain += 2*low_analysis_taps[n];
  low_scale = 1.0F / gain;
  for (n=-low_analysis_L; n <= low_analysis_L; n++)
    low_analysis_taps[n] *= low_scale;
  for (n=-low_synthesis_L; n <= low_synthesis_L; n++)
    low_synthesis_taps[n] *= gain;

  gain = high_analysis_taps[0];
  for (sign_flip=-1.0F, n=1; n <= high_analysis_L; n++, sign_flip=-sign_flip)
    gain += 2*sign_flip*high_analysis_taps[n];
  high_scale = 1.0F / gain;
  for (n=-high_analysis_L; n <= high_analysis_L; n++)
    high_analysis_taps[n] *= high_scale;
  for (n=-high_synthesis_L; n <= high_synthesis_L; n++)
    high_synthesis_taps[n] *= gain;
}

/*****************************************************************************/
/*                          kdu_kernels::~kdu_kernels                        */
/*****************************************************************************/

kdu_kernels::~kdu_kernels()
{
  if (downshifts != NULL)
    delete[] downshifts;
  delete[] lifting_factors;
  delete[] (low_analysis_taps-low_analysis_L);
  delete[] (high_analysis_taps-high_analysis_L);
  delete[] (low_synthesis_taps-low_synthesis_L);
  delete[] (high_synthesis_taps-high_synthesis_L);
  delete[] (work1-work_L);
  delete[] (work2-work_L);
  delete[] bibo_step_gains;
}

/*****************************************************************************/
/*                      kdu_kernels::get_lifting_factors                     */
/*****************************************************************************/

float *
  kdu_kernels::get_lifting_factors(int &num_steps,
                                   float &low_scale, float &high_scale)
{
  num_steps = this->num_steps;
  low_scale = this->low_scale;
  high_scale = this->high_scale;
  return lifting_factors;
}

/*****************************************************************************/
/*                     kdu_kernels::get_impulse_response                     */
/*****************************************************************************/

float *
  kdu_kernels::get_impulse_response(kdu_kernel_type which, int &half_length)
{
  switch (which) {
    case KDU_ANALYSIS_LOW:
      half_length = low_analysis_L;
      return low_analysis_taps;
    case KDU_ANALYSIS_HIGH:
      half_length = high_analysis_L;
      return high_analysis_taps;
    case KDU_SYNTHESIS_LOW:
      half_length = low_synthesis_L;
      return low_synthesis_taps;
    case KDU_SYNTHESIS_HIGH:
      half_length = high_synthesis_L;
      return high_synthesis_taps;
    default:
      assert(0);
    }
  return NULL;
}

/*****************************************************************************/
/*                        kdu_kernels::get_energy_gain                       */
/*****************************************************************************/

double
  kdu_kernels::get_energy_gain(kdu_kernel_type which, int level_idx)
{
  if (level_idx == 0)
    return (which==KDU_SYNTHESIS_LOW)?1.0:0.0;
  int extra_levels = level_idx - max_expansion_levels;
  if (extra_levels < 0)
    extra_levels = 0; 
  else
    level_idx -= extra_levels;

  int L, n, k;

  if (which == KDU_SYNTHESIS_LOW)
    {
      L = low_synthesis_L;
      for (n=-L; n <= L; n++)
        work1[n] = low_synthesis_taps[n];
    }
  else if (which == KDU_SYNTHESIS_HIGH)
    {
      L = high_synthesis_L;
      for (n=-L; n <= L; n++)
        work1[n] = high_synthesis_taps[n];
    }
  else
    assert(0); // Function only computes synthesis energy gains.
  for (level_idx--; level_idx > 0; level_idx--)
    {
      float *tbuf=work1; work1=work2; work2=tbuf;
      int new_L = 2*L + low_synthesis_L;
      assert(new_L <= work_L);
      for (n=-new_L; n <= new_L; n++)
        work1[n] = 0.0F;
      for (n=-L; n <= L; n++)
        for (k=-low_synthesis_L; k <= low_synthesis_L; k++)
          work1[2*n+k] += work2[n]*low_synthesis_taps[k];
      L = new_L;
    }

  double val, energy = 0.0;
  for (n=-L; n <= L; n++)
    {
      val = work1[n];
      energy += val*val;
    }
  while (extra_levels--)
    energy *= 2.0;

  return energy;
}

/*****************************************************************************/
/*                         kdu_kernels::get_bibo_gains                       */
/*****************************************************************************/

double *
  kdu_kernels::get_bibo_gains(int level_idx,
                              double &low_gain, double &high_gain)
{
  if (level_idx == 0)
    {
      low_gain = 1.0;
      high_gain = 0.0;
      return NULL;
    }

  if (level_idx > max_expansion_levels)
    level_idx = max_expansion_levels;
  float *work_low=work1, *work_high=work2;

  // In the seqel, `work_low' will hold the analysis kernel used to compute
  // the even sub-sequence entry at location 0, while `work_high' will hold the
  // analysis kernels used to compute the odd sub-sequence entry at location
  // 1.  The lifting procedure is followed to alternately update these
  // analysis kernels.

  int k, lev, low_L, high_L, gap;

  // Initialize analysis vectors and gains for a 1 level lazy wavelet
  for (k=-work_L; k <= work_L; k++)
    work_low[k] = work_high[k] = 0.0F;
  work_low[0] = 1.0F;
  low_L = high_L = 0;
  low_gain = high_gain = 1.0;

  for (gap=1, lev=1; lev <= level_idx; lev++, gap<<=1)
    { // Work through the levels
      /* Copy the low analysis vector from the last level to the high analysis
         vector for the current level. */
      for (k=0; k <= low_L; k++)
        work_high[k] = work_high[-k] = work_low[k];
      for (; k <= high_L; k++)
        work_high[k] = work_high[-k] = 0.0F;
      high_L = low_L;
      high_gain = low_gain;

      for (int step=0; step < num_steps; step+=2)
        { // Work through the lifting steps in this level
          float factor;
          // Start by updating the odd sub-sequence analysis kernel
          factor = lifting_factors[step];
          assert((low_L+gap) <= work_L);
          for (k=-low_L; k <= low_L; k++)
            {
              work_high[k-gap] += work_low[k]*factor;
              work_high[k+gap] += work_low[k]*factor;
            }
          high_L = ((low_L+gap) > high_L)?(low_L+gap):high_L;
          for (high_gain=0.0, k=-high_L; k <= high_L; k++)
            high_gain += fabs(work_high[k]);
          bibo_step_gains[step] = high_gain;
          
          // Now update the even sub-sequence analysis kernel
          if ((step+1) < num_steps)
            {
              factor = lifting_factors[step+1];
              assert((high_L+gap) <= work_L);
              for (k=-high_L; k <= high_L; k++)
                {
                  work_low[k-gap] += work_high[k]*factor;
                  work_low[k+gap] += work_high[k]*factor;
                }
              low_L = ((high_L+gap) > low_L)?(high_L+gap):low_L;
              for (low_gain=0.0, k=-low_L; k <= low_L; k++)
                low_gain += fabs(work_low[k]);
              bibo_step_gains[step+1] = low_gain;
            }
        }

      // Now incorporate the subband scaling factors
      for (k=-high_L; k <= high_L; k++)
        work_high[k] *= high_scale;
      high_gain *= high_scale;
      for (k=-low_L; k <= low_L; k++)
        work_low[k] *= low_scale;
      low_gain *= low_scale;
    }

  return bibo_step_gains;
}