hplx_hor_analysis_by_convolution.c

JPEG2000实现的源码

/*****************************************************************************/
/* Copyright 1998, Hewlett-Packard Company                                   */
/* All rights reserved                                                       */
/* File: "hplx_hor_analysis_by_convolution.c"                                */
/* Description: Implementation of the `hplx_analysis_stage' object for       */
/*              horizontal analysis via convolution.                         */
/* Author: David Taubman                                                     */
/* Affiliation: Hewlett-Packard and                                          */
/*              The University of New South Wales, Australia                 */
/* Acknowledgements: Partly developed in collaboration with                  */
/*                   Christos Chrysafix of HP Labs                           */
/* Version: VM9.0                                                            */
/* Last Revised: 20 April, 2001                                              */
/*****************************************************************************/

/*****************************************************************************/
/* 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.                                                  */
/*****************************************************************************/

/*****************************************************************************/
/* 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 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.                                   */
/*****************************************************************************/

/*****************************************************************************/
/* Modified by Brendt Wohlberg (Los Alamos National Laboratory) to support   */
/* use of even length filters on tiles of arbitrary dimension and offset     */
/* (as per proposal WG1N1842 submitted at WG1 Rochester meeting). All        */
/* modifications are Copyright 2000 University of California.                */
/*****************************************************************************/

/*****************************************************************************/
/* Modified by Jianxin Wei, Australian Defence Force Academy (ADFA)          */
/* to implement the Odd Tile-Size Low-Pass First Convention (OTLPF_CONVENTION)*/
/* Copyright 2000 University of New South Wales, Australia.                  */
/* All rights reserved for modified parts.                                   */
/*****************************************************************************/

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

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

/*****************************************************************************/
/*                     hplx_hor_analysis_by_convolution_obj                  */
/*****************************************************************************/

typedef
  struct hplx_hor_analysis_by_convolution_obj {
    hplx_analysis_stage_obj base;
    int tile_width, first_idx;
    int branch_neg_supports[2];
    int branch_pos_supports[2];
    float *branch_taps[2];
    int even_taps;
    int extend;
    float *in_buf;
    float *branch_bufs[2];
  } hplx_hor_analysis_by_convolution_obj,
    *hplx_hor_analysis_by_convolution_ref;
    /* Extends the basic `hplx_analysis_stage' object to implement
       horizontal analysis by convolution.
           The `tile_width' field identifies the width of the current
       tile at the input to this function.
           The `first_idx' field identifies the index of the first sample
       on each row of the current tile, relative to the absolute
       canvas coordinate system.
           The negative and positive supports for each kernel conform to
       the policy described at length in the "ifc.h" header file.  Also,
       note that the kernel taps for each branch (0 for low, 1 for high)
       appear in inner-product order (not impulse response order) and are
       identified by pointers to the centres of their respective regions
       of supports, not the beginning.
           The `extend' field hold the maximum of the low- and high-pass
       `neg_support' and `pos_support' values.  This is the distances by which
       the input row buffer must be symmetrically extended to the left and
       right.  We use only one value for both ends mainly for simplicity;
       there really is no harm in extending too far.
           The `even_taps' field identifies whether or not the low- and
       high-pass kernels have an even number of taps.  This affects the
       boundary extension policy, as outlined carefully in "ifc.h".
           The `in_buf' buffer is large enough to receive a copy of the
       input sample values, while providing extra storage at the left
       boundary (negative indices) and right boundary (indices beyond
       the end of the row) to simplify the implementation of symmetric
       extension.
           The `branch_bufs' array holds pointers to the low- and high-pass
       subband rows generated from the input buffer. */

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

/* LANL Even Length Begin */
/*****************************************************************************/
/* STATIC                      analysis_end_filter                           */
/*****************************************************************************/

/* Implementation of (L,L) -> (L',H') postprocessing */
static void analysis_end_filter(float* sp, float* dpL, float* dpH) {
  float low, high;

  low = (sp[0] + sp[1])/2.0f;
  high = (sp[0] - sp[1]);
  *dpL = low;
  *dpH = high;
}
/* LANL Even Length End */

/*****************************************************************************/
/* STATIC                     check_kernel_symmetry                          */
/*****************************************************************************/

static void
  check_kernel_symmetry(int neg_support, int pos_support, float *taps,
                        int even_taps, int high_pass)

 /* Checks that the supplied kernel has the relevant symmetry (even or odd,
    as identified by the `even_taps' argument).  Note that we can be
    quite sure that the kernel has been normalized to have a unit
    DC or Nyquist gain, as appropriate, which allows us to establish
    reliable tollerances.  If the test fails, an appropriate error
    message is generated. */

{
  double tolerance;
  int failed;
  int n;

  tolerance = 0.0001;
  failed = 0;
  if (even_taps)
    {
      if ((neg_support-high_pass) != (pos_support+high_pass-1))
        failed = 1;
      if (high_pass && !failed)
        {
          for (n=1; n <= neg_support; n++)
            if (fabs(taps[-n]+taps[n-1]) > tolerance)
              failed = 1;
        }
      if ((!high_pass) && !failed)
        {
          for (n=1; n <= pos_support; n++)
            if (fabs(taps[n]-taps[-n+1]) > tolerance)
              failed = 1;
        }
    }
  else
    {
      if (neg_support != pos_support)
        failed = 1;
      if (!failed)
        for (n=1; n <= pos_support; n++)
          if (fabs(taps[n]-taps[-n]) > tolerance)
            failed = 1;
    }
  if (failed)
    local_error("Wavelet kernels may be implemented by convolution "
                "only if they satisfy the symmetry requirements imposed "
                "by the boundary extension policy.  All other kernels "
                "must be implemented by lifting -- i.e. the identifier "
                "supplied via the `-Fkernels' argument must be preceded "
                "with a capital `L'.");
}

/*****************************************************************************/
/*                              perform_1d_analysis                          */
/*****************************************************************************/

static void
  perform_1d_analysis(hplx_hor_analysis_by_convolution_ref self,
                      float *in_buf, float *branch_bufs[], int width,
                      int first_sample_odd)

 /* This function does most of the work of `__push_line_float' and
    `__push_line_fixed'.  It symmetrically extends the supplied input
    buffer as necessary and writes the two sets of subband samples into the
    two buffers pointed to by the entries of the `branch_bufs' array. */

{
  hplx_analysis_stage_ref base = (hplx_analysis_stage_ref)self;
  float *sp, *dp;
  int n, b;

  if (width == 1)
    { /* Generally speaking, the transforms only work for signals of length
         2 or more. */
      *(branch_bufs[first_sample_odd]) = in_buf[0];
      return;
    }

  /* Perform symmetric extension.  This is almost completely trivial, but
     we must beware of the possibility that the entire row is smaller than
     the kernel support, even down to a length of 1.  As a consequence,
     the best way to go about applying the extension is to extend each
     boundary one sample at a time. */

  sp = in_buf; dp = sp + width - 1;
  if (base->usePsExtension)
    {
      float sp0 = sp[0];
      float dp0 = dp[0];

      if (first_sample_odd)
        local_error("Cannot start image, tile or frame with an odd indexed "
                    "sample when using the PSE!  Use regular symmetric "
                    "extension!");
      for (n=1; n <= self->extend; n++)
        {
          if (self->even_taps)
            local_error("Cannot use PSE with even length filters!");
          else
            { /* Odd symmetry for odd length filters. */
              sp[-n] = 2 * sp0 - sp[n];
              dp[n]  = 2 * dp0 - dp[-n];
            }
        }
    }
  else if (self->even_taps)
    {
      /* LANL Even Length Begin */
      /* Disable odd indexing even filter error message in modified code */
      /* LANL Even Length End */

      for (n=1; n <= self->extend; n++)
        {
          sp[-n] = sp[n-1];
          dp[n] = dp[-n+1];
        }
    }
  else
    { /* The usual case of odd-length filters is the most general and also
         the simplest. */
      for (n=1; n <= self->extend; n++)
        {
          sp[-n] = sp[n];
          dp[n] = dp[-n];
        }
    }

  /* Apply the two kernels one at a time. */

  for (b=0; b < base->horiz_limit; b++)
    {
      int neg, pos, k;
      float *taps, sum;

      /* LANL Even Length Begin */
      /* Set even length filter centres at (1/2, 1/2) for even
         offset inputs */
      int offset = (self->even_taps)?b:(b^first_sample_odd);
      if (self->even_taps && !first_sample_odd) {
        if (b)
          offset++;
        else
          offset--;
      }
      /* LANL Even Length End */

      dp = branch_bufs[b];
      neg = self->branch_neg_supports[b];
      pos = self->branch_pos_supports[b];