hplx_hor_synthesis_by_lifting.c

JPEG2000实现的源码

/*****************************************************************************/
/* Copyright 1998, Hewlett-Packard Company                                   */
/* All rights reserved                                                       */
/* File: "hplx_hor_synthesis_by_lifting.c"                                   */
/* Description: Implementation of the `hplx_synthesis_stage' object for      */
/*              horizontal analysis via lifting.                             */
/* Author: David Taubman                                                     */
/* Affiliation: Hewlett-Packard and                                          */
/*              The University of New South Wales, Australia                 */
/* Acknowledgements: Partly developed in collaboration with                  */
/*                   Christos Chrysafis of HP Labs                           */
/* Version: VM8.6                                                            */
/* Last Revised: 17 October, 2000                                            */
/*****************************************************************************/

/*****************************************************************************/
/* Modifications tagged with comment lines or delimiters involving the       */
/* string, "David T mod" or "mod by David T", have been made by David        */
/* Taubman; they are copyrighted by Hewlett-Packard Company with all rights  */
/* reserved for the modified parts.                                          */
/*****************************************************************************/

/*****************************************************************************/
/* Modifications tagged with comment lines or delimiters involving the       */
/* string "BBDWT mod" have been made to implement the single-sample          */
/* overlap Wavelet transform within the frames paradigm -- this is the       */
/* BBDWT transform proposed by Canon Research France and accepted at the     */
/* meeting in Seoul, Korea.                                                  */
/*****************************************************************************/

/*****************************************************************************/
/* Modifications marked "Christos Chrysafis (HPL)" made to ensure that       */
/* the boundary extension policy for lifting is identical to that used for   */
/* convolution so that convolution vs. lifting is purely an implementation   */
/* decision.  Mods might cause problems for non-symmetric kernels, but that  */
/* should not be relevant to the standard.  The modifications are            */
/* copyrighted by Hewlett-Packard Company with all rights reserved for the   */
/* modified parts.                                                           */
/*****************************************************************************/

/*****************************************************************************/
/* Modified by David Taubman to support arbitrary reference points for the   */
/* transform and the various regular partitions, so as to facilitate         */
/* cropping and geometric transformations in the compressed domain and to    */
/* enable full support for CRF's single-sample overlap Wavelet transform,    */
/* as originally documented in Seoul.  Changes are too numerous to flag      */
/* individually within the code.  Changes copyrighted by HP with all rights  */
/* reserved for the modified parts.                                          */
/*****************************************************************************/

/*****************************************************************************/
/* Modified by David Taubman to support interface modifications, arbitrary   */
/* changes in coding parameters from component to component and from tile    */
/* to tile, to support the full generality of PART-1 of the JPEG2000         */
/* standard, and to support most anticipated generality of PART-2.  Changes  */
/* are too numerous to flag individually within the code, which in some      */
/* places has been completely rewritten.  All changes copyrighted by HP with */
/* all rights reserved for the modified parts.                               */
/*****************************************************************************/

/*****************************************************************************/
/* Modified by David Taubman to allow SSO-DWT style extension at tile        */
/* boundaries.  Changes copyrighted by HP with all rights reserved for the   */
/* modified parts.                                                           */
/*****************************************************************************/

/*****************************************************************************/
/* Modified for general wavelet decompositions.                              */
/* Copyright 2000 Science Applications International Corporation (SAIC).     */
/* Copyright 1995 University of Arizona, Arizona Board of Regents.           */
/* All Rights Reserved for modified parts.                                   */
/*****************************************************************************/

#include <local_services.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#include <ifc.h>
#include "hplx_synthesis_local.h"

/* ========================================================================= */
/* ---------------------- Definition of Extended Object -------------------- */
/* ========================================================================= */

/*****************************************************************************/
/*                      hplx_hor_synthesis_by_lifting_obj                    */
/*****************************************************************************/

typedef
  struct hplx_hor_synthesis_by_lifting_obj {
    hplx_synthesis_stage_obj base;
    int tile_width, first_idx;
    int kernel_type;
    int num_steps;
    int neg_supports[INFO__MAX_LIFTING_STEPS];
    int pos_supports[INFO__MAX_LIFTING_STEPS];
    float *float_taps[INFO__MAX_LIFTING_STEPS]; /* Only for float-lifting. */
    int *int_taps[INFO__MAX_LIFTING_STEPS]; /* Only for int-lifting. */
    int int_downshifts[INFO__MAX_LIFTING_STEPS]; /* Only for int-lifting. */
    float step_dc_gains[INFO__MAX_LIFTING_STEPS]; /* BBDWT mod by David T */
    float branch_scale[2]; /* Only for float-lifting. */
    int neg_extend, pos_extend;
    float *branch_buf_float[2];
    ifc_int *branch_buf_int[2];
    int num_frames, frame_cols, frame_ref, first_frame_cols;
    float *frame_work_buf_float[2];
    ifc_int *frame_work_buf_int[2];
  } hplx_hor_synthesis_by_lifting_obj, *hplx_hor_synthesis_by_lifting_ref;
  /* The `neg_extend' and `pos_extend' values are the same for all tiles;
     they are the maximum negative and positive support extents, respectively,
     over all steps and are, therefore, safe extension distances for the
     sub-sequences subjected to lifting steps.  Each of the even and odd
     sub-sequence buffers (i.e. the buffers referenced by the entries of
     `branch_buf_float' or `branch_buf_int') has an extra `neg_extend'
     entries which may be indexed below 0 and `pos_extend' entries which
     may be indexed beyond the nominal end of the array (i.e. after
     `branch_cols'[0] or `branch_cols'[1], as appropriate).
         The `num_frames' field holds the number of frames.  The SSO-DWT
     is applied if and only if `num_frames' is not equal to 1.  Even
     though the SSO-DWT is an overlapping transform, the number of frames
     is computed based upon a disjoint partition.  The overlap samples
     are then determined subsequently.  This is important, since it can
     happen that the overlap sample is an entire frame.
         The `frame_cols' field holds the nominal horizontal frame
     dimension; it is meaningful only if `num_frames' > 1.
         The `frame_ref' field holds the horizontal anchor coordinate
     for the frame partition.
         `first_frame_cols' holds the width of the first frame in the current
     tile. */


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

/*****************************************************************************/
/* STATIC                         undo_float_step                            */
/*****************************************************************************/

static void
  undo_float_step(int the_stage,
                  int neg_support, int pos_support, float *taps,
                  float *src, float *dst, int src_cols, int dst_cols,
                  int src_offset, int src_on_left_boundary,
                  int src_on_right_boundary, int skip_first, int skip_last,
                  float stage_dc_gain)

 /* Undoes a single lifting step, removing a contribution due to the
    `src' sequence from the `dst' sequence.  The elements of the `taps'
    array have indices running from -`neg_support' to +`pos_support'.  These
    supports refer to the nominal alignment of the `src' and `dst'
    sub-sequences and are independent of particular boundary locations.
    It can happen that one of the sub-sequences is truncated by a boundary
    while the other is not.  In this case, `src_offset' identifies the
    amount by which the `src' sub-sequence is displaced to the right of
    the `dst' sub-sequence.  It is always a value in the range -1 to 1.  The
    value of `src_offset' should simply be subtracted from `pos_support' and
    added to `neg_support'.  The `src_on_left_boundary' flag indicates
    whether or not the first sample of the `src' buffer lies on the left
    hand boundary of the transform block.  Similarly, `src_on_right_boundary'
    indicates whether or not the last sample of the `src' buffer lies on
    the right hand boundary of the transform block.  This information tells us
    how to perform symmetric extension on the lifting steps.  The
    `skip_first' flag indicates whether or not the lifting step which
    updates the first sample in `dst' should be skipped.  Similarly, the
    `skip_last' flag indicates whether or not the lifting step which updates
    the last sample in `dst' should be skipped.  Skipped lifting steps should
    be replaced by the removal of the `stage_dc_gain' factor from the
    relevant sample in the `dst' array (i.e. division).  Lifting steps must
    sometimes be skipped in this way to perform the SSO-DWT. */

{
  int even_taps, n, p;
  int left_ext, right_ext;
  float sum;
  float *sp, *dp;

  assert(src_cols && dst_cols);
  even_taps = (pos_support+neg_support+1)&1;

  /* Determine extension types.  Note that the extensions for even
     length filters might not replicate the effects of symmetric extension
     in convolution implementations in all cases. */
  if (even_taps)
    { /* Extension for even tap filters. */
      if (src_offset)
        local_error("Cannot start image, tile or frame with an odd indexed "
                    "sample when using even length filters!  Use odd length "
                    "filters!");
      left_ext = symmetric_extension_type(1,(src_cols+dst_cols)&1,
                                          the_stage&1,0);
      right_ext = symmetric_extension_type(1,(src_cols+dst_cols)&1,
                                           the_stage&1,1);
    }
  else
    { /* Extension for odd tap filters. */
      pos_support-=src_offset;
      neg_support+=src_offset;
      taps += src_offset;
      left_ext = (src_on_left_boundary)?0:1;
      right_ext = (src_on_right_boundary)?0:1;
    }

  /* Perform extensions, alternating between left and right extension
     until we are done. */
  sp = src; dp = sp + src_cols - 1;
  pos_support += dst_cols - src_cols;
  for (n=p=1; (n <= neg_support) || (p <= pos_support); )
    {
      int did_something;

      did_something = 0;
      if (n <= neg_support)
        { /* Still need extension on the left. */
          if (left_ext)
            { /* Half-point extension. */
              if (n <= (src_cols+p-1))
                { sp[-n] = sp[n-1]*left_ext; n++; did_something = 1; }
            }
          else
            { /* Nice case for odd length filters; should always work. */
              if (n < (src_cols+p-1))
                { sp[-n] = sp[n]; n++; did_something = 1; }
            }
        }
      if (p <= pos_support)
        { /* Still need extension on the right. */
          if (right_ext)
            { /* Half-point extension. */
              if (p <= (src_cols+n-1))
                { dp[p] = dp[-p+1]*right_ext; p++; did_something = 1; }
            }
          else
            { /* Whole-point extension. */
              if (p < (src_cols+n-1))
                { dp[p] = dp[-p]; p++; did_something = 1; }
            }
        }
      assert(did_something);
    }
  pos_support -= dst_cols - src_cols;

  /* Apply lifting step.  Note that `taps' points to "centre" of support. */

  if (skip_first)
    { *dst /= stage_dc_gain; dst_cols--; dst++; src++; }
  if (skip_last)
    { dst[dst_cols-1] /= stage_dc_gain; dst_cols--; }
  for (; dst_cols > 0; dst_cols--, src++, dst++)
    {
      for (sum=0.0F, n=-neg_support; n <= pos_support; n++)
        sum += src[n] * taps[n];
      *dst -= sum;
    }
}

/*****************************************************************************/
/* STATIC                          undo_int_step                             */
/*****************************************************************************/

static void
  undo_int_step(int the_stage,
                int neg_support, int pos_support, int *taps, int downshift,
                ifc_int *src, ifc_int *dst, int src_cols, int dst_cols,
                int src_offset, int src_on_left_boundary,
                int src_on_right_boundary, int skip_first, int skip_last)

 /* Undoes a single lifting step, removing a contribution due to the
    `src' sequence from the `dst' sequence.  The elements of the `taps'
    array have indices running from -`neg_support' to +`pos_support'.  These
    supports refer to the nominal alignment of the `src' and `dst'
    sub-sequences and are independent of particular boundary locations.
    It can happen that one of the sub-sequences is truncated by a boundary
    while the other is not.  In this case, `src_offset' identifies the
    amount by which the `src' sub-sequence is displaced to the right of
    the `dst' sub-sequence.  It is always a value in the range -1 to 1.  The
    value of `src_offset' should simply be subtracted from `pos_support' and
    added to `neg_support'.  The `src_on_left_boundary' flag indicates
    whether or not the first sample of the `src' buffer lies on the left
    hand boundary of the transform block.  Similarly, `src_on_right_boundary'
    indicates whether or not the last sample of the `src' buffer lies on
    the right hand boundary of the transform block.  This information tells us
    how to perform symmetric extension on the lifting steps.  The
    `skip_first' flag indicates whether or not the lifting step which
    updates the first sample in `dst' should be skipped.  Similarly, the
    `skip_last' flag indicates whether or not the lifting step which updates
    the last sample in `dst' should be skipped.  Lifting steps must sometimes
    be skipped to perform the SSO-DWT. */


{
  int even_taps, n, p;
  int left_ext, right_ext;
  std_int sum, offset;
  ifc_int *sp,*dp;

  assert(src_cols && dst_cols);
  even_taps = (pos_support+neg_support+1)&1;

  /* Determine extension types.  Note that the extensions for even
     length filters might not replicate the effects of symmetric extension
     in convolution implementations in all cases. */
  if (even_taps)
    { /* Extension for even tap filters. */
      if (src_offset)
        local_error("Cannot start image, tile or frame with an odd indexed "
                    "sample when using even length filters!  Use odd length "
                    "filters!");
      left_ext = symmetric_extension_type(1,(src_cols+dst_cols)&1,
                                          the_stage&1,0);