HEAVY CHANGE. Further cleaning of the matrix classes of VCG.
Get rid of all the unused stuff. internally use Eigen for computing Inverse. Removed the stupid incomprehensible method Invert() that changed the matrix itself. Nobody was using it in a reasonable way.
This commit is contained in:
parent
f156a5a82c
commit
2e65cae10e
|
|
@ -20,75 +20,6 @@
|
|||
* for more details. *
|
||||
* *
|
||||
****************************************************************************/
|
||||
/****************************************************************************
|
||||
History
|
||||
|
||||
$Log: not supported by cvs2svn $
|
||||
Revision 1.18 2007/04/19 14:30:26 pietroni
|
||||
added RotationMatrix method to calculate rotation matrix along an axis
|
||||
|
||||
Revision 1.17 2007/04/07 23:06:47 pietroni
|
||||
Added function RotationMatrix
|
||||
|
||||
Revision 1.16 2007/01/29 00:20:25 pietroni
|
||||
-Used scalar type passed as template argument istead of double to prevent warnings.. in Rotate function
|
||||
|
||||
Revision 1.15 2006/09/25 23:05:29 ganovelli
|
||||
added constructor from matrix44 excluding a row and colum
|
||||
|
||||
Revision 1.14 2006/06/22 08:00:05 ganovelli
|
||||
bug in operator + with MatrixxDig
|
||||
|
||||
Revision 1.13 2006/01/20 16:41:44 pietroni
|
||||
added operators:
|
||||
operator -= ( const Matrix33Diag<S> &p )
|
||||
Matrix33 operator - ( const Matrix33Diag<S> &p )
|
||||
Matrix33 operator + ( const Matrix33 &m )
|
||||
Matrix33 operator + ( const Matrix33Diag<S> &p )
|
||||
|
||||
Revision 1.12 2005/11/14 10:28:25 cignoni
|
||||
Changed Invert -> FastInvert for the function based on the maple expansion
|
||||
|
||||
Revision 1.11 2005/10/13 15:45:23 ponchio
|
||||
Changed a Zero in SetZero in WeightedCrossCovariance() (again)
|
||||
|
||||
Revision 1.10 2005/10/05 17:06:12 pietroni
|
||||
corrected sintax error on singular value decomposition
|
||||
|
||||
Revision 1.9 2005/09/29 09:53:58 ganovelli
|
||||
added inverse by SVD
|
||||
|
||||
Revision 1.8 2005/06/10 14:51:54 cignoni
|
||||
Changed a Zero in SetZero in WeightedCrossCovariance()
|
||||
|
||||
Revision 1.7 2005/06/10 11:46:49 pietroni
|
||||
Added Norm Function
|
||||
|
||||
Revision 1.6 2005/06/07 14:29:56 ganovelli
|
||||
changed from Matrix33Ide to MatrixeeDiag
|
||||
|
||||
Revision 1.5 2005/05/23 15:05:26 ganovelli
|
||||
Matrix33Diag Added: it implements diagonal matrix. Added only operator += in Matrix33
|
||||
|
||||
Revision 1.4 2005/04/11 14:11:22 pietroni
|
||||
changed swap to math::Swap in Traspose Function
|
||||
|
||||
Revision 1.3 2004/10/18 15:03:02 fiorin
|
||||
Updated interface: all Matrix classes have now the same interface
|
||||
|
||||
Revision 1.2 2004/07/13 06:48:26 cignoni
|
||||
removed uppercase references in include
|
||||
|
||||
Revision 1.1 2004/05/28 13:09:05 ganovelli
|
||||
created
|
||||
|
||||
Revision 1.1 2004/05/28 13:00:39 ganovelli
|
||||
created
|
||||
|
||||
|
||||
****************************************************************************/
|
||||
|
||||
|
||||
#ifndef __VCGLIB_MATRIX33_H
|
||||
#define __VCGLIB_MATRIX33_H
|
||||
|
||||
|
|
@ -99,24 +30,11 @@ created
|
|||
|
||||
namespace vcg {
|
||||
|
||||
template <class S>
|
||||
class Matrix33Diag:public Point3<S>{
|
||||
public:
|
||||
|
||||
/** @name Matrix33
|
||||
Class Matrix33Diag.
|
||||
This is the class for definition of a diagonal matrix 3x3.
|
||||
@param S (Templete Parameter) Specifies the ScalarType field.
|
||||
*/
|
||||
Matrix33Diag(const S & p0,const S & p1,const S & p2):Point3<S>(p0,p1,p2){};
|
||||
Matrix33Diag(const Point3<S>&p ):Point3<S>(p){};
|
||||
};
|
||||
|
||||
template<class S>
|
||||
/** @name Matrix33
|
||||
Class Matrix33.
|
||||
This is the class for definition of a matrix 3x3.
|
||||
@param S (Templete Parameter) Specifies the ScalarType field.
|
||||
@param S (Template Parameter) Specifies the ScalarType field.
|
||||
*/
|
||||
class Matrix33
|
||||
{
|
||||
|
|
@ -206,14 +124,6 @@ public:
|
|||
return *this;
|
||||
}
|
||||
|
||||
/// Modificatore somma per matrici 3x3
|
||||
Matrix33 & operator += ( const Matrix33Diag<S> &p )
|
||||
{
|
||||
a[0] += p[0];
|
||||
a[4] += p[1];
|
||||
a[8] += p[2];
|
||||
return *this;
|
||||
}
|
||||
|
||||
/// Modificatore sottrazione per matrici 3x3
|
||||
Matrix33 & operator -= ( const Matrix33 &m )
|
||||
|
|
@ -223,15 +133,6 @@ public:
|
|||
return *this;
|
||||
}
|
||||
|
||||
/// Modificatore somma per matrici 3x3
|
||||
Matrix33 & operator -= ( const Matrix33Diag<S> &p )
|
||||
{
|
||||
a[0] -= p[0];
|
||||
a[4] -= p[1];
|
||||
a[8] -= p[2];
|
||||
return *this;
|
||||
}
|
||||
|
||||
/// Modificatore divisione per scalare
|
||||
Matrix33 & operator /= ( const S &s )
|
||||
{
|
||||
|
|
@ -265,30 +166,6 @@ public:
|
|||
for(i=0;i<9;++i) this->a[i] = r.a[i];
|
||||
}
|
||||
|
||||
/// Dot product with a diagonal matrix
|
||||
Matrix33 operator * ( const Matrix33Diag< S> & t ) const
|
||||
{
|
||||
Matrix33<S> r;
|
||||
|
||||
int i,j;
|
||||
for(i=0;i<3;++i)
|
||||
for(j=0;j<3;++j)
|
||||
r[i][j] = (*this)[i][j]*t[j];
|
||||
|
||||
return r;
|
||||
}
|
||||
|
||||
/// Dot product modifier with a diagonal matrix
|
||||
void operator *=( const Matrix33Diag< S> & t )
|
||||
{
|
||||
Matrix33<S> r;
|
||||
int i,j;
|
||||
for(i=0;i<3;++i)
|
||||
for(j=0;j<3;++j)
|
||||
r[i][j] = (*this)[i][j]*t[j];
|
||||
for(i=0;i<9;++i) this->a[i] = r.a[i];
|
||||
}
|
||||
|
||||
/// Modificatore prodotto per costante
|
||||
Matrix33 & operator *= ( const S t )
|
||||
{
|
||||
|
|
@ -316,17 +193,6 @@ public:
|
|||
|
||||
return r;
|
||||
}
|
||||
|
||||
/// Operatore sottrazione di matrici 3x3 con matrici diagonali
|
||||
Matrix33 operator - ( const Matrix33Diag<S> &p ) const
|
||||
{
|
||||
Matrix33<S> r=a;
|
||||
r[0][0] -= p[0];
|
||||
r[1][1] -= p[1];
|
||||
r[2][2] -= p[2];
|
||||
return r;
|
||||
}
|
||||
|
||||
/// Operatore sottrazione per matrici 3x3
|
||||
Matrix33 operator + ( const Matrix33 &m ) const
|
||||
{
|
||||
|
|
@ -336,17 +202,6 @@ public:
|
|||
|
||||
return r;
|
||||
}
|
||||
|
||||
/// Operatore addizione di matrici 3x3 con matrici diagonali
|
||||
Matrix33 operator + ( const Matrix33Diag<S> &p ) const
|
||||
{
|
||||
Matrix33<S> r=(*this);
|
||||
r[0][0] += p[0];
|
||||
r[1][1] += p[1];
|
||||
r[2][2] += p[2];
|
||||
return r;
|
||||
}
|
||||
|
||||
/** Operatore per il prodotto matrice-vettore.
|
||||
@param v A point in $R^{3}$
|
||||
@return Il vettore risultante in $R^{3}$
|
||||
|
|
@ -486,42 +341,6 @@ public:
|
|||
a[2]*(a[3]*a[7]-a[4]*a[6]) ;
|
||||
}
|
||||
|
||||
// Warning, this Inversion code can be HIGHLY NUMERICALLY UNSTABLE!
|
||||
// In most case you are advised to use the Invert() method based on SVD decomposition.
|
||||
|
||||
Matrix33 & FastInvert()
|
||||
{
|
||||
// Maple produsse:
|
||||
S t4 = a[0]*a[4];
|
||||
S t6 = a[0]*a[5];
|
||||
S t8 = a[1]*a[3];
|
||||
S t10 = a[2]*a[3];
|
||||
S t12 = a[1]*a[6];
|
||||
S t14 = a[2]*a[6];
|
||||
S t17 = 1/(t4*a[8]-t6*a[7]-t8*a[8]+t10*a[7]+t12*a[5]-t14*a[4]);
|
||||
S a0 = a[0];
|
||||
S a1 = a[1];
|
||||
S a3 = a[3];
|
||||
S a4 = a[4];
|
||||
a[0] = (a[4]*a[8]-a[5]*a[7])*t17;
|
||||
a[1] = -(a[1]*a[8]-a[2]*a[7])*t17;
|
||||
a[2] = (a1 *a[5]-a[2]*a[4])*t17;
|
||||
a[3] = -(a[3]*a[8]-a[5]*a[6])*t17;
|
||||
a[4] = (a0 *a[8]-t14 )*t17;
|
||||
a[5] = -(t6 - t10)*t17;
|
||||
a[6] = (a3 *a[7]-a[4]*a[6])*t17;
|
||||
a[7] = -(a[0]*a[7]-t12)*t17;
|
||||
a[8] = (t4-t8)*t17;
|
||||
|
||||
return *this;
|
||||
}
|
||||
|
||||
void show(FILE * /*fp*/)
|
||||
{
|
||||
for(int i=0;i<3;++i)
|
||||
printf("| %g \t%g \t%g |\n",a[3*i+0],a[3*i+1],a[3*i+2]);
|
||||
}
|
||||
|
||||
// return the Trace of the matrix i.e. the sum of the diagonal elements
|
||||
S Trace() const
|
||||
{
|
||||
|
|
@ -663,27 +482,14 @@ vcg::Matrix33<S> TransformationMatrix(const vcg::Point3<S> dirX,
|
|||
/////then find the inverse
|
||||
return (Trans);
|
||||
}
|
||||
|
||||
template <class S>
|
||||
void Invert(Matrix33<S> &m)
|
||||
{
|
||||
Matrix33<S> v;
|
||||
Point3<typename Matrix33<S>::ScalarType> e;
|
||||
SingularValueDecomposition(m,&e[0],v);
|
||||
e[0]=1/e[0];e[1]=1/e[1];e[2]=1/e[2];
|
||||
m.Transpose();
|
||||
m = v * Matrix33Diag<S>(e) * m;
|
||||
}
|
||||
|
||||
template <class S>
|
||||
Matrix33<S> Inverse(const Matrix33<S>&m)
|
||||
{
|
||||
Matrix33<S> v,m_copy=m;
|
||||
Point3<S> e;
|
||||
SingularValueDecomposition(m_copy,&e[0],v);
|
||||
m_copy.Transpose();
|
||||
e[0]=1/e[0];e[1]=1/e[1];e[2]=1/e[2];
|
||||
return v * Matrix33Diag<S>(e) * m_copy;
|
||||
Eigen::Matrix3d mm,mmi;
|
||||
m.ToEigenMatrix(mm);
|
||||
mmi=mm.inverse();
|
||||
Matrix33<S> res;
|
||||
res.FromEigenMatrix(mmi);
|
||||
}
|
||||
|
||||
///given 2 vector centered into origin calculate the rotation matrix from first to the second
|
||||
|
|
|
|||
|
|
@ -19,84 +19,6 @@
|
|||
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
|
||||
* for more details. *
|
||||
* *
|
||||
****************************************************************************/
|
||||
/****************************************************************************
|
||||
History
|
||||
|
||||
$Log: not supported by cvs2svn $
|
||||
Revision 1.34 2007/07/12 06:42:01 cignoni
|
||||
added the missing static Construct() member
|
||||
|
||||
Revision 1.33 2007/07/03 16:06:48 corsini
|
||||
add DCM to Euler Angles conversion
|
||||
|
||||
Revision 1.32 2007/03/08 14:39:27 corsini
|
||||
final fix to euler angles transformation
|
||||
|
||||
Revision 1.31 2007/02/06 09:57:40 corsini
|
||||
fix euler angles computation
|
||||
|
||||
Revision 1.30 2007/02/05 14:16:33 corsini
|
||||
add from euler angles to rotation matrix conversion
|
||||
|
||||
Revision 1.29 2005/12/02 09:46:49 croccia
|
||||
Corrected bug in == and != Matrix44 operators
|
||||
|
||||
Revision 1.28 2005/06/28 17:42:47 ganovelli
|
||||
added Matrix44Diag
|
||||
|
||||
Revision 1.27 2005/06/17 05:28:47 cignoni
|
||||
Completed Shear Matrix code and comments,
|
||||
Added use of swap inside Transpose
|
||||
Added more complete comments on the usage of Decompose
|
||||
|
||||
Revision 1.26 2005/06/10 15:04:12 cignoni
|
||||
Added Various missing functions: SetShearXY, SetShearXZ, SetShearYZ, SetScale for point3 and Decompose
|
||||
Completed *=(scalar); made uniform GetRow and GetColumn
|
||||
|
||||
Revision 1.25 2005/04/14 11:35:09 ponchio
|
||||
*** empty log message ***
|
||||
|
||||
Revision 1.24 2005/03/18 00:14:39 cignoni
|
||||
removed small gcc compiling issues
|
||||
|
||||
Revision 1.23 2005/03/15 11:40:56 cignoni
|
||||
Added operator*=( std::vector<PointType> ...) to apply a matrix to a vector of vertexes (replacement of the old style mesh.Apply(tr).
|
||||
|
||||
Revision 1.22 2004/12/15 18:45:50 tommyfranken
|
||||
*** empty log message ***
|
||||
|
||||
Revision 1.21 2004/10/22 14:41:30 ponchio
|
||||
return in operator+ added.
|
||||
|
||||
Revision 1.20 2004/10/18 15:03:14 fiorin
|
||||
Updated interface: all Matrix classes have now the same interface
|
||||
|
||||
Revision 1.19 2004/10/07 14:23:57 ganovelli
|
||||
added function to take rows and comlumns. Added toMatrix and fromMatrix to comply
|
||||
RotationTYpe prototype in Similarity.h
|
||||
|
||||
Revision 1.18 2004/05/28 13:01:50 ganovelli
|
||||
changed scalar to ScalarType
|
||||
|
||||
Revision 1.17 2004/05/26 15:09:32 cignoni
|
||||
better comments in set rotate
|
||||
|
||||
Revision 1.16 2004/05/07 10:05:50 cignoni
|
||||
Corrected abuse of for index variable scope
|
||||
|
||||
Revision 1.15 2004/05/04 23:19:41 cignoni
|
||||
Clarified initial comment, removed vector*matrix operator (confusing!)
|
||||
Corrected translate and Rotate, removed gl stuff.
|
||||
|
||||
Revision 1.14 2004/05/04 02:34:03 ganovelli
|
||||
wrong use of operator [] corrected
|
||||
|
||||
Revision 1.13 2004/04/07 10:45:54 cignoni
|
||||
Added: [i][j] access, V() for the raw float values, constructor from T[16]
|
||||
|
||||
Revision 1.12 2004/03/25 14:57:49 ponchio
|
||||
|
||||
****************************************************************************/
|
||||
|
||||
#ifndef __VCGLIB_MATRIX44
|
||||
|
|
@ -108,7 +30,8 @@ Revision 1.12 2004/03/25 14:57:49 ponchio
|
|||
#include <vcg/space/point4.h>
|
||||
#include <vector>
|
||||
#include <iostream>
|
||||
|
||||
#include <eigenlib/Eigen/Core>
|
||||
#include <eigenlib/Eigen/LU>
|
||||
|
||||
namespace vcg {
|
||||
|
||||
|
|
@ -146,19 +69,6 @@ for 'column' vectors.
|
|||
|
||||
*/
|
||||
|
||||
template <class S>
|
||||
class Matrix44Diag:public Point4<S>{
|
||||
public:
|
||||
/** @name Matrix33
|
||||
Class Matrix33Diag.
|
||||
This is the class for definition of a diagonal matrix 4x4.
|
||||
@param S (Templete Parameter) Specifies the ScalarType field.
|
||||
*/
|
||||
Matrix44Diag(const S & p0,const S & p1,const S & p2,const S & p3):Point4<S>(p0,p1,p2,p3){};
|
||||
Matrix44Diag( const Point4<S> & p ):Point4<S>(p){};
|
||||
};
|
||||
|
||||
|
||||
/** This class represent a 4x4 matrix. T is the kind of element in the matrix.
|
||||
*/
|
||||
template <class T> class Matrix44 {
|
||||
|
|
@ -173,23 +83,11 @@ public:
|
|||
/** $name Constructors
|
||||
* No automatic casting and default constructor is empty
|
||||
*/
|
||||
Matrix44() {};
|
||||
~Matrix44() {};
|
||||
Matrix44() {}
|
||||
~Matrix44() {}
|
||||
Matrix44(const Matrix44 &m);
|
||||
Matrix44(const T v[]);
|
||||
|
||||
/// Number of columns
|
||||
inline unsigned int ColumnsNumber() const
|
||||
{
|
||||
return 4;
|
||||
};
|
||||
|
||||
/// Number of rows
|
||||
inline unsigned int RowsNumber() const
|
||||
{
|
||||
return 4;
|
||||
};
|
||||
|
||||
T &ElementAt(const int row, const int col);
|
||||
T ElementAt(const int row, const int col) const;
|
||||
//T &operator[](const int i);
|
||||
|
|
@ -227,7 +125,6 @@ public:
|
|||
Matrix44 operator+(const Matrix44 &m) const;
|
||||
Matrix44 operator-(const Matrix44 &m) const;
|
||||
Matrix44 operator*(const Matrix44 &m) const;
|
||||
Matrix44 operator*(const Matrix44Diag<T> &m) const;
|
||||
Point4<T> operator*(const Point4<T> &v) const;
|
||||
|
||||
bool operator==(const Matrix44 &m) const;
|
||||
|
|
@ -248,6 +145,13 @@ public:
|
|||
template <class Matrix44Type>
|
||||
void FromMatrix(const Matrix44Type & m){for(int i = 0; i < 16; i++) V()[i]=m.V()[i];}
|
||||
|
||||
template <class EigenMatrix44Type>
|
||||
void ToEigenMatrix(EigenMatrix44Type & m) const {
|
||||
for(int i = 0; i < 4; i++)
|
||||
for(int j = 0; j < 4; j++)
|
||||
m(i,j)=(*this)[i][j];
|
||||
}
|
||||
|
||||
template <class EigenMatrix44Type>
|
||||
void FromEigenMatrix(const EigenMatrix44Type & m){
|
||||
for(int i = 0; i < 4; i++)
|
||||
|
|
@ -340,7 +244,6 @@ template <class T> Point3<T> operator*(const Matrix44<T> &m, const Point3<T> &p)
|
|||
|
||||
template <class T> Matrix44<T> &Transpose(Matrix44<T> &m);
|
||||
//return NULL matrix if not invertible
|
||||
template <class T> Matrix44<T> &Invert(Matrix44<T> &m);
|
||||
template <class T> Matrix44<T> Inverse(const Matrix44<T> &m);
|
||||
|
||||
typedef Matrix44<short> Matrix44s;
|
||||
|
|
@ -418,14 +321,6 @@ template <class T> Matrix44<T> Matrix44<T>::operator*(const Matrix44 &m) const {
|
|||
return ret;
|
||||
}
|
||||
|
||||
template <class T> Matrix44<T> Matrix44<T>::operator*(const Matrix44Diag<T> &m) const {
|
||||
Matrix44 ret = (*this);
|
||||
for(int i = 0; i < 4; ++i)
|
||||
for(int j = 0; j < 4; ++j)
|
||||
ret[i][j]*=m[i];
|
||||
return ret;
|
||||
}
|
||||
|
||||
template <class T> Point4<T> Matrix44<T>::operator*(const Point4<T> &v) const {
|
||||
Point4<T> ret;
|
||||
for(int i = 0; i < 4; i++){
|
||||
|
|
@ -477,17 +372,6 @@ template <class T> void Matrix44<T>::operator-=(const Matrix44 &m) {
|
|||
}
|
||||
template <class T> void Matrix44<T>::operator*=( const Matrix44 & m ) {
|
||||
*this = *this *m;
|
||||
|
||||
/*for(int i = 0; i < 4; i++) { //sbagliato
|
||||
Point4<T> t(0, 0, 0, 0);
|
||||
for(int k = 0; k < 4; k++) {
|
||||
for(int j = 0; j < 4; j++) {
|
||||
t[k] += ElementAt(i, k) * m.ElementAt(k, j);
|
||||
}
|
||||
}
|
||||
for(int l = 0; l < 4; l++)
|
||||
ElementAt(i, l) = t[l];
|
||||
} */
|
||||
}
|
||||
|
||||
template < class PointType , class T > void operator*=( std::vector<PointType> &vert, const Matrix44<T> & m ) {
|
||||
|
|
@ -625,44 +509,6 @@ template <class T> Matrix44<T> &Matrix44<T>::SetRotateRad(T AngleRad, const Poin
|
|||
return *this;
|
||||
}
|
||||
|
||||
/* Shear Matrixes
|
||||
XY
|
||||
1 k 0 0 x x+ky
|
||||
0 1 0 0 y y
|
||||
0 0 1 0 z z
|
||||
0 0 0 1 1 1
|
||||
|
||||
1 0 k 0 x x+kz
|
||||
0 1 0 0 y y
|
||||
0 0 1 0 z z
|
||||
0 0 0 1 1 1
|
||||
|
||||
1 1 0 0 x x
|
||||
0 1 k 0 y y+kz
|
||||
0 0 1 0 z z
|
||||
0 0 0 1 1 1
|
||||
|
||||
*/
|
||||
|
||||
template <class T> Matrix44<T> & Matrix44<T>:: SetShearXY( const T sh) {// shear the X coordinate as the Y coordinate change
|
||||
SetIdentity();
|
||||
ElementAt(0,1) = sh;
|
||||
return *this;
|
||||
}
|
||||
|
||||
template <class T> Matrix44<T> & Matrix44<T>:: SetShearXZ( const T sh) {// shear the X coordinate as the Z coordinate change
|
||||
SetIdentity();
|
||||
ElementAt(0,2) = sh;
|
||||
return *this;
|
||||
}
|
||||
|
||||
template <class T> Matrix44<T> &Matrix44<T>:: SetShearYZ( const T sh) {// shear the Y coordinate as the Z coordinate change
|
||||
SetIdentity();
|
||||
ElementAt(1,2) = sh;
|
||||
return *this;
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
Given a non singular, non projective matrix (e.g. with the last row equal to [0,0,0,1] )
|
||||
This procedure decompose it in a sequence of
|
||||
|
|
@ -795,8 +641,9 @@ bool Decompose(Matrix44<T> &M, Point3<T> &ScaleV, Point3<T> &ShearV, Point3<T> &
|
|||
|
||||
|
||||
template <class T> T Matrix44<T>::Determinant() const {
|
||||
LinearSolve<T> solve(*this);
|
||||
return solve.Determinant();
|
||||
Eigen::Matrix4d mm;
|
||||
this->ToEigenMatrix(mm);
|
||||
return mm.determinant();
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -811,17 +658,6 @@ template <class T> Point3<T> operator*(const Matrix44<T> &m, const Point3<T> &p)
|
|||
return s;
|
||||
}
|
||||
|
||||
//template <class T> Point3<T> operator*(const Point3<T> &p, const Matrix44<T> &m) {
|
||||
// T w;
|
||||
// Point3<T> s;
|
||||
// s[0] = m.ElementAt(0, 0)*p[0] + m.ElementAt(1, 0)*p[1] + m.ElementAt(2, 0)*p[2] + m.ElementAt(3, 0);
|
||||
// s[1] = m.ElementAt(0, 1)*p[0] + m.ElementAt(1, 1)*p[1] + m.ElementAt(2, 1)*p[2] + m.ElementAt(3, 1);
|
||||
// s[2] = m.ElementAt(0, 2)*p[0] + m.ElementAt(1, 2)*p[1] + m.ElementAt(2, 2)*p[2] + m.ElementAt(3, 2);
|
||||
// w = m.ElementAt(0, 3)*p[0] + m.ElementAt(1, 3)*p[1] + m.ElementAt(2, 3)*p[2] + m.ElementAt(3, 3);
|
||||
// if(w != 0) s /= w;
|
||||
// return s;
|
||||
//}
|
||||
|
||||
template <class T> Matrix44<T> &Transpose(Matrix44<T> &m) {
|
||||
for(int i = 1; i < 4; i++)
|
||||
for(int j = 0; j < i; j++) {
|
||||
|
|
@ -830,202 +666,12 @@ template <class T> Matrix44<T> &Transpose(Matrix44<T> &m) {
|
|||
return m;
|
||||
}
|
||||
|
||||
/*
|
||||
To invert a matrix you can
|
||||
either invert the matrix inplace calling
|
||||
|
||||
vcg::Invert(yourMatrix);
|
||||
|
||||
or get the inverse matrix of a given matrix without touching it:
|
||||
|
||||
invertedMatrix = vcg::Inverse(untouchedMatrix);
|
||||
|
||||
*/
|
||||
template <class T> Matrix44<T> & Invert(Matrix44<T> &m) {
|
||||
LinearSolve<T> solve(m);
|
||||
|
||||
for(int j = 0; j < 4; j++) { //Find inverse by columns.
|
||||
Point4<T> col(0, 0, 0, 0);
|
||||
col[j] = 1.0;
|
||||
col = solve.Solve(col);
|
||||
for(int i = 0; i < 4; i++)
|
||||
m.ElementAt(i, j) = col[i];
|
||||
}
|
||||
return m;
|
||||
}
|
||||
|
||||
template <class T> Matrix44<T> Inverse(const Matrix44<T> &m) {
|
||||
LinearSolve<T> solve(m);
|
||||
Eigen::Matrix4d mm,mmi;
|
||||
m.ToEigenMatrix(mm);
|
||||
mmi=mm.inverse();
|
||||
Matrix44<T> res;
|
||||
for(int j = 0; j < 4; j++) { //Find inverse by columns.
|
||||
Point4<T> col(0, 0, 0, 0);
|
||||
col[j] = 1.0;
|
||||
col = solve.Solve(col);
|
||||
for(int i = 0; i < 4; i++)
|
||||
res.ElementAt(i, j) = col[i];
|
||||
}
|
||||
return res;
|
||||
}
|
||||
|
||||
|
||||
|
||||
/* LINEAR SOLVE IMPLEMENTATION */
|
||||
|
||||
template <class T> LinearSolve<T>::LinearSolve(const Matrix44<T> &m): Matrix44<T>(m) {
|
||||
if(!Decompose()) {
|
||||
for(int i = 0; i < 4; i++)
|
||||
index[i] = i;
|
||||
Matrix44<T>::SetZero();
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template <class T> T LinearSolve<T>::Determinant() const {
|
||||
T det = d;
|
||||
for(int j = 0; j < 4; j++)
|
||||
det *= this-> ElementAt(j, j);
|
||||
return det;
|
||||
}
|
||||
|
||||
|
||||
/*replaces a matrix by its LU decomposition of a rowwise permutation.
|
||||
d is +or -1 depeneing of row permutation even or odd.*/
|
||||
#define TINY 1e-100
|
||||
|
||||
template <class T> bool LinearSolve<T>::Decompose() {
|
||||
|
||||
/* Matrix44<T> A;
|
||||
for(int i = 0; i < 16; i++)
|
||||
A[i] = operator[](i);
|
||||
SetIdentity();
|
||||
Point4<T> scale;
|
||||
// Set scale factor, scale(i) = max( |a(i,j)| ), for each row
|
||||
for(int i = 0; i < 4; i++ ) {
|
||||
index[i] = i; // Initialize row index list
|
||||
T scalemax = (T)0.0;
|
||||
for(int j = 0; j < 4; j++)
|
||||
scalemax = (scalemax > math::Abs(A.ElementAt(i,j))) ? scalemax : math::Abs(A.ElementAt(i,j));
|
||||
scale[i] = scalemax;
|
||||
}
|
||||
|
||||
// Loop over rows k = 1, ..., (N-1)
|
||||
int signDet = 1;
|
||||
for(int k = 0; k < 3; k++ ) {
|
||||
// Select pivot row from max( |a(j,k)/s(j)| )
|
||||
T ratiomax = (T)0.0;
|
||||
int jPivot = k;
|
||||
for(int i = k; i < 4; i++ ) {
|
||||
T ratio = math::Abs(A.ElementAt(index[i], k))/scale[index[i]];
|
||||
if(ratio > ratiomax) {
|
||||
jPivot = i;
|
||||
ratiomax = ratio;
|
||||
}
|
||||
}
|
||||
// Perform pivoting using row index list
|
||||
int indexJ = index[k];
|
||||
if( jPivot != k ) { // Pivot
|
||||
indexJ = index[jPivot];
|
||||
index[jPivot] = index[k]; // Swap index jPivot and k
|
||||
index[k] = indexJ;
|
||||
signDet *= -1; // Flip sign of determinant
|
||||
}
|
||||
// Perform forward elimination
|
||||
for(int i=k+1; i < 4; i++ ) {
|
||||
T coeff = A.ElementAt(index[i],k)/A.ElementAt(indexJ,k);
|
||||
for(int j=k+1; j < 4; j++ )
|
||||
A.ElementAt(index[i],j) -= coeff*A.ElementAt(indexJ,j);
|
||||
A.ElementAt(index[i],k) = coeff;
|
||||
//for( j=1; j< 4; j++ )
|
||||
// ElementAt(index[i],j) -= A.ElementAt(index[i], k)*ElementAt(indexJ, j);
|
||||
}
|
||||
}
|
||||
for(int i = 0; i < 16; i++)
|
||||
operator[](i) = A[i];
|
||||
|
||||
d = signDet;
|
||||
// this = A;
|
||||
return true; */
|
||||
|
||||
d = 1; //no permutation still
|
||||
|
||||
T scaling[4];
|
||||
int i,j,k;
|
||||
//Saving the scvaling information per row
|
||||
for(i = 0; i < 4; i++) {
|
||||
T largest = 0.0;
|
||||
for(j = 0; j < 4; j++) {
|
||||
T t = math::Abs(this->ElementAt(i, j));
|
||||
if (t > largest) largest = t;
|
||||
}
|
||||
|
||||
if (largest == 0.0) { //oooppps there is a zero row!
|
||||
return false;
|
||||
}
|
||||
scaling[i] = (T)1.0 / largest;
|
||||
}
|
||||
|
||||
int imax = 0;
|
||||
for(j = 0; j < 4; j++) {
|
||||
for(i = 0; i < j; i++) {
|
||||
T sum = this->ElementAt(i,j);
|
||||
for(int k = 0; k < i; k++)
|
||||
sum -= this->ElementAt(i,k)*this->ElementAt(k,j);
|
||||
this->ElementAt(i,j) = sum;
|
||||
}
|
||||
T largest = 0.0;
|
||||
for(i = j; i < 4; i++) {
|
||||
T sum = this->ElementAt(i,j);
|
||||
for(k = 0; k < j; k++)
|
||||
sum -= this->ElementAt(i,k)*this->ElementAt(k,j);
|
||||
this->ElementAt(i,j) = sum;
|
||||
T t = scaling[i] * math::Abs(sum);
|
||||
if(t >= largest) {
|
||||
largest = t;
|
||||
imax = i;
|
||||
}
|
||||
}
|
||||
if (j != imax) {
|
||||
for (int k = 0; k < 4; k++) {
|
||||
T dum = this->ElementAt(imax,k);
|
||||
this->ElementAt(imax,k) = this->ElementAt(j,k);
|
||||
this->ElementAt(j,k) = dum;
|
||||
}
|
||||
d = -(d);
|
||||
scaling[imax] = scaling[j];
|
||||
}
|
||||
index[j]=imax;
|
||||
if (this->ElementAt(j,j) == 0.0) this->ElementAt(j,j) = (T)TINY;
|
||||
if (j != 3) {
|
||||
T dum = (T)1.0 / (this->ElementAt(j,j));
|
||||
for (i = j+1; i < 4; i++)
|
||||
this->ElementAt(i,j) *= dum;
|
||||
}
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
template <class T> Point4<T> LinearSolve<T>::Solve(const Point4<T> &b) {
|
||||
Point4<T> x(b);
|
||||
int first = -1, ip;
|
||||
for(int i = 0; i < 4; i++) {
|
||||
ip = index[i];
|
||||
T sum = x[ip];
|
||||
x[ip] = x[i];
|
||||
if(first!= -1)
|
||||
for(int j = first; j <= i-1; j++)
|
||||
sum -= this->ElementAt(i,j) * x[j];
|
||||
else
|
||||
if(sum) first = i;
|
||||
x[i] = sum;
|
||||
}
|
||||
for (int i = 3; i >= 0; i--) {
|
||||
T sum = x[i];
|
||||
for (int j = i+1; j < 4; j++)
|
||||
sum -= this->ElementAt(i, j) * x[j];
|
||||
x[i] = sum / this->ElementAt(i, i);
|
||||
}
|
||||
return x;
|
||||
res.FromEigenMatrix(mmi);
|
||||
}
|
||||
|
||||
} //namespace
|
||||
|
|
|
|||
|
|
@ -142,8 +142,7 @@ template <class T> void View<T>::GetView() {
|
|||
glGetIntegerv(GL_VIEWPORT, (GLint*)viewport);
|
||||
|
||||
matrix = proj*model;
|
||||
inverse = matrix;
|
||||
Invert(inverse);
|
||||
inverse = vcg::Inverse(matrix);
|
||||
}
|
||||
|
||||
template <class T> void View<T>::SetView(const float *_proj,
|
||||
|
|
@ -162,16 +161,13 @@ template <class T> void View<T>::SetView(const float *_proj,
|
|||
}
|
||||
|
||||
template <class T> Point3<T> View<T>::ViewPoint() const {
|
||||
Matrix44<T> mi=model;
|
||||
Invert(mi);
|
||||
return mi* Point3<T>(0, 0, 0);
|
||||
return vcg::Inverse(model)* Point3<T>(0, 0, 0);
|
||||
}
|
||||
// Note that p it is assumed to be in model coordinate.
|
||||
template <class T> Plane3<T> View<T>::ViewPlaneFromModel(const Point3<T> &p)
|
||||
{
|
||||
//compute normal, pointing away from view.
|
||||
Matrix44<T> imodel = model;
|
||||
Invert(imodel);
|
||||
Matrix44<T> imodel = vcg::Inverse(model);
|
||||
Point3<T> vp=ViewPoint();
|
||||
vcg::Point3f n = imodel * vcg::Point3<T>(0.0f, 0, -1.0f) - vp;
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue