rdwt.c

spiht for linux this is used to decod and encode vedio i wich all enjoy

/*
 * 
 * QccPack: Quantization, compression, and coding libraries
 * Copyright (C) 1997-2009  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 QccWAVWaveletRedundantAnalysis1D(const QccVector input_signal,
                                            QccVector output_signal_low,
                                            QccVector output_signal_high,
                                            int signal_length,
                                            const QccWAVWavelet *wavelet)
{
  int return_value;

  switch (wavelet->implementation)
    {
    case QCCWAVWAVELET_IMPLEMENTATION_FILTERBANK:
      if (QccWAVFilterBankRedundantAnalysis(input_signal,
                                            output_signal_low,
                                            output_signal_high,
                                            signal_length,
                                            &wavelet->filter_bank,
                                            wavelet->boundary))
        {
          QccErrorAddMessage("(QccWAVWaveletRedundantAnalysis1D): Error calling QccWAVFilterBankRedundantAnalysis()");
          goto Error;
        }
      break;
    case QCCWAVWAVELET_IMPLEMENTATION_LIFTED:
      if (QccWAVLiftingRedundantAnalysis(input_signal,
                                         output_signal_low,
                                         output_signal_high,
                                         signal_length,
                                         &wavelet->lifting_scheme,
                                         wavelet->boundary))
        {
          QccErrorAddMessage("(QccWAVWaveletRedundantAnalysis1D): Error calling QccWAVLiftingRedundantAnalysis()");
          goto Error;
        }
      break;
    default:
      QccErrorAddMessage("(QccWAVWaveletRedundantAnalysis1D): Undefined implementation (%d)",
                         wavelet->implementation);
      goto Error;
    }

  return_value = 0;
  goto Return;
 Error:
  return_value = 1;
 Return:
  return(return_value);
}


static int QccWAVWaveletRedundantSynthesis1D(const QccVector input_signal_low,
                                             const QccVector input_signal_high,
                                             QccVector output_signal,
                                             int signal_length,
                                             const QccWAVWavelet *wavelet)
{
  int return_value;

  switch (wavelet->implementation)
    {
    case QCCWAVWAVELET_IMPLEMENTATION_FILTERBANK:
      if (QccWAVFilterBankRedundantSynthesis(input_signal_low,
                                             input_signal_high,
                                             output_signal,
                                             signal_length,
                                             &wavelet->filter_bank,
                                             wavelet->boundary))
        {
          QccErrorAddMessage("(QccWAVWaveletRedundantSynthesis1D): Error calling QccWAVFilterBankRedundantSynthesis()");
          goto Error;
        }
      break;
    case QCCWAVWAVELET_IMPLEMENTATION_LIFTED:
      if (QccWAVLiftingRedundantSynthesis(input_signal_low,
                                          input_signal_high,
                                          output_signal,
                                          signal_length,
                                          &wavelet->lifting_scheme,
                                          wavelet->boundary))
        {
          QccErrorAddMessage("(QccWAVWaveletRedundantSynthesis1D): Error calling QccWAVLiftingRedundantSynthesis()");
          goto Error;
        }
      break;
    default:
      QccErrorAddMessage("(QccWAVWaveletRedundantSynthesis1D): Undefined implementation (%d)",
                         wavelet->implementation);
      goto Error;
    }

  return_value = 0;
  goto Return;
 Error:
  return_value = 1;
 Return:
  return(return_value);
}


int QccWAVWaveletRedundantDWT1D(const QccVector input_signal,
                                QccMatrix output_signals,
                                int signal_length,
                                int num_scales,
                                const QccWAVWavelet *wavelet)
{
  int return_value = 0;
  QccVector temp_signal = NULL;
  QccMatrix output_signals2 = NULL;
  int even_length;
  int odd_length;
  int scale;
  int max_scales;
  
  if (input_signal == NULL)
    return(0);
  if (output_signals == NULL)
    return(0);
  if (wavelet == NULL)
    return(0);
  
  if (num_scales <= 0)
    return(0);
  
  if (signal_length <= 0)
    return(0);
  
  max_scales = (int)floor(QccMathLog2(signal_length));
  
  if (max_scales < num_scales)
    {
      QccErrorAddMessage("(QccWAVWaveletRedundantDWT1D): %d transform scales not supported for length %d signal - use %d or fewer scales",
                         num_scales,
                         signal_length,
                         max_scales);
      goto Error;
    }
  
  if (QccWAVWaveletRedundantAnalysis1D(input_signal,
                                       output_signals[0],
                                       output_signals[num_scales],
                                       signal_length,
                                       wavelet))
    {
      QccErrorAddMessage("(QccWAVWaveletRedundantDWT1D): Error calling QccWAVWaveletRedundantAnalysis1D()");
      goto Error;
    }
  
  if (num_scales == 1)
    goto Return;
  
  if ((temp_signal = QccVectorAlloc(signal_length)) == NULL)
    {
      QccErrorAddMessage("(QccWAVWaveletRedundantDWT1D): Error calling QccVectorAlloc()");
      goto Error;
    }
  
  even_length = QccWAVWaveletDWTSubbandLength(signal_length,
                                              1, 0, 0, 0);
  odd_length = QccWAVWaveletDWTSubbandLength(signal_length,
                                             1, 1, 0, 0);
  
  if ((output_signals2 =
       (QccMatrix)malloc(sizeof(QccVector) * num_scales)) == NULL)
    {
      QccErrorAddMessage("(QccWAVWaveletRedundantDWT1D): Error allocating memory");
      goto Error;
    }
  
  for (scale = 0; scale < num_scales; scale++)
    output_signals2[scale] = &(output_signals[scale][even_length]);
  
  if (QccWAVWaveletLWT(output_signals[0],
                       temp_signal,
                       signal_length,
                       0, 0))
    {
      QccErrorAddMessage("(QccWAVWaveletRedundantDWT1D): Error calling QccWAVWaveletLWT()");
      goto Error;
    }
  
  if (QccWAVWaveletRedundantDWT1D(temp_signal,
                                  output_signals,
                                  even_length,
                                  num_scales - 1,
                                  wavelet))
    {
      QccErrorAddMessage("(QccWAVWaveletRedundantDWT1D): Error calling QccWAVWaveletRedundantDWT1D()");
      goto Error;
    }
  
  if (QccWAVWaveletRedundantDWT1D(&(temp_signal[even_length]),
                                  output_signals2,
                                  odd_length,
                                  num_scales - 1,
                                  wavelet))
    {
      QccErrorAddMessage("(QccWAVWaveletRedundantDWT1D): Error calling QccWAVWaveletRedundantDWT1D()");
      goto Error;
    }
  
  for (scale = 0; scale < num_scales; scale++)
    {
      if (QccWAVWaveletInverseLWT(output_signals[scale],
                                  temp_signal,
                                  signal_length,
                                  0, 0))
        {
          QccErrorAddMessage("(QccWAVWaveletRedundantDWT1D): Error calling QccWAVWaveletInverseLWT()");
          goto Error;
        }
      if (QccVectorCopy(output_signals[scale], temp_signal, signal_length))
        {
          QccErrorAddMessage("(QccWAVWaveletRedundantDWT1D): Error calling QccVector()");
          goto Error;
        }
    }
  
  return_value = 0;
  goto Return;
 Error:
  return_value = 1;
 Return:
  QccVectorFree(temp_signal);
  if (output_signals2 != NULL)
    QccFree(output_signals2);
  return(return_value);
}


void QccWAVWaveletRedundantDWT1DSubsampleGetPattern(int subband,
                                                    int num_scales,
                                                    int pattern,
                                                    int *index_skip,
                                                    int *start_index,
                                                    const QccWAVWavelet
                                                    *wavelet)
{
  int scale = (!subband) ? num_scales : num_scales - (subband - 1);

  *index_skip = 1 << scale;

  *start_index = pattern;

  if (QccWAVWaveletBiorthogonal(wavelet) && subband)
    *start_index += (1 << (scale - 1));

  *start_index = QccMathModulus(*start_index, *index_skip);
}


int QccWAVWaveletRedundantDWT1DSubsample(const QccMatrix input_signals,
                                         QccVector output_signal,
                                         int signal_length,
                                         int num_scales,
                                         int subsample_pattern,
                                         const QccWAVWavelet *wavelet)
{
  int start_index;
  int index_skip;
  int subband;
  int index1;
  int index2;
  int index3;
  int max_scales;
  int baseband_phase_group_size;
  int num_subbands;
  int scale;

  if (input_signals == NULL)
    return(0);
  if (output_signal == NULL)
    return(0);
  if (wavelet == NULL)
    return(0);
  if (!signal_length)
    return(0);
  
  max_scales = (int)floor(QccMathLog2(signal_length));
  
  if (max_scales < num_scales)
    {
      QccErrorAddMessage("(QccWAVWaveletRedundantDWT1DSubsample): %d transform scales not supported for length %d signal - use %d or fewer scales",
                         num_scales,
                         signal_length,
                         max_scales);
      return(1);
    }
  
  index1 = 0;

  num_subbands = num_scales + 1;

  baseband_phase_group_size = (1 << num_scales);

  if (subsample_pattern > (baseband_phase_group_size - 1))
    {
      QccErrorAddMessage("(QccWAVWaveletRedundantDWT1DSubsample): subband_pattern can take on only values 0 to %d for a transform of %d scales",
                         baseband_phase_group_size - 1,
                         num_scales);
      return(1);
    }

  index_skip = (1 << num_scales);
  
  for (subband = 0; subband < num_subbands; subband++)
    {
      QccWAVWaveletRedundantDWT1DSubsampleGetPattern(subband,
                                                     num_scales,
                                                     subsample_pattern,
                                                     &index_skip,
                                                     &start_index,
                                                     wavelet);

      scale = (!subband) ? num_scales : num_scales - (subband - 1);

      for (index2 = 0, index3 = 0;
           index3 < QccWAVWaveletDWTSubbandLength(signal_length,
                                                  scale,
                                                  subband,
                                                  0, 0);
           index1++, index2 += index_skip, index3++)
        if (wavelet->boundary == QCCWAVWAVELET_BOUNDARY_SYMMETRIC_EXTENSION)
          output_signal[index1] = 
            input_signals[subband]
            [QccFilterCalcSymmetricExtension(index2 + start_index,
                                             signal_length,
                                             QCCFILTER_SYMMETRICWHOLE)];
        else
          output_signal[index1] = 
            input_signals[subband][QccMathModulus(index2 + start_index,
                                                  signal_length)];
    }
  
  return(0);
}


static
int QccWAVWaveletInverseRedundantDWT1DRecursion(QccMatrix input_signals,
                                                QccVector output_signal,
                                                int signal_length,
                                                int num_scales,
                                                const QccWAVWavelet *wavelet)
{
  int return_value;
  int max_scales;
  QccVector temp_signal = NULL;
  QccMatrix input_signals2 = NULL;
  int scale;
  int even_length;
  int odd_length;
  
  if (input_signals == NULL)
    return(0);
  if (output_signal == NULL)
    return(0);
  if (wavelet == NULL)
    return(0);
  if (signal_length <= 0)
    return(0);
  if (num_scales <= 0)
    return(0);
  
  max_scales = (int)floor(QccMathLog2(signal_length));
  
  if (max_scales < num_scales)
    {
      QccErrorAddMessage("(QccWAVWaveletInverseRedundantDWT1DRecursion): %d transform scales not supported for length %d signal - use %d or fewer scales",
                         num_scales,
                         signal_length,
                         max_scales);
      goto Error;
    }
  
  if (num_scales != 1)
    {
      if ((temp_signal = QccVectorAlloc(signal_length)) == NULL)
        {
          QccErrorAddMessage("(QccWAVWaveletInverseRedundantDWT1DRecursion): Error calling QccVectorAlloc()");
          goto Error;
        }
      
      for (scale = 0; scale < num_scales; scale++)
        {
          if (QccWAVWaveletLWT(input_signals[scale],
                               temp_signal,
                               signal_length,
                               0, 0))
            {
              QccErrorAddMessage("(QccWAVWaveletInverseRedundantDWT1DRecursion): Error calling QccWAVWaveletLWT()");
              goto Error;
            }
          if (QccVectorCopy(input_signals[scale], temp_signal, signal_length))
            {
              QccErrorAddMessage("(QccWAVWaveletInverseRedundantDWT1DRecursion): Error calling QccVector()");
              goto Error;
            }