/* *  * 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. *  *//* *  This code was written by Joe Boettcher <jbb15@msstate.edu> based on *  the originally developed algorithm and code by Suxia Cui. */#include "libQccPack.h"#ifdef LIBQCCPACKSPIHT_H#define QCCVIDRDWTBLOCK_WINDOWSIZE 15typedef struct{  int blocksize;  QccMatrix block;} QccVIDRDWTBlockBlock;typedef struct{  int intraframe;  int motion_vector_bits;  int intraframe_bits;} QccVIDRDWTBlockStatistics;static void QccVIDRDWTBlockBlockInitialize(QccVIDRDWTBlockBlock *block){  block->blocksize = 0;  block->block = NULL;}static int QccVIDRDWTBlockBlockAlloc(QccVIDRDWTBlockBlock *block){  if ((block->block =       QccMatrixAlloc(block->blocksize, block->blocksize)) == NULL)    {      QccErrorAddMessage("(QccVIDRDWTBlockBlockAlloc): Error calling QccMatrixAlloc()");      return(1);    }  return(0);}static void QccVIDRDWTBlockBlockFree(QccVIDRDWTBlockBlock *block){  QccMatrixFree(block->block, block->blocksize);}static void QccVIDRDWTBlockPrintStatistics(const QccVIDRDWTBlockStatistics                                           *statistics){  int total_bits;  int field_width;  total_bits =    statistics->motion_vector_bits + statistics->intraframe_bits;  field_width = (int)ceil(log10((double)total_bits));  printf("            Frame type: ");  if (statistics->intraframe)    printf("I\n");  else    printf("P\n");  printf("    Motion-vector bits: %*d (%5.1f%%)\n",         field_width,         statistics->motion_vector_bits,         QccMathPercent(statistics->motion_vector_bits, total_bits));  printf("       Intraframe bits: %*d (%5.1f%%)\n",         field_width,         statistics->intraframe_bits,         QccMathPercent(statistics->intraframe_bits, total_bits));  printf("\n");}static double QccVIDRDWTBlockMAE(QccVIDRDWTBlockBlock *block1,                                 QccVIDRDWTBlockBlock *block2,                                 int num_levels){  int row, col;  double mae = 0;  int level;  int index_skip = 1;  double scale_factor;  for (level = 1; level <= num_levels; level++)    {      index_skip = (1 << level);      scale_factor = 1.0 / index_skip;      for (row = 0; row < (block1->blocksize / index_skip);           row++)        for (col = 0; col < (block1->blocksize / index_skip);             col++)          {            //HL            mae += scale_factor *              fabs(block1->block                   [row][block1->blocksize / index_skip + col] -                   block2->block                   [row][block1->blocksize / index_skip + col]);            //LH            mae += scale_factor *              fabs(block1->block                   [block1->blocksize / index_skip + row][col] -                   block2->block                   [block1->blocksize / index_skip + row][col]);             //HH            mae += scale_factor *              fabs(block1->block                   [block1->blocksize / index_skip + row]                   [block1->blocksize / index_skip + col] -                   block2->block                   [block1->blocksize / index_skip + row]                   [block1->blocksize / index_skip + col]);          }    }    //LL    scale_factor = 1.0 / index_skip;  for (row = 0; row < (block1->blocksize / index_skip); row++)    for (col = 0; col < (block1->blocksize / index_skip); col++)      mae += scale_factor *        fabs(block1->block[row][col] - block2->block[row][col]);    return(mae);}static double QccVIDRDWTBlockExtractCoefficient(QccWAVSubbandPyramid                                                *subband_pyramid,                                                int row,                                                int col){  if (row < 0)    row = 0;  else    if (row >= subband_pyramid->num_rows)      row = subband_pyramid->num_rows - 1;  if (col < 0)    col = 0;  else    if (col >= subband_pyramid->num_cols)      col = subband_pyramid->num_cols - 1;  return(subband_pyramid->matrix[row][col]);}static int QccVIDRDWTBlockExtractBlock(QccWAVSubbandPyramid *subband_pyramid,                                       QccVIDRDWTBlockBlock *wavelet_block,                                       int ref_offset_row,                                        int ref_offset_col,                                       int num_rows,                                        int num_cols,                                       int num_levels){  int row, col;  int offset_row, offset_col;  int level;  int index_skip;  int pattern;  int row_pattern, col_pattern;  int v, h;    index_skip = (1 << num_levels);  offset_row = ref_offset_row / index_skip;  offset_col = ref_offset_col / index_skip;  row_pattern = QccMathModulus(ref_offset_row, index_skip);  col_pattern = QccMathModulus(ref_offset_col, index_skip);    for (level = 1; level <= num_levels; level++)    {      index_skip = (1 << level);      offset_row = ref_offset_row / index_skip;      offset_col = ref_offset_col / index_skip;      pattern = (1 << (level - 1));      v = (row_pattern & pattern) / pattern;      h = (col_pattern & pattern) / pattern;            for (row = 0; row < (wavelet_block->blocksize / index_skip);           row++)	for (col = 0; col < (wavelet_block->blocksize / index_skip);             col++)          {	    //HL	    wavelet_block->block              [row][wavelet_block->blocksize / index_skip + col] =              QccVIDRDWTBlockExtractCoefficient(subband_pyramid,                                                offset_row + row,                                                num_cols / index_skip +                                                offset_col + col + h);	    //LH	    wavelet_block->block              [wavelet_block->blocksize / index_skip + row][col] =	      QccVIDRDWTBlockExtractCoefficient(subband_pyramid,                                                num_rows / index_skip +                                                offset_row + row + v,                                                offset_col + col);	    //HH	    wavelet_block->block              [wavelet_block->blocksize / index_skip + row]              [wavelet_block->blocksize / index_skip + col] =	      QccVIDRDWTBlockExtractCoefficient(subband_pyramid,                                                num_rows / index_skip +                                                offset_row + row + v,                                                num_cols / index_skip +                                                offset_col + col + h);	  }    }    //LL  for (row = 0; row < (wavelet_block->blocksize / index_skip); row++)    for (col = 0; col < (wavelet_block->blocksize / index_skip); col++)      wavelet_block->block[row][col] =        QccVIDRDWTBlockExtractCoefficient(subband_pyramid,                                          offset_row + row,                                          offset_col + col);    return(0);}static int QccVIDRDWTBlockInsertBlock(QccWAVSubbandPyramid *subband_pyramid,                                      QccVIDRDWTBlockBlock *wavelet_block,                                      int ref_offset_row,                                       int ref_offset_col,                                      int num_rows,                                      int num_cols,                                      int num_levels){  int row, col;  int offset_row, offset_col;  int level;  int index_skip;    index_skip = (1 << num_levels);  offset_row = ref_offset_row / index_skip;  offset_col = ref_offset_col / index_skip;  for (level = 1; level <= num_levels; level++)    {      index_skip = (1 << level);      offset_row = ref_offset_row / index_skip;      offset_col = ref_offset_col / index_skip;            for (row = 0; row < (wavelet_block->blocksize / index_skip);           row++)	for (col = 0; col < (wavelet_block->blocksize / index_skip);             col++)          {	    //HL	    subband_pyramid->matrix              [offset_row + row][num_cols / index_skip + offset_col + col] =	      wavelet_block->block              [row][wavelet_block->blocksize / index_skip + col];	    //LH	    subband_pyramid->matrix              [num_rows / index_skip + offset_row + row][offset_col + col] = 	      wavelet_block->block              [wavelet_block->blocksize / index_skip + row][col];	    //HH	    subband_pyramid->matrix              [num_rows / index_skip + offset_row + row]              [num_cols / index_skip + offset_col + col] =	      wavelet_block->block              [wavelet_block->blocksize / index_skip + row]              [wavelet_block->blocksize / index_skip + col];	  }    }    //LL  for (row = 0; row < (wavelet_block->blocksize / index_skip); row++)    for (col = 0; col < (wavelet_block->blocksize / index_skip); col++)      subband_pyramid->matrix[offset_row + row][offset_col + col] =        wavelet_block->block[row][col];    return(0);}static QccWAVSubbandPyramid ****QccVIDRDWTBlockRDWTPhasesAlloc(int num_levels,                                                               int num_rows,                                                               int num_cols,                                                               int subpixel_accuracy){  QccWAVSubbandPyramid ****rdwt_phases = NULL;  int num_phases;  int phase_row, phase_col;  int num_subpixels;  int subpixel_row, subpixel_col;    switch (subpixel_accuracy)    {    case QCCVID_ME_FULLPIXEL:      num_subpixels = 1;      break;    case QCCVID_ME_HALFPIXEL:      num_subpixels = 2;      break;    case QCCVID_ME_QUARTERPIXEL:      num_subpixels = 4;      break;    case QCCVID_ME_EIGHTHPIXEL:      num_subpixels = 8;      break;    default:      QccErrorAddMessage("(QccVIDRDWTBlockRDWTPhasesAlloc): Unrecognized subpixel accuracy");      return(NULL);    }  num_phases = (1 << num_levels);    if ((rdwt_phases =       (QccWAVSubbandPyramid ****)calloc(num_subpixels,                                         sizeof(QccWAVSubbandPyramid ***))) ==      NULL)    {      QccErrorAddMessage("(QccVIDRDWTBlockRDWTPhasesAlloc): Error allocating memory");      return(NULL);    }  for (subpixel_row = 0; subpixel_row < num_subpixels; subpixel_row++)    {      if ((rdwt_phases[subpixel_row] =           (QccWAVSubbandPyramid ***)calloc(num_subpixels,                                            sizeof(QccWAVSubbandPyramid **)))          == NULL)        {          QccErrorAddMessage("(QccVIDRDWTBlockRDWTPhasesAlloc): Error allocating memory");          return(NULL);        }            for (subpixel_col = 0; subpixel_col < num_subpixels; subpixel_col++)        {          if ((rdwt_phases[subpixel_row][subpixel_col] =