rdwtblock.c

QccPack implementation in C

/*
 * 
 * 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 <[e-mail]> based on
 *  the originally developed algorithm and code by Suxia Cui.
 */

#include "libQccPack.h"


#ifdef LIBQCCPACKSPIHT_H


#define QCCVIDRDWTBLOCK_WINDOWSIZE 15


typedef 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] =