hplx_vert_convolution.c

JPEG2000实现的源码

/*****************************************************************************/
/* Copyright 1998, Hewlett-Packard Company                                   */
/* All rights reserved                                                       */
/* File: "hplx_vert_convolution.c"                                           */
/* Description: Common buffer management tools for vertical synthesis        */
/*              implementations by 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: VM8.5                                                            */
/* Last Revised: 11 September, 2000                                          */
/*****************************************************************************/

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

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

#ifndef MIN
#  define MIN(a,b) ((a<b)?a:b)
#endif
#ifndef MAX
#  define MAX(a,b) ((a>b)?a:b)
#endif

/*****************************************************************************/
/* EXTERN                  hplx_rolling_buffer__advance                      */
/*****************************************************************************/

float *
  hplx_rolling_buffer__advance(hplx_rolling_buffer_ptr window)
{
  int num_buffers, n;

  num_buffers = window->next_row_idx - window->min_row_idx;
  if (num_buffers < window->window_rows)
    num_buffers++;
  else
    {
      float *tmp;

      window->min_row_idx++;
      tmp = window->buffer[0];
      for (n=1; n < num_buffers; n++)
        window->buffer[n-1] = window->buffer[n];
      window->buffer[n-1] = tmp;
    }
  assert(window->next_row_idx < window->rows);
  window->next_row_idx++;
  return(window->buffer[num_buffers-1]);
}

/*****************************************************************************/
/* EXTERN                  hplx_rolling_buffer__get_rows                     */
/*****************************************************************************/

float **
  hplx_rolling_buffer__get_rows(hplx_rolling_buffer_ptr window,
                                int min_idx, int num_rows,
                                int sign_flips[], int wb,
                                int isAnalysis, int usePsExtension,
                                float* lowEndValues)
{
  float **result;
  int avail_min, avail_max;
  int desired_min, desired_max;
  int known_min, known_max;
  int n, direction;

  assert(num_rows <= window->window_rows);
  avail_min = window->min_row_idx;
  avail_max = window->next_row_idx - 1;
  desired_min = min_idx;
  desired_max = min_idx + num_rows - 1;
  if ((desired_max > avail_max) && (avail_max < (window->rows-1)))
    return(NULL); /* Caller must advance window before call can be satisfied */

  /* LANL Even Length Begin */
  /* This step should only be performed when the bottom extension is
     "special" */
  if (window->bottom_extension_special &&
    (desired_max > avail_max) && !window->have_zero_row)
  /* LANL Even Length End */

    { /* Advance window past last row and then fill in the last row with
         all zeros, for use in `special' extensions, i.e. odd boundary
         extensions with even-length filters in the high-pass band. */
      float *zp;

      window->have_zero_row = 1;
      window->rows++; /* Temporarily lie to the `advance' function. */
      zp = hplx_rolling_buffer__advance(window);
      window->rows--;
      window->next_row_idx--;
      for (n=window->cols; n > 0; n--)
        *(zp++) = 0.0F;
      avail_min = window->min_row_idx; /* Has probably changed. */
    }
  if ((desired_min < avail_min) && (avail_min > 0))
    assert(0); /* We will never come back to this row. */
  if (window->have_zero_row && window->bottom_extension_special)
    {
      assert(window->bottom_extension_odd);
      avail_max++; /* The zero row at the end is treated as available and
                      symmetric extension occurs about the zero row.  Note
                      that this only happens when synthesizing an odd
                      height image with even-length filters, and then only
                      in the high-pass subband. */
    }

  /* Point symmetric extension currently doesn't handle this case */
  if ((usePsExtension) && (desired_min<avail_min) && (desired_max>avail_max))
      local_error("Point symmetric extension currently doesn't handle "
                  "such a high number of levels of decomposition with "
                  "this tile size.");

  known_min = MAX(desired_min,avail_min);
  known_max = MIN(desired_max,avail_max);
  result = window->receiving_buf;

  /* LANL Even Length Begin */
  /* The general algorithm doesn't correctly handle extending a single
     row with "non-special" extensions at the top and bottom */
  if (window->rows == 1 &&
      !window->top_extension_special &&
      !window->bottom_extension_special) {
    for (n = 0; n < window->window_rows; n++) {
      result[n] = window->buffer[0];
      if (sign_flips != NULL)
        sign_flips[n] = 0;
    }
    return (result);
  }
  /* LANL Even Length End */

  for (n=known_min; n <= known_max; n++)
    {
      result[n-desired_min] = window->buffer[n-avail_min];
      if (sign_flips != NULL)
        sign_flips[n-desired_min] = 0;
    }
  if (desired_min < known_min)
    {
      direction = -1;
      n = known_min;
      while (desired_min < known_min)
        {
          n += direction;
          if (n < avail_min)
            {
              assert((n == (avail_min-1)) && (direction == -1));
              n = avail_min + window->top_extension_odd;
              direction = 1;
            }
          if (n > avail_max)
            {
              assert((n == (avail_max+1)) && (direction == 1));
              n = avail_max - window->bottom_extension_odd;
              direction = -1;
            }
          if (n < avail_min)
            { /* Looks redundant, but necessary for buffers with height of 1 */
              assert((n == (avail_min-1)) && (direction == -1));
              n = avail_min + window->top_extension_odd;
              direction = 1;
            }
          known_min--;
          result[known_min-desired_min] = window->buffer[n-avail_min];
          if (sign_flips != NULL)
              sign_flips[known_min-desired_min] = (direction == 1);
          /* Begin PSE stuff */
          if (usePsExtension)
            {
              int y = known_min - desired_min;
              int x;

              /* Highpass synthesis uses 0 as the point of symmetry */
              if ((!isAnalysis) && (wb == 1))
                for (x = 0; x < window->cols; x++)
                  window->ps_buf[wb][y][x] = 0 - result[y][x];
              else
                for (x = 0; x < window->cols; x++)
                  window->ps_buf[wb][y][x] = 2*window->buffer[0][x] - result[y][x];

              result[y] = window->ps_buf[wb][y];

              if (sign_flips != NULL)
                sign_flips[y] = 0;
            }
          /* End PSE stuff */
        }
    }
  if (desired_max > known_max)
    {
      direction = 1;
      n = known_max;
      while (desired_max > known_max)
        {
          n += direction;
          if (n < avail_min)
            {
              assert((n == (avail_min-1)) && (direction == -1));
              n = avail_min + window->top_extension_odd;
              direction = 1;
            }
          if (n > avail_max)
            {
              assert((n == (avail_max+1)) && (direction == 1));
              n = avail_max - window->bottom_extension_odd;
              direction = -1;
            }
          if (n < avail_min)
            { /* Looks redundant, but necessary for buffers with height of 1 */
              assert((n == (avail_min-1)) && (direction == -1));
              n = avail_min + window->top_extension_odd;
              direction = 1;
            }
          known_max++;
          result[known_max-desired_min] = window->buffer[n-avail_min];
          if (sign_flips != NULL)
            sign_flips[known_max-desired_min] = (direction == -1);
          /* Begin PSE stuff */
          if (usePsExtension)
            {
              int y = known_max - desired_min;
              int end = avail_max - avail_min;
              int x;

              /* Highpass synthesis uses 0 as the point of symmetry */
              if ((!isAnalysis) && (wb == 1))
                for (x = 0; x < window->cols; x++)
                  window->ps_buf[wb][y][x] = 0 - result[y][x];
              else if (lowEndValues == NULL)
                for (x = 0; x < window->cols; x++)
                  window->ps_buf[wb][y][x] = 2*window->buffer[end][x] - result[y][x];
              else
                for (x = 0; x < window->cols; x++)
                {

                  window->ps_buf[wb][y][x] = 2*lowEndValues[x];
                  window->ps_buf[wb][y][x] -= result[y][x];
                }

              result[y] = window->ps_buf[wb][y];
            }
          /* End PSE stuff */
        }
    }
  return(result);
}