nagmatrix.h

VC视频对象的跟踪提取原代码（vc视频监控源码）

// 
// NagMatrix.h
// 
// A class for manipulating (2D) realno matrices
// (uses NAG or BLAS/LAPACK routines)

#ifndef __NAG_MATRIX_H__
#define __NAG_MATRIX_H__

#include <cstdlib>
#include <cstdio>
#include <cassert>

#include "NagVector.h"

namespace ReadingPeopleTracker
{

#ifdef DEBUG
#define MATRIX_DEBUG
#endif

class NagMatrix 
{
    friend ostream &operator<<(ostream&, const NagMatrix&);
    friend istream &operator>>(istream&, NagMatrix&);
//    friend class NagVector;
    
protected:
    unsigned int rows;
    unsigned int columns;
    realno *data;		// array of values
    bool own_memory;	        // if true then delete data when finished
    
public:

    // constructors, destructor
    
    NagMatrix(unsigned int n, unsigned int m, realno *the_data) 
	{
	    rows = n;
	    columns = m;
	    data = the_data;
	    own_memory = false;
#ifdef USE_LAPACK
// 	    best_work = 0;  // optimal workspace size yet unknown
#else

#endif
	}
    
    NagMatrix(unsigned int n, unsigned int m);
    NagMatrix(unsigned int n, unsigned int m, realno initial_value);
    NagMatrix(unsigned int n, unsigned int m, realno (*func) (unsigned int, unsigned int));
    NagMatrix() 
	{
	    rows = columns = 0;
	    data = NULL;
	    own_memory = false;
	}
    NagMatrix(const NagMatrix &m)
	{
	    // first indicate that we have no data yet
	    rows = columns = 0;
	    data = NULL;
	    own_memory = false;

	    // then use our operator=
	    (*this) = m;
	}

    NagMatrix(unsigned int n, unsigned int m, const NagVector &v);
    
    ~NagMatrix();
    
    //data access
    inline unsigned int no_rows() const { return rows; }
    inline unsigned int no_columns() const { return columns; }
    inline bool has_own_memory() const { return own_memory; }

    // careful with this one!
    inline realno *get_data() { return data; }
    const realno *get_data_const() const { return data; }    
    
    inline void set (unsigned int i, unsigned int j, realno value)
	{ 
	    assert (i < rows);
	    assert (j < columns);
	    
	    data[i + j * rows] = value;
	}
    inline realno read (unsigned int i, unsigned int j) const
	{
	    assert (i < rows);
	    assert (j < columns);
	    
	    return data[i + j * rows];
	}
    inline realno *get (unsigned int i, unsigned int j)
	{ 
	    assert (i < rows);
	    assert (j < columns);
	    
	    return &(data[i + j * rows]);
	}
    
    inline realno &operator[] (unsigned int i)
	{
	    return data[i];
	}
    
    void operator=(const NagMatrix &m2);
    void reconstruct(const unsigned int rows, const unsigned int columns);
    void reconstruct(const unsigned int rows, const unsigned int columns, const NagVector &v);
    void reset()
	{
	    rows = columns  = 0;
	    data = NULL;
	    own_memory = false;
	}
    void destruct()
	{
	    if ((data != NULL) && (own_memory))
		delete [] data;
	    reset();
	}
    
#ifdef _SVID_SOURCE
    // the following uses SVID 48-bit random numbers...
    
    // set each element randomly
    void set_unif_random(realno min = -1, realno max = 1.0)
	{
	    NagVector tmp1(rows * columns, data);
	    tmp1.set_unif_random(min, max);
	    tmp1.reset();
	}
    
    void set_gauss_random(realno mean = 0, realno sd = 1.0)
	{
	    NagVector tmp1(rows * columns, data);
	    tmp1.set_gauss_random(mean, sd);
	    tmp1.reset();
	}
#endif // ifdef _SVID_SOURCE
    
    // mathematical functions
    
    void transpose(NagMatrix &r) const;
    void flip_horizontally();
    void invert(NagMatrix &r) const;
    void pseudo_invert(NagMatrix &r) const;
    
// get_eigenvectors(result, evals))
//  calculates eigenvalues, eigenvectors of real symmetric matrix
//
//  on entry:
//    (*this) real sysmetric matrix
//  on exit:
//     result = (v_1 ... v_n) eigenvectors
//      evals = (e_1 ... e_n) eigenvalues in ascending order
    
    
    void get_eigenvectors(NagMatrix &result, NagVector &evals) const;
    
// get_eigensystem
// returns complex eigenvalues, eigenvectors 
// of a real non-symmetric matrix
    
    void get_eigensystem(NagMatrix &res_r, NagMatrix &res_i,
			 NagVector &evals_r, NagVector &evals_i) const;
    
// QR factorise a real matrix
    
    void QR_factorise(NagMatrix &Q) const;
    
// msqrt
// returns the positive definite square root
// of a real symmetric matrix
    
    void msqrt(NagMatrix &res);
    
// sv_decompose
// Singular Value Decomposition
// A =  Q D P'
// where Q , P are orthogonal matrices
    
    void sv_decompose(NagMatrix &Q, NagVector &sv, NagMatrix &P_t) const;
    
    void scale(const realno s, NagMatrix &r) const;
    
    void add(const NagMatrix &m, NagMatrix &r) const;
    void subtract(const NagMatrix &m, NagMatrix &r) const;
    void multiply(const NagMatrix &m, NagMatrix &r) const;
    void multiply(const NagVector &m, NagVector &r) const;
    void multiply(NagVector &r) const
	{
	    NagVector tmp;
	    r.copy(tmp);
	    multiply(tmp, r);
	}
    void multiply(NagMatrix &m) const
	{
	    NagMatrix tmp;
	    tmp = m;
	    multiply(tmp, m);
	}
    void solve_equations(NagMatrix B, NagMatrix &res) const;
    
    // result = alpha * op(A) * op(B)
    // where op(A) = A   when transa = 'N'
    // and   op(A) = A^T when transa = 'T'
    void multiply(const NagMatrix &m2, NagMatrix &result, 
		  char transa, char transb, realno alpha = 1.0) const;
    void multiply(const NagVector &m, NagVector &r, char trans,
		  realno alpha = 1.0);
    // misc

    void clear(realno s = 0);
    
    void output() const;
    
    // column -- returns a NagVector referencing the i'th column
    NagVector column(unsigned int i) {NagVector res(rows, get(0,i)); return res; }
    
    // get_column -- copy contents of column into result
    void get_column(unsigned int col, NagVector &r) const;
    void get_row(unsigned int row, NagVector &r) const;
    void set_column(unsigned int col, const NagVector &column);
    void set_row(unsigned int row, const NagVector &row_v);
    
    void set_block(unsigned int row_off, unsigned int col_off, const NagMatrix &block_m);
    void get_block(unsigned int row_top, unsigned int row_bottom,
		   unsigned int col_left, unsigned int col_right,
		   NagMatrix &block) const;
    
    
    // pack a symmetric  into an ordered vector
    // (lower triangle used)
    void pack_to_triangle(NagVector &res) const;
    // unpack a vector into a symmetric matrix
    void unpack_triangle(const NagVector &data, bool symmetric = true);
    
    void copy(NagMatrix &res) const;
    
    // ortho transform
    // calculates Qt A Q
    // (where Qt is Q transform)
    void ortho_transform(NagMatrix &Q, NagMatrix &res, bool flag = true) const;
    void square_entries(NagMatrix &r) const;
    void scale_diagonal(realno fac = 0.5);
    void read_diagonal(NagVector &diag);
    void set_diagonal(NagVector &diag);
    
    realno trace();
    
    // X_full = *this + (*this)^T
    void lower_to_full(NagMatrix &X_full);
    
public:
    // public helpers
    static realno identity_fn(unsigned int i, unsigned int j);
    
private:    
    // private helpers
    
    // calculate sum (i,j) { A_ij B_ji }
    realno trace_A_B(NagMatrix &A, NagMatrix &B);
    // calculate sum (i,j) { A_ij B_ij }
    realno trace_A_Bt(NagMatrix &A, NagMatrix &B);
    void matrix_error(const char *message) const;
    
#ifdef USE_LAPACK
//     int best_work;   // optimal workspace size
#else
    // default value for IFAIL variable, determines what NAG should do on failure
    static const int DEF_IFAIL; // show error message but continue
#endif
    
};

ostream &operator<< (ostream&out, const NagMatrix&m);
istream &operator>> (istream&in, NagMatrix&m);  

inline NagMatrix operator*(const NagMatrix &m1, const NagMatrix &m2)
{
    NagMatrix res;
    m1.multiply(m2, res);
    return res;
}

inline NagMatrix operator+(const NagMatrix &m1, const NagMatrix &m2)
{
    NagMatrix res;
    m1.add(m2, res);
    return res;
}

inline NagMatrix operator-(const NagMatrix &m1, const NagMatrix &m2)
{
    NagMatrix res;
    m1.subtract(m2, res);
    return res;
}

inline NagVector operator*(const NagMatrix &m, const NagVector &v)
{
    NagVector res;
    m.multiply(v, res);
    return res;
}

inline NagMatrix operator*(const double k, const NagMatrix &m)
{
    NagMatrix res;
    m.scale(k, res);
    return res;
}

inline NagMatrix operator*(const NagMatrix &m, const double k)
{
    NagMatrix res;
    m.scale(k, res);
    return res;
}

inline NagMatrix operator/(const NagMatrix &m, const double k)
{
    NagMatrix res;
    m.scale(1.0 / k, res);
    return res;
}

} // namespace ReadingPeopleTracker

#endif