synthesis.cpp

JPEG2000的C++实现代码

/*****************************************************************************/
// File: synthesis.cpp [scope = CORESYS/TRANSFORMS]
// 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 reverse DWT (subband/wavelet synthesis).  The implementation
uses lifting to reduce memory and processing, while keeping as much of the
implementation as possible common to both the reversible and the irreversible
processing paths.  The implementation is generic to the extent that it
supports any odd-length symmetric wavelet kernels -- although only 3 are
currently accepted by the "kdu_kernels" object.
******************************************************************************/

#include <assert.h>
#include <string.h>
#include <math.h>
#include "kdu_messaging.h"
#include "kdu_compressed.h"
#include "kdu_sample_processing.h"
#include "kdu_kernels.h"
#include "synthesis_local.h"

// Set things up for the inclusion of assembler optimized routines
// for specific architectures.  The reason for this is to exploit
// the availability of SIMD type instructions on many modern processors.

#if defined KDU_PENTIUM_MSVC
#  define KDU_SIMD_OPTIMIZATIONS
#  include "msvc_dwt_mmx_local.h" // Header contains all asm commands in-line
static int simd_exists = msvc_dwt_mmx_exists();
#elif defined KDU_PENTIUM_GCC
#  define KDU_SIMD_OPTIMIZATIONS
#  include "gcc_dwt_mmx_local.h" // Header declares functs in separate .s file
static int simd_exists = gcc_dwt_mmx_exists();
#endif // KDU_PENTIUM_GCC

/* ========================================================================= */
/*                               kdu_synthesis                               */
/* ========================================================================= */

/*****************************************************************************/
/*                        kdu_synthesis::kdu_synthesis                       */
/*****************************************************************************/

kdu_synthesis::kdu_synthesis(kdu_resolution resolution,
                             kdu_sample_allocator *allocator,
                             bool use_shorts, float normalization)
  // In the future, we may create separate, optimized objects for each kernel.
{
  state = new kd_synthesis(resolution,allocator,use_shorts,normalization);
}

/* ========================================================================= */
/*                              kd_synthesis                                 */
/* ========================================================================= */

/*****************************************************************************/
/*                       kd_synthesis::kd_synthesis                          */
/*****************************************************************************/

kd_synthesis::kd_synthesis(kdu_resolution resolution,
                           kdu_sample_allocator *allocator,
                           bool use_shorts, float normalization)
{
  reversible = resolution.get_reversible();
  this->use_shorts = use_shorts;
  int kernel_id = resolution.get_kernel_id();
  kdu_kernels kernels(kernel_id,reversible);
  float low_gain, high_gain;
  float *factors =
    kernels.get_lifting_factors(L_max,low_gain,high_gain);
  int n;

  assert(L_max <= 4); // We have statically sized the array to improve locality
  for (n=0; n < L_max; n++)
    {
      steps[n].augend_parity = (n+1) & 1; // Step 0 updates odd locations
      steps[n].lambda = factors[n];
      if (kernels.get_lifting_downshift(n,steps[n].downshift))
        { // Reversible case
          steps[n].i_lambda = (kdu_int32)
            floor(0.5 + steps[n].lambda*(1<<steps[n].downshift));
        }
      else
        { // Irreversible case
          steps[n].i_lambda = steps[n].downshift = 0;
          kdu_int32 fix_lambda = (kdu_int32) floor(0.5 + factors[n]*(1<<16));
          steps[n].fixpoint.fix_lambda = fix_lambda;
          steps[n].fixpoint.i_lambda = 0;
          while (fix_lambda >= (1<<15))
            { steps[n].fixpoint.i_lambda++; fix_lambda -= (1<<16); }
          while (fix_lambda < -(1<<15))
            { steps[n].fixpoint.i_lambda--; fix_lambda += (1<<16); }
          steps[n].fixpoint.remainder = (kdu_int16) fix_lambda;
          steps[n].fixpoint.pre_offset = (kdu_int16)
            floor(0.5 + ((double)(1<<15)) / ((double) fix_lambda));
        }
    }

  kdu_dims dims;
  kdu_coords min, max;

  // Get output dimensions.

  resolution.get_dims(dims);
  min = dims.pos; max = min + dims.size; max.x--; max.y--;
  y_out_next = min.y;
  y_out_max = max.y;
  x_out_min = min.x;
  x_out_max = max.x;
  
  empty = !dims;
  if (empty)
    return;

  // Get input dimensions.  Note that the following code contains assert
  // statements which may fail if the low- and high-pass filter lengths do not
  // differ by 2 (as reported by `kdu_kernels' and used in `kdu_resolution' to
  // compute the effects of kernel support on regions of interest.

  resolution.access_next().get_dims(dims); // low-pass dimensions.
  low_width = dims.size.x;
  min = dims.pos; max = min + dims.size; max.x--; max.y--;
  min.x+=min.x; min.y+=min.y; max.x+=max.x; max.y+=max.y;
  y_in_next = min.y;
  y_in_max = max.y;
  x_in_min = min.x;
  x_in_max = max.x;

  resolution.access_subband(HH_BAND).get_dims(dims); // high-pass dimensions.
  high_width = dims.size.x;
  min = dims.pos; max = min + dims.size; max.x--; max.y--;
  min.x+=min.x+1; min.y+=min.y+1; max.x+=max.x+1; max.y+=max.y+1;

  assert(min.y <= y_in_next+1);
  if (min.y < y_in_next)
    { y_in_next--; assert(min.y==y_in_next); }
  assert(max.y >= y_in_max-1);
  if (max.y > y_in_max)
    { y_in_max++; assert(max.y==y_in_max); }

  assert(min.x <= x_in_min+1);
  if (min.x < x_in_min)
    { x_in_min--; assert(min.x==x_in_min); }
  assert(max.x >= x_in_max-1);
  if (max.x > x_in_max)
    { x_in_max++; assert(max.x==x_in_max); }

  assert((y_in_next <= y_out_next) && (y_in_max >= y_out_max));
  assert((x_in_min <= x_out_min) && (x_in_max >= x_out_max));
  assert((y_in_next <= y_in_max) && (x_in_min <= x_in_max));
  
  unit_height = (y_in_next==y_in_max);
  unit_width = (x_in_min==x_in_max);

  // Pre-allocate the line buffers.

  augend.pre_create(allocator,low_width,high_width,reversible,use_shorts);
  new_state.pre_create(allocator,low_width,high_width,reversible,use_shorts);
  for (n=0; n < L_max; n++)
    steps[n].state.pre_create(allocator,low_width,high_width,
                              reversible,use_shorts);
  initialized = false; // Finalize creation in the first `push' call.

  // Now determine the normalizing upshift and subband nominal ranges.

  float LL_range, HL_range, LH_range, HH_range;
  LL_range = HL_range = LH_range = HH_range = normalization;
  normalizing_upshift = 0;
  if (!reversible)
    {
      int lev_idx = resolution.get_dwt_level(); assert(lev_idx > 0);
      double bibo_low, bibo_high, bibo_prev;
      kernels.get_bibo_gains(lev_idx-1,bibo_prev,bibo_high);
      double *bibo_steps = kernels.get_bibo_gains(lev_idx,bibo_low,bibo_high);
      double bibo_max = 0.0;

      // Find BIBO and nominal ranges for the vertical analysis transform.
      if (unit_height)
        bibo_max = normalization;
      else
        {
          LL_range /= low_gain;  HL_range /= low_gain;
          LH_range /= high_gain; HH_range /= high_gain;
          bibo_prev *= normalization; // BIBO horizontal range at stage input
          for (n=0; n < L_max; n++)
            if ((bibo_prev * bibo_steps[n]) > bibo_max)
              bibo_max = bibo_prev * bibo_steps[n];
        }
      // Find BIBO gains for horizontal analysis
      if (!unit_width)
        {
          LL_range /= low_gain;  LH_range /= low_gain;
          HL_range /= high_gain; HH_range /= high_gain;
          bibo_prev = bibo_low / low_gain; // If bounded by vertical low band
          if ((bibo_high / high_gain) > bibo_prev)
            bibo_prev = bibo_high / high_gain; // Bounded by vertical high band
          bibo_prev *= normalization; // BIBO vertical range at horiz. input
          for (n=0; n < L_max; n++)
            if ((bibo_prev * bibo_steps[n]) > bibo_max)
              bibo_max = bibo_prev * bibo_steps[n];
        }
      double overflow_limit = 1.0 * (double)(1<<(16-KDU_FIX_POINT));
          // This is the largest numeric range which can be represented in
          // our signed 16-bit fixed-point representation without overflow.
      while (bibo_max > 0.75*overflow_limit)
        { // Leave some extra headroom to allow for the effects of
          // dequantization artefacts.
          normalizing_upshift++;
          LL_range*=0.5F; LH_range*=0.5F; HL_range*=0.5F; HH_range*=0.5F;
          bibo_max *= 0.5;
        }
    }

  // Finally, create the subband interfaces.

  assert(resolution.which() > 0);
  if (resolution.which() == 1)
    hor_low[0] = kdu_decoder(resolution.access_next().access_subband(LL_BAND),
                             allocator,use_shorts,LL_range);
  else
    hor_low[0] = kdu_synthesis(resolution.access_next(),
                               allocator,use_shorts,LL_range);
  hor_high[0]  = kdu_decoder(resolution.access_subband(HL_BAND),
                             allocator,use_shorts,HL_range);
  hor_low[1]   = kdu_decoder(resolution.access_subband(LH_BAND),
                             allocator,use_shorts,LH_range);
  hor_high[1]  = kdu_decoder(resolution.access_subband(HH_BAND),
                             allocator,use_shorts,HH_range);
}

/*****************************************************************************/
/*                      kd_synthesis::~kd_synthesis                          */
/*****************************************************************************/

kd_synthesis::~kd_synthesis()
{
  hor_low[0].destroy();
  hor_low[1].destroy();
  hor_high[0].destroy();
  hor_high[1].destroy();
}

/*****************************************************************************/
/*                           kd_synthesis::pull                              */
/*****************************************************************************/

void
  kd_synthesis::pull(kdu_line_buf &line, bool allow_exchange)
{
  if (empty)
    return;

  if (!initialized)
    { // Finish creating all the buffers.
      augend.create(); augend.deactivate();
      new_state.create(); new_state.deactivate();
      for (int n=0; n < L_max; n++)
        { steps[n].state.create(); steps[n].state.deactivate(); }
      initialized = true;
    }

  int c, k;
  kd_line_cosets *out = (y_out_next & 1)?(&new_state):(&augend);

  assert(y_out_next <= y_out_max);
  if (unit_height)
    { // No transform performed in this special case.
      horizontal_synthesis(*out);
      if (reversible && (y_out_next & 1))
        { // Need to halve integer sample values.
          if (!use_shorts)
            { // Working with 32-bit data
              kdu_sample32 *dp;
              for (c=0; c < 2; c++)
                for (dp=out->cosets[c].get_buf32(),
                     k=out->cosets[c].get_width(); k--; dp++)
                  dp->ival >>= 1;
            }
          else
            { // Working with 16-bit data
              kdu_sample16 *dp;
              for (c=0; c < 2; c++)
                for (dp=out->cosets[c].get_buf16(),
                     k=out->cosets[c].get_width(); k--; dp++)
                  dp->ival >>= 1;
            }
        }
    }