component_mix.c

JPEG2000实现的源码

/*****************************************************************************/
/* Copyright 1998, Hewlett-Packard Company                                   */
/* All rights reserved                                                       */
/* File: "component_mix.c"                                                   */
/* Description: Reference implementation of the `component_mix' 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.   */
/*****************************************************************************/

#include <local_services.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <assert.h>
#include <ifc.h>
#include <markers.h>
#include "component_mix_local.h"
#include "xccxct.h"

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

int std_load_kernel(char *id, char *kernels_dir, std_kernel_info_ptr kernel);
std_user_defined_ptr std_user_defined_summary(std_kernel_info_ptr kernel);

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

/*****************************************************************************/
/* STATIC                       apply_ycbcr_xform                            */
/*****************************************************************************/

static void
apply_ycbcr_xform(ifc_int *buf1, ifc_int *buf2, ifc_int *buf3, int width)
{
  float r, g, b;
  
  for(;width>0;width--){
    r = (float) *buf1; g = (float) *buf2; b = (float) *buf3;
    *(buf1++) = (int) floor(0.5 + 0.29900*r + 0.58700*g + 0.11400*b);
    *(buf2++) = (int) floor(0.5 - 0.16875*r - 0.33126*g + 0.50000*b);
    *(buf3++) = (int) floor(0.5 + 0.50000*r - 0.41869*g - 0.08131*b);
  }
}

/*****************************************************************************/
/* STATIC                        apply_rct_xform                             */
/*****************************************************************************/

static void
apply_rct_xform(ifc_int *buf1, ifc_int *buf2, ifc_int *buf3, int width)
{
  ifc_int r, g, b;
  
  for(;width>0;width--){
    /* Fix for bug reported by Michael Adams */
    r = *buf1; g = *buf2; b = *buf3;
    *(buf1++) = (r+b+g+g)>>2;
    *(buf2++) = b-g;
    *(buf3++) = r-g;
  }
}

/* UofA Begin */
/*****************************************************************************/
/* STATIC                       get_nominal_gain                             */
/*****************************************************************************/

static float
get_nominal_gain(std_tap_info_ptr taps)
     
     /* This function estimates the expansion in nominal range associated with
    a single one-dimensional stage in the Wavelet decomposition, by computing
    the absolute gain of the analysis kernel at DC and at Nyquist and
    returning the maximum of these two. */
     
{
  int n, num_taps;
  double dc, nyq, val;
  
  num_taps = taps->neg_support + taps->pos_support + 1;
  for(dc=nyq=0.0,n=0;n<num_taps;n++){
    val = taps->taps[n];
    dc += val;
    nyq += (n&1)?val:(-val);
  }
  dc = fabs(dc);
  nyq = fabs(nyq);
  return((float)((dc>nyq)?dc:nyq));
}

/*****************************************************************************/
/* STATIC                   compute_waveform_l2_norm                         */
/*****************************************************************************/

static float
compute_waveform_l2_norm(std_tap_info_ptr taps)
{
  double sum, val;
  int n, num_taps;
  
  num_taps = taps->neg_support + taps->pos_support + 1;
  for(sum=0.0,n=0;n<num_taps;n++){
    val = taps->taps[n];
    sum += val*val;
  }
  return((float) sqrt(sum));
}

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

static void
convolve_waveforms(std_tap_info_ptr synth, std_tap_info_ptr in,
           std_tap_info_ptr out, int interp_factor)
     /* This function is key to the determination of synthesis waveforms and their
    L2-norms.  Essentially, the function takes the synthesis waveform provided
    via the `synth' argument, interpolates by `interp_factor', inserting
    `interp_factor'-1 zeros between the entries, and convolves with the
    filter represented by `in', writing the result into `out'. */
{
  int out_neg, out_pos, in_taps, out_taps, synth_taps, i, k;
  float *sp, *op, *ip, *opp, val;
  
  out_neg = synth->neg_support*interp_factor + in->neg_support;
  out_pos = synth->pos_support*interp_factor + in->pos_support;
  out->neg_support = out_neg;
  out->pos_support = out_pos;
  out_taps = out_neg+out_pos+1;
  out->taps = (float *)
    local_realloc(out->taps,sizeof(float)*(size_t) out_taps);
  
  /* Perform interpolated convolution, as described above. */
  
  for (i=0; i < out_taps; i++)
    out->taps[i] = 0.0F; /* Ready to accumulate. */
  synth_taps = synth->neg_support + synth->pos_support + 1;
  in_taps = in->neg_support + in->pos_support + 1;
  op = NULL;
  for(sp=synth->taps,opp=out->taps,i=0;i<synth_taps;i++,opp+=interp_factor){
    val = *(sp++);
    for (k=0,ip=in->taps,op=opp;k<in_taps;k++)
      *(op++) += val * *(ip++);
  }
  assert(op == (out->taps + out_taps));
}

/*****************************************************************************/
/* STATIC                     create_unit_waveform                           */
/*****************************************************************************/

static void
create_unit_waveform(std_tap_info_ptr taps)
{
  taps->neg_support = taps->pos_support = 0;
  taps->taps = (float *)
    local_malloc(MIXING_MEM_KEY,sizeof(float));
  taps->taps[0] = 1.0F;
  taps->int_rdx = -1;
  taps->int_taps = NULL;
}

/*****************************************************************************/
/* STATIC                      install_third_d_kernels                       */
/*****************************************************************************/

static void
install_third_d_kernels(component_tile_kernel_ptr tile_kernel, 
            char *id, char *kernels_dir)
{
  char *id_ptr, *dp;
  int id_len;
  
  id_ptr = id;
  id_len = strlen(id_ptr);
  
  tile_kernel->third_d_identifier = dp = (char *)
    local_malloc(MIXING_MEM_KEY,id_len+1);
  while (id_len--)
    *(dp++) = *(id_ptr++);
  *dp = '\0';
  
  /* Now we are ready to load the kernel specification. */
  
  tile_kernel->third_d_kernel_user_defined =
    std_load_kernel(tile_kernel->third_d_identifier,kernels_dir,
            &(tile_kernel->third_d_kernel));
}

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

static int
get_third_d_kernel_type(the_component_mix_ref self,int which,
            char *default_filter,int *low_neg_support,
            int *low_pos_support,int *high_neg_support,
            int *high_pos_support)
{
  component_tile_kernel_ptr  tile_kernel;
  int tap_idx = 0, kernel_type = 0;
  std_tap_info_ptr low = 0, high = 0;
  
  if (self->current_tile_idx >= self->num_tiles)
    local_error("Attempting to access 3-D DWT information beyond "
        "the last image tile!");
  tile_kernel = self->tile_kernels + self->current_tile_idx;
  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;
  
  if (strcmp(tile_kernel->third_d_identifier, default_filter)) {
    char *data;
    int data_size;