---

/home/slang/XVision2/src/Tools/XVMatrix.h

// *** BEGIN_XVISION2_COPYRIGHT_NOTICE ***
// *** END_XVISION2_COPYRIGHT_NOTICE ***

//---------------------------------------------------------------------
//*
//* XVMatrix.h
//*
//* 02/20/95 Gregory D. Hager
//* 
//* Header file for a XVMatrix package for small systems
//* A simple matrix package for small applications.  Matrices are stored in
//* row-major order using a ragged-array data structure.  Submatrices are
//* also supported as are row and column vectors.
//*
//* Bounds checking and size checking are also included by default.
//* Defining the NO_BOUNDS_CHECK turns off bounds checking.
//* Defining check_size turns off size checking on matrix operations.
//* Check in XVMatrix.cc at the start of the file to see whether
//* NO_BOUNDS_CHECK or CHECKSIZE are defined as you would like them to be.
//*
//* The SUBSCRIPT_START_WITH_1 preprocessor variable, if defined,
//* allows you to index these structures from 1 to n, rather than from
//* 0 to n-1.
//*
//* Matrices should be templates, but for some reason templates aren't
//* working so I'll just hardcode the type for now.
//*
//* This header file contains the definitions for the classes
//* _rowvec, XVRowVector, XVColumnVector, XVMatrix, and XVXVDiagonalMatrix
//*
//---------------------------------------------------------------------

//#define SUBSCRIPT_START_WITH_1

#ifndef _XVMATRIX_H_
#define _XVMATRIX_H_

#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <string>

#include <XVImageScalar.h>

// Next two #defines used to be in XVMatrix.cc, but should be here.
#define NO_BOUNDS_CHECK
//#define CHECKSIZE(m,msg) /* do nothing */

#define FR_EPSFCN 1.0E-16       // FrReal is the type of the matrices
#define FrReal double           // FR_EPSFCN is the minimum available value

extern void _panic(char *);     // For simple exception handling.

typedef struct {                // Data structure for reference counts
  FrReal *data;
  int refCount;
} RefCounter;

#ifndef CHECKSIZE
#define CHECKSIZE(m,mess) \
{if ((m.rowNum != rowNum) || (m.colNum != colNum)) _panic(mess);}
#endif

class _rowvec
{
private:
  FrReal *data;
  int colNum;

public:

  _rowvec(FrReal *datain,int ncols) {data = datain; colNum = ncols;}

  FrReal&
  operator[] (int n) {
    if ((n < 0) || (n >= colNum))
      _panic("XVMatrix second argument out of range");
    return data[n];
  }

  const FrReal&
  operator[] (int n) const {
    if ((n < 0) || (n >= colNum))
      _panic("XVMatrix second argument out of range");
    return data[n];
  }

};

class XVRowVector;
class XVColVector;
class XVMatrix;
class XVDiagonalMatrix;


// The XVMatrix class definition starts here
//---------------------------------------------------------------------
//*
//* XVMatrix.h
//*
//* 02/20/95 Gregory D. Hager
//* 
//* Header file for a XVMatrix package for small systems
//* A simple matrix package for small applications.  Matrices are stored in
//* row-major order using a ragged-array data structure.  Submatrices are
//* also supported as are row and column vectors.
//*
//* Bounds checking and size checking are also included by default.
//* Defining the NO_BOUNDS_CHECK turns off bounds checking.
//* Defining check_size turns off size checking on matrix operations.
//* Check in XVMatrix.cc at the start of the file to see whether
//* NO_BOUNDS_CHECK or CHECKSIZE are defined as you would like them to be.
//*
//* The SUBSCRIPT_START_WITH_1 preprocessor variable, if defined,
//* allows you to index these structures from 1 to n, rather than from
//* 0 to n-1.
//*
//* Matrices should be templates, but for some reason templates aren't
//* working so I'll just hardcode the type for now.
//*
//* This header file contains the definitions for the classes
//* _rowvec, XVRowVector, XVColumnVector, XVMatrix, and XVDiagonalMatrix
//*
//---------------------------------------------------------------------

class XVMatrix
{
  friend class XVRowVector;
  friend class XVColVector;
  friend class XVDiagonalMatrix;
  friend XVMatrix add_column(XVMatrix &x, XVColVector &c);

private:
  static FrReal deriv_from_pts(FrReal y1, FrReal y2, FrReal y3,
                      FrReal dstep, FrReal x);

  static void deriv_n(XVColVector (*fn)(const XVColVector&, const XVColVector&),
      const XVColVector& val, int n, const XVColVector& extra, XVColVector& deriv);

  int init(XVMatrix &m,int startr, int startc, int nrows, int ncols);

  // initialise all fields to 0
  void init_empty(void);

protected:
  RefCounter *refPtr;
  int rowNum;
  int colNum;

  FrReal *data;
  FrReal **rowPtrs;

  int csize;
  int dsize;
  int trsize;

  RefCounter *ref();
  void unref();

  int dataShared();

  XVMatrix(XVMatrix &m, int startr, int startc, int nrows, int ncols);


public:
  XVMatrix();
  XVMatrix(const XVMatrix& m);
  XVMatrix(int rr,int cc);
  XVMatrix(int rr,int cc,FrReal* x);
  virtual ~XVMatrix ();
  int n_of_rows() const;
  int n_of_cols() const;

  int resize(int nrows, int ncols);

  XVMatrix& operator=(FrReal x);
  XVMatrix& operator=(const XVMatrix &m);
  XVMatrix& operator<<(const XVMatrix &m);
  XVMatrix& operator<<(FrReal*);
  XVMatrix& operator+=(FrReal x);
  XVMatrix& operator+=(const XVMatrix &mat);
  XVMatrix& operator-=(FrReal x);
  XVMatrix& operator-=(const XVMatrix &mat);
  XVMatrix& operator*=(FrReal x);
  XVMatrix& operator/=(FrReal x);
  XVMatrix operator*(FrReal x) const;
  XVMatrix operator*(const XVMatrix &mat) const;
  XVColVector operator*(XVColVector &mat) const;
  XVRowVector operator*(XVRowVector &mat) const;
  XVMatrix operator/(FrReal x) const;
  XVMatrix operator+(const XVMatrix &mat) const;
  XVMatrix operator-(const XVMatrix &mat) const;
  XVMatrix operator-() const;
  FrReal SumSquare() const;
  FrReal Sum() const;

#ifndef NO_BOUNDS_CHECK
  _rowvec operator[](int n);
#else
  FrReal* operator[](int n);
#endif

  FrReal* operator[](int n) const;
  XVRowVector row(int i);
  XVRowVector Row(int i);
  XVColVector col(int j);
  XVColVector Column(int j);

#ifdef SUBSCRIPT_START_WITH_1  

  XVMatrix operator()(int sr, int lr, int sc, int lc);
  XVMatrix Rows(int first_row, int last_row);
  XVMatrix Columns(int first_col, int last_col);
        
#else

  XVMatrix operator()(int sr, int lr, int sc, int lc);
  XVMatrix Rows(int first_row, int last_row);
  XVMatrix Columns(int first_col, int last_col);

#endif // SUBSCRIPT_START_WITH_1
 
  XVMatrix i() const;
  XVMatrix t() const;
  XVMatrix ip() const;
  XVMatrix op() const;

  void LUDcmp(int* perm, int& d);
  void LUBksb(int* perm, XVColVector& b);
  void solveByLUD(const XVColVector& b, XVColVector& x);
  XVColVector LUDsolve(const XVColVector& b);
  void LDLtDcmp();
  XVMatrix sqrt();
  void SVDcmp(XVColVector& w, XVMatrix& v);
  void SVBksb(const XVColVector& w,const XVMatrix& v,
           const XVColVector& b, XVColVector& x);
  void solveBySVD(const XVColVector& b, XVColVector& x);
  XVColVector SVDsolve(const XVColVector& b);
  static XVMatrix jacobian(XVColVector (*fn)(const XVColVector&, const XVColVector&),
                int odim, const XVColVector& where, const XVColVector& extra);
  friend ostream &operator << (ostream &s,const XVMatrix &m);
  XVMatrix Map(FrReal (*fn)(FrReal)) const;

  XVMatrix add_column(const XVColVector &c) const;

};



class XVRowVector : public XVMatrix
{
  friend class XVMatrix;

protected:
  XVRowVector(XVMatrix &m, int i);

public:
  XVRowVector();
  XVRowVector(int nn);
  XVRowVector(const XVRowVector &v);
#ifdef __xvimagescalar_h
  XVRowVector(const XVImageScalar<double> &x);
#endif

  void resize(int i);
  FrReal &operator [](int n);
  const FrReal &operator [](int n) const;
  XVRowVector &operator=(const XVRowVector &v);
  XVRowVector &operator<<(const XVRowVector &v);
  XVRowVector &operator=(const XVMatrix &m);
  XVRowVector &operator=(FrReal x);
  XVColVector t() const;
};



class XVColVector : public XVMatrix
{
  friend class XVMatrix;

protected:
  XVColVector (XVMatrix &m, int j);
  XVColVector (XVColVector &m, int startr, int nrows);

public:

  XVColVector();
  XVColVector(int nn);
  XVColVector (const XVColVector &v);
  ~XVColVector();

  void resize(int i);
  FrReal &operator [](int n);
  const FrReal &operator [](int n) const;
  XVColVector &operator=(const XVColVector &v);
  XVColVector &operator<<(const XVColVector &v);
  XVColVector &operator<<(FrReal *);
  XVColVector &operator=(const XVMatrix &m);
  XVColVector &operator=(FrReal x);

  // J.W. added the below Aug 15th
  
  XVColVector operator+(const XVColVector &mat) const;
  XVColVector operator-(const XVColVector &mat) const;
  XVColVector operator*(FrReal x) const;
  XVColVector operator-();

  XVColVector Rows(int first_row, int last_row);
  XVRowVector t() const;
  FrReal ip() const;
  FrReal ip(const XVColVector &) const;
};

extern void SVD(XVMatrix& , XVMatrix&, XVMatrix&, XVMatrix& );

//#ifdef __xvimagescalar_h
template <class T> XVColVector  &operator <<(XVColVector &x,const XVImageScalar<T> &y);
template <class T> XVRowVector  &operator <<(XVRowVector &x,const XVImageScalar<T> &y);
template <class T> XVImageScalar<T> &operator >>(const XVColVector &x,const XVImageScalar<T> &y);
template <class T> XVImageScalar<T> &operator >>(const XVRowVector &x,const XVImageScalar<T> &y);

//#endif





class XVDiagonalMatrix
{
private:  
  XVColVector t;

public:
  XVDiagonalMatrix(int n= 0);
  XVDiagonalMatrix(const XVColVector &tin);

  int n_of_cols() const;
  int n_of_rows() const;
  XVColVector& diagonal();
  void resize(int n);

  FrReal &operator() (int n);
  FrReal operator() (int n) const;
  XVDiagonalMatrix &operator = (const XVDiagonalMatrix &m);
  XVDiagonalMatrix &operator = (FrReal x);

  friend XVMatrix operator*(const XVDiagonalMatrix &, const XVMatrix &);
  friend XVRowVector operator*(const XVDiagonalMatrix &, const XVRowVector &);
  friend XVMatrix operator*(const XVMatrix&, const XVDiagonalMatrix &);
  friend XVDiagonalMatrix operator*(const XVDiagonalMatrix &x,const XVDiagonalMatrix &y);
  friend ostream &operator << (ostream &,const XVDiagonalMatrix &);
};


//----------------------------------------------------------------
// utilities
//----------------------------------------------------------------

FrReal invs(FrReal x);

template <class T>
XVMatrix add_column(XVMatrix &x, XVImageScalar<T> &I);


#include "XVMatrix.icc"

#endif

---