/vol/vipdata/irtk/geometry++/include/irtkMatrix.h

/*=========================================================================

  Library   : Image Registration Toolkit (IRTK)
  Module    : $Id: irtkMatrix.h 6 2009-01-23 13:44:31Z DR $
  Copyright : Imperial College, Department of Computing
              Visual Information Processing (VIP), 2008 onwards
  Date      : $Date: 2009-01-23 13:44:31 +0000 (Fri, 23 Jan 2009) $
  Version   : $Revision: 6 $
  Changes   : $Author: DR $

=========================================================================*/

#ifndef _IRTKMATRIX_H

#define _IRTKMATRIX_H

#ifdef USE_VXL
#include <vnl/vnl_matrix.h>
#endif

class irtkMatrix : public irtkObject
{

protected:

  int _rows;

  int _cols;

  double **_matrix;

public:

  irtkMatrix();

  irtkMatrix(int, int);

  irtkMatrix(const irtkMatrix &);

  ~irtkMatrix();

  void Initialize(int, int);

  //
  // Matrix access functions
  //

  int Rows() const;

  int Cols() const;

  void   Put(int, int, double);

  double Get(int, int) const;

  //
  // Operators for matrix access
  //

  double& operator()(int, int);

  double  operator()(int, int) const;

  irtkMatrix  operator()(int, int, int, int);

  void    operator()(irtkMatrix &, int, int);

  //
  // Matrix operators for doubles
  //

  irtkMatrix& operator-=(const double&);

  irtkMatrix& operator+=(const double&);

  irtkMatrix& operator*=(const double&);

  irtkMatrix& operator/=(const double&);

  irtkMatrix  operator- (const double&);

  irtkMatrix  operator+ (const double&);

  irtkMatrix  operator* (const double&);

  irtkMatrix  operator/ (const double&);

  //
  // Matrix operators for matrices
  //

  irtkMatrix& operator =(const irtkMatrix&);

  irtkMatrix& operator-=(const irtkMatrix&);

  irtkMatrix& operator+=(const irtkMatrix&);

  irtkMatrix& operator*=(const irtkMatrix&);

  irtkMatrix  operator- (const irtkMatrix&);

  irtkMatrix  operator+ (const irtkMatrix&);

  irtkMatrix  operator* (const irtkMatrix&);

  irtkMatrix  operator! (void);

  irtkMatrix  operator~ (void);

  int operator==(const irtkMatrix &);

  // Matrix exponential via Pade approximation.
  // See Golub and Van Loan, Matrix Computations, Algorithm 11.3-1.
  friend irtkMatrix expm(irtkMatrix);

  // Matrix logarithm.
  friend irtkMatrix logm(irtkMatrix);

  // Matrix square root.
  friend irtkMatrix sqrtm(irtkMatrix);

  friend irtkMatrix FrechetMean(irtkMatrix *, int, int = 10);

  friend irtkMatrix FrechetMean(irtkMatrix *, double *, int, int = 10);

#ifndef USE_STL
  int operator!=(const irtkMatrix &);
#endif

  //
  // Matrix functions
  //

  double Norm(void) const;

  // The infinity norm is the maximum of the absolute value row sums.
  double InfinityNorm(void) const;

  double Det() const;

  void   SVD(irtkMatrix &, irtkVector &, irtkMatrix &) const;

  void   Ident();

  bool IsIdentity() const;

  void   Invert();

  void   Transpose();

  void   Eigenvalues(irtkMatrix &, irtkVector &);

  void   LeastSquaresFit(const irtkVector &, irtkVector &);

  //
  // Matrix in- and output
  //

  friend ostream& operator<< (ostream&, const irtkMatrix&);

  friend istream& operator>> (istream&, irtkMatrix&);

  void Print();

  void Read (char *);

  void Write(char *);

  void Import (char *, int, int);

#ifdef USE_VXL

  template <class T> void Matrix2Vnl(vnl_matrix<T> *m) const;

  template <class T> void Vnl2Matrix(vnl_matrix<T> *m);

#else

  void Matrix2NR(float **) const;
  void Matrix2NR(double **) const;

  void NR2Matrix(float **);
  void NR2Matrix(double **);

  //
  // Matrix operators for vectors
  //

  irtkVector  operator* (const irtkVector&);

#endif
};

//
// Access operators
//

inline int irtkMatrix::Rows() const
{
  return _rows;
}

inline int irtkMatrix::Cols() const
{
  return _cols;
}

inline void irtkMatrix::Put(int rows, int cols, double matrix)
{
#ifdef NO_BOUNDS
  _matrix[cols][rows] = matrix;
#else
  if ((rows >= 0) && (rows < _rows) && (cols >= 0) && (cols < _cols)) {
    _matrix[cols][rows] = matrix;
  } else {
    cout << "irtkMatrix::Put: parameter out of range\n";
  }
#endif
}

inline double irtkMatrix::Get(int rows, int cols) const
{
#ifdef NO_BOUNDS
  return _matrix[cols][rows];
#else
  if ((rows >= 0) && (rows < _rows) && (cols >= 0) && (cols < _cols)) {
    return _matrix[cols][rows];
  } else {
    cout << "irtkMatrix::Get: parameter out of range\n";
    return 0;
  }
#endif
}

inline double &irtkMatrix::operator()(int rows, int cols)
{
#ifdef NO_BOUNDS
  return _matrix[cols][rows];
#else
  if ((rows >= 0) && (rows < _rows) && (cols >= 0) && (cols < _cols)) {
    return _matrix[cols][rows];
  } else {
    cout << "irtkMatrix::operator(): parameter out of range\n";
    return _matrix[0][0];
  }
#endif
}

inline double irtkMatrix::operator()(int rows, int cols) const
{
#ifdef NO_BOUNDS
  return _matrix[cols][rows];
#else
  if ((rows >= 0) && (rows < _rows) && (cols >= 0) && (cols < _cols)) {
    return _matrix[cols][rows];
  } else {
    cout << "irtkMatrix::operator(): parameter out of range\n";
    return 0;
  }
#endif
}

//
// Matrix operators for doubles
//

inline irtkMatrix& irtkMatrix::operator-=(const double &x)
{
  int i, j;

  for (j = 0; j < _cols; j++) {
    for (i = 0; i < _rows; i++) {
      _matrix[j][i] -= x;
    }
  }
  return *this;
}

inline irtkMatrix& irtkMatrix::operator+=(const double &x)
{
  int i, j;

  for (j = 0; j < _cols; j++) {
    for (i = 0; i < _rows; i++) {
      _matrix[j][i] += x;
    }
  }
  return *this;
}

inline irtkMatrix& irtkMatrix::operator*=(const double &x)
{
  int i, j;

  for (j = 0; j < _cols; j++) {
    for (i = 0; i < _rows; i++) {
      _matrix[j][i] *= x;
    }
  }
  return *this;
}

inline irtkMatrix& irtkMatrix::operator/=(const double &x)
{
  int i, j;

  for (j = 0; j < _cols; j++) {
    for (i = 0; i < _rows; i++) {
      _matrix[j][i] /= x;
    }
  }
  return *this;
}

inline irtkMatrix irtkMatrix::operator-(const double &x)
{
  int i, j;
  irtkMatrix m;

  m = *this;
  for (j = 0; j < _cols; j++) {
    for (i = 0; i < m._rows; i++) {
      m._matrix[j][i] = _matrix[j][i] - x;
    }
  }
  return m;
}

inline irtkMatrix irtkMatrix::operator+(const double &x)
{
  int i, j;
  irtkMatrix m;

  m = *this;
  for (j = 0; j < _cols; j++) {
    for (i = 0; i < m._rows; i++) {
      m._matrix[j][i] = _matrix[j][i] + x;
    }
  }
  return m;
}

inline irtkMatrix irtkMatrix::operator*(const double &x)
{
  int i, j;
  irtkMatrix m;

  m = *this;
  for (j = 0; j < _cols; j++) {
    for (i = 0; i < m._rows; i++) {
      m._matrix[j][i] = _matrix[j][i] * x;
    }
  }
  return m;
}

inline irtkMatrix irtkMatrix::operator/(const double &x)
{
  int i, j;
  irtkMatrix m;

  m = *this;
  for (j = 0; j < _cols; j++) {
    for (i = 0; i < m._rows; i++) {
      m._matrix[j][i] = _matrix[j][i] / x;
    }
  }
  return m;
}

//
// Matrix operators for matrices
//

inline int irtkMatrix::operator==(const irtkMatrix& m)
{
  int i, j;

  if ((m._rows != _rows) || (m._cols != _cols)) return 0;
  for (j = 0; j < _cols; j++) {
    for (i = 0; i < m._rows; i++) {
      if (m._matrix[j][i] != _matrix[j][i]) return 0;
    }
  }
  return 1;
}

#ifndef USE_STL
inline int irtkMatrix::operator!=(const irtkMatrix& m)
{
  return !(*this == m);
}
#endif

inline double irtkMatrix::Norm(void) const
{
  int i, j;
  double norm = 0;

  // The norm of a matrix M is defined as trace(M * M~)
  for (j = 0; j < _cols; j++) {
    for (i = 0; i < _rows; i++) {
      norm += _matrix[j][i]*_matrix[j][i];
    }
  }
  return sqrt(norm);
}

inline double irtkMatrix::InfinityNorm(void) const
{
  int i, j;
  double normInf = -1.0 * DBL_MAX;
  double sum;

  for (i = 0; i < _rows; ++i) {
    sum = 0;
    for (j = 0; j < _cols; ++j) {
      sum += abs(_matrix[j][i]);
    }
    if (sum > normInf)
      normInf = sum;
  }
  return normInf;
}


#endif

Generated on Mon Mar 23 16:55:52 2009 for IRTK by Doxygen 1.5.4