/* *  * 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. *  */#include "libQccPack.h"int QccWAVsfqWaveletAnalysis(const QccIMGImageComponent *input_image,                             QccWAVSubbandPyramid *subband_pyramid,                             QccWAVZerotree *zerotree,                             const QccSQScalarQuantizer *mean_quantizer,                             const QccWAVWavelet *wavelet,                             const QccWAVPerceptualWeights *perceptual_weights,                             double lambda){  int return_value;  int num_levels;  if ((subband_pyramid == NULL) ||      (mean_quantizer == NULL) ||      (wavelet == NULL))    return(0);  num_levels = zerotree->num_levels;  if (QccMatrixCopy(subband_pyramid->matrix,                    input_image->image,                    input_image->num_rows,                    input_image->num_cols))    {      QccErrorAddMessage("(QccWAVsfqWaveletAnalysis): Error calling QccMatrixCopy()");      goto QccError;    }    if (QccWAVSubbandPyramidSubtractMean(subband_pyramid,                                       &(zerotree->image_mean),                                       mean_quantizer))    {      QccErrorAddMessage("(QccWAVsfqWaveletAnalysis): Error calling QccWAVSubbandPyramidSubtractMean()");      goto QccError;    }    if (num_levels)    if (QccWAVSubbandPyramidDWT(subband_pyramid,                                num_levels,                                wavelet))      {        QccErrorAddMessage("(QccWAVsfqWaveletAnalysis): Error calling QccWAVSubbandPyramidDWT()");        goto QccError;      }  if (perceptual_weights != NULL)    if (QccWAVPerceptualWeightsApply(subband_pyramid,                                     perceptual_weights))      {        QccErrorAddMessage("(QccWAVsfqWaveletAnalysis): Error calling QccWAVPerceptualWeightsApply()");        QccErrorExit();      }  return_value = 0;  goto QccReturn; QccError:  return_value = 1; QccReturn:  return(return_value);}static int QccWAVsfqBasebandEncodeProcess(QccIMGImageComponent *baseband_image,                                          QccSQScalarQuantizer                                           *baseband_quantizer,                                          QccChannel *channels,                                          QccIMGImageComponent                                           *reconstructed_baseband,                                          double *distortion,                                           double *rate){  int return_value;  QccVector distortion_vector = NULL;  if ((baseband_image == NULL) ||      (baseband_quantizer == NULL) ||      (channels == NULL))    return(0);  if ((distortion_vector = QccVectorAlloc(baseband_image->num_rows *                                          baseband_image->num_cols)) ==       NULL)    {      QccErrorAddMessage("(QccWAVsfqBasebandEncodeProcess): Error calling QccVectorAlloc()");      goto QccError;    }  if (QccIMGImageComponentScalarQuantize(baseband_image,                                         baseband_quantizer,                                         distortion_vector,                                         &(channels[0])))    {      QccErrorAddMessage("(QccWAVsfqBasebandEncodeProcess): Error calling QccIMGImageComponentScalarQuantize()");      goto QccError;    }  if (QccIMGImageComponentInverseScalarQuantize(&(channels[0]),                                                baseband_quantizer,                                                reconstructed_baseband))    {      QccErrorAddMessage("(QccWAVsfqBasebandEncodeProcess): Error calling QccIMGImageComponentInverseScalarQuantize()");      goto QccError;    }  if (QccChannelNormalize(&(channels[0])))    {      QccErrorAddMessage("(QccWAVsfqBasebandEncodeProcess): Error calling QccChannelNormalize()");      goto QccError;    }  if (distortion != NULL)    *distortion =       QccVectorMean(distortion_vector,                    baseband_image->num_rows *                    baseband_image->num_cols);  if (rate != NULL)    *rate = QccChannelEntropy(&(channels[0]), 1);  return_value = 0;  goto QccReturn; QccError:  return_value = 1; QccReturn:  QccVectorFree(distortion_vector);  return(return_value);}int QccWAVsfqBasebandEncode(const QccWAVSubbandPyramid *subband_pyramid,                            QccWAVZerotree *zerotree,                            QccSQScalarQuantizer *baseband_quantizer,                            QccChannel *channels,                            double lambda,                            QccIMGImageComponent *reconstructed_baseband){  int return_value;  QccIMGImageComponent baseband_image;  double distortion;  double rate;  double baseband_quantizer_start_stepsize;  double baseband_quantizer_end_stepsize;  int done = 0;  int last_time = 0;  QccVector J = NULL;  int pass;  int num_passes;  double min_J = MAXDOUBLE;  double winner = 0;  double max_coefficient = -MAXDOUBLE;  int row, col;  if ((subband_pyramid == NULL) ||      (zerotree == NULL) ||      (baseband_quantizer == NULL) ||      (channels == NULL))    return(0);  if (baseband_quantizer->stepsize <= 0.0)    {      baseband_quantizer_start_stepsize =        QCCWAVSFQ_BASEBANDQUANTIZER_START_STEPSIZE;      baseband_quantizer_end_stepsize =        QCCWAVSFQ_BASEBANDQUANTIZER_END_STEPSIZE;    }  else    {      baseband_quantizer_start_stepsize =        baseband_quantizer_end_stepsize =        baseband_quantizer->stepsize;      last_time = 1;    }    num_passes =    ceil((baseband_quantizer_end_stepsize -          baseband_quantizer_start_stepsize) /         QCCWAVSFQ_BASEBANDQUANTIZER_INCREMENT) + 1;  baseband_image.num_rows = zerotree->num_rows[0];  baseband_image.num_cols = zerotree->num_cols[0];  baseband_image.image = subband_pyramid->matrix;  for (row = 0; row < baseband_image.num_rows; row++)    for (col = 0; col < baseband_image.num_cols; col++)      if (fabs(baseband_image.image[row][col]) > max_coefficient)        max_coefficient = fabs(baseband_image.image[row][col]);  if ((J = QccVectorAlloc(num_passes)) == NULL)    {      QccErrorAddMessage("(QccWAVsfqBasebandEncode): Error calling QccVectorAlloc()");      goto QccError;    }  for (pass = 0; pass < num_passes; pass++)    J[pass] = MAXDOUBLE;  baseband_quantizer->stepsize =     baseband_quantizer_start_stepsize;    pass = 0;  do    {      if (last_time)        done = 1;            if (QccSQScalarQuantization(max_coefficient,                                  baseband_quantizer,                                  NULL,                                  &(channels[0].alphabet_size)))        {          QccErrorAddMessage("(QccWAVsfqBasebandEncode): Error calling QccSQScalarQuantization()");          goto QccError;        }      channels[0].alphabet_size =        channels[0].alphabet_size * 2 + 1;      baseband_quantizer->num_levels = channels[0].alphabet_size;      if (QccWAVsfqBasebandEncodeProcess(&baseband_image,                                         baseband_quantizer,                                         channels,                                         reconstructed_baseband,                                         &distortion, &rate))        {          QccErrorAddMessage("(QccWAVsfqBasebandEncode): Error calling QccWAVsfqBasebandEncodeProcess()");          goto QccError;        }      if (!done)        {          J[pass] =            distortion + lambda * rate;          /****/          /*            printf("    step: %f\n", baseband_quantizer->stepsize);            printf("     MSE: %f\n", distortion);            printf("    rate: %f\n", rate);            printf("       J: %f\n\n", J[pass]);          */          if (pass == num_passes - 1)            {              for (pass = 0;                   pass < num_passes;                   pass++)                if (J[pass] < min_J)                  {                    min_J = J[pass];                    winner = baseband_quantizer_start_stepsize +                      QCCWAVSFQ_BASEBANDQUANTIZER_INCREMENT * pass;                  }                            baseband_quantizer->stepsize = winner;              last_time = 1;            }          else            {              baseband_quantizer->stepsize +=                QCCWAVSFQ_BASEBANDQUANTIZER_INCREMENT;              pass++;            }        }    }  while (!done);  return_value = 0;  goto QccReturn; QccError:  return_value = 1; QccReturn:  QccVectorFree(J);  return(return_value);}static int QccWAVsfqCalcCodewordLengths(QccVector codeword_lengths,                                        int num_symbols,                                        QccChannel *channels,                                        QccWAVZerotree *zerotree){  int subband, row, col;  int symbol;  int channel_index;  int num_nonnull_symbols = 0;    for (symbol = 0; symbol < num_symbols; symbol++)    codeword_lengths[symbol] = 0;    for (subband = 1;       subband < zerotree->num_subbands; subband++)    for (row = 0, channel_index = 0;          row < zerotree->num_rows[subband]; row++)      for (col = 0; col < zerotree->num_cols[subband];            col++, channel_index++)        if (!QccWAVZerotreeNullSymbol(zerotree->zerotree[subband][row][col]))          {            codeword_lengths[channels[subband].channel_symbols                            [channel_index]]++;            num_nonnull_symbols++;          }    for (symbol = 0; symbol < num_symbols; symbol++)    codeword_lengths[symbol] =      -QccMathLog2(codeword_lengths[symbol] / num_nonnull_symbols);    return(0);}static double QccWAVsfqCalcChildrenCost(QccVector *costs,                                        int child_subband,                                        int parent_row,                                        int parent_col,                                        int num_child_cols){  int child_row, child_col;  int num_children;  double cost = 0;    num_children = (child_subband < 4) ? 1 : 2;    for (child_row = num_children*parent_row;       child_row < num_children*(parent_row + 1);       child_row++)    for (child_col = num_children*parent_col;         child_col < num_children*(parent_col + 1);         child_col++)      cost +=        costs[child_subband][child_row*num_child_cols + child_col];    return(cost);  }static double QccWAVsfqSelectMinimumCost(QccWAVZerotree *zerotree,                                          int subband,                                          int row, int col, int cost_index,                                         QccVector *J,                                          QccVector *residue_tree_J,                                         QccVector *delta_J,                                         QccVector *squared_subband_images,                                         int *zerotree_pruned){  int level;  double cost1, cost2;  double minimum_cost;  int child_subband;  int child_subband_start, child_subband_end;    level = QccWAVSubbandPyramidCalcLevelFromSubband(subband, zerotree->num_levels);    if (subband)    child_subband_start = child_subband_end =      subband + 3;