component_demix.c

JPEG2000实现的源码

/*****************************************************************************/
/* Copyright 1998, Hewlett-Packard Company                                   */
/* All rights reserved                                                       */
/* File: "component_demix.c"                                                 */
/* Description: Reference implementation of the `component_demix' object.    */
/* Author: David Taubman                                                     */
/* Affiliation: Hewlett-Packard and                                          */
/*              The University of New South Wales, Australia                 */
/* Version: VM9.0                                                            */
/* Last Revised: 18 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 by David Taubman to implement changes in Tokyo between CD and    */
/* FCD for Part-1 of the standard.  Copyrighted by HP with all rights        */
/* reserved for the modified parts.                                          */
/*****************************************************************************/

/*****************************************************************************/
/* Modified by Ali Bilgin to support wavelet transforms across the component */
/* direction. Copyright 2000 University of Arizona, Arizona Board of Regents.*/
/* All rights reserved for the modified parts.                               */
/*****************************************************************************/

/*****************************************************************************/
/* Modified to incorporate linear component transform.  Material identified  */
/* by "Kodak/SAIC Linear Transform" comments has been added by Austin Lan of */
/* Kodak and/or Tom Flohr of SAIC.                                           */
/*****************************************************************************/

/*****************************************************************************/
/* Modified by Brendt Wohlberg (Los Alamos National Laboratory) to include   */
/* component-direction wavelet transforms based on those defined in Part 2   */
/* Annex G. All modifications are Copyright 2001 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 <markers.h>
#include "component_demix_local.h"
#include <stdio.h>

/* LANL begin */
#ifdef LANL_CMPNT_WVLT
#define DEMIX
#include "VMInterface.h"
#define DEMIX
#endif
/* LANL end */

void std_translate_kernel(char *id, std_user_defined_ptr user,
                          std_kernel_info_ptr kernel);

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

/*****************************************************************************/
/* STATIC                      invert_ycbcr_xform                            */
/*****************************************************************************/

static void
  invert_ycbcr_xform(ifc_int *buf1, ifc_int *buf2, ifc_int *buf3, int width)
{
  float y, cb, cr;

  for (; width > 0; width--)
    {
      y = (float) *buf1; cb = (float) *buf2; cr = (float) *buf3;
      *(buf1++) = (int) floor(0.5 + y - 0.00001*cb + 1.40200*cr);
      *(buf2++) = (int) floor(0.5 + y - 0.34413*cb - 0.71414*cr);
      *(buf3++) = (int) floor(0.5 + y + 1.77200*cb + 0.00004*cr);
    }
}

/*****************************************************************************/
/* STATIC                        invert_rct_xform                            */
/*****************************************************************************/

static void
  invert_rct_xform(ifc_int *buf1, ifc_int *buf2, ifc_int *buf3, int width)
{
  ifc_int y, cb, cr, g;

  for (; width > 0; width--)
    {
      /* Fix for bug reported by Michael Adams */
      y = *buf1; cb = *buf2; cr = *buf3;
      g = y - ((cb+cr)>>2);
      *(buf1++) = g + cr;
      *(buf2++) = g;
      *(buf3++) = g + cb;
    }
}

/* UofA Begin */
/*****************************************************************************/
/* STATIC                        destroy_taps                                */
/*****************************************************************************/

static void
  destroy_taps(std_tap_info_ptr taps)
{
  if (taps->taps != NULL)
    local_free(taps->taps);
  if (taps->int_taps != NULL)
    local_free(taps->int_taps);
}

/*****************************************************************************/
/* STATIC                       destroy_lifting                              */
/*****************************************************************************/

static void
  destroy_lifting(std_lifting_info_ptr lifting)
{
  int n;

  if (lifting->steps != NULL)
    {
      for (n=0; n < lifting->num_steps; n++)
        destroy_taps(lifting->steps + n);
      local_free(lifting->steps);
    }
}

/*****************************************************************************/
/* STATIC                        destroy_kernel                              */
/*****************************************************************************/

static void
  destroy_kernel(std_kernel_info_ptr kernel)
{
  int i;

  destroy_lifting(&(kernel->lifting));
  for (i=0; i < 2; i++)
    {
      destroy_taps(kernel->low+i);
      destroy_taps(kernel->high+i);
    }
}

/*****************************************************************************/
/* STATIC                   get_third_d_kernel_type                          */
/*****************************************************************************/

static int
  get_third_d_kernel_type(component_tile_kernel_ptr tile_kernel, int which,
                          int *low_neg_support, int *low_pos_support,
                          int *high_neg_support, int *high_pos_support)
{
  int tap_idx = 0, kernel_type = 0;
  std_tap_info_ptr low = 0, high = 0;

  if (which == INFO__ANALYSIS)
    tap_idx = 0;
  else if (which == INFO__SYNTHESIS)
    tap_idx = 1;
  else
    local_error("Illegal kernel selector supplied to function, "
                "`get_kernel_type'!!  Must be one of INFO__ANALYSIS or "
                "INFO__SYNTHESIS.");

  low = tile_kernel->third_d_kernel.low;
  high = tile_kernel->third_d_kernel.high;
  kernel_type = tile_kernel->third_d_kernel.kernel_type;

  if (low_neg_support != NULL)
    *low_neg_support = low[tap_idx].neg_support;
  if (low_pos_support != NULL)
    *low_pos_support = low[tap_idx].pos_support;
  if (high_neg_support != NULL)
    *high_neg_support = high[tap_idx].neg_support;
  if (high_pos_support != NULL)
    *high_pos_support = high[tap_idx].pos_support;
  return(kernel_type);
}

/*****************************************************************************/
/* STATIC                get_third_d_convolution_taps                        */
/*****************************************************************************/

static void
  get_third_d_convolution_taps(component_tile_kernel_ptr tile_kernel,
                               float **low_taps, float **high_taps)
{
  std_tap_info_ptr low = 0, high = 0;
  int kernel_type = 0;

  low = tile_kernel->third_d_kernel.low;
  high = tile_kernel->third_d_kernel.high;
  kernel_type = tile_kernel->third_d_kernel.kernel_type;

  if (kernel_type != INFO__CONVOLUTION)
    local_error("Attempting to retrieve convolution taps for a "
                "Wavelet kernel which is to be implemented by lifting "
                "rather than convolution!!");
  if (low_taps != NULL)
    *low_taps = low[1].taps + low[1].neg_support;
  if (high_taps != NULL)
    *high_taps = high[1].taps + high[1].neg_support;
}

/*****************************************************************************/
/* STATIC                    get_third_d_int_steps                           */
/*****************************************************************************/

static int
  get_third_d_int_steps(component_tile_kernel_ptr tile_kernel, 
                        int neg_support[], int pos_support[], int *taps[], 
                        int downshift[])
{
  std_lifting_info_ptr info = 0;
  int n, kernel_type = 0;

  info = &(tile_kernel->third_d_kernel.lifting);
  kernel_type = tile_kernel->third_d_kernel.kernel_type;

  if (kernel_type != INFO__INT_LIFTING)
    local_error("Attempting to retrieve reversible lifting steps "
                "for a Wavelet kernel which is to be implemented using "
                "either non-reversible lifting steps or convolution!");
  assert(info->num_steps <= INFO__MAX_LIFTING_STEPS);
  for (n=0; n < info->num_steps; n++)
    {
      neg_support[n] = info->steps[n].neg_support;
      pos_support[n] = info->steps[n].pos_support;
      taps[n] = info->steps[n].int_taps + info->steps[n].neg_support;
      downshift[n] = info->steps[n].int_rdx;
    }
  return(info->num_steps);
}

/*****************************************************************************/
/* 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)
{
  int even_taps, n;
  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;
  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;
    }
  sp = src; dp = sp + src_cols - 1;
  /* extension to the left*/
  for (n=1; n <= neg_support; n++)
    {
      /* Half point extension*/
      if(left_ext)
        {
          assert( n  <= src_cols );
          sp[-n] = sp[n-1]*left_ext;
        }
      /* whole point extension*/
      else
        {
          assert( n <=src_cols);
          /* TJF: above was formerly assert( n < src_cols -1) */
          sp[-n] = sp[n];
        }
    }
  /* Extension to the right*/
  /* Extension to the right*/
  for (n=1; n <= pos_support+dst_cols-src_cols; n++)