/*****************************************************************************/// File: synthesis.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 reverse DWT (subband/wavelet synthesis).  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 "synthesis_local.h"// Set things up for the inclusion of assembler optimized routines// for specific architectures.  The reason for this is to exploit// the availability of SIMD type instructions on many modern processors.#if defined KDU_PENTIUM_MSVC#  define KDU_SIMD_OPTIMIZATIONS#  include "msvc_dwt_mmx_local.h" // Header contains all asm commands in-linestatic int simd_exists = msvc_dwt_mmx_exists();#elif defined KDU_PENTIUM_GCC#  define KDU_SIMD_OPTIMIZATIONS#  include "gcc_dwt_mmx_local.h" // Header declares functs in separate .s filestatic int simd_exists = gcc_dwt_mmx_exists();#endif // KDU_PENTIUM_GCC/* ========================================================================= *//*                               kdu_synthesis                               *//* ========================================================================= *//*****************************************************************************//*                        kdu_synthesis::kdu_synthesis                       *//*****************************************************************************/kdu_synthesis::kdu_synthesis(kdu_resolution resolution,                             kdu_sample_allocator *allocator,                             bool use_shorts, float normalization)  // In the future, we may create separate, optimized objects for each kernel.{  state = new kd_synthesis(resolution,allocator,use_shorts,normalization);}/* ========================================================================= *//*                              kd_synthesis                                 *//* ========================================================================= *//*****************************************************************************//*                       kd_synthesis::kd_synthesis                          *//*****************************************************************************/kd_synthesis::kd_synthesis(kdu_resolution resolution,                           kdu_sample_allocator *allocator,                           bool use_shorts, float normalization){  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;  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_out_next = min.y;  y_out_max = max.y;  x_out_min = min.x;  x_out_max = max.x;    empty = !dims;  if (empty)    return;  // Get input dimensions.  Note that the following code contains assert  // statements which may fail if the low- and high-pass filter lengths do not  // differ by 2 (as reported by `kdu_kernels' and used in `kdu_resolution' to  // compute the effects of kernel support on regions of interest.  resolution.access_next().get_dims(dims); // low-pass dimensions.  low_width = dims.size.x;  min = dims.pos; max = min + dims.size; max.x--; max.y--;  min.x+=min.x; min.y+=min.y; max.x+=max.x; max.y+=max.y;  y_in_next = min.y;  y_in_max = max.y;  x_in_min = min.x;  x_in_max = max.x;  resolution.access_subband(HH_BAND).get_dims(dims); // high-pass dimensions.  high_width = dims.size.x;  min = dims.pos; max = min + dims.size; max.x--; max.y--;  min.x+=min.x+1; min.y+=min.y+1; max.x+=max.x+1; max.y+=max.y+1;  assert(min.y <= y_in_next+1);  if (min.y < y_in_next)    { y_in_next--; assert(min.y==y_in_next); }  assert(max.y >= y_in_max-1);  if (max.y > y_in_max)    { y_in_max++; assert(max.y==y_in_max); }  assert(min.x <= x_in_min+1);  if (min.x < x_in_min)    { x_in_min--; assert(min.x==x_in_min); }  assert(max.x >= x_in_max-1);  if (max.x > x_in_max)    { x_in_max++; assert(max.x==x_in_max); }  assert((y_in_next <= y_out_next) && (y_in_max >= y_out_max));  assert((x_in_min <= x_out_min) && (x_in_max >= x_out_max));  assert((y_in_next <= y_in_max) && (x_in_min <= x_in_max));    unit_height = (y_in_next==y_in_max);  unit_width = (x_in_min==x_in_max);  // 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 upshift and subband nominal ranges.  float LL_range, HL_range, LH_range, HH_range;  LL_range = HL_range = LH_range = HH_range = normalization;  normalizing_upshift = 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.75*overflow_limit)        { // Leave some extra headroom to allow for the effects of          // dequantization artefacts.          normalizing_upshift++;          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.  assert(resolution.which() > 0);  if (resolution.which() == 1)    hor_low[0] = kdu_decoder(resolution.access_next().access_subband(LL_BAND),                             allocator,use_shorts,LL_range);  else    hor_low[0] = kdu_synthesis(resolution.access_next(),                               allocator,use_shorts,LL_range);  hor_high[0]  = kdu_decoder(resolution.access_subband(HL_BAND),                             allocator,use_shorts,HL_range);  hor_low[1]   = kdu_decoder(resolution.access_subband(LH_BAND),                             allocator,use_shorts,LH_range);  hor_high[1]  = kdu_decoder(resolution.access_subband(HH_BAND),                             allocator,use_shorts,HH_range);}/*****************************************************************************//*                      kd_synthesis::~kd_synthesis                          *//*****************************************************************************/kd_synthesis::~kd_synthesis(){  hor_low[0].destroy();  hor_low[1].destroy();  hor_high[0].destroy();  hor_high[1].destroy();}/*****************************************************************************//*                           kd_synthesis::pull                              *//*****************************************************************************/void  kd_synthesis::pull(kdu_line_buf &line, bool allow_exchange){  if (empty)    return;  if (!initialized)    { // Finish creating all the buffers.      augend.create(); augend.deactivate();      new_state.create(); new_state.deactivate();      for (int n=0; n < L_max; n++)        { steps[n].state.create(); steps[n].state.deactivate(); }      initialized = true;    }  int c, k;  kd_line_cosets *out = (y_out_next & 1)?(&new_state):(&augend);  assert(y_out_next <= y_out_max);  if (unit_height)    { // No transform performed in this special case.      horizontal_synthesis(*out);      if (reversible && (y_out_next & 1))        { // Need to halve integer sample values.          if (!use_shorts)            { // Working with 32-bit data              kdu_sample32 *dp;              for (c=0; c < 2; c++)                for (dp=out->cosets[c].get_buf32(),                     k=out->cosets[c].get_width(); k--; dp++)                  dp->ival >>= 1;            }          else            { // Working with 16-bit data              kdu_sample16 *dp;              for (c=0; c < 2; c++)                for (dp=out->cosets[c].get_buf16(),                     k=out->cosets[c].get_width(); k--; dp++)                  dp->ival >>= 1;            }        }    }