rwmh.c

QccPack-0.54-1 released (2007-04-30) is being developed and tested on Fedora Core Linux. QccPack pro

/*
 * 
 * QccPack: Quantization, compression, and coding libraries
 * Copyright (C) 1997-2007  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 QCCVIDRWMH_BLOCKSIZE 16

#define QCCVIDRWMH_WINDOWSIZE 15


typedef struct
{
  int intraframe;
  int motion_vector_bits;
  int intraframe_bits;
} QccVIDRWMHStatistics;


typedef struct
{
  int blocksize;
  int num_levels;
  QccMatrix *block;
} QccVIDRWMHBlock;


static void QccVIDRWMHBlockInitialize(QccVIDRWMHBlock *block)
{
  block->blocksize = 0;
  block->num_levels = 0;
  block->block = NULL;
}


static int QccVIDRWMHBlockAlloc(QccVIDRWMHBlock *block)
{
  if ((block->block =
       QccWAVWaveletRedundantDWT2DAlloc(block->blocksize,
                                        block->blocksize,
                                        block->num_levels)) ==
      NULL)
    {
      QccErrorAddMessage("(QccVIDRWMHBlockAlloc): Error calling QccWAVWaveletRedundantDWT2DAlloc()");
      return(1);
    }

  return(0);
}


static void QccVIDRWMHBlockFree(QccVIDRWMHBlock *block)
{
  QccWAVWaveletRedundantDWT2DFree(block->block,
                                  block->blocksize,
                                  block->num_levels);
}


static int QccVIDRWMHBlockExtract(QccVIDRWMHBlock *block,
                                  QccMatrix *frame_rdwt,
                                  int num_rows,
                                  int num_cols,
                                  int num_levels,
                                  double row,
                                  double col,
                                  int subpixel_accuracy)
{
  int num_subbands;
  int subband;
  int step;
  int row2, col2;
  int block_row, block_col;
  int extraction_row, extraction_col;

  switch (subpixel_accuracy)
    {
    case QCCVID_ME_FULLPIXEL:
      step = 1;
      row2 = (int)row;
      col2 = (int)col;
      break;
    case QCCVID_ME_HALFPIXEL:
      step = 2;
      row2 = (int)(row * 2);
      col2 = (int)(col * 2);
      break;
    case QCCVID_ME_QUARTERPIXEL:
      step = 4;
      row2 = (int)(row * 4);
      col2 = (int)(col * 4);
      break;
    case QCCVID_ME_EIGHTHPIXEL:
      step = 8;
      row2 = (int)(row * 8);
      col2 = (int)(col * 8);
      break;
    default:
      QccErrorAddMessage("(QccVIDRWMHBlockExtract): Unrecognized subpixel accuracy");
      return(1);
    }

  num_subbands = QccWAVSubbandPyramidNumLevelsToNumSubbands(num_levels);

  for (subband = 0; subband < num_subbands; subband++)
    for (block_row = 0; block_row < block->blocksize; block_row++)
      for (block_col = 0; block_col < block->blocksize; block_col++)
        {
          extraction_row = row2 + block_row * step;
          extraction_col = col2 + block_col * step;
          
          if (extraction_row < 0)
            extraction_row = 0;
          else
            if (extraction_row >= num_rows)
              extraction_row = num_rows - 1;
          
          if (extraction_col < 0)
            extraction_col = 0;
          else
            if (extraction_col >= num_cols)
              extraction_col = num_cols - 1;

          block->block[subband][block_row][block_col] =
            frame_rdwt[subband][extraction_row][extraction_col];
        }

  return(0);
}


static int QccVIDRWMHBlockInsert(QccVIDRWMHBlock *block,
                                 QccMatrix *frame_rdwt,
                                 int num_rows,
                                 int num_cols,
                                 int num_levels,
                                 int row,
                                 int col)
{
  int num_subbands;
  int subband;
  int block_row, block_col;

  num_subbands = QccWAVSubbandPyramidNumLevelsToNumSubbands(num_levels);

  for (subband = 0; subband < num_subbands; subband++)
    for (block_row = 0; block_row < block->blocksize; block_row++)
      for (block_col = 0; block_col < block->blocksize; block_col++)
        frame_rdwt[subband][block_row + row][block_col + col] =
          block->block[subband][block_row][block_col];

  return(0);
}


static double QccVIDRWMHBlockMAE(QccVIDRWMHBlock *block1,
                                 QccVIDRWMHBlock *block2,
                                 int num_levels)
{
  int num_subbands;
  int subband;
  int block_row, block_col;
  double mae = 0;
  double mae1;
  int level;
  int scale_factor;
  int end_subband;

  num_subbands = QccWAVSubbandPyramidNumLevelsToNumSubbands(num_levels);

  end_subband = 0;

  for (subband = num_subbands - 1; subband >= end_subband; subband--)
    {
      level =
        QccWAVSubbandPyramidCalcLevelFromSubband(subband,
                                                 num_levels);
      scale_factor = (1 << (level * 2 + 1));
      
      mae1 = 0.0;

      for (block_row = 0; block_row < block1->blocksize; block_row++)
        for (block_col = 0; block_col < block1->blocksize; block_col++)
          mae1 +=
            fabs(block1->block[subband][block_row][block_col] -
                 block2->block[subband][block_row][block_col]);

      mae += mae1 / scale_factor;
    }
  
  mae /=
    (double)(block1->blocksize * block1->blocksize * num_subbands);
  
  return(mae);
}


static void QccVIDRWMHPrintStatistics(const QccVIDRWMHStatistics
                                      *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 int QccVIDRWMHCreateReferenceFrame(QccMatrix reconstructed_frame,
                                          QccMatrix *reconstructed_frame_rdwt,
                                          QccMatrix *reference_frame_rdwt,
                                          int num_rows,
                                          int num_cols,
                                          int reference_num_rows,
                                          int reference_num_cols,
                                          int num_levels,
                                          int subpixel_accuracy,
                                          const QccFilter *filter1,
                                          const QccFilter *filter2,
                                          const QccFilter *filter3,
                                          const QccWAVWavelet *wavelet)
{
  int subband;
  int num_subbands;
  QccIMGImageComponent reconstructed_subband;
  QccIMGImageComponent reference_subband;

  QccIMGImageComponentInitialize(&reconstructed_subband);
  QccIMGImageComponentInitialize(&reference_subband);

  if (QccWAVWaveletRedundantDWT2D(reconstructed_frame,
                                  reconstructed_frame_rdwt,
                                  num_rows,
                                  num_cols,
                                  num_levels,
                                  wavelet))
    {
      QccErrorAddMessage("(QccVIDRWMHCreateReferenceFrame): Error calling QccWAVWaveletRedundantDWT2D()");
      return(1);
    }
  
  num_subbands = QccWAVSubbandPyramidNumLevelsToNumSubbands(num_levels);

  reconstructed_subband.num_rows = num_rows;
  reconstructed_subband.num_cols = num_cols;
  reference_subband.num_rows = reference_num_rows;
  reference_subband.num_cols = reference_num_cols;

  for (subband = 0; subband < num_subbands; subband++)
    {
      reconstructed_subband.image =
        (QccIMGImageArray)reconstructed_frame_rdwt[subband];
      reference_subband.image =
        (QccIMGImageArray)reference_frame_rdwt[subband];

      if (QccVIDMotionEstimationCreateReferenceFrame(&reconstructed_subband,
                                                     &reference_subband,
                                                     subpixel_accuracy,
                                                     filter1,
                                                     filter2,
                                                     filter3))
        {
          QccErrorAddMessage("(QccVIDRWMHCreateReferenceFrame): Error calling QccVIDMotionEstimationCreateReferenceFrame()");
          return(1);
        }
    }

  return(0);
}


static int QccVIDRWMHMotionEstimationSearch(QccMatrix *reference_frame_rdwt,
                                            int reference_num_rows,
                                            int reference_num_cols,
                                            int num_levels,
                                            QccVIDRWMHBlock *current_block,
                                            QccVIDRWMHBlock *reference_block,
                                            double search_row,
                                            double search_col,
                                            double window_size,
                                            double search_step,
                                            int subpixel_accuracy,
                                            double *mv_horizontal,
                                            double *mv_vertical)
{
  double u, v;
  double reference_frame_row, reference_frame_col;
  double current_mae;
  double min_mae = MAXDOUBLE;

  for (v = -window_size; v <= window_size; v += search_step)
    for (u = -window_size; u <= window_size; u += search_step)
      {
        reference_frame_row = search_row + v;
        reference_frame_col = search_col + u;

        if (QccVIDRWMHBlockExtract(reference_block,
                                   reference_frame_rdwt,
                                   reference_num_rows,
                                   reference_num_cols,
                                   num_levels,
                                   reference_frame_row,
                                   reference_frame_col,
                                   subpixel_accuracy))
          {
            QccErrorAddMessage("(QccVIDRWMHMotionEstimationSearch): Error calling QccVIDRWMHBlockExtract()");
            return(1);
          }
        
        current_mae =
          QccVIDRWMHBlockMAE(current_block,
                             reference_block,
                             num_levels);
        
        if (current_mae < min_mae)
          {
            min_mae = current_mae;
            *mv_horizontal = u;
            *mv_vertical = v;
          }
      }

  return(0);
}


static int QccVIDRWMHMotionEstimation(QccMatrix *current_frame_rdwt,
                                      QccMatrix *reference_frame_rdwt,
                                      int num_rows,
                                      int num_cols,
                                      int reference_num_rows,
                                      int reference_num_cols,
                                      int blocksize,
                                      int num_levels,
                                      int subpixel_accuracy,
                                      QccIMGImageComponent *horizontal_motion,
                                      QccIMGImageComponent *vertical_motion)
{
  int return_value;
  int row, col;
  int num_subbands;
  double u = 0;
  double v = 0;
  QccVIDRWMHBlock current_block;
  QccVIDRWMHBlock reference_block;
  double search_row, search_col;
  double current_window_size;
  double final_search_step;
  double current_search_step;
  int mv_row, mv_col;

  QccVIDRWMHBlockInitialize(&current_block);
  QccVIDRWMHBlockInitialize(&reference_block);

  switch (subpixel_accuracy)
    {
    case QCCVID_ME_FULLPIXEL:
      final_search_step = 1.0;
      break;
    case QCCVID_ME_HALFPIXEL:
      final_search_step = 0.5;
      break;
    case QCCVID_ME_QUARTERPIXEL:
      final_search_step = 0.25;
      break;
    case QCCVID_ME_EIGHTHPIXEL:
      final_search_step = 0.125;
      break;
    default:
      QccErrorAddMessage("(QccVIDRWMHMotionEstimation): Unrecognized subpixel accuracy");
      goto Error;
    }

  current_block.blocksize = blocksize;
  current_block.num_levels = num_levels;
  if (QccVIDRWMHBlockAlloc(&current_block))
    {
      QccErrorAddMessage("(QccVIDRWMHMotionEstimation): Error calling QccVIDRWMHBlockAlloc()");
      goto Error;
    }
  reference_block.blocksize = blocksize;
  reference_block.num_levels = num_levels;
  if (QccVIDRWMHBlockAlloc(&reference_block))
    {
      QccErrorAddMessage("(QccVIDRWMHMotionEstimation): Error calling QccVIDRWMHBlockAlloc()");
      goto Error;
    }

  num_subbands = QccWAVSubbandPyramidNumLevelsToNumSubbands(num_levels);