hplx_vert_synthesis_by_convolution.c

JPEG2000实现的源码

/*****************************************************************************/
/* Copyright 1998, Hewlett-Packard Company                                   */
/* All rights reserved                                                       */
/* File: "hplx_vert_synthesis_by_convolution.c"                              */
/* Description: Implementation of the `hplx_synthesis_stage' object for      */
/*              vertical synthesis 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_vert_convolution.h"
#include "hplx_synthesis_local.h"

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

/*****************************************************************************/
/*                            hplx_polyphase_kernel                          */
/*****************************************************************************/

typedef
  struct hplx_polyphase_kernel {
    int neg_support, pos_support;
    float *taps;
  } hplx_polyphase_kernel, *hplx_polyphase_kernel_ptr;
  /* This structure represents a single polyphase operator, of which four
     are required to completely synthesize a 1-D signal from its low- and
     high-pass subbands.  The kernel represents the contribution of one
     of the input subbands to either the even or odd sub-sequence of the
     output.  The kernel taps appear in inner product fashion and have the
     designated region of support.  Note that the `taps' array is dynamically
     allocated and must be destroyed upon termination of the containing
     object. */

/*****************************************************************************/
/*                  hplx_vert_synthesis_by_convolution_obj                   */
/*****************************************************************************/

typedef
  struct hplx_vert_synthesis_by_convolution_obj {
    hplx_synthesis_stage_obj base;
    int cols;
    int next_row_idx;
    hplx_polyphase_kernel branch_to_even[2];
    hplx_polyphase_kernel branch_to_odd[2];
    int even_taps;
    hplx_rolling_buffer branch_windows[2];
    float *conversion_buf;
    float *taps_copy;
    int *sign_flips_buffer;
    float *ps_prev_buf;
    float *ps_low_buf;
    float *ps_high_buf;
    int first_row_idx;
    /* LANL Even Length Begin */
    /* Flag indicating whether the final lowpass row has been generated
       by the preprocessor */
    int extra_right_low_ready;
    /* Index of the current position in the lookahead buffer for highpass
       rows */
    int next_high_lookahead;
    /* Number of rows in the highpass lookahead buffer */
    int high_lookahead_rows;
    /* LANL Even Length End */
  } hplx_vert_synthesis_by_convolution_obj,
    *hplx_vert_synthesis_by_convolution_ref;
  /* Extends the basic `hplx_synthesis_stage' object for the purpose of
     vertical synthesis by convolution.
         The `cols' field identifies the width of rows processed by this
     synthesis stage.
         The `next_row_idx' field holds the zero-based index of the next
     row to be synthesized.
         The `branch_to_even' array holds two polyphase kernels,
     corresponding to the operators which map the two input branches
     (low- and high-pass subbands) into the even indexed output rows,
     i.e. rows 0,2,4, etc.
         The `branch_to_odd' array holds two polyphase kernels,
     corresponding to the operators which map the two input branches
     (low- and high-pass subbands) into the odd indexed output rows,
     i.e. 1,3,5, etc.
         `even_taps' indicates whether the filter kernels have even or
     odd length.  This affects the type of symmetric extension which
     must be applied.
         The `branch_windows' array manages a rolling window of input
     rows for each of the two input subbands.
         The `conversion_buf' field points to a buffer which is large
     enough to hold an entire reconstructed row of the image.  It is
     only required by calls to `pull_line_fixed'.
         The `taps_copy' and `sign_flips_buffer' arrays must both be
     sufficiently large to cover the support of any of the high-pass
     polyphase kernels (i.e. those in `branch_to_even[1]' and
     `branch_to_odd[1]').  They are used only with even-length filters
     in order to mimic the sign reversal which is required for symmetric
     extension of the high-pass subbands.
         `next_frame_idx' is a counter which identifies the index of the next
     frame to be processed.  It is reset to 0 at the start of each tile.
         The `first_row_idx' field holds the absolute index of the
     first row, with respect to the absolute canvas coordinate system. */


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

/* LANL Even Length Begin */
/*****************************************************************************/
/* STATIC                      synthesis_end_filter                          */
/*****************************************************************************/

/* Implementation of (L',H') -> (L,L) preprocessing */
static void synthesis_end_filter(float* spL, float* spH,
                                 float* dp0, float* dp1, int ncols) {
  int n;
  float d0, d1;

  for (n = 0; n < ncols; n++) {
    d0 = (spL[n] + spH[n]/2.0F);
    d1 = (spL[n] - spH[n]/2.0F);

    dp0[n] = d0;
    dp1[n] = d1;
    spH[n] = 0.0;
 
  }
}
/* LANL Even Length End */

/*****************************************************************************/
/* STATIC                    create_polyphase_kernel                         */
/*****************************************************************************/

static void
  create_polyphase_kernel(hplx_polyphase_kernel_ptr kernel,
                          int input_odd, int output_odd,
                          int neg_support, int pos_support, float *taps)

 /* This function fills out the contents of the polyphase kernel structure
    referenced by `kernel', by reversing and sub-sampling the impulse
    response supplied via `taps', which may be indexed as `taps'[k] where
    k ranges from -`neg_support' to +`pos_support'.  The `input_odd'
    argument identifies whether this impulse response is notionally
    centred on even (low-pass) or odd (high-pass) samples, while the
    `output_odd' argument identifies whether the polyphase kernel is to
    generate the even or odd sub-sequence of the reconstructed signal. */

{
  int displacement, n;

  if (output_odd)
    output_odd = 1; /* Make sure. */
  if (input_odd)
    input_odd = 1; /* Make sure. */
  displacement = output_odd - input_odd;
  kernel->pos_support = (displacement + neg_support)>>1;
  kernel->neg_support = (pos_support - displacement)>>1;
  if ((kernel->neg_support < 0) || (kernel->pos_support < 0))
    return;
  kernel->taps = (float *)
    local_malloc(XFORM_MEM_KEY,
                 sizeof(float)*(kernel->neg_support+kernel->pos_support+1));
  kernel->taps += kernel->neg_support;
  for (n=-kernel->neg_support; n <= kernel->pos_support; n++)
    kernel->taps[n] = taps[displacement - 2*n];
}

/*****************************************************************************/
/* STATIC                    complete_pull_line                              */
/*****************************************************************************/

static void
  complete_pull_line(hplx_synthesis_stage_ref base, float *line_buf)

 /* This function does almost all the work of `__pull_line_float'. */

{
  hplx_vert_synthesis_by_convolution_ref self =
    (hplx_vert_synthesis_by_convolution_ref) base;
  hplx_polyphase_kernel_ptr kernel;
  hplx_rolling_buffer_ptr window=NULL;  /* TJF: =NULL to avoid warnings */
  float *taps, **source[2], *dp, sum;
  int min_row_idx, neg, pos, support;
  int b, n, k;
  int *sign_flips[2];

  /* LANL Even Length Begin */
  /* Set flags to distinguish between cases based on input length
     and offset respectively */
  int oddevencase = 0;
  int oddoddcase = 0;
  int evenevencase = 0;

  if (self->even_taps) {
    if ((self->branch_windows[0].rows+self->branch_windows[1].rows)%2) {
      if (self->first_row_idx%2)
        oddoddcase = 1;
      else
        oddevencase = 1;
    } else {
      if (!(self->first_row_idx%2))
        evenevencase = 1;
    }
  }
  /* LANL Even Length End */

  if ((self->branch_windows[0].rows+self->branch_windows[1].rows) == 1)
    { /* Signals of length (height) < 2 must be treated specially. */

      /* OTLPF_CONVENTION begin; JX Wei ADFA, WJ Zeng Sharp */
      int first_row_odd;
      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 */

      base->branches[first_row_odd]->pull_line_float(
        base->branches[first_row_odd],line_buf,self->cols);
      return;
    }

    /* OTLPF_CONVENTION begin; JX Wei ADFA, WJ Zeng Sharp */
  if(!base->otlpf_convention)
  {
        /* LANL Even Length Begin */
        /* For even length filters, switch kernel order so that filter centers
           are at (1/2, 1/2) for even input offset and (-1/2, -1/2) for odd input
           offset. */
        if (self->even_taps) {
                if (self->first_row_idx%2) {
                if (self->next_row_idx & 1)
                        kernel = self->branch_to_even;
                else
                        kernel = self->branch_to_odd;
                } else {
                if (self->next_row_idx & 1)
                        kernel = self->branch_to_even;
                else
                        kernel = self->branch_to_odd;
                }
        }       else {
        if (self->next_row_idx & 1)
                kernel = self->branch_to_odd;
        else
                kernel = self->branch_to_even;
        }
  /* LANL Even Length End */
  }
  else
  {
        if ((self->next_row_idx-self->first_row_idx) & 1)
    kernel = self->branch_to_odd;
  else
    kernel = self->branch_to_even;
  }