/******************************************************************************Description:   Implements the embedded block coding algorithm, including distortionestimation and R-D covex hull analysis, in addition to the codingpasses themselves.  The low level services offered by the MQ arithmetic coderappear in "mq_encoder.cpp" and "mq_encoder.h".******************************************************************************/#include <math.h>#include <string.h>#include "kdu_messaging.h"#include "kdu_block_coding.h"#include "block_coding_common.h"#include "mq_encoder.h"static kdu_byte *significance_luts[4] =  {lh_sig_lut, hl_sig_lut, lh_sig_lut, hh_sig_lut};#define DISTORTION_LSBS 5#define SIGNIFICANCE_DISTORTIONS (1<<DISTORTION_LSBS)#define REFINEMENT_DISTORTIONS (1<<(DISTORTION_LSBS+1))static kdu_int32 significance_distortion_lut[SIGNIFICANCE_DISTORTIONS];static kdu_int32 significance_distortion_lut_lossless[SIGNIFICANCE_DISTORTIONS];static kdu_int32 refinement_distortion_lut[REFINEMENT_DISTORTIONS];static kdu_int32 refinement_distortion_lut_lossless[REFINEMENT_DISTORTIONS];#define EXTRA_ENCODE_CWORDS 3 // Number of extra context-words between stripes.#define MAX_POSSIBLE_PASSES (31*3-2)/* ========================================================================= *//*                   Local Class and Structure Definitions                   *//* ========================================================================= *//*****************************************************************************//*                             kd_block_encoder                              *//*****************************************************************************/class kd_block_encoder : public kdu_block_encoder_base {  /* Although we can supply a constructor and a virtual destructor in the     future, we have no need for these for the moment. */  protected:    void encode(kdu_block *block, bool reversible, double msb_wmse,                kdu_uint16 estimated_slope_threshold);  private: // Internal implementation    void reset_states()      { // See Table 12.1 in the book by Taubman and Marcellin        for (int n=0; n < 18; n++)          states[n].init(0,0);        states[KAPPA_SIG_BASE].init(4,0);        states[KAPPA_RUN_BASE].init(3,0);      }  private: // Data    mqe_state states[18];  };/* ========================================================================= *//*            Initialization of Distortion Estimation Tables                 *//* ========================================================================= */static void initialize_significance_distortion_luts();static void initialize_refinement_distortion_luts();static class encoder_local_init {    public: encoder_local_init()              { initialize_significance_distortion_luts();                initialize_refinement_distortion_luts(); }  } _do_it;/*****************************************************************************//* STATIC          initialize_significance_distortion_luts                   *//*****************************************************************************/static void  initialize_significance_distortion_luts(){  double fp_scale = (double)(1<<16);  for (kdu_int32 n=0; n < SIGNIFICANCE_DISTORTIONS; n++)    {      kdu_int32 idx = n | (1<<DISTORTION_LSBS);      double v_tilde = ((double) idx) / ((double)(1<<DISTORTION_LSBS));      assert((v_tilde >= 1.0) && (v_tilde < 2.0));      double sqe_before = v_tilde*v_tilde;      double sqe_after = (v_tilde-1.5)*(v_tilde-1.5);      significance_distortion_lut[n] = (int)        floor(0.5 + fp_scale*(sqe_before-sqe_after));      significance_distortion_lut_lossless[n] = (int)        floor(0.5 + fp_scale*sqe_before);    }}/*****************************************************************************//* STATIC            initialize_refinement_distortion_luts                   *//*****************************************************************************/static void  initialize_refinement_distortion_luts(){  double fp_scale = (double)(1<<16);  for (kdu_int32 n=0; n < REFINEMENT_DISTORTIONS; n++)    {      double v_tilde = ((double) n) / ((double)(1<<DISTORTION_LSBS));      assert(v_tilde < 2.0);      double sqe_before = (v_tilde-1.0)*(v_tilde-1.0);      v_tilde = (n >> DISTORTION_LSBS)?(v_tilde-1.0):v_tilde;      assert((v_tilde >= 0.0) && (v_tilde < 1.0));      double sqe_after = (v_tilde-0.5)*(v_tilde-0.5);      refinement_distortion_lut[n] = (int)        floor(0.5 + fp_scale*(sqe_before-sqe_after));      refinement_distortion_lut_lossless[n] = (int)        floor(0.5 + fp_scale*sqe_before);    }}/* ========================================================================= *//*             Binding of MQ and Raw Symbol Coding Services                  *//* ========================================================================= */#define USE_FAST_MACROS // Comment this out if you want functions instead.#ifdef USE_FAST_MACROS#  define _mq_check_out_(coder)                                     \     register kdu_int32 A; register kdu_int32 C; register kdu_int32 t; \     kdu_int32 temp; kdu_byte *store;                                 \     coder.check_out(A,C,t,temp,store)#  define _mq_check_in_(coder)                                      \     coder.check_in(A,C,t,temp,store)#  define _mq_enc_(coder,symbol,state)                              \     _mq_encode_(symbol,state,A,C,t,temp,store)#  define _mq_enc_run_(coder,run)                                   \     _mq_encode_run_(run,A,C,t,temp,store)#  define _raw_check_out_(coder)                                    \     register kdu_int32 t; register kdu_int32 temp; kdu_byte *store;   \     coder.check_out(t,temp,store)#  define _raw_check_in_(coder)                                     \     coder.check_in(t,temp,store)#  define _raw_enc_(coder,symbol)                                   \     _raw_encode_(symbol,t,temp,store)#else // Do not use fast macros#  define _mq_check_out_(coder)#  define _mq_check_in_(coder)#  define _mq_enc_(coder,symbol,state) coder.mq_encode(symbol,state)#  define _mq_enc_run_(coder,run) coder.mq_encode_run(run)#  define _raw_check_out_(coder)#  define _raw_check_in_(coder)#  define _raw_enc_(coder,symbol) coder.raw_encode(symbol)#endif // USE_FAST_MACROS /* The coding pass functions defined below all return a 32-bit integer,    which represents the normalized reduction in MSE associated with the    coded symbols.  Specifically, the MSE whose reduction is returned is    equal to 2^16 * sum_i (x_i/2^p - x_i_hat/2^p)^2 where x_i denotes the    integer sample values in the `samples' array and x_i_hat denotes the    quantized representation available from the current coding pass and all    previous coding passes, assuming a mid-point reconstruction rule.       The mid-point reconstruction rule satisfies x_i_hat = (q_i+1/2)*Delta    where q_i denotes the quantization indices and Delta is the quantization    step size.  This rule is modified only if the `lossless_pass' argument is    true, which is permitted only when symbols coded in the coding pass    result in a lossless representation of the corresponding subband samples.    Of course, this can only happen in the last bit-plane when the reversible    compression path is being used.  In this case, the function uses the fact    that all coded symbols have 0 distortion.       It should be noted that the MSE reduction can be negative, meaning    that the coding of symbols actually increases distortion. *//* ========================================================================= *//*                           Coding pass functions                           *//* ========================================================================= *//*****************************************************************************//* STATIC                   encode_sig_prop_pass_raw                         *//*****************************************************************************/static kdu_int32  encode_sig_prop_pass_raw(mq_encoder &coder, int p, bool causal,                           kdu_int32 *samples, kdu_int32 *contexts,                           int width, int num_stripes, int context_row_gap,                           bool lossless_pass){  /* Ideally, register storage is available for 9 32-bit integers. Two     are declared inside the "_raw_check_out_" macro.  The order of priority     for these registers corresponds roughly to the order in which their     declarations appear below.  Unfortunately, none of these register     requests are likely to be honored by the register-starved X86 family     of processors, but the register declarations may prove useful to     compilers for other architectures or for hand optimizations of     assembly code. */  register kdu_int32 *cp = contexts;  register int c;  register kdu_int32 cword;  _raw_check_out_(coder); // Declares t and temp as registers.  register kdu_int32 sym;  register kdu_int32 val;  register kdu_int32 *sp = samples;  register kdu_int32 shift = 31-p; assert(shift > 0);  int r, width_by2=width+width, width_by3=width_by2+width;  kdu_int32 distortion_change = 0;  kdu_int32 *distortion_lut = significance_distortion_lut;  if (lossless_pass)    distortion_lut = significance_distortion_lut_lossless;  assert((context_row_gap - width) == EXTRA_ENCODE_CWORDS);  for (r=num_stripes; r > 0; r--, cp += EXTRA_ENCODE_CWORDS, sp += width_by3)    for (c=width; c > 0; c--, sp++, cp++)      {        if (*cp == 0)          continue;        cword = *cp;        if ((cword & (NBRHD_MASK<<0)) && !(cword & (SIG_PROP_MEMBER_MASK<<0)))          { // Process first row of stripe column (row 0)            val = sp[0]<<shift; // Move bit p to sign bit.            sym = (kdu_int32)(((kdu_uint32) val)>>31); // Move bit into LSB            _raw_enc_(coder,sym);            if (val >= 0) // New magnitude bit was 0, so still insignificant              { cword |= (PI_BIT<<0); goto row_1; }            // Compute distortion change            val =  (val>>(31-DISTORTION_LSBS)) & (SIGNIFICANCE_DISTORTIONS-1);            distortion_change += distortion_lut[val];            // Encode sign bit            sym = sp[0];            sym = (kdu_int32)(((kdu_uint32) sym)>>31); // Move sign into LSB            _raw_enc_(coder,sym);            // Broadcast neighbourhood context changes            if (!causal)              {                cp[-context_row_gap-1] |=(SIGMA_BR_BIT<<9);                cp[-context_row_gap  ] |=(SIGMA_BC_BIT<<9)|(sym<<NEXT_CHI_POS);                cp[-context_row_gap+1] |=(SIGMA_BL_BIT<<9);              }            cp[-1] |= (SIGMA_CR_BIT<<0);            cp[1]  |= (SIGMA_CL_BIT<<0);            cword |= (SIGMA_CC_BIT<<0) | (PI_BIT<<0) | (sym<<CHI_POS);          }row_1:        if ((cword & (NBRHD_MASK<<3)) && !(cword & (SIG_PROP_MEMBER_MASK<<3)))          { // Process second row of stripe column (row 1)            val = sp[width]<<shift; // Move bit p to sign bit.            sym = (kdu_int32)(((kdu_uint32) val)>>31); // Move bit into LSB            _raw_enc_(coder,sym);            if (val >= 0) // New magnitude bit was 0, so still insignificant              { cword |= (PI_BIT<<3); goto row_2; }            // Compute distortion change            val =  (val>>(31-DISTORTION_LSBS)) & (SIGNIFICANCE_DISTORTIONS-1);            distortion_change += distortion_lut[val];            // Encode sign bit            sym = sp[width];            sym = (kdu_int32)(((kdu_uint32) sym)>>31); // Move sign into LSB            _raw_enc_(coder,sym);            // Broadcast neighbourhood context changes            cp[-1] |= (SIGMA_CR_BIT<<3);            cp[1]  |= (SIGMA_CL_BIT<<3);            cword |= (SIGMA_CC_BIT<<3) | (PI_BIT<<3) | (sym<<(CHI_POS+3));          }row_2:        if ((cword & (NBRHD_MASK<<6)) && !(cword & (SIG_PROP_MEMBER_MASK<<6)))          { // Process third row of stripe column (row 2)            val = sp[width_by2]<<shift; // Move bit p to sign bit.            sym = (kdu_int32)(((kdu_uint32) val)>>31); // Move bit into LSB            _raw_enc_(coder,sym);            if (val >= 0) // New magnitude bit was 0, so still insignificant              { cword |= (PI_BIT<<6); goto row_3; }            // Compute distortion change            val =  (val>>(31-DISTORTION_LSBS)) & (SIGNIFICANCE_DISTORTIONS-1);            distortion_change += distortion_lut[val];            // Encode sign bit            sym = sp[width_by2];