/*****************************************************************************/// File: analysis.cpp [scope = CORESYS/TRANSFORMS]// Version: Kakadu, V2.2// Author: David Taubman// Last Revised: 20 June, 2001/*****************************************************************************/// Copyright 2001, David Taubman, The University of New South Wales (UNSW)// The copyright owner is Unisearch Ltd, Australia (commercial arm of UNSW)// Neither this copyright statement, nor the licensing details below// may be removed from this file or dissociated from its contents./*****************************************************************************/// Licensee: Book Owner// License number: 99999// The Licensee has been granted a NON-COMMERCIAL license to the contents of// this source file, said Licensee being the owner of a copy of the book,// "JPEG2000: Image Compression Fundamentals, Standards and Practice," by// Taubman and Marcellin (Kluwer Academic Publishers, 2001).  A brief summary// of the license appears below.  This summary is not to be relied upon in// preference to the full text of the license agreement, which was accepted// upon breaking the seal of the compact disc accompanying the above-mentioned// book.// 1. The Licensee has the right to Non-Commercial Use of the Kakadu software,//    Version 2.2, including distribution of one or more Applications built//    using the software, provided such distribution is not for financial//    return.// 2. The Licensee has the right to personal use of the Kakadu software,//    Version 2.2.// 3. The Licensee has the right to distribute Reusable Code (including//    source code and dynamically or statically linked libraries) to a Third//    Party, provided the Third Party possesses a license to use the Kakadu//    software, Version 2.2, and provided such distribution is not for//    financial return./******************************************************************************Description:   Implements the forward DWT (subband/wavelet analysis).  The implementationuses lifting to reduce memory and processing, while keeping as much of theimplementation as possible common to both the reversible and the irreversibleprocessing paths.  The implementation is generic to the extent that itsupports any odd-length symmetric wavelet kernels -- although only 3 arecurrently accepted by the "kdu_kernels" object.******************************************************************************/#include <assert.h>#include <string.h>#include <math.h>#include "kdu_messaging.h"#include "kdu_compressed.h"#include "kdu_sample_processing.h"#include "kdu_kernels.h"#include "analysis_local.h"/* ========================================================================= *//*                               kdu_analysis                                *//* ========================================================================= *//*****************************************************************************//*                         kdu_analysis::kdu_analysis                        *//*****************************************************************************/kdu_analysis::kdu_analysis(kdu_resolution resolution,                           kdu_sample_allocator *allocator,                           bool use_shorts, float normalization,                           kdu_roi_node *roi)  // In the future, we may create separate, optimized objects for each kernel.{  state = new kd_analysis(resolution,allocator,use_shorts,normalization,roi);}/* ========================================================================= *//*                              kd_analysis                                  *//* ========================================================================= *//*****************************************************************************//*                        kd_analysis::kd_analysis                           *//*****************************************************************************/kd_analysis::kd_analysis(kdu_resolution resolution,                         kdu_sample_allocator *allocator,                         bool use_shorts, float normalization,                         kdu_roi_node *roi){  reversible = resolution.get_reversible();  this->use_shorts = use_shorts;  int kernel_id = resolution.get_kernel_id();  kdu_kernels kernels(kernel_id,reversible);  float low_gain, high_gain;  float *factors =    kernels.get_lifting_factors(L_max,low_gain,high_gain);  int n;  if ((roi != NULL) && !resolution.propagate_roi())    {      roi->release();      roi = NULL;    }  if (roi != NULL)    roi_level.create(resolution,roi);  assert(L_max <= 4); // We have statically sized the array to improve locality  for (n=0; n < L_max; n++)    {      steps[n].augend_parity = (n+1) & 1; // Step 0 updates odd locations      steps[n].lambda = factors[n];      if (kernels.get_lifting_downshift(n,steps[n].downshift))        { // Reversible case          steps[n].i_lambda = (kdu_int32)            floor(0.5 + steps[n].lambda*(1<<steps[n].downshift));        }      else        { // Irreversible case          steps[n].i_lambda = steps[n].downshift = 0;          kdu_int32 fix_lambda = (kdu_int32) floor(0.5 + factors[n]*(1<<16));          steps[n].fixpoint.fix_lambda = fix_lambda;          steps[n].fixpoint.i_lambda = 0;          while (fix_lambda >= (1<<15))            { steps[n].fixpoint.i_lambda++; fix_lambda -= (1<<16); }          while (fix_lambda < -(1<<15))            { steps[n].fixpoint.i_lambda--; fix_lambda += (1<<16); }          steps[n].fixpoint.remainder = (kdu_int16) fix_lambda;          steps[n].fixpoint.pre_offset = (kdu_int16)            floor(0.5 + ((double)(1<<15)) / ((double) fix_lambda));        }    }  kdu_dims dims;  kdu_coords min, max;  // Get output dimensions.  resolution.get_dims(dims);  min = dims.pos; max = min + dims.size; max.x--; max.y--;  y_next = min.y;  y_max = max.y;  x_min = min.x;  x_max = max.x;    empty = !dims;  unit_height = (y_next==y_max);  unit_width = (x_min==x_max);  low_width = ((max.x+2)>>1) - ((min.x+1)>>1);  high_width = ((max.x+1)>>1) - (min.x>>1);  output_rows_remaining = y_max+1-y_next;  if (empty)    return;  // Pre-allocate the line buffers.  augend.pre_create(allocator,low_width,high_width,reversible,use_shorts);  new_state.pre_create(allocator,low_width,high_width,reversible,use_shorts);  for (n=0; n < L_max; n++)    steps[n].state.pre_create(allocator,low_width,high_width,                              reversible,use_shorts);  initialized = false; // Finalize creation in the first `push' call.  // Now determine the normalizing downshift and subband nominal ranges.  float LL_range, HL_range, LH_range, HH_range;  LL_range = HL_range = LH_range = HH_range = normalization;  normalizing_downshift = 0;  if (!reversible)    {      int lev_idx = resolution.get_dwt_level(); assert(lev_idx > 0);      double bibo_low, bibo_high, bibo_prev;      kernels.get_bibo_gains(lev_idx-1,bibo_prev,bibo_high);      double *bibo_steps = kernels.get_bibo_gains(lev_idx,bibo_low,bibo_high);      double bibo_max = 0.0;      // Find BIBO and nominal ranges for the vertical analysis transform.      if (unit_height)        bibo_max = normalization;      else        {          LL_range /= low_gain;  HL_range /= low_gain;          LH_range /= high_gain; HH_range /= high_gain;          bibo_prev *= normalization; // BIBO horizontal range at stage input          for (n=0; n < L_max; n++)            if ((bibo_prev * bibo_steps[n]) > bibo_max)              bibo_max = bibo_prev * bibo_steps[n];        }      // Find BIBO gains for horizontal analysis      if (!unit_width)        {          LL_range /= low_gain;  LH_range /= low_gain;          HL_range /= high_gain; HH_range /= high_gain;          bibo_prev = bibo_low / low_gain; // If bounded by vertical low band          if ((bibo_high / high_gain) > bibo_prev)            bibo_prev = bibo_high / high_gain; // Bounded by vertical high band          bibo_prev *= normalization; // BIBO vertical range at horiz. input          for (n=0; n < L_max; n++)            if ((bibo_prev * bibo_steps[n]) > bibo_max)              bibo_max = bibo_prev * bibo_steps[n];        }      double overflow_limit = 1.0 * (double)(1<<(16-KDU_FIX_POINT));          // This is the largest numeric range which can be represented in          // our signed 16-bit fixed-point representation without overflow.      while (bibo_max > 0.95*overflow_limit)        { // Leave a little extra headroom to allow for approximations in          // the numerical BIBO gain calculations.          normalizing_downshift++;          LL_range*=0.5F; LH_range*=0.5F; HL_range*=0.5F; HH_range*=0.5F;          bibo_max *= 0.5;        }    }  // Finally, create the subband interfaces.  kdu_roi_node *LL_node=NULL, *HL_node=NULL, *LH_node=NULL, *HH_node=NULL;  if (roi != NULL)    {      LL_node = roi_level.acquire_node(LL_BAND);      HL_node = roi_level.acquire_node(HL_BAND);      LH_node = roi_level.acquire_node(LH_BAND);      HH_node = roi_level.acquire_node(HH_BAND);    }  assert(resolution.which() > 0);  if (resolution.which() == 1)    hor_low[0] = kdu_encoder(resolution.access_next().access_subband(LL_BAND),                             allocator,use_shorts,LL_range,LL_node);  else    hor_low[0] = kdu_analysis(resolution.access_next(),                              allocator,use_shorts,LL_range,LL_node);  hor_high[0] = kdu_encoder(resolution.access_subband(HL_BAND),                            allocator,use_shorts,HL_range,HL_node);  hor_low[1] = kdu_encoder(resolution.access_subband(LH_BAND),                           allocator,use_shorts,LH_range,LH_node);  hor_high[1] = kdu_encoder(resolution.access_subband(HH_BAND),                            allocator,use_shorts,HH_range,HH_node);}/*****************************************************************************//*                       kd_analysis::~kd_analysis                           *//*****************************************************************************/kd_analysis::~kd_analysis(){  hor_low[0].destroy();  hor_low[1].destroy();  hor_high[0].destroy();  hor_high[1].destroy();  if (roi_level.exists())    roi_level.destroy(); // Important to do this last, giving descendants a                         // chance to call the `release' function on their                         // `roi_node' interfaces.}/*****************************************************************************//*                            kd_analysis::push                              *//*****************************************************************************/void  kd_analysis::push(kdu_line_buf &line, bool allow_exchange){  assert(y_next <= y_max);  assert(reversible == line.is_absolute());  if (empty)    {      y_next++;      output_rows_remaining--;      return;    }  int k, c;  if (!initialized)    { // Finish creating all the buffers.      augend.create(); augend.deactivate();      new_state.create(); new_state.deactivate();      for (k=0; k < L_max; k++)        { steps[k].state.create(); steps[k].state.deactivate(); }      initialized = true;    }  // Determine the appropriate input line.  kd_line_cosets *in = (y_next & 1)?(&augend):(&new_state);  if (!in->is_active())    in->activate();  in->lnum = y_next++;  // Copy the samples from `line', de-interleaving even and odd cosets  assert(line.get_width() == (low_width+high_width));  c = x_min & 1; // Index of first coset to be de-interleaved.  k = (line.get_width()+1)>>1; // May move one extra sample.  if (!use_shorts)    { // Working with 32-bit data      kdu_sample32 *sp = line.get_buf32();      kdu_sample32 *dp1 = in->cosets[c].get_buf32();      kdu_sample32 *dp2 = in->cosets[1-c].get_buf32();      if (normalizing_downshift == 0)        for (; k--; sp+=2, dp1++, dp2++)          { *dp1 = sp[0]; *dp2 = sp[1]; }      else        {          float scale = 1.0F / (float)(1<<normalizing_downshift);          for (; k--; sp+=2, dp1++, dp2++)            {