hplx_vert_analysis_by_convolution.c

JPEG2000实现的源码

/*****************************************************************************/
/* Copyright 1998, Hewlett-Packard Company                                   */
/* All rights reserved                                                       */
/* File: "hplx_vert_analysis_by_convolution.c"                               */
/* Description: Implementation of the `hplx_analysis_stage' object for       */
/*              vertical 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                                              */
/*****************************************************************************/

/*****************************************************************************/
/* 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_vert_convolution.h"
#include "hplx_analysis_local.h"

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

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

/*****************************************************************************/
/*                      hplx_vert_analysis_by_lifting_obj                    */
/*****************************************************************************/

typedef
  struct hplx_vert_analysis_by_convolution_obj {
    hplx_analysis_stage_obj base;
    int cols;
    int branch_rows[2];
    int branch_next_row_idx[2];
    int branch_neg_supports[2];
    int branch_pos_supports[2];
    float *branch_taps[2];
    int even_taps;
    hplx_rolling_buffer window;
    float *out_buf[2];
    float *ps_buf;
    int first_row_idx;
    /* LANL Even Length Begin */
    /* Flag indicating whether the 2nd last lowpass row has been produced.
       Required for postprocessing at the bottom */
    int second_last_low_seen;
    /* Flag indicating that the last highpass row has been generated
       by the postprocessing stage */
    int final_high_ready;
    /* LANL Even Length End */
  } hplx_vert_analysis_by_convolution_obj,
    *hplx_vert_analysis_by_convolution_ref;
  /* Extends the basic `hplx_analysis_stage' object for the purpose of
     vertical analysis by convolution.
         The `cols' field identifies the width of rows processed by this
     analysis stage.
         The `branch_rows' array holds the height of the low-pass and
     high-pass subbands, respectively, in the current tile.
         The `branch_next_row_idx' array holds the zero-based indices of
     the next row to be produced for each of the two subbands.
         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 `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 `window' structure manages a rolling buffer with a window
     into the source image.
         The `out_buf' array points to a working buffer which has
     at least `cols' entries.  This buffer collects the output samples for
     a newly generated low- or high-pass subband row before they are sent
     on to the next analysis stage or the quantizer, as appropriate.
         The `first_row_idx' field holds the absolute index of the
     first row, in the canvas coordinate system. */

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

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

/* Implementation of (L,L) -> (L',H') postprocessing */
static void analysis_end_filter(float* sp0, float* sp1,
                                float* dpL, float* dpH,
                                int ncols) {
  int n;
  float low, high;

  for (n = 0; n < ncols; n++) {
    low = (sp0[n] + sp1[n])/2.0f;
    high = (sp0[n] - sp1[n]);

    dpL[n] = low;
    dpH[n] = 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'.");
}

/*****************************************************************************/
/* STATIC                      complete_push_line                            */
/*****************************************************************************/

static void
  complete_push_line(hplx_analysis_stage_ref base)

 /* This function does all the work of `__push_line_float' or
    `__push_line_fixed', after the new row of samples has been entered
    into the rolling buffer.  It iteratively tries to construct low- and
    high-pass output rows and push them downstream, until nothing more
    can be generated. */

{
  hplx_vert_analysis_by_convolution_ref self =
    (hplx_vert_analysis_by_convolution_ref) base;
  hplx_rolling_buffer_ptr window;
  float *taps, **source, *dp, sum;
  int b, n, k, done, neg, pos, support, row_idx;
  int first_row_odd;

  /* LANL Even Length Begin */
  /* Flags indicating input length and offset parity. Required for
     even length filters */
  int oddoddcase = 0;
  int evenevencase = 0;
  int oddevencase = 0;
  /* LANL Even Length End */

  /* OTLPF_CONVENTION begin; JX Wei ADFA, WJ Zeng Sharp */
  if(!base->otlpf_convention)
    first_row_odd = self->first_row_idx & 1;
  else
    first_row_odd =0;
  /** always treat it as if it were even. **/
  /* OTLPF_CONVENTION end; JX Wei ADFA, WJ Zeng Sharp */

  window = &(self->window);

  /* LANL Even Length Begin */
  /* Flags indicating input length and offset parity. Required for
     even length filters */
  oddoddcase = self->even_taps && window->rows%2 && first_row_odd;
  evenevencase = self->even_taps && !(window->rows%2) && !first_row_odd;
  oddevencase = self->even_taps && window->rows%2 && !first_row_odd;
  /* LANL Even Length End */

  /* LANL Even Length Begin */
  /* Disable odd indexing even filter error message in modified code */
  if (first_row_odd && base->usePsExtension)
    local_error("Cannot start image, tile or frame with an odd indexed "
                "sample when using the PSE!  Use "
                "regular symmetric extension with odd length filters!");
  /* Ensure PSE isn't used with even length filters */
  if (self->even_taps && base->usePsExtension)
    local_error("Cannot use PSE with even length filters");
  /* LANL Even Length End */

  if (window->rows == 1)
    { /* Signals of length 1 must be treated specially. */
      base->branches[first_row_odd]->push_line_float(
            base->branches[first_row_odd],window->buffer[0],self->cols);
      return;
    }
  do {
    int continued[2];
    float *bottom = 0;

    done = 1;
    for (b=0; b < base->vert_limit; b++)
      {
        row_idx = self->branch_next_row_idx[b];
        /* LANL Even Length Begin */
        /* Allow for slightly different `row_idx' behaviour in new
           effective filter alignments */
        if (row_idx >= self->branch_rows[b])
        /* LANL Even Length End */
          {
            continued[b] = 1;
            continue;
          }
        /* LANL Even Length Begin */
        /* Correct the row index calculations for even length filters */
        if (self->even_taps)
          row_idx = (row_idx << 1) + b;
        else
          row_idx = (row_idx << 1) + (b ^ first_row_odd);

        /* Kludge to get correct effective filter alignment for even
         length filters with filter centres at (1/2,1/2) */
        if (oddevencase || evenevencase) {
          if (b)
            row_idx++;
          else
            row_idx--;
        }
        /* LANL Even Length End */
        neg = self->branch_neg_supports[b];
        pos = self->branch_pos_supports[b];
        support = neg + pos + 1;
        source =
          hplx_rolling_buffer__get_rows(window,row_idx-neg,support,NULL,b,1,
                                        base->usePsExtension,NULL);
        if (source == NULL)
          {
            continued[b] = 1;
            continue;
          }
        done = 0;
        taps = self->branch_taps[b] - neg;
        for (dp=self->out_buf[b], n=0; n<self->cols; n++)
          {
            for (sum=0.0F, k=0; k < support; k++)
              sum += taps[k] * source[k][n];
            *(dp++) = sum;
          }

        if (b == 0)
          bottom = source[neg+1];

        continued[b] = 0;
      }

    /* Begin PSE stuff */
    if (base->usePsExtension)
      {