Initial commit
This commit is contained in:
parent
94895e32a6
commit
546d4d88e1
vcg/math
|
|
@ -0,0 +1,178 @@
|
||||||
|
#include <vcg/math/matrix44.h>
|
||||||
|
|
||||||
|
namespace vcg
|
||||||
|
{
|
||||||
|
/** \addtogroup math */
|
||||||
|
/* @{ */
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Computes all eigenvalues and eigenvectors of a real symmetric matrix .
|
||||||
|
* On output, elements of the input matrix above the diagonal are destroyed.
|
||||||
|
* \param d returns the eigenvalues of a.
|
||||||
|
* \param v is a matrix whose columns contain, the normalized eigenvectors
|
||||||
|
* \param nrot returns the number of Jacobi rotations that were required.
|
||||||
|
*/
|
||||||
|
template <typename TYPE>
|
||||||
|
static void Jacobi(Matrix44<TYPE> &w, Point4<TYPE> &d, Matrix44<TYPE> &v, int &nrot)
|
||||||
|
{
|
||||||
|
int j,iq,ip,i;
|
||||||
|
//assert(w.IsSymmetric());
|
||||||
|
TYPE tresh, theta, tau, t, sm, s, h, g, c;
|
||||||
|
Point4<TYPE> b, z;
|
||||||
|
|
||||||
|
v.SetIdentity();
|
||||||
|
|
||||||
|
for (ip=0;ip<4;++ip) //Initialize b and d to the diagonal of a.
|
||||||
|
{
|
||||||
|
b[ip]=d[ip]=w[ip][ip];
|
||||||
|
z[ip]=0.0; //This vector will accumulate terms of the form tapq as in equation (11.1.14).
|
||||||
|
}
|
||||||
|
nrot=0;
|
||||||
|
for (i=0;i<50;i++)
|
||||||
|
{
|
||||||
|
sm=0.0;
|
||||||
|
for (ip=0;ip<3;++ip) // Sum off diagonal elements
|
||||||
|
{
|
||||||
|
for (iq=ip+1;iq<4;++iq)
|
||||||
|
sm += fabs(w[ip][iq]);
|
||||||
|
}
|
||||||
|
if (sm == 0.0) //The normal return, which relies on quadratic convergence to machine underflow.
|
||||||
|
{
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
if (i < 4)
|
||||||
|
tresh=0.2*sm/(4*4); //...on the first three sweeps.
|
||||||
|
else
|
||||||
|
tresh=0.0; //...thereafter.
|
||||||
|
for (ip=0;ip<4-1;++ip)
|
||||||
|
{
|
||||||
|
for (iq=ip+1;iq<4;iq++)
|
||||||
|
{
|
||||||
|
g=100.0*fabs(w[ip][iq]);
|
||||||
|
//After four sweeps, skip the rotation if the off-diagonal element is small.
|
||||||
|
if(i>4 && (float)(fabs(d[ip])+g) == (float)fabs(d[ip]) && (float)(fabs(d[iq])+g) == (float)fabs(d[iq]))
|
||||||
|
w[ip][iq]=0.0;
|
||||||
|
else if (fabs(w[ip][iq]) > tresh)
|
||||||
|
{
|
||||||
|
h=d[iq]-d[ip];
|
||||||
|
if ((float)(fabs(h)+g) == (float)fabs(h))
|
||||||
|
t=(w[ip][iq])/h; //t =1/(2#)
|
||||||
|
else
|
||||||
|
{
|
||||||
|
theta=0.5*h/(w[ip][iq]); //Equation (11.1.10).
|
||||||
|
t=1.0/(fabs(theta)+sqrt(1.0+theta*theta));
|
||||||
|
if (theta < 0.0) t = -t;
|
||||||
|
}
|
||||||
|
c=1.0/sqrt(1+t*t);
|
||||||
|
s=t*c;
|
||||||
|
tau=s/(1.0+c);
|
||||||
|
h=t*w[ip][iq];
|
||||||
|
z[ip] -= h;
|
||||||
|
z[iq] += h;
|
||||||
|
d[ip] -= h;
|
||||||
|
d[iq] += h;
|
||||||
|
w[ip][iq]=0.0;
|
||||||
|
for (j=0;j<=ip-1;j++) { //Case of rotations 1 <= j < p.
|
||||||
|
JacobiRotate<TYPE>(w,s,tau,j,ip,j,iq) ;
|
||||||
|
}
|
||||||
|
for (j=ip+1;j<=iq-1;j++) { //Case of rotations p < j < q.
|
||||||
|
JacobiRotate<TYPE>(w,s,tau,ip,j,j,iq);
|
||||||
|
}
|
||||||
|
for (j=iq+1;j<4;j++) { //Case of rotations q< j <= n.
|
||||||
|
JacobiRotate<TYPE>(w,s,tau,ip,j,iq,j);
|
||||||
|
}
|
||||||
|
for (j=0;j<4;j++) {
|
||||||
|
JacobiRotate<TYPE>(w,s,tau,j,ip,j,iq);
|
||||||
|
}
|
||||||
|
++nrot;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
for (ip=0;ip<4;ip++)
|
||||||
|
{
|
||||||
|
b[ip] += z[ip];
|
||||||
|
d[ip]=b[ip]; //Update d with the sum of ta_pq ,
|
||||||
|
z[ip]=0.0; //and reinitialize z.
|
||||||
|
}
|
||||||
|
}
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
*
|
||||||
|
*/
|
||||||
|
template< typename TYPE >
|
||||||
|
void JacobiRotate(Matrix44<TYPE> &A, TYPE s, TYPE tau, int i,int j,int k,int l)
|
||||||
|
{
|
||||||
|
TYPE g=A[i][j];
|
||||||
|
TYPE h=A[k][l];
|
||||||
|
A[i][j]=g-s*(h+g*tau);
|
||||||
|
A[k][l]=h+s*(g-h*tau);
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Given a matrix <I>A<SUB>m×n</SUB></I>, this routine computes its singular value decomposition,
|
||||||
|
* i.e. <I>A=U·W·V<SUP>T</SUP></I>. The matrix <I>A</I> will be destroyed!
|
||||||
|
* \param A ...
|
||||||
|
* \param W the diagonal matrix of singular values <I>W</I>, stored as a vector <I>W[1...N]</I>
|
||||||
|
* \param V the matrix <I>V</I> (not the transpose <I>V<SUP>T</SUP></I>)
|
||||||
|
*/
|
||||||
|
template <typename MATRIX_TYPE>
|
||||||
|
static void SingularValueDecomposition(MATRIX_TYPE &A, typename MATRIX_TYPE::ScalarType *W, MATRIX_TYPE &V)
|
||||||
|
{
|
||||||
|
typedef typename MATRIX_TYPE::ScalarType ScalarType;
|
||||||
|
|
||||||
|
};
|
||||||
|
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Solves A·X = B for a vector X, where A is specified by the matrices <I>U<SUB>m×n</SUB></I>,
|
||||||
|
* <I>W<SUB>n×1</SUB></I> and <I>V<SUB>n×n</SUB></I> as returned by <CODE>SingularValueDecomposition</CODE>.
|
||||||
|
* No input quantities are destroyed, so the routine may be called sequentially with different b’s.
|
||||||
|
* \param x is the output solution vector (<I>x<SUB>n×1</SUB></I>)
|
||||||
|
* \param b is the input right-hand side (<I>b<SUB>n×1</SUB></I>)
|
||||||
|
*/
|
||||||
|
template <typename MATRIX_TYPE>
|
||||||
|
static void SingularValueBacksubstitution(const MATRIX_TYPE &U,
|
||||||
|
const typename MATRIX_TYPE::ScalarType *W,
|
||||||
|
const MATRIX_TYPE &V,
|
||||||
|
typename MATRIX_TYPE::ScalarType *x,
|
||||||
|
const typename MATRIX_TYPE::ScalarType *b)
|
||||||
|
{
|
||||||
|
unsigned int jj, j, i;
|
||||||
|
ScalarType s;
|
||||||
|
ScalarType tmp = new ScalarType[U._columns];
|
||||||
|
for (j=0; j<U._columns; j++) //Calculate U^T * B.
|
||||||
|
{
|
||||||
|
s = 0;
|
||||||
|
if (W[j]!=0) //Nonzero result only if wj is nonzero.
|
||||||
|
{
|
||||||
|
for (i=0; i<U._rows; i++)
|
||||||
|
s += U[i][j]*b[i];
|
||||||
|
s /= w[j]; //This is the divide by wj .
|
||||||
|
}
|
||||||
|
tmp[j]=s;
|
||||||
|
}
|
||||||
|
for (j=0;j<U._columns;j++) //Matrix multiply by V to get answer.
|
||||||
|
{
|
||||||
|
s = 0;
|
||||||
|
for (jj=0; jj<U._columns; jj++)
|
||||||
|
s += V[j][jj]*tmp[jj];
|
||||||
|
x[j]=s;
|
||||||
|
}
|
||||||
|
delete []tmp;
|
||||||
|
};
|
||||||
|
|
||||||
|
// Computes (a^2 + b^2)^(1/2) without destructive underflow or overflow.
|
||||||
|
template <typename TYPE>
|
||||||
|
inline TYPE pythagora(TYPE a, TYPE b)
|
||||||
|
{
|
||||||
|
TYPE abs_a = fabs(a);
|
||||||
|
TYPE abs_b = fabs(b);
|
||||||
|
if (abs_a > abs_b)
|
||||||
|
return abs_a*sqrt(1.0+sqr(abs_b/abs_a));
|
||||||
|
else
|
||||||
|
return (abs_b == 0.0 ? 0.0 : abs_b*sqrt(1.0+sqr(abs_a/abs_b)));
|
||||||
|
};
|
||||||
|
|
||||||
|
/*! @} */
|
||||||
|
}; // end of namespace
|
||||||
|
|
@ -0,0 +1,557 @@
|
||||||
|
/****************************************************************************
|
||||||
|
* VCGLib o o *
|
||||||
|
* Visual and Computer Graphics Library o o *
|
||||||
|
* _ O _ *
|
||||||
|
* Copyright(C) 2004 \/)\/ *
|
||||||
|
* Visual Computing Lab /\/| *
|
||||||
|
* ISTI - Italian National Research Council | *
|
||||||
|
* \ *
|
||||||
|
* All rights reserved. *
|
||||||
|
* *
|
||||||
|
* This program is free software; you can redistribute it and/or modify *
|
||||||
|
* it under the terms of the GNU General Public License as published by *
|
||||||
|
* the Free Software Foundation; either version 2 of the License, or *
|
||||||
|
* (at your option) any later version. *
|
||||||
|
* *
|
||||||
|
* This program 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 General Public License (http://www.gnu.org/licenses/gpl.txt) *
|
||||||
|
* for more details. *
|
||||||
|
* *
|
||||||
|
****************************************************************************/
|
||||||
|
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <math.h>
|
||||||
|
#include <memory.h>
|
||||||
|
#include <assert.h>
|
||||||
|
#include <algorithm>
|
||||||
|
|
||||||
|
|
||||||
|
namespace vcg
|
||||||
|
{
|
||||||
|
/** \addtogroup math */
|
||||||
|
/* @{ */
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* This class represent a generic <I>m</I>×<I>n</I> matrix. The class is templated over the scalar type field.
|
||||||
|
* @param TYPE (Templete Parameter) Specifies the ScalarType field.
|
||||||
|
*/
|
||||||
|
template<class TYPE>
|
||||||
|
class Matrix
|
||||||
|
{
|
||||||
|
|
||||||
|
public:
|
||||||
|
typedef TYPE ScalarType;
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Default constructor
|
||||||
|
* All the elements are initialized to zero.
|
||||||
|
* \param m the number of matrix rows
|
||||||
|
* \param n the number of matrix columns
|
||||||
|
*/
|
||||||
|
Matrix(unsigned int m, unsigned int n)
|
||||||
|
{
|
||||||
|
_rows = m;
|
||||||
|
_columns = n;
|
||||||
|
_data = new ScalarType[m*n];
|
||||||
|
memset(_data, ScalarType(0.0), m*n*sizeof(ScalarType));
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Constructor
|
||||||
|
* The matrix elements are initialized with the values of the elements in \i values.
|
||||||
|
* \param m the number of matrix rows
|
||||||
|
* \param n the number of matrix columns
|
||||||
|
* \param values the values of the matrix elements
|
||||||
|
*/
|
||||||
|
Matrix(unsigned int m, unsigned int n, TYPE *values)
|
||||||
|
{
|
||||||
|
_rows = m;
|
||||||
|
_columns = n;
|
||||||
|
_data = new ScalarType[m*n];
|
||||||
|
unsigned int i;
|
||||||
|
for (i=0; i<_rows*_columns; i++)
|
||||||
|
_data[i] = values[i];
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Copy constructor
|
||||||
|
* The matrix elements are initialized with the value of the corresponding element in \i m
|
||||||
|
* \param m the matrix to be copied
|
||||||
|
*/
|
||||||
|
Matrix(const Matrix<TYPE> &m)
|
||||||
|
{
|
||||||
|
_rows = m._rows;
|
||||||
|
_columns = m._columns;
|
||||||
|
_data = new ScalarType[_rows*_columns];
|
||||||
|
for (unsigned int i=0; i<_rows*_columns; i++)
|
||||||
|
_data[i] = m._data[i];
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Default destructor
|
||||||
|
*/
|
||||||
|
~Matrix()
|
||||||
|
{
|
||||||
|
delete []_data;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Number of columns
|
||||||
|
*/
|
||||||
|
inline unsigned int ColumnsNumber()
|
||||||
|
{
|
||||||
|
return _columns;
|
||||||
|
};
|
||||||
|
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Number of rows
|
||||||
|
*/
|
||||||
|
inline unsigned int RowsNumber()
|
||||||
|
{
|
||||||
|
return _rows;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Equality operator.
|
||||||
|
* \param m
|
||||||
|
* \return true iff the matrices have same size and its elements have same values.
|
||||||
|
*/
|
||||||
|
bool operator==(const Matrix<TYPE> &m) const
|
||||||
|
{
|
||||||
|
if (_rows==m._rows && _columns==m._columns)
|
||||||
|
{
|
||||||
|
bool result = true;
|
||||||
|
for (unsigned int i=0; i<_rows*_columns && result; i++)
|
||||||
|
result = (_data[i]==m._data[i]);
|
||||||
|
return result;
|
||||||
|
}
|
||||||
|
return false;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Inequality operator
|
||||||
|
* \param m
|
||||||
|
* \return true iff the matrices have different size or if their elements have different values.
|
||||||
|
*/
|
||||||
|
bool operator!=(const Matrix<TYPE> &m) const
|
||||||
|
{
|
||||||
|
if (_rows==m._rows && _columns==m._columns)
|
||||||
|
{
|
||||||
|
bool result = false;
|
||||||
|
for (unsigned int i=0; i<_rows*_columns && !result; i++)
|
||||||
|
result = (_data[i]!=m._data[i]);
|
||||||
|
return result;
|
||||||
|
}
|
||||||
|
return true;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Return the element stored in the <I>i</I>-th rows at the <I>j</I>-th column
|
||||||
|
* \param i the row index
|
||||||
|
* \param j the column index
|
||||||
|
* \return the element
|
||||||
|
*/
|
||||||
|
inline TYPE ElementAt(unsigned int i, unsigned int j)
|
||||||
|
{
|
||||||
|
assert(i>=0 && i<_rows);
|
||||||
|
assert(j>=0 && j<_columns);
|
||||||
|
return _data[i*_columns+j];
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Calculate and return the matrix determinant (Laplace)
|
||||||
|
* \return the matrix determinant
|
||||||
|
*/
|
||||||
|
TYPE Determinant() const
|
||||||
|
{
|
||||||
|
assert(_rows == _columns);
|
||||||
|
switch (_rows)
|
||||||
|
{
|
||||||
|
case 2:
|
||||||
|
{
|
||||||
|
return _data[0]*_data[3]-_data[1]*_data[2];
|
||||||
|
break;
|
||||||
|
};
|
||||||
|
case 3:
|
||||||
|
{
|
||||||
|
return _data[0]*(_data[4]*_data[8]-_data[5]*_data[7]) -
|
||||||
|
_data[1]*(_data[3]*_data[8]-_data[5]*_data[6]) +
|
||||||
|
_data[2]*(_data[3]*_data[7]-_data[4]*_data[6]) ;
|
||||||
|
break;
|
||||||
|
};
|
||||||
|
default:
|
||||||
|
{
|
||||||
|
// da migliorare: si puo' cercare la riga/colonna con maggior numero di zeri
|
||||||
|
ScalarType det = 0;
|
||||||
|
for (unsigned int j=0; j<_columns; j++)
|
||||||
|
if (_data[j]!=0)
|
||||||
|
det += _data[j]*this->Cofactor(0, j);
|
||||||
|
|
||||||
|
return det;
|
||||||
|
}
|
||||||
|
};
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Return the cofactor <I>A<SUB>i,j</SUB></I> of the <I>a<SUB>i,j</SUB></I> element
|
||||||
|
* \return ...
|
||||||
|
*/
|
||||||
|
TYPE Cofactor(unsigned int i, unsigned int j) const
|
||||||
|
{
|
||||||
|
assert(_rows == _columns);
|
||||||
|
assert(_rows>2);
|
||||||
|
TYPE *values = new TYPE[(_rows-1)*(_columns-1)];
|
||||||
|
unsigned int u, v, p, q, s, t;
|
||||||
|
for (u=0, p=0, s=0, t=0; u<_rows; u++, t+=_rows)
|
||||||
|
{
|
||||||
|
if (i==u)
|
||||||
|
continue;
|
||||||
|
|
||||||
|
for (v=0, q=0; v<_columns; v++)
|
||||||
|
{
|
||||||
|
if (j==v)
|
||||||
|
continue;
|
||||||
|
values[s+q] = _data[t+v];
|
||||||
|
q++;
|
||||||
|
}
|
||||||
|
p++;
|
||||||
|
s+=(_rows-1);
|
||||||
|
}
|
||||||
|
Matrix<TYPE> temp(_rows-1, _columns-1, values);
|
||||||
|
return (pow(-1, i+j)*temp.Determinant());
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Subscript operator:
|
||||||
|
* \param i the index of the row
|
||||||
|
* \return a reference to the <I>i</I>-th matrix row
|
||||||
|
*/
|
||||||
|
inline TYPE* operator[](const unsigned int i)
|
||||||
|
{
|
||||||
|
assert(i>=0 && i<_columns);
|
||||||
|
return _data + i*_columns;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Const subscript operator
|
||||||
|
* \param i the index of the row
|
||||||
|
* \return a reference to the <I>i</I>-th matrix row
|
||||||
|
*/
|
||||||
|
inline const TYPE* operator[](const unsigned int i) const
|
||||||
|
{
|
||||||
|
assert(i>=0 && i<_columns);
|
||||||
|
return _data + i*_columns;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Get the <I>j</I>-th column on the matrix.
|
||||||
|
* \param j the column index.
|
||||||
|
* \return the reference to the column elements.
|
||||||
|
*/
|
||||||
|
TYPE* GetColumn(const unsigned int j)
|
||||||
|
{
|
||||||
|
assert(j>=0 && j<_columns);
|
||||||
|
ScalarType *v = new ScalarType[_columns];
|
||||||
|
unsigned int i, p;
|
||||||
|
for (i=0, p=j; i<_rows; i++, p+=_columns)
|
||||||
|
v[i] = _data[p];
|
||||||
|
return v;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Get the <I>i</I>-th row on the matrix.
|
||||||
|
* \param i the column index.
|
||||||
|
* \return the reference to the row elements.
|
||||||
|
*/
|
||||||
|
TYPE* GetRow(const unsigned int i)
|
||||||
|
{
|
||||||
|
assert(i>=0 && i<_rows);
|
||||||
|
ScalarType *v = new ScalarType[_rows];
|
||||||
|
unsigned int j, p;
|
||||||
|
for (j=0, p=i*_columns; j<_columns; j++, p++)
|
||||||
|
v[j] = _data[p];
|
||||||
|
return v;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Assignment operator
|
||||||
|
* \param m ...
|
||||||
|
*/
|
||||||
|
Matrix<TYPE>& operator=(const Matrix<TYPE> &m)
|
||||||
|
{
|
||||||
|
if (this != &m)
|
||||||
|
{
|
||||||
|
assert(_rows == m._rows);
|
||||||
|
assert(_columns == m._columns);
|
||||||
|
for (unsigned int i=0; i<_rows*_columns; i++)
|
||||||
|
_data[i] = m._data[i];
|
||||||
|
}
|
||||||
|
return *this;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Adds a matrix <I>m</I> to this matrix.
|
||||||
|
* \param m reference to matrix to add to this
|
||||||
|
* \return the matrix sum.
|
||||||
|
*/
|
||||||
|
Matrix<TYPE>& operator+=(const Matrix<TYPE> &m)
|
||||||
|
{
|
||||||
|
assert(_rows == m._rows);
|
||||||
|
assert(_columns == m._columns);
|
||||||
|
for (unsigned int i=0; i<_rows*_columns; i++)
|
||||||
|
_data[i] += m._data[i];
|
||||||
|
return *this;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Subtracts a matrix <I>m</I> to this matrix.
|
||||||
|
* \param m reference to matrix to subtract
|
||||||
|
* \return the matrix difference.
|
||||||
|
*/
|
||||||
|
Matrix<TYPE>& operator-=(const Matrix<TYPE> &m)
|
||||||
|
{
|
||||||
|
assert(_rows == m._rows);
|
||||||
|
assert(_columns == m._columns);
|
||||||
|
for (unsigned int i=0; i<_rows*_columns; i++)
|
||||||
|
_data[i] -= m._data[i];
|
||||||
|
return *this;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* (Modifier) Add to each element of this matrix the scalar constant <I>k</I>.
|
||||||
|
* \param k the scalar constant
|
||||||
|
* \return the modified matrix
|
||||||
|
*/
|
||||||
|
Matrix<TYPE>& operator+=(const TYPE k)
|
||||||
|
{
|
||||||
|
for (unsigned int i=0; i<_rows*_columns; i++)
|
||||||
|
_data[i] += k;
|
||||||
|
return *this;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* (Modifier) Subtract from each element of this matrix the scalar constant <I>k</I>.
|
||||||
|
* \param k the scalar constant
|
||||||
|
* \return the modified matrix
|
||||||
|
*/
|
||||||
|
Matrix<TYPE>& operator-=(const TYPE k)
|
||||||
|
{
|
||||||
|
for (unsigned int i=0; i<_rows*_columns; i++)
|
||||||
|
_data[i] -= k;
|
||||||
|
return *this;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* (Modifier) Multiplies each element of this matrix by the scalar constant <I>k</I>.
|
||||||
|
* \param k the scalar constant
|
||||||
|
* \return the modified matrix
|
||||||
|
*/
|
||||||
|
Matrix<TYPE>& operator*=(const TYPE k)
|
||||||
|
{
|
||||||
|
for (unsigned int i=0; i<_rows*_columns; i++)
|
||||||
|
_data[i] *= k;
|
||||||
|
return *this;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* (Modifier) Divides each element of this matrix by the scalar constant <I>k</I>.
|
||||||
|
* \param k the scalar constant
|
||||||
|
* \return the modified matrix
|
||||||
|
*/
|
||||||
|
Matrix<TYPE>& operator/=(const TYPE k)
|
||||||
|
{
|
||||||
|
assert(k!=0);
|
||||||
|
for (unsigned int i=0; i<_rows*_columns; i++)
|
||||||
|
_data[i] /= k;
|
||||||
|
return *this;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Matrix multiplication: calculates the cross product.
|
||||||
|
* \param reference to the matrix to multiply by
|
||||||
|
* \result the matrix product
|
||||||
|
*/
|
||||||
|
Matrix<TYPE> operator*(const Matrix<TYPE> &m)
|
||||||
|
{
|
||||||
|
assert(_columns == m._rows);
|
||||||
|
Matrix<TYPE> result(_rows, m._columns);
|
||||||
|
unsigned int i, j, k, p, q, r;
|
||||||
|
for (i=0, p=0, r=0; i<result._rows; i++, p+=_columns, r+=result._columns)
|
||||||
|
for (j=0; j<result._columns; j++)
|
||||||
|
{
|
||||||
|
ScalarType temp = 0;
|
||||||
|
for (k=0, q=0; k<_columns; k++, q+=m._columns)
|
||||||
|
temp+=(_data[p+k]*m._data[q+j]);
|
||||||
|
result._data[r+j] = temp;
|
||||||
|
}
|
||||||
|
|
||||||
|
return result;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Scalar sum.
|
||||||
|
* \param k
|
||||||
|
* \return the resultant matrix
|
||||||
|
*/
|
||||||
|
Matrix<TYPE> operator+(const TYPE k)
|
||||||
|
{
|
||||||
|
Matrix<TYPE> result(_rows, _columns);
|
||||||
|
for (unsigned int i=0; i<_rows*_columns; i++)
|
||||||
|
result._data[i] = _data[i]+k;
|
||||||
|
return result;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Scalar difference.
|
||||||
|
* \param k
|
||||||
|
* \return the resultant matrix
|
||||||
|
*/
|
||||||
|
Matrix<TYPE> operator-(const TYPE k)
|
||||||
|
{
|
||||||
|
Matrix<TYPE> result(_rows, _columns);
|
||||||
|
for (unsigned int i=0; i<_rows*_columns; i++)
|
||||||
|
results._data[i] = _data[i]-k;
|
||||||
|
return result;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Negate all matrix elements
|
||||||
|
* \return ...
|
||||||
|
*/
|
||||||
|
Matrix<TYPE> operator-() const
|
||||||
|
{
|
||||||
|
Matrix<TYPE> result(_rows, _columns, _data);
|
||||||
|
for (unsigned int i=0; i<_columns*_rows; i++)
|
||||||
|
result._data[i] = -1*_data[i];
|
||||||
|
return result;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Scalar multiplication.
|
||||||
|
* \param k value to multiply every member by
|
||||||
|
* \return the resultant matrix
|
||||||
|
*/
|
||||||
|
Matrix<TYPE> operator*(const TYPE k)
|
||||||
|
{
|
||||||
|
Matrix<TYPE> result(_rows, _columns);
|
||||||
|
for (unsigned int i=0; i<_rows*_columns; i++)
|
||||||
|
results._data[i] = _data[i]*k;
|
||||||
|
return result;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Scalar division.
|
||||||
|
* \param k value to divide every member by
|
||||||
|
* \return the resultant matrix
|
||||||
|
*/
|
||||||
|
Matrix<TYPE> operator/(const TYPE k)
|
||||||
|
{
|
||||||
|
Matrix<TYPE> result(_rows, _columns);
|
||||||
|
for (unsigned int i=0; i<_rows*_columns; i++)
|
||||||
|
results._data[i] = _data[i]/k;
|
||||||
|
return result;
|
||||||
|
};
|
||||||
|
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Set all the matrix elements to zero.
|
||||||
|
*/
|
||||||
|
void SetZero()
|
||||||
|
{
|
||||||
|
for (unsigned int i=0; i<_rows*_columns; i++)
|
||||||
|
_data[i] = ScalarType(0.0);
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Set the values of <I>j</I>-th column to v[j]
|
||||||
|
* \param j the column index
|
||||||
|
* \param v ...
|
||||||
|
*/
|
||||||
|
void SetColumn(const unsigned int j, TYPE* v)
|
||||||
|
{
|
||||||
|
assert(j>=0 && j<_columns);
|
||||||
|
unsigned int i, p;
|
||||||
|
for (i=0, p=0; i<_rows; i++, p+=_columns)
|
||||||
|
_data[p] = v[i];
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Set the elements of the <I>i</I>-th row to v[j]
|
||||||
|
* \param i the row index
|
||||||
|
* \param v ...
|
||||||
|
*/
|
||||||
|
void SetRow(const unsigned int i, TYPE* v)
|
||||||
|
{
|
||||||
|
assert(i>=0 && i<_rows);
|
||||||
|
unsigned int j, p;
|
||||||
|
for (j=0, p=i*_rows; j<_columns; j++, p++)
|
||||||
|
_data[p] = v[j];
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Resize the current matrix.
|
||||||
|
* \param m the number of matrix rows.
|
||||||
|
* \param n the number of matrix columns.
|
||||||
|
*/
|
||||||
|
void Resize(const unsigned int m, const unsigned int n)
|
||||||
|
{
|
||||||
|
assert(m>=2);
|
||||||
|
assert(n>=2);
|
||||||
|
_rows = m;
|
||||||
|
_columns = n;
|
||||||
|
delete []_data;
|
||||||
|
_data = new ScalarType[m*n];
|
||||||
|
for (unsigned int i=0; i<m*n; i++)
|
||||||
|
_data[i] = 0;
|
||||||
|
};
|
||||||
|
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Matrix transposition operation: set the current matrix to its transpose
|
||||||
|
*/
|
||||||
|
void Transpose()
|
||||||
|
{
|
||||||
|
ScalarType *temp = new ScalarType[_rows*_columns];
|
||||||
|
unsigned int i, j, p, q;
|
||||||
|
for (i=0, p=0; i<_rows; i++, p+=_columns)
|
||||||
|
for (j=0, q=0; j<_columns; j++, q+=_rows)
|
||||||
|
temp[q+i] = _data[p+j];
|
||||||
|
|
||||||
|
std::swap(_columns, _rows);
|
||||||
|
std::swap(_data, temp);
|
||||||
|
delete []temp;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*!
|
||||||
|
* Print all matrix elements
|
||||||
|
*/
|
||||||
|
void Dump()
|
||||||
|
{
|
||||||
|
unsigned int i, j, p;
|
||||||
|
for (i=0, p=0; i<_rows; i++, p+=_columns)
|
||||||
|
{
|
||||||
|
printf("[\t");
|
||||||
|
for (j=0; j<_columns; j++)
|
||||||
|
printf("%g\t", _data[p+j]);
|
||||||
|
|
||||||
|
printf("]\n");
|
||||||
|
}
|
||||||
|
printf("\n");
|
||||||
|
};
|
||||||
|
|
||||||
|
protected:
|
||||||
|
/// the number of matrix rows
|
||||||
|
unsigned int _rows;
|
||||||
|
|
||||||
|
/// the number of matrix rows
|
||||||
|
unsigned int _columns;
|
||||||
|
|
||||||
|
/// the matrix elements
|
||||||
|
ScalarType *_data;
|
||||||
|
};
|
||||||
|
|
||||||
|
/*! @} */
|
||||||
|
}; // end of namespace
|
||||||
Loading…
Reference in New Issue