matrix.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.
 * 
 */


#include "libQccPack.h"


QccMatrix QccMatrixAlloc(int num_rows, int num_cols)
{
  QccMatrix matrix;
  int row;
  
  if ((num_rows <= 0) || (num_cols <= 0))
    return(NULL);

  if ((matrix = 
       (QccMatrix)malloc(sizeof(QccVector)*num_rows)) == NULL)
    {
      QccErrorAddMessage("(QccMatrixAlloc): Error allocating memory");
      return(NULL);
    }
  
  for (row = 0; row < num_rows; row++)
    if ((matrix[row] = QccVectorAlloc(num_cols)) == NULL)
      {
        QccErrorAddMessage("(QccMatrixAlloc): Error allocating memory");
        QccFree(matrix);
        return(NULL);
      }
  
  return(matrix);
}


void QccMatrixFree(QccMatrix matrix, int num_rows)
{
  int row;
  
  if (matrix != NULL)
    {
      for (row = 0; row < num_rows; row++)
        QccVectorFree(matrix[row]);

      QccFree(matrix);
    }
}


int QccMatrixZero(QccMatrix matrix, int num_rows, int num_cols)
{
  int row, col;

  if (matrix == NULL)
    return(0);

  for (row = 0; row < num_rows; row++)
    for (col = 0; col < num_cols; col++)
      matrix[row][col] = 0;

  return(0);
}


int QccMatrixCopy(QccMatrix matrix1, const QccMatrix matrix2,
                  int num_rows, int num_cols)
{
  int row, col;

  if ((matrix1 == NULL) || 
      (matrix2 == NULL) ||
      (num_rows <= 0) ||
      (num_cols <= 0))
    return(0);

  for (row = 0; row < num_rows; row++)
    for (col = 0; col < num_cols; col++)
      matrix1[row][col] = matrix2[row][col];

  return(0);
}


double QccMatrixMaxValue(const QccMatrix matrix,
                         int num_rows, int num_cols)
{
  int row;
  double max_value = -MAXDOUBLE;
  double tmp;
  
  if (matrix != NULL)
    for (row = 0; row < num_rows; row++)
      {
        tmp = QccVectorMaxValue(matrix[row], num_cols, NULL);
        if (tmp > max_value)
          max_value = tmp;
      }

  return(max_value);
}


double QccMatrixMinValue(const QccMatrix matrix,
                         int num_rows, int num_cols)
{
  int row;
  double min_value = MAXDOUBLE;
  double tmp;
  
  if (matrix != NULL)
    for (row = 0; row < num_rows; row++)
      {
        tmp = QccVectorMinValue(matrix[row], num_cols, NULL);
        if (tmp < min_value)
          min_value = tmp;
      }

  return(min_value);
}


int QccMatrixPrint(const QccMatrix matrix, int num_rows, int num_cols)
{
  int row, col;

  if (matrix == NULL)
    return(0);

  for (row = 0; row < num_rows; row++)
    {
      for (col = 0; col < num_cols; col++)
        printf("%g ", matrix[row][col]);
      printf("\n");
    }

  return(0);
}


int QccMatrixTranspose(const QccMatrix matrix1, QccMatrix matrix2,
                       int num_rows, int num_cols)
{
  int row, col;

  if (matrix1 == NULL)
    return(0);
  if (matrix2 == NULL)
    return(0);

  for (row = 0; row < num_rows; row++)
    for (col = 0; col < num_cols; col++)
      matrix2[col][row] = matrix1[row][col];

  return(0);
}


int QccMatrixAdd(QccMatrix matrix1, const QccMatrix matrix2,
                 int num_rows, int num_cols)
{
  int row;

  if (matrix1 == NULL)
    return(0);
  if (matrix2 == NULL)
    return(0);

  for (row = 0; row < num_rows; row++)
    if (QccVectorAdd(matrix1[row], matrix2[row], num_cols))
      {
        QccErrorAddMessage("(QccMatrixAdd): Error calling QccVectorAdd()");
        return(1);
      }

  return(0);
}


int QccMatrixSubtract(QccMatrix matrix1, const QccMatrix matrix2,
                      int num_rows, int num_cols)
{
  int row;

  if (matrix1 == NULL)
    return(0);
  if (matrix2 == NULL)
    return(0);

  for (row = 0; row < num_rows; row++)
    if (QccVectorSubtract(matrix1[row], matrix2[row], num_cols))
      {
        QccErrorAddMessage("(QccMatrixSubtract): Error calling QccVectorSubtract()");
        return(1);
      }

  return(0);
}


double QccMatrixMean(const QccMatrix matrix,
                     int num_rows, int num_cols)
{
  double sum = 0.0;
  int row, col;
  
  for (row = 0; row < num_rows; row++)
    for (col = 0; col < num_cols; col++)
      sum += matrix[row][col];
  
  sum /= num_rows * num_cols;
  
  return(sum);
}


double QccMatrixVariance(const QccMatrix matrix,
                         int num_rows, int num_cols)
{
  double variance = 0.0;
  double mean;
  int row, col;

  mean = QccMatrixMean(matrix, num_rows, num_cols);

  for (row = 0; row < num_rows; row++)
    for (col = 0; col < num_cols; col++)
      variance +=
        (matrix[row][col] - mean) *
        (matrix[row][col] - mean);

  variance /= (num_rows * num_cols);

  return(variance);
}


double QccMatrixMaxSignalPower(const QccMatrix matrix,
                               int num_rows, int num_cols)
{
  int current_row, component;
  double max_power = -1;
  double power;
  
  for (current_row = 0; current_row < num_rows; current_row++)
    {
      power = 0;
      for (component = 0; component < num_cols; component++)
        power += matrix[current_row][component]*
          matrix[current_row][component];
      if (power > max_power)
        max_power = power;
    }
  
  return(max_power);
  
}


int QccMatrixVectorMultiply(const QccMatrix matrix,
                            const QccVector vector1,
                            QccVector vector2,
                            int num_rows, int num_cols)
{
  int row;

  if (matrix == NULL)
    return(0);
  if (vector1 == NULL)
    return(0);
  if (vector2 == NULL)
    return(0);

  for (row = 0; row < num_rows; row++)
    vector2[row] =
      QccVectorDotProduct(matrix[row], vector1, num_cols);

  return(0);
}


int QccMatrixDCT(const QccMatrix input_block,
                 QccMatrix output_block,
                 int num_rows, int num_cols)
{
  int return_value;
  QccVector col_vector1 = NULL;
  QccVector col_vector2 = NULL;
  int row, col;
  
  if ((input_block == NULL) ||
      (output_block == NULL))
    return(0);
  if ((num_rows <= 0) || (num_cols <= 0))
    return(0);
  
  if ((col_vector1 = QccVectorAlloc(num_rows)) == NULL)
    {
      QccErrorAddMessage("(QccMatrixDCT): Error calling QccVectorAlloc()");
      goto Error;
    }
  if ((col_vector2 = QccVectorAlloc(num_rows)) == NULL)
    {
      QccErrorAddMessage("(QccMatrixDCT): Error calling QccVectorAlloc()");
      goto Error;
    }
  
  for (row = 0; row < num_rows; row++)
    if (QccVectorDCT(input_block[row],
                     output_block[row],
                     num_cols))
      {
        QccErrorAddMessage("(QccMatrixDCT): Error calling QccVectorDCT()");
        goto Error;
      }
  
  for (col = 0; col < num_cols; col++)
    {
      for (row = 0; row < num_rows; row++)
        col_vector1[row] = output_block[row][col];
      if (QccVectorDCT(col_vector1,
                       col_vector2,
                       num_rows))
        {
          QccErrorAddMessage("(QccMatrixDCT): Error calling QccVectorDCT()");
          goto Error;
        }
      for (row = 0; row < num_rows; row++)
        output_block[row][col] = col_vector2[row];
    }
  
  return_value = 0;
  goto Return;
 Error:
  return_value = 1;
 Return:
  QccVectorFree(col_vector1);
  QccVectorFree(col_vector2);
  return(return_value);
}


int QccMatrixInverseDCT(const QccMatrix input_block,
                        QccMatrix output_block,
                        int num_rows, int num_cols)
{
  int return_value;
  QccVector col_vector1 = NULL;
  QccVector col_vector2 = NULL;
  int row, col;
  
  if ((input_block == NULL) ||
      (output_block == NULL))
    return(0);
  if ((num_rows <= 0) || (num_cols <= 0))
    return(0);
  
  if ((col_vector1 = QccVectorAlloc(num_rows)) == NULL)
    {
      QccErrorAddMessage("(QccMatrixInverseDCT): Error calling QccVectorAlloc()");
      goto Error;
    }
  if ((col_vector2 = QccVectorAlloc(num_rows)) == NULL)
    {
      QccErrorAddMessage("(QccMatrixInverseDCT): Error calling QccVectorAlloc()");
      goto Error;
    }
  
  for (row = 0; row < num_rows; row++)
    if (QccVectorInverseDCT(input_block[row],
                            output_block[row],
                            num_cols))
      {
        QccErrorAddMessage("(QccMatrixInverseDCT): Error calling QccVectorInverseDCT()");
        goto Error;
      }
  
  for (col = 0; col < num_cols; col++)
    {
      for (row = 0; row < num_rows; row++)
        col_vector1[row] = output_block[row][col];
      if (QccVectorInverseDCT(col_vector1,
                              col_vector2,
                              num_rows))
        {
          QccErrorAddMessage("(QccMatrixInverseDCT): Error calling QccVectorInverseDCT()");
          goto Error;
        }
      for (row = 0; row < num_rows; row++)
        output_block[row][col] = col_vector2[row];
    }
  
  return_value = 0;
  goto Return;
 Error:
  return_value = 1;
 Return:
  QccVectorFree(col_vector1);
  QccVectorFree(col_vector2);
  return(return_value);
}


int QccMatrixAddNoiseToRegion(QccMatrix matrix,
                              int num_rows, int num_cols,
                              int start_row, int start_col,
                              int region_num_rows,
                              int region_num_cols,
                              double noise_variance)
{
  int stop_row, stop_col;
  int row, col;
  double val;
  
  if (matrix == NULL)
    return(0);
  
  if ((region_num_rows <= 0) ||
      (region_num_cols <= 0))
    return(0);
  
  if (start_row < 0)
    start_row = 0;
  if (start_col < 0)
    start_col = 0;
  
  stop_row = start_row + region_num_rows - 1;
  stop_col = start_col + region_num_cols - 1;
  
  if (stop_row >= num_rows)
    stop_row = num_rows - 1;
  if (stop_col >= num_cols)
    stop_col = num_cols - 1;
  
  for (row = start_row; row <= stop_row; row++)
    for (col = start_col; col <= stop_col; col++)
      {
        /*
         * The following uses the fact that the variance of
         * a (zero-mean) uniform random variable on the range (-c, c) is
         *     var = (c^2)/3
         */
        
        /*  Uniform(0.0, 1.0), mean = 0.5  */
        val = QccMathRand();
        /*  Uniform(-1.0, 1.0), variance = 1/3, mean = 0.0  */
        val = (val - 0.5)*2.0;
        /*  Uniform(-c, c), Variance = noise_variance  */
        val *= sqrt(3.0*noise_variance);
        
        matrix[row][col] += val;
      }
  
  return(0);
}