/*
 * 
 * QccPack: Quantization, compression, and coding libraries
 * Copyright (C) 1997-2005  James E. Fowler
 * 
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 * 
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Library General Public License for more details.
 * 
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 675 Mass Ave, Cambridge,
 * MA 02139, USA.
 * 
 */

/* 
 *
 * Written by
 *
 * Chao Tian, at Cornell University, 2003 
 *
 */


#include "libQccPack.h"


#define QCCTCE_BOUNDARY_VALUE (1e-6)

//definition used in context information
#define QCCTCE_Z 0 //zero, AKA, insignificant
#define QCCTCE_NZN_NEW 2 //non-zero-neighbor
#define QCCTCE_S 3 //significant
#define QCCTCE_S_NEW 4
#define QCCTCE_NZN 5

// more refinement can be made if direction of Tarp filter is considered
// but the improvement is minor
#define QCCWAVTCE_ALPHA 0.4 
#define QCCWAVTCE_ALPHA_HIGH 0.4
#define QCCWAVTCE_ALPHA_LOW 0.4

//ugly fix, use 1-D IIR filter to provide PMF estimate
#define QCCWAVTCE_ALPHA_1D 0.995 
#define QCCWAVTCE_ALPHA_1D_O 0.005

// this should be 0.5, however, 0.3 seems to work a little better (minor)
// reason unknown
#define QCCWAVTCE_REFINE_HOLDER 0.3 

// threshold for cross-scale prediction...
#define QCCWAVTCE_PREDICT_THRESHOLD 0.05  
//weight factor for cross-scale and current scale
#define QCCWAVTCE_CURRENT_SCALE 0.7 
#define QCCWAVTCE_PARENT_SCALE 0.3


static int QccWAVtceUpdateModel(QccENTArithmeticModel *model, double prob)
{
  double probabilities[2]; 
  probabilities[1] = prob;
  probabilities[0] = 1 - probabilities[1];
  
  if (QccENTArithmeticSetModelProbabilities(model,probabilities,0))
    {
      QccErrorAddMessage("(QccWAVtceUpdateModel): Error calling QccENTArithmeticSetModelProbabilities()");
      return(1);
    } 
  return(0);
}


static int QccWAVtceEncodeDWT(QccWAVSubbandPyramid *image_subband_pyramid,
                              char **sign_array,
                              const QccIMGImageComponent *image,
                              int num_levels,
                              double *image_mean,
                              int *max_coefficient_bits,
                              double stepsize,
                              const QccWAVWavelet *wavelet)
{
  double coefficient_magnitude;
  double max_coefficient;
  int num_subbands,subband;
  int row, col;
  int subband_origin_row;
  int subband_origin_col;
  int subband_num_rows;
  int subband_num_cols;
  
  num_subbands =
    QccWAVSubbandPyramidNumLevelsToNumSubbands(num_levels);
  
  if (QccMatrixCopy(image_subband_pyramid->matrix,
                    image->image,
                    image->num_rows,
                    image->num_cols))
    {
      QccErrorAddMessage("(QccWAVtceEncodeDWT): Error calling QccMatrixCopy()");
      return(1);
    }
  
  if (QccWAVSubbandPyramidSubtractMean(image_subband_pyramid,
                                       image_mean,
                                       NULL))
    {
      QccErrorAddMessage("(QccWAVtceEncodeDWT): Error calling QccWAVSubbandPyramidSubtractMean()");
      return(1);
    }
  
  if (QccWAVSubbandPyramidDWT(image_subband_pyramid,
                              num_levels,
                              wavelet))
    {
      QccErrorAddMessage("(QccWAVtceEncodeDWT): Error calling QccWAVSubbandPyramidDWT()");
      return(1);
    }
  
  for (subband = 0; subband < num_subbands; subband++)
    {
      if (QccWAVSubbandPyramidSubbandSize(image_subband_pyramid,
                                          subband,
                                          &subband_num_rows,
                                          &subband_num_cols))
        {
          QccErrorAddMessage("(QccWAVtceEncodeDWT): Error calling QccWAVSubbandPyramidSubbandSize()");
          return(1);
        }

      if (QccWAVSubbandPyramidSubbandOffsets(image_subband_pyramid,
                                             subband,
                                             &subband_origin_row,
                                             &subband_origin_col))
        {
          QccErrorAddMessage("(QccWAVtceEncodeDWT): Error calling QccWAVSubbandPyramidSubbandOffsets()");
          return(1);
        }
      max_coefficient = -MAXFLOAT;
      
      if (stepsize > 0)
        for (row = 0; row < subband_num_rows; row++)
          for (col = 0; col < subband_num_cols; col++)
            {
              coefficient_magnitude =
                floor(fabs(image_subband_pyramid->matrix
                           [subband_origin_row + row]
                           [subband_origin_col + col] / stepsize));
              
              if (image_subband_pyramid->matrix
                  [subband_origin_row + row]
                  [subband_origin_col + col] < 0)
                sign_array
                  [subband_origin_row + row]
                  [subband_origin_col + col] = 1;
              else
                sign_array
                  [subband_origin_row + row]
                  [subband_origin_col + col] = 0;
              
              image_subband_pyramid->matrix
                [subband_origin_row + row]
                [subband_origin_col + col] = coefficient_magnitude;     
              
              if (coefficient_magnitude > max_coefficient)
                max_coefficient = coefficient_magnitude;     
            }
      else
        for (row = 0; row < subband_num_rows; row++)
          for (col = 0; col < subband_num_cols; col++)
            {
              coefficient_magnitude =
                fabs(image_subband_pyramid->matrix
                     [subband_origin_row + row]
                     [subband_origin_col + col]);
              
              if (image_subband_pyramid->matrix
                  [subband_origin_row + row]
                  [subband_origin_col + col] < 0)
                sign_array
                  [subband_origin_row + row]
                  [subband_origin_col + col] = 1;
              else
                sign_array
                  [subband_origin_row + row]
                  [subband_origin_col + col] = 0;
              
              image_subband_pyramid->matrix
                [subband_origin_row + row]
                [subband_origin_col + col] = coefficient_magnitude;     

              if (coefficient_magnitude > max_coefficient)
                max_coefficient = coefficient_magnitude;
            }

      max_coefficient_bits[subband] = 
        (int)floor(QccMathMax(0, QccMathLog2(max_coefficient + 0.0001)));  

    }   

  return(0);
}


static int QccWAVtceEncodeBitPlaneInfo(QccBitBuffer *output_buffer,
                                       int num_subbands,
                                       int *max_coefficient_bits)
{
  int subband;
  int num_zeros = 0;
  
  for (subband = 1; subband < num_subbands; subband++)
    {
      for (num_zeros = max_coefficient_bits[0] - max_coefficient_bits[subband];
           num_zeros > 0;
           num_zeros--)
        if (QccBitBufferPutBit(output_buffer, 0))
          return(1);

      if (QccBitBufferPutBit(output_buffer, 1))
        return(1);
    }

  return(0);
}


static int QccWAVtceDecodeBitPlaneInfo(QccBitBuffer *input_buffer,
                                       int num_subbands,
                                       int *max_bits,
                                       int max_coefficient_bits)
{
  int subband;
  int bit_value;
  
  max_bits[0] = max_coefficient_bits;
  
  for (subband = 1; subband < num_subbands; subband++)
    {
      max_bits[subband] = max_coefficient_bits;
      do
        {
          if (QccBitBufferGetBit(input_buffer,&bit_value))
            return(1);
          max_bits[subband]--;
        }
      while (bit_value == 0);
      max_bits[subband]++;  
    }
  
  return(0);
}


int QccWAVtceEncodeHeader(QccBitBuffer *output_buffer, 
                          int num_levels, 
                          int num_rows, int num_cols,
                          double image_mean,
                          double stepsize,
                          int max_coefficient_bits)
{
  int return_value;
  
  if (QccBitBufferPutChar(output_buffer, (unsigned char)num_levels))
    {
      QccErrorAddMessage("(QccWAVtceEncodeHeader): Error calling QccBitBufferPuChar()");
      goto Error;
    }
  
  if (QccBitBufferPutInt(output_buffer, num_rows))
    {
      QccErrorAddMessage("(QccWAVtceEncodeHeader): Error calling QccBitBufferPutInt()");
      goto Error;
    }
  
  if (QccBitBufferPutInt(output_buffer, num_cols))
    {
      QccErrorAddMessage("(QccWAVtceEncodeHeader): Error calling QccBitBufferPutInt()");
      goto Error;
    }
  
  if (QccBitBufferPutDouble(output_buffer, image_mean))
    {
      QccErrorAddMessage("(QccWAVtceEncodeHeader): Error calling QccBitBufferPutDouble()");
      goto Error;
    }
  
  if (QccBitBufferPutDouble(output_buffer, stepsize))
    {
      QccErrorAddMessage("(QccWAVtceEncodeHeader): Error calling QccBitBufferPutDouble()");
      goto Error;
    }
  
  if (QccBitBufferPutInt(output_buffer, max_coefficient_bits))
    {
      QccErrorAddMessage("(QccWAVtceEncodeHeader): Error calling QccBitBufferPutInt()");
      goto Error;
    }
  
  return_value = 0;
  goto Return;
 Error:
  return_value = 1;
 Return:
  return(return_value);
}


static int QccWAVtceRevEst(QccWAVSubbandPyramid *coefficients,
                           char **significance_map,
                           int subband,
                           double *subband_significance,
                           double **p)
{
  
  int return_value;
  int subband_origin_row;
  int subband_origin_col;
  int subband_num_rows;
  int subband_num_cols;
  int row, col;
  int current_row, current_col;
  double p1, p3;
  QccVector p2 = NULL;
  char *p_char;
  int v;  
  double filter_coef;
  double alpha_v, alpha_h;
  
  if (subband % 3==1)
    { 
      //LH horizental more important?
      alpha_h = QCCWAVTCE_ALPHA_HIGH;
      alpha_v = QCCWAVTCE_ALPHA_LOW;
    }
  else
    {
      if (subband % 3==2)
        {
          alpha_h = QCCWAVTCE_ALPHA_LOW;
          alpha_v = QCCWAVTCE_ALPHA_HIGH;   
        }
      else
        {
          alpha_h = QCCWAVTCE_ALPHA;
          alpha_v = QCCWAVTCE_ALPHA;
        }
    }
  
  filter_coef = (1 - alpha_h) * (1 - alpha_v) / (2*alpha_h + 2*alpha_v);
  
  subband_significance[subband] = 0.0;
  
  if (QccWAVSubbandPyramidSubbandSize(coefficients,
                                      subband,
                                      &subband_num_rows,
                                      &subband_num_cols))
    {
      QccErrorAddMessage("(QccWAVtceRevEst): Error calling QccWAVSubbandPyramidSubbandSize()");
      goto Error;
    }
  if (QccWAVSubbandPyramidSubbandOffsets(coefficients,
                                         subband,
                                         &subband_origin_row,
                                         &subband_origin_col))
    {
      QccErrorAddMessage("(QccWAVtceRevEst): Error calling QccWAVSubbandPyramidSubbandOffsets()");
      goto Error;
    }