dwt.c

来自「QccPack implementation in C」· C语言 代码 · 共 879 行 · 第 1/2 页

/*
 * 
 * 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"


static int QccWAVWaveletDWTSubbandPhase(int level, int signal_origin,
                                        int subsample_pattern)
{
  int phase = -(signal_origin % (1 << level));
  
  phase &= (1 << (level - 1));

  phase ^= (subsample_pattern & (1 << (level - 1))) << 1;

  if (phase)
    return(QCCWAVWAVELET_PHASE_ODD);
  else
    return(QCCWAVWAVELET_PHASE_EVEN);
}


int QccWAVWaveletDWTSubbandLength(int original_length, int level,
                                  int highband, int signal_origin,
                                  int subsample_pattern)
{
  int scale;
  int length;
  int phase;

  length = original_length;

  if (!level)
    {
      length = (highband) ? 0 : original_length;
      return(length);
    }

  for (scale = 0; scale < level - 1; scale++)
    {
      phase =
        QccWAVWaveletDWTSubbandPhase(scale + 1, signal_origin,
                                     subsample_pattern);
      if (length % 2)
        length =
          (phase == QCCWAVWAVELET_PHASE_ODD) ?
          (length >> 1) : (length >> 1) + 1;
      else
        length >>= 1;
      phase >>= 1;
    }

  phase =
    QccWAVWaveletDWTSubbandPhase(level, signal_origin, subsample_pattern);

  if (length % 2)
    length =
      (((phase == QCCWAVWAVELET_PHASE_EVEN) && highband) ||
       ((phase == QCCWAVWAVELET_PHASE_ODD) && !highband)) ?
      (length >> 1) : (length >> 1) + 1;
  else
    length >>= 1;

  return(length);
}


int QccWAVWaveletDWT1D(const QccVector input_signal,
                       QccVector output_signal,
                       int signal_length,
                       int signal_origin,
                       int subsample_pattern,
                       int num_scales,
                       const QccWAVWavelet *wavelet)
{
  int return_value;
  int scale;
  QccVector temp_vector = NULL;
  int current_length;
  int phase;

  if (input_signal == NULL)
    return(0);
  if (output_signal == NULL)
    return(0);
  if (wavelet == NULL)
    return(0);

  if (num_scales < 0)
    return(0);
  if (num_scales == 0)
    {
      QccVectorCopy(output_signal, input_signal, signal_length);
      return(0);
    }

  if ((temp_vector = QccVectorAlloc(signal_length)) == NULL)
    {
      QccErrorAddMessage("(QccWAVWaveletDWT1D): Error calling QccVectorAlloc()");
      goto QccError;
    }

  if (QccVectorCopy(temp_vector, input_signal, signal_length))
    {
      QccErrorAddMessage("(QccWAVWaveletDWT1D): Error calling QccVectorCopy()");
      goto QccError;
    }
    
  current_length = signal_length;

  for (scale = 0; scale < num_scales; scale++)
    {
      phase =
        QccWAVWaveletDWTSubbandPhase(scale + 1, signal_origin,
                                     subsample_pattern);

      if (QccWAVWaveletAnalysis1D(temp_vector,
                                  output_signal,
                                  current_length,
                                  phase,
                                  wavelet))
        {
          QccErrorAddMessage("(QccWAVWaveletDWT1D): Error calling QccWAVWaveletAnalysis1D()");
          goto QccError;
        }

      current_length =
        QccWAVWaveletDWTSubbandLength(signal_length, scale + 1, 0,
                                      signal_origin, subsample_pattern);

      if (QccVectorCopy(temp_vector, output_signal, current_length))
        {
          QccErrorAddMessage("(QccWAVWaveletDWT1D): Error calling QccVectorCopy()");
          goto QccError;
        }
    }

  return_value = 0;
  goto QccReturn;
 QccError:
  return_value = 1;
 QccReturn:
  QccVectorFree(temp_vector);
  return(return_value);
}


int QccWAVWaveletInverseDWT1D(const QccVector input_signal,
                              QccVector output_signal,
                              int signal_length,
                              int signal_origin,
                              int subsample_pattern,
                              int num_scales,
                              const QccWAVWavelet *wavelet)
{
  int return_value;
  int scale;
  QccVector temp_vector = NULL;
  int current_length;
  int phase;

  if (input_signal == NULL)
    return(0);
  if (output_signal == NULL)
    return(0);
  if (wavelet == NULL)
    return(0);

  if (num_scales < 0)
    return(0);
  if (num_scales == 0)
    {
      QccVectorCopy(output_signal, input_signal, signal_length);
      return(0);
    }

  if ((temp_vector = QccVectorAlloc(signal_length)) == NULL)
    {
      QccErrorAddMessage("(QccWAVWaveletInverseDWT1D): Error calling QccVectorAlloc()");
      goto QccError;
    }

  if (QccVectorCopy(temp_vector, input_signal, signal_length))
    {
      QccErrorAddMessage("(QccWAVWaveletInverseDWT1D): Error calling QccVectorCopy()");
      goto QccError;
    }
    
  for (scale = num_scales - 1; scale >= 0; scale--)
    {
      current_length = 
        QccWAVWaveletDWTSubbandLength(signal_length, scale,
                                      0, signal_origin,
                                      subsample_pattern);

      phase =
        QccWAVWaveletDWTSubbandPhase(scale + 1, signal_origin,
                                     subsample_pattern);

      if (QccWAVWaveletSynthesis1D(temp_vector,
                                   output_signal,
                                   current_length,
                                   phase,
                                   wavelet))
        {
          QccErrorAddMessage("(QccWAVWaveletInverseDWT1D): Error calling QccWAVWaveletSynthesis1D()");
          goto QccError;
        }

      if (QccVectorCopy(temp_vector, output_signal, current_length))
        {
          QccErrorAddMessage("(QccWAVWaveletInverseDWT1D): Error calling QccVectorCopy()");
          goto QccError;
        }
    }

  return_value = 0;
  goto QccReturn;
 QccError:
  return_value = 1;
 QccReturn:
  QccVectorFree(temp_vector);
  return(return_value);
}


int QccWAVWaveletDWT2D(const QccMatrix input_matrix,
                       QccMatrix output_matrix,
                       int num_rows,
                       int num_cols,
                       int origin_row,
                       int origin_col,
                       int subsample_pattern_row,
                       int subsample_pattern_col,
                       int num_scales,
                       const QccWAVWavelet *wavelet)
{
  int return_value;
  int scale;
  QccMatrix temp_matrix = NULL;
  int original_num_rows = 0;
  int original_num_cols = 0;
  int phase_row;
  int phase_col;

  if (input_matrix == NULL)
    return(0);
  if (output_matrix == NULL)
    return(0);
  if (wavelet == NULL)
    return(0);

  if (num_scales < 0)
    return(0);
  if (num_scales == 0)
    {
      QccMatrixCopy(output_matrix, input_matrix, 
                    num_rows, num_cols);
      return(0);
    }

  original_num_rows = num_rows;
  original_num_cols = num_cols;

  if ((temp_matrix = QccMatrixAlloc(num_rows, num_cols)) == NULL)
    {
      QccErrorAddMessage("(QccWAVWaveletDWT2D): Error calling QccMatrixAlloc()");
      goto QccError;
    }

  if (QccMatrixCopy(temp_matrix, input_matrix, num_rows, num_cols))
    {
      QccErrorAddMessage("(QccWAVWaveletDWT2D): Error calling QccMatrixCopy()");
      goto QccError;
    }
    
  for (scale = 0; scale < num_scales; scale++)
    {
      phase_row = QccWAVWaveletDWTSubbandPhase(scale + 1, origin_row,
                                               subsample_pattern_row);
      phase_col = QccWAVWaveletDWTSubbandPhase(scale + 1, origin_col,
                                               subsample_pattern_col);

      if (QccWAVWaveletAnalysis2D(temp_matrix,
                                  output_matrix,
                                  num_rows,
                                  num_cols,
                                  phase_row,
                                  phase_col,
                                  wavelet))
        {
          QccErrorAddMessage("(QccWAVWaveletDWT2D): Error calling QccWAVWaveletAnalysis2D()");
          goto QccError;
        }

      num_rows =
        QccWAVWaveletDWTSubbandLength(original_num_rows, scale + 1,
                                      0, origin_row, subsample_pattern_row);
      num_cols =
        QccWAVWaveletDWTSubbandLength(original_num_cols, scale + 1,
                                      0, origin_col, subsample_pattern_col);

      if (QccMatrixCopy(temp_matrix, output_matrix, num_rows, num_cols))
        {
          QccErrorAddMessage("(QccWAVWaveletDWT2D): Error calling QccMatrixCopy()");
          goto QccError;
        }
    }

  return_value = 0;
  goto QccReturn;
 QccError:
  return_value = 1;
 QccReturn:
  QccMatrixFree(temp_matrix, original_num_rows);
  return(return_value);
}


int QccWAVWaveletInverseDWT2D(const QccMatrix input_matrix,
                              QccMatrix output_matrix,
                              int num_rows,
                              int num_cols,
                              int origin_row,
                              int origin_col,
                              int subsample_pattern_row,
                              int subsample_pattern_col,
                              int num_scales,
                              const QccWAVWavelet *wavelet)
{
  int return_value;
  int scale;
  QccMatrix temp_matrix = NULL;
  int original_num_rows = 0;
  int original_num_cols = 0;
  int phase_row;
  int phase_col;

  if (input_matrix == NULL)
    return(0);
  if (output_matrix == NULL)
    return(0);
  if (wavelet == NULL)
    return(0);

  if (num_scales < 0)
    return(0);
  if (num_scales == 0)
    {
      QccMatrixCopy(output_matrix, input_matrix, num_rows, num_cols);
      return(0);
    }

  original_num_rows = num_rows;
  original_num_cols = num_cols;

  if ((temp_matrix = QccMatrixAlloc(num_rows, num_cols)) == NULL)
    {
      QccErrorAddMessage("(QccWAVWaveletInverseDWT2D): Error calling QccMatrixAlloc()");
      goto QccError;
    }

  if (QccMatrixCopy(temp_matrix, input_matrix, num_rows, num_cols))
    {
      QccErrorAddMessage("(QccWAVWaveletInverseDWT2D): Error calling QccMatrixCopy()");
      goto QccError;
    }
    
  for (scale = num_scales - 1; scale >= 0; scale--)
    {
      num_rows =
        QccWAVWaveletDWTSubbandLength(original_num_rows, scale,
                                      0, origin_row, subsample_pattern_row);
      num_cols =
        QccWAVWaveletDWTSubbandLength(original_num_cols, scale,
                                      0, origin_col, subsample_pattern_col);

      phase_row = QccWAVWaveletDWTSubbandPhase(scale + 1, origin_row,
                                               subsample_pattern_row);
      phase_col = QccWAVWaveletDWTSubbandPhase(scale + 1, origin_col,
                                               subsample_pattern_col);

      if (QccWAVWaveletSynthesis2D(temp_matrix,
                                   output_matrix,
                                   num_rows,
                                   num_cols,
                                   phase_row,
                                   phase_col,
                                   wavelet))
        {
          QccErrorAddMessage("(QccWAVWaveletInverseDWT2D): Error calling QccWAVWaveletSynthesis2D()");
          goto QccError;
        }

      if (QccMatrixCopy(temp_matrix, output_matrix, num_rows, num_cols))
        {
          QccErrorAddMessage("(QccWAVWaveletInverseDWT2D): Error calling QccMatrixCopy()");
          goto QccError;
        }
    }

  return_value = 0;
  goto QccReturn;
 QccError:
  return_value = 1;
 QccReturn:
  QccMatrixFree(temp_matrix, original_num_rows);
  return(return_value);
}


int QccWAVWaveletDyadicDWT3D(const QccVolume input_volume,
                             QccVolume output_volume,
                             int num_frames,
                             int num_rows,
                             int num_cols,
                             int origin_frame,
                             int origin_row,
                             int origin_col,
                             int subsample_pattern_frame,