(just fixed a warning-producing redundant assert)

This commit is contained in:
mtarini 2013-06-05 11:08:55 +00:00
parent 2c3d20ca40
commit 90cdbb6214
1 changed files with 308 additions and 308 deletions

View File

@ -35,8 +35,8 @@
namespace vcg { namespace vcg {
/* /*
Annotations: Annotations:
Opengl stores matrix in column-major order. That is, the matrix is stored as: Opengl stores matrix in column-major order. That is, the matrix is stored as:
a0 a4 a8 a12 a0 a4 a8 a12
@ -69,153 +69,153 @@ for 'column' vectors.
*/ */
/** This class represent a 4x4 matrix. T is the kind of element in the matrix. /** This class represent a 4x4 matrix. T is the kind of element in the matrix.
*/ */
template <class T> class Matrix44 { template <class T> class Matrix44 {
protected: protected:
T _a[16]; T _a[16];
public: public:
typedef T ScalarType; typedef T ScalarType;
///@{ ///@{
/** $name Constructors /** $name Constructors
* No automatic casting and default constructor is empty * No automatic casting and default constructor is empty
*/ */
Matrix44() {} Matrix44() {}
~Matrix44() {} ~Matrix44() {}
Matrix44(const Matrix44 &m); Matrix44(const Matrix44 &m);
Matrix44(const T v[]); Matrix44(const T v[]);
T &ElementAt(const int row, const int col); T &ElementAt(const int row, const int col);
T ElementAt(const int row, const int col) const; T ElementAt(const int row, const int col) const;
//T &operator[](const int i); //T &operator[](const int i);
//const T &operator[](const int i) const; //const T &operator[](const int i) const;
T *V(); T *V();
const T *V() const ; const T *V() const ;
T *operator[](const int i); T *operator[](const int i);
const T *operator[](const int i) const; const T *operator[](const int i) const;
// return a copy of the i-th column // return a copy of the i-th column
Point4<T> GetColumn4(const int& i)const{ Point4<T> GetColumn4(const int& i)const{
assert(i>=0 && i<4); assert(i>=0 && i<4);
return Point4<T>(ElementAt(0,i),ElementAt(1,i),ElementAt(2,i),ElementAt(3,i)); return Point4<T>(ElementAt(0,i),ElementAt(1,i),ElementAt(2,i),ElementAt(3,i));
//return Point4<T>(_a[i],_a[i+4],_a[i+8],_a[i+12]); //return Point4<T>(_a[i],_a[i+4],_a[i+8],_a[i+12]);
} }
Point3<T> GetColumn3(const int& i)const{ Point3<T> GetColumn3(const int& i)const{
assert(i>=0 && i<4); assert(i>=0 && i<4);
return Point3<T>(ElementAt(0,i),ElementAt(1,i),ElementAt(2,i)); return Point3<T>(ElementAt(0,i),ElementAt(1,i),ElementAt(2,i));
} }
Point4<T> GetRow4(const int& i)const{ Point4<T> GetRow4(const int& i)const{
assert(i>=0 && i<4); assert(i>=0 && i<4);
return Point4<T>(ElementAt(i,0),ElementAt(i,1),ElementAt(i,2),ElementAt(i,3)); return Point4<T>(ElementAt(i,0),ElementAt(i,1),ElementAt(i,2),ElementAt(i,3));
// return *((Point4<T>*)(&_a[i<<2])); alternativa forse + efficiente // return *((Point4<T>*)(&_a[i<<2])); alternativa forse + efficiente
} }
Point3<T> GetRow3(const int& i)const{ Point3<T> GetRow3(const int& i)const{
assert(i>=0 && i<4); assert(i>=0 && i<4);
return Point3<T>(ElementAt(i,0),ElementAt(i,1),ElementAt(i,2)); return Point3<T>(ElementAt(i,0),ElementAt(i,1),ElementAt(i,2));
// return *((Point4<T>*)(&_a[i<<2])); alternativa forse + efficiente // return *((Point4<T>*)(&_a[i<<2])); alternativa forse + efficiente
} }
Matrix44 operator+(const Matrix44 &m) const; Matrix44 operator+(const Matrix44 &m) const;
Matrix44 operator-(const Matrix44 &m) const; Matrix44 operator-(const Matrix44 &m) const;
Matrix44 operator*(const Matrix44 &m) const; Matrix44 operator*(const Matrix44 &m) const;
Point4<T> operator*(const Point4<T> &v) const; Point4<T> operator*(const Point4<T> &v) const;
bool operator==(const Matrix44 &m) const; bool operator==(const Matrix44 &m) const;
bool operator!= (const Matrix44 &m) const; bool operator!= (const Matrix44 &m) const;
Matrix44 operator-() const; Matrix44 operator-() const;
Matrix44 operator*(const T k) const; Matrix44 operator*(const T k) const;
void operator+=(const Matrix44 &m); void operator+=(const Matrix44 &m);
void operator-=(const Matrix44 &m); void operator-=(const Matrix44 &m);
void operator*=( const Matrix44 & m ); void operator*=( const Matrix44 & m );
void operator*=( const T k ); void operator*=( const T k );
template <class Matrix44Type> template <class Matrix44Type>
void ToMatrix(Matrix44Type & m) const {for(int i = 0; i < 16; i++) m.V()[i]=V()[i];} void ToMatrix(Matrix44Type & m) const {for(int i = 0; i < 16; i++) m.V()[i]=V()[i];}
void ToEulerAngles(T &alpha, T &beta, T &gamma); void ToEulerAngles(T &alpha, T &beta, T &gamma);
template <class Matrix44Type> template <class Matrix44Type>
void FromMatrix(const Matrix44Type & m){for(int i = 0; i < 16; i++) V()[i]=m.V()[i];} void FromMatrix(const Matrix44Type & m){for(int i = 0; i < 16; i++) V()[i]=m.V()[i];}
template <class EigenMatrix44Type> template <class EigenMatrix44Type>
void ToEigenMatrix(EigenMatrix44Type & m) const { void ToEigenMatrix(EigenMatrix44Type & m) const {
for(int i = 0; i < 4; i++) for(int i = 0; i < 4; i++)
for(int j = 0; j < 4; j++) for(int j = 0; j < 4; j++)
m(i,j)=(*this)[i][j]; m(i,j)=(*this)[i][j];
} }
template <class EigenMatrix44Type> template <class EigenMatrix44Type>
void FromEigenMatrix(const EigenMatrix44Type & m){ void FromEigenMatrix(const EigenMatrix44Type & m){
for(int i = 0; i < 4; i++) for(int i = 0; i < 4; i++)
for(int j = 0; j < 4; j++) for(int j = 0; j < 4; j++)
ElementAt(i,j)=m(i,j); ElementAt(i,j)=m(i,j);
} }
void FromEulerAngles(T alpha, T beta, T gamma); void FromEulerAngles(T alpha, T beta, T gamma);
void SetZero(); void SetZero();
void SetIdentity(); void SetIdentity();
void SetDiagonal(const T k); void SetDiagonal(const T k);
Matrix44 &SetScale(const T sx, const T sy, const T sz); Matrix44 &SetScale(const T sx, const T sy, const T sz);
Matrix44 &SetScale(const Point3<T> &t); Matrix44 &SetScale(const Point3<T> &t);
Matrix44<T>& SetColumn(const unsigned int ii,const Point4<T> &t); Matrix44<T>& SetColumn(const unsigned int ii,const Point4<T> &t);
Matrix44<T>& SetColumn(const unsigned int ii,const Point3<T> &t); Matrix44<T>& SetColumn(const unsigned int ii,const Point3<T> &t);
Matrix44 &SetTranslate(const Point3<T> &t); Matrix44 &SetTranslate(const Point3<T> &t);
Matrix44 &SetTranslate(const T sx, const T sy, const T sz); Matrix44 &SetTranslate(const T sx, const T sy, const T sz);
Matrix44 &SetShearXY(const T sz); Matrix44 &SetShearXY(const T sz);
Matrix44 &SetShearXZ(const T sy); Matrix44 &SetShearXZ(const T sy);
Matrix44 &SetShearYZ(const T sx); Matrix44 &SetShearYZ(const T sx);
///use radiants for angle. ///use radiants for angle.
Matrix44 &SetRotateDeg(T AngleDeg, const Point3<T> & axis); Matrix44 &SetRotateDeg(T AngleDeg, const Point3<T> & axis);
Matrix44 &SetRotateRad(T AngleRad, const Point3<T> & axis); Matrix44 &SetRotateRad(T AngleRad, const Point3<T> & axis);
T Determinant() const; T Determinant() const;
template <class Q> void Import(const Matrix44<Q> &m) { template <class Q> void Import(const Matrix44<Q> &m) {
for(int i = 0; i < 16; i++) for(int i = 0; i < 16; i++)
_a[i] = (T)(m.V()[i]); _a[i] = (T)(m.V()[i]);
} }
template <class Q> template <class Q>
static inline Matrix44 Construct( const Matrix44<Q> & b ) static inline Matrix44 Construct( const Matrix44<Q> & b )
{ {
Matrix44<T> tmp; tmp.FromMatrix(b); Matrix44<T> tmp; tmp.FromMatrix(b);
return tmp; return tmp;
} }
static inline const Matrix44 &Identity( ) static inline const Matrix44 &Identity( )
{ {
static Matrix44<T> tmp; tmp.SetIdentity(); static Matrix44<T> tmp; tmp.SetIdentity();
return tmp; return tmp;
} }
// for the transistion to eigen // for the transistion to eigen
Matrix44 transpose() const Matrix44 transpose() const
{ {
Matrix44 res = *this; Matrix44 res = *this;
Transpose(res); Transpose(res);
return res; return res;
} }
void transposeInPlace() { Transpose(*this); } void transposeInPlace() { Transpose(*this); }
void print() { void print() {
unsigned int i, j, p; unsigned int i, j, p;
for (i=0, p=0; i<4; i++, p+=4) for (i=0, p=0; i<4; i++, p+=4)
{ {
std::cout << "[\t"; std::cout << "[\t";
for (j=0; j<4; j++) for (j=0; j<4; j++)
std::cout << _a[p+j] << "\t"; std::cout << _a[p+j] << "\t";
std::cout << "]\n"; std::cout << "]\n";
} }
std::cout << "\n"; std::cout << "\n";
} }
}; };
@ -224,16 +224,16 @@ public:
/** Class for solving A * x = b. */ /** Class for solving A * x = b. */
template <class T> class LinearSolve: public Matrix44<T> { template <class T> class LinearSolve: public Matrix44<T> {
public: public:
LinearSolve(const Matrix44<T> &m); LinearSolve(const Matrix44<T> &m);
Point4<T> Solve(const Point4<T> &b); // solve A <20> x = b Point4<T> Solve(const Point4<T> &b); // solve A <20> x = b
///If you need to solve some equation you can use this function instead of Matrix44 one for speed. ///If you need to solve some equation you can use this function instead of Matrix44 one for speed.
T Determinant() const; T Determinant() const;
protected: protected:
///Holds row permutation. ///Holds row permutation.
int index[4]; //hold permutation int index[4]; //hold permutation
///Hold sign of row permutation (used for determinant sign) ///Hold sign of row permutation (used for determinant sign)
T d; T d;
bool Decompose(); bool Decompose();
}; };
/*** Postmultiply */ /*** Postmultiply */
@ -254,135 +254,135 @@ typedef Matrix44<double> Matrix44d;
template <class T> Matrix44<T>::Matrix44(const Matrix44<T> &m) { template <class T> Matrix44<T>::Matrix44(const Matrix44<T> &m) {
memcpy((T *)_a, (T *)m._a, 16 * sizeof(T)); memcpy((T *)_a, (T *)m._a, 16 * sizeof(T));
} }
template <class T> Matrix44<T>::Matrix44(const T v[]) { template <class T> Matrix44<T>::Matrix44(const T v[]) {
memcpy((T *)_a, v, 16 * sizeof(T)); memcpy((T *)_a, v, 16 * sizeof(T));
} }
template <class T> T &Matrix44<T>::ElementAt(const int row, const int col) { template <class T> T &Matrix44<T>::ElementAt(const int row, const int col) {
assert(row >= 0 && row < 4); assert(row >= 0 && row < 4);
assert(col >= 0 && col < 4); assert(col >= 0 && col < 4);
return _a[(row<<2) + col]; return _a[(row<<2) + col];
} }
template <class T> T Matrix44<T>::ElementAt(const int row, const int col) const { template <class T> T Matrix44<T>::ElementAt(const int row, const int col) const {
assert(row >= 0 && row < 4); assert(row >= 0 && row < 4);
assert(col >= 0 && col < 4); assert(col >= 0 && col < 4);
return _a[(row<<2) + col]; return _a[(row<<2) + col];
} }
//template <class T> T &Matrix44<T>::operator[](const int i) { //template <class T> T &Matrix44<T>::operator[](const int i) {
// assert(i >= 0 && i < 16); // assert(i >= 0 && i < 16);
// return ((T *)_a)[i]; // return ((T *)_a)[i];
//} //}
// //
//template <class T> const T &Matrix44<T>::operator[](const int i) const { //template <class T> const T &Matrix44<T>::operator[](const int i) const {
// assert(i >= 0 && i < 16); // assert(i >= 0 && i < 16);
// return ((T *)_a)[i]; // return ((T *)_a)[i];
//} //}
template <class T> T *Matrix44<T>::operator[](const int i) { template <class T> T *Matrix44<T>::operator[](const int i) {
assert(i >= 0 && i < 4); assert(i >= 0 && i < 4);
return _a+i*4; return _a+i*4;
} }
template <class T> const T *Matrix44<T>::operator[](const int i) const { template <class T> const T *Matrix44<T>::operator[](const int i) const {
assert(i >= 0 && i < 4); assert(i >= 0 && i < 4);
return _a+i*4; return _a+i*4;
} }
template <class T> T *Matrix44<T>::V() { return _a;} template <class T> T *Matrix44<T>::V() { return _a;}
template <class T> const T *Matrix44<T>::V() const { return _a;} template <class T> const T *Matrix44<T>::V() const { return _a;}
template <class T> Matrix44<T> Matrix44<T>::operator+(const Matrix44 &m) const { template <class T> Matrix44<T> Matrix44<T>::operator+(const Matrix44 &m) const {
Matrix44<T> ret; Matrix44<T> ret;
for(int i = 0; i < 16; i++) for(int i = 0; i < 16; i++)
ret.V()[i] = V()[i] + m.V()[i]; ret.V()[i] = V()[i] + m.V()[i];
return ret; return ret;
} }
template <class T> Matrix44<T> Matrix44<T>::operator-(const Matrix44 &m) const { template <class T> Matrix44<T> Matrix44<T>::operator-(const Matrix44 &m) const {
Matrix44<T> ret; Matrix44<T> ret;
for(int i = 0; i < 16; i++) for(int i = 0; i < 16; i++)
ret.V()[i] = V()[i] - m.V()[i]; ret.V()[i] = V()[i] - m.V()[i];
return ret; return ret;
} }
template <class T> Matrix44<T> Matrix44<T>::operator*(const Matrix44 &m) const { template <class T> Matrix44<T> Matrix44<T>::operator*(const Matrix44 &m) const {
Matrix44 ret; Matrix44 ret;
for(int i = 0; i < 4; i++) for(int i = 0; i < 4; i++)
for(int j = 0; j < 4; j++) { for(int j = 0; j < 4; j++) {
T t = 0.0; T t = 0.0;
for(int k = 0; k < 4; k++) for(int k = 0; k < 4; k++)
t += ElementAt(i, k) * m.ElementAt(k, j); t += ElementAt(i, k) * m.ElementAt(k, j);
ret.ElementAt(i, j) = t; ret.ElementAt(i, j) = t;
} }
return ret; return ret;
} }
template <class T> Point4<T> Matrix44<T>::operator*(const Point4<T> &v) const { template <class T> Point4<T> Matrix44<T>::operator*(const Point4<T> &v) const {
Point4<T> ret; Point4<T> ret;
for(int i = 0; i < 4; i++){ for(int i = 0; i < 4; i++){
T t = 0.0; T t = 0.0;
for(int k = 0; k < 4; k++) for(int k = 0; k < 4; k++)
t += ElementAt(i,k) * v[k]; t += ElementAt(i,k) * v[k];
ret[i] = t; ret[i] = t;
} }
return ret; return ret;
} }
template <class T> bool Matrix44<T>::operator==(const Matrix44 &m) const { template <class T> bool Matrix44<T>::operator==(const Matrix44 &m) const {
for(int i = 0; i < 4; ++i) for(int i = 0; i < 4; ++i)
for(int j = 0; j < 4; ++j) for(int j = 0; j < 4; ++j)
if(ElementAt(i,j) != m.ElementAt(i,j)) if(ElementAt(i,j) != m.ElementAt(i,j))
return false; return false;
return true; return true;
} }
template <class T> bool Matrix44<T>::operator!=(const Matrix44 &m) const { template <class T> bool Matrix44<T>::operator!=(const Matrix44 &m) const {
for(int i = 0; i < 4; ++i) for(int i = 0; i < 4; ++i)
for(int j = 0; j < 4; ++j) for(int j = 0; j < 4; ++j)
if(ElementAt(i,j) != m.ElementAt(i,j)) if(ElementAt(i,j) != m.ElementAt(i,j))
return true; return true;
return false; return false;
} }
template <class T> Matrix44<T> Matrix44<T>::operator-() const { template <class T> Matrix44<T> Matrix44<T>::operator-() const {
Matrix44<T> res; Matrix44<T> res;
for(int i = 0; i < 16; i++) for(int i = 0; i < 16; i++)
res.V()[i] = -V()[i]; res.V()[i] = -V()[i];
return res; return res;
} }
template <class T> Matrix44<T> Matrix44<T>::operator*(const T k) const { template <class T> Matrix44<T> Matrix44<T>::operator*(const T k) const {
Matrix44<T> res; Matrix44<T> res;
for(int i = 0; i < 16; i++) for(int i = 0; i < 16; i++)
res.V()[i] =V()[i] * k; res.V()[i] =V()[i] * k;
return res; return res;
} }
template <class T> void Matrix44<T>::operator+=(const Matrix44 &m) { template <class T> void Matrix44<T>::operator+=(const Matrix44 &m) {
for(int i = 0; i < 16; i++) for(int i = 0; i < 16; i++)
V()[i] += m.V()[i]; V()[i] += m.V()[i];
} }
template <class T> void Matrix44<T>::operator-=(const Matrix44 &m) { template <class T> void Matrix44<T>::operator-=(const Matrix44 &m) {
for(int i = 0; i < 16; i++) for(int i = 0; i < 16; i++)
V()[i] -= m.V()[i]; V()[i] -= m.V()[i];
} }
template <class T> void Matrix44<T>::operator*=( const Matrix44 & m ) { template <class T> void Matrix44<T>::operator*=( const Matrix44 & m ) {
*this = *this *m; *this = *this *m;
} }
template < class PointType , class T > void operator*=( std::vector<PointType> &vert, const Matrix44<T> & m ) { template < class PointType , class T > void operator*=( std::vector<PointType> &vert, const Matrix44<T> & m ) {
typename std::vector<PointType>::iterator ii; typename std::vector<PointType>::iterator ii;
for(ii=vert.begin();ii!=vert.end();++ii) for(ii=vert.begin();ii!=vert.end();++ii)
(*ii).P()=m * (*ii).P(); (*ii).P()=m * (*ii).P();
} }
template <class T> void Matrix44<T>::operator*=( const T k ) { template <class T> void Matrix44<T>::operator*=( const T k ) {
for(int i = 0; i < 16; i++) for(int i = 0; i < 16; i++)
_a[i] *= k; _a[i] *= k;
} }
template <class T> template <class T>
@ -390,7 +390,7 @@ void Matrix44<T>::ToEulerAngles(T &alpha, T &beta, T &gamma)
{ {
alpha = atan2(ElementAt(1,2), ElementAt(2,2)); alpha = atan2(ElementAt(1,2), ElementAt(2,2));
beta = asin(-ElementAt(0,2)); beta = asin(-ElementAt(0,2));
gamma = atan2(ElementAt(0,1), ElementAt(0,0)); gamma = atan2(ElementAt(0,1), ElementAt(0,0));
} }
template <class T> template <class T>
@ -421,48 +421,48 @@ void Matrix44<T>::FromEulerAngles(T alpha, T beta, T gamma)
} }
template <class T> void Matrix44<T>::SetZero() { template <class T> void Matrix44<T>::SetZero() {
memset((T *)_a, 0, 16 * sizeof(T)); memset((T *)_a, 0, 16 * sizeof(T));
} }
template <class T> void Matrix44<T>::SetIdentity() { template <class T> void Matrix44<T>::SetIdentity() {
SetDiagonal(1); SetDiagonal(1);
} }
template <class T> void Matrix44<T>::SetDiagonal(const T k) { template <class T> void Matrix44<T>::SetDiagonal(const T k) {
SetZero(); SetZero();
ElementAt(0, 0) = k; ElementAt(0, 0) = k;
ElementAt(1, 1) = k; ElementAt(1, 1) = k;
ElementAt(2, 2) = k; ElementAt(2, 2) = k;
ElementAt(3, 3) = 1; ElementAt(3, 3) = 1;
} }
template <class T> Matrix44<T> &Matrix44<T>::SetScale(const Point3<T> &t) { template <class T> Matrix44<T> &Matrix44<T>::SetScale(const Point3<T> &t) {
SetScale(t[0], t[1], t[2]); SetScale(t[0], t[1], t[2]);
return *this; return *this;
} }
template <class T> Matrix44<T> &Matrix44<T>::SetScale(const T sx, const T sy, const T sz) { template <class T> Matrix44<T> &Matrix44<T>::SetScale(const T sx, const T sy, const T sz) {
SetZero(); SetZero();
ElementAt(0, 0) = sx; ElementAt(0, 0) = sx;
ElementAt(1, 1) = sy; ElementAt(1, 1) = sy;
ElementAt(2, 2) = sz; ElementAt(2, 2) = sz;
ElementAt(3, 3) = 1; ElementAt(3, 3) = 1;
return *this; return *this;
} }
template <class T> Matrix44<T> &Matrix44<T>::SetTranslate(const Point3<T> &t) { template <class T> Matrix44<T> &Matrix44<T>::SetTranslate(const Point3<T> &t) {
SetTranslate(t[0], t[1], t[2]); SetTranslate(t[0], t[1], t[2]);
return *this; return *this;
} }
template <class T> Matrix44<T> &Matrix44<T>::SetTranslate(const T tx, const T ty, const T tz) { template <class T> Matrix44<T> &Matrix44<T>::SetTranslate(const T tx, const T ty, const T tz) {
SetIdentity(); SetIdentity();
ElementAt(0, 3) = tx; ElementAt(0, 3) = tx;
ElementAt(1, 3) = ty; ElementAt(1, 3) = ty;
ElementAt(2, 3) = tz; ElementAt(2, 3) = tz;
return *this; return *this;
} }
template <class T> Matrix44<T> &Matrix44<T>::SetColumn(const unsigned int ii,const Point3<T> &t) { template <class T> Matrix44<T> &Matrix44<T>::SetColumn(const unsigned int ii,const Point3<T> &t) {
assert((ii >= 0) && (ii < 4)); assert( ii < 4 );
ElementAt(0, ii) = t.X(); ElementAt(0, ii) = t.X();
ElementAt(1, ii) = t.Y(); ElementAt(1, ii) = t.Y();
ElementAt(2, ii) = t.Z(); ElementAt(2, ii) = t.Z();
@ -470,53 +470,53 @@ template <class T> Matrix44<T> &Matrix44<T>::SetColumn(const unsigned int ii,con
} }
template <class T> Matrix44<T> &Matrix44<T>::SetColumn(const unsigned int ii,const Point4<T> &t) { template <class T> Matrix44<T> &Matrix44<T>::SetColumn(const unsigned int ii,const Point4<T> &t) {
assert((ii < 4)); assert( ii < 4 );
ElementAt(0, ii) = t[0]; ElementAt(0, ii) = t[0];
ElementAt(1, ii) = t[1]; ElementAt(1, ii) = t[1];
ElementAt(2, ii) = t[2]; ElementAt(2, ii) = t[2];
ElementAt(3, ii) = t[3]; ElementAt(3, ii) = t[3];
return *this; return *this;
} }
template <class T> Matrix44<T> &Matrix44<T>::SetRotateDeg(T AngleDeg, const Point3<T> & axis) { template <class T> Matrix44<T> &Matrix44<T>::SetRotateDeg(T AngleDeg, const Point3<T> & axis) {
return SetRotateRad(math::ToRad(AngleDeg),axis); return SetRotateRad(math::ToRad(AngleDeg),axis);
} }
template <class T> Matrix44<T> &Matrix44<T>::SetRotateRad(T AngleRad, const Point3<T> & axis) { template <class T> Matrix44<T> &Matrix44<T>::SetRotateRad(T AngleRad, const Point3<T> & axis) {
//angle = angle*(T)3.14159265358979323846/180; e' in radianti! //angle = angle*(T)3.14159265358979323846/180; e' in radianti!
T c = math::Cos(AngleRad); T c = math::Cos(AngleRad);
T s = math::Sin(AngleRad); T s = math::Sin(AngleRad);
T q = 1-c; T q = 1-c;
Point3<T> t = axis; Point3<T> t = axis;
t.Normalize(); t.Normalize();
ElementAt(0,0) = t[0]*t[0]*q + c; ElementAt(0,0) = t[0]*t[0]*q + c;
ElementAt(0,1) = t[0]*t[1]*q - t[2]*s; ElementAt(0,1) = t[0]*t[1]*q - t[2]*s;
ElementAt(0,2) = t[0]*t[2]*q + t[1]*s; ElementAt(0,2) = t[0]*t[2]*q + t[1]*s;
ElementAt(0,3) = 0; ElementAt(0,3) = 0;
ElementAt(1,0) = t[1]*t[0]*q + t[2]*s; ElementAt(1,0) = t[1]*t[0]*q + t[2]*s;
ElementAt(1,1) = t[1]*t[1]*q + c; ElementAt(1,1) = t[1]*t[1]*q + c;
ElementAt(1,2) = t[1]*t[2]*q - t[0]*s; ElementAt(1,2) = t[1]*t[2]*q - t[0]*s;
ElementAt(1,3) = 0; ElementAt(1,3) = 0;
ElementAt(2,0) = t[2]*t[0]*q -t[1]*s; ElementAt(2,0) = t[2]*t[0]*q -t[1]*s;
ElementAt(2,1) = t[2]*t[1]*q +t[0]*s; ElementAt(2,1) = t[2]*t[1]*q +t[0]*s;
ElementAt(2,2) = t[2]*t[2]*q +c; ElementAt(2,2) = t[2]*t[2]*q +c;
ElementAt(2,3) = 0; ElementAt(2,3) = 0;
ElementAt(3,0) = 0; ElementAt(3,0) = 0;
ElementAt(3,1) = 0; ElementAt(3,1) = 0;
ElementAt(3,2) = 0; ElementAt(3,2) = 0;
ElementAt(3,3) = 1; ElementAt(3,3) = 1;
return *this; return *this;
} }
/* /*
Given a non singular, non projective matrix (e.g. with the last row equal to [0,0,0,1] ) 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 This procedure decompose it in a sequence of
Scale,Shear,Rotation e Translation - Scale,Shear,Rotation e Translation
- ScaleV and Tranv are obiviously scaling and translation. - ScaleV and Tranv are obiviously scaling and translation.
- ShearV contains three scalars with, respectively - ShearV contains three scalars with, respectively,
ShearXY, ShearXZ e ShearYZ ShearXY, ShearXZ and ShearYZ
- RotateV contains the rotations (in degree!) around the x,y,z axis - RotateV contains the rotations (in degree!) around the x,y,z axis
The input matrix is modified leaving inside it a simple roto translation. The input matrix is modified leaving inside it a simple roto translation.
@ -535,14 +535,14 @@ double srv() { return (double(rand()%40)-20)/2.0; } // small random value
Matrix44d Scl; Scl.SetScale(ScV); Matrix44d Scl; Scl.SetScale(ScV);
Matrix44d Sxy; Sxy.SetShearXY(ShV[0]); Matrix44d Sxy; Sxy.SetShearXY(ShV[0]);
Matrix44d Sxz; Sxz.SetShearXZ(ShV[1]); Matrix44d Sxz; Sxz.SetShearXZ(ShV[1]);
Matrix44d Syz; Syz.SetShearYZ(ShV[2]); Matrix44d Syz; Syz.SetShearYZ(ShV[2]);
Matrix44d Rtx; Rtx.SetRotate(math::ToRad(RtV[0]),Point3d(1,0,0)); Matrix44d Rtx; Rtx.SetRotate(math::ToRad(RtV[0]),Point3d(1,0,0));
Matrix44d Rty; Rty.SetRotate(math::ToRad(RtV[1]),Point3d(0,1,0)); Matrix44d Rty; Rty.SetRotate(math::ToRad(RtV[1]),Point3d(0,1,0));
Matrix44d Rtz; Rtz.SetRotate(math::ToRad(RtV[2]),Point3d(0,0,1)); Matrix44d Rtz; Rtz.SetRotate(math::ToRad(RtV[2]),Point3d(0,0,1));
Matrix44d Trn; Trn.SetTranslate(TrV); Matrix44d Trn; Trn.SetTranslate(TrV);
Matrix44d StartM = Trn * Rtx*Rty*Rtz * Syz*Sxz*Sxy *Scl; Matrix44d StartM = Trn * Rtx*Rty*Rtz * Syz*Sxz*Sxy *Scl;
Matrix44d ResultM=StartM; Matrix44d ResultM=StartM;
Decompose(ResultM,ScVOut,ShVOut,RtVOut,TrVOut); Decompose(ResultM,ScVOut,ShVOut,RtVOut,TrVOut);
@ -556,8 +556,9 @@ double srv() { return (double(rand()%40)-20)/2.0; } // small random value
Trn.SetTranslate(TrVOut); Trn.SetTranslate(TrVOut);
// Now Rebuild is equal to StartM // Now Rebuild is equal to StartM
Matrix44d RebuildM = Trn * Rtx*Rty*Rtz * Syz*Sxz*Sxy * Scl ; Matrix44d RebuildM = Trn * Rtx*Rty*Rtz * Syz*Sxz*Sxy * Scl ;
*/ */
template <class T> template <class T>
bool Decompose(Matrix44<T> &M, Point3<T> &ScaleV, Point3<T> &ShearV, Point3<T> &RotV,Point3<T> &TranV) bool Decompose(Matrix44<T> &M, Point3<T> &ScaleV, Point3<T> &ShearV, Point3<T> &RotV,Point3<T> &TranV)
{ {
@ -565,9 +566,8 @@ bool Decompose(Matrix44<T> &M, Point3<T> &ScaleV, Point3<T> &ShearV, Point3<T> &
return false; return false;
if(math::Abs(M.Determinant())<1e-10) return false; // matrix should be at least invertible... if(math::Abs(M.Determinant())<1e-10) return false; // matrix should be at least invertible...
// First Step recover the traslation // First Step recover the traslation
TranV=M.GetColumn3(3); TranV=M.GetColumn3(3);
// Second Step Recover Scale and Shearing interleaved // Second Step Recover Scale and Shearing interleaved
ScaleV[0]=Norm(M.GetColumn3(0)); ScaleV[0]=Norm(M.GetColumn3(0));
@ -577,10 +577,10 @@ bool Decompose(Matrix44<T> &M, Point3<T> &ScaleV, Point3<T> &ShearV, Point3<T> &
ShearV[0]=R[0]*M.GetColumn3(1); // xy shearing ShearV[0]=R[0]*M.GetColumn3(1); // xy shearing
R[1]= M.GetColumn3(1)-R[0]*ShearV[0]; R[1]= M.GetColumn3(1)-R[0]*ShearV[0];
assert(math::Abs(R[1]*R[0])<1e-10); assert(math::Abs(R[1]*R[0])<1e-10);
ScaleV[1]=Norm(R[1]); // y scaling ScaleV[1]=Norm(R[1]); // y scaling
R[1]=R[1]/ScaleV[1]; R[1]=R[1]/ScaleV[1];
ShearV[0]=ShearV[0]/ScaleV[1]; ShearV[0]=ShearV[0]/ScaleV[1];
ShearV[1]=R[0]*M.GetColumn3(2); // xz shearing ShearV[1]=R[0]*M.GetColumn3(2); // xz shearing
R[2]= M.GetColumn3(2)-R[0]*ShearV[1]; R[2]= M.GetColumn3(2)-R[0]*ShearV[1];
@ -598,81 +598,81 @@ bool Decompose(Matrix44<T> &M, Point3<T> &ScaleV, Point3<T> &ShearV, Point3<T> &
ShearV[2]=R[1]*M.GetColumn3(2); // yz shearing ShearV[2]=R[1]*M.GetColumn3(2); // yz shearing
ShearV[2]=ShearV[2]/ScaleV[2]; ShearV[2]=ShearV[2]/ScaleV[2];
int i,j; int i,j;
for(i=0;i<3;++i) for(i=0;i<3;++i)
for(j=0;j<3;++j) for(j=0;j<3;++j)
M[i][j]=R[j][i]; M[i][j]=R[j][i];
// Third and last step: Recover the rotation // Third and last step: Recover the rotation
//now the matrix should be a pure rotation matrix so its determinant is +-1 //now the matrix should be a pure rotation matrix so its determinant is +-1
double det=M.Determinant(); double det=M.Determinant();
if(math::Abs(det)<1e-10) return false; // matrix should be at least invertible... if(math::Abs(det)<1e-10) return false; // matrix should be at least invertible...
assert(math::Abs(math::Abs(det)-1.0)<1e-10); // it should be +-1... assert(math::Abs(math::Abs(det)-1.0)<1e-10); // it should be +-1...
if(det<0) { if(det<0) {
ScaleV *= -1; ScaleV *= -1;
M *= -1; M *= -1;
} }
double alpha,beta,gamma; // rotations around the x,y and z axis double alpha,beta,gamma; // rotations around the x,y and z axis
beta=asin( M[0][2]); beta=asin( M[0][2]);
double cosbeta=cos(beta); double cosbeta=cos(beta);
if(math::Abs(cosbeta) > 1e-5) if(math::Abs(cosbeta) > 1e-5)
{ {
alpha=asin(-M[1][2]/cosbeta); alpha=asin(-M[1][2]/cosbeta);
if((M[2][2]/cosbeta) < 0 ) alpha=M_PI-alpha; if((M[2][2]/cosbeta) < 0 ) alpha=M_PI-alpha;
gamma=asin(-M[0][1]/cosbeta); gamma=asin(-M[0][1]/cosbeta);
if((M[0][0]/cosbeta)<0) gamma = M_PI-gamma; if((M[0][0]/cosbeta)<0) gamma = M_PI-gamma;
} }
else else
{ {
alpha=asin(-M[1][0]); alpha=asin(-M[1][0]);
if(M[1][1]<0) alpha=M_PI-alpha; if(M[1][1]<0) alpha=M_PI-alpha;
gamma=0; gamma=0;
} }
RotV[0]=math::ToDeg(alpha); RotV[0]=math::ToDeg(alpha);
RotV[1]=math::ToDeg(beta); RotV[1]=math::ToDeg(beta);
RotV[2]=math::ToDeg(gamma); RotV[2]=math::ToDeg(gamma);
return true; return true;
} }
template <class T> T Matrix44<T>::Determinant() const { template <class T> T Matrix44<T>::Determinant() const {
Eigen::Matrix4d mm; Eigen::Matrix4d mm;
this->ToEigenMatrix(mm); this->ToEigenMatrix(mm);
return mm.determinant(); return mm.determinant();
} }
template <class T> Point3<T> operator*(const Matrix44<T> &m, const Point3<T> &p) { template <class T> Point3<T> operator*(const Matrix44<T> &m, const Point3<T> &p) {
T w; T w;
Point3<T> s; Point3<T> s;
s[0] = m.ElementAt(0, 0)*p[0] + m.ElementAt(0, 1)*p[1] + m.ElementAt(0, 2)*p[2] + m.ElementAt(0, 3); s[0] = m.ElementAt(0, 0)*p[0] + m.ElementAt(0, 1)*p[1] + m.ElementAt(0, 2)*p[2] + m.ElementAt(0, 3);
s[1] = m.ElementAt(1, 0)*p[0] + m.ElementAt(1, 1)*p[1] + m.ElementAt(1, 2)*p[2] + m.ElementAt(1, 3); s[1] = m.ElementAt(1, 0)*p[0] + m.ElementAt(1, 1)*p[1] + m.ElementAt(1, 2)*p[2] + m.ElementAt(1, 3);
s[2] = m.ElementAt(2, 0)*p[0] + m.ElementAt(2, 1)*p[1] + m.ElementAt(2, 2)*p[2] + m.ElementAt(2, 3); s[2] = m.ElementAt(2, 0)*p[0] + m.ElementAt(2, 1)*p[1] + m.ElementAt(2, 2)*p[2] + m.ElementAt(2, 3);
w = m.ElementAt(3, 0)*p[0] + m.ElementAt(3, 1)*p[1] + m.ElementAt(3, 2)*p[2] + m.ElementAt(3, 3); w = m.ElementAt(3, 0)*p[0] + m.ElementAt(3, 1)*p[1] + m.ElementAt(3, 2)*p[2] + m.ElementAt(3, 3);
if(w!= 0) s /= w; if(w!= 0) s /= w;
return s; return s;
} }
template <class T> Matrix44<T> &Transpose(Matrix44<T> &m) { template <class T> Matrix44<T> &Transpose(Matrix44<T> &m) {
for(int i = 1; i < 4; i++) for(int i = 1; i < 4; i++)
for(int j = 0; j < i; j++) { for(int j = 0; j < i; j++) {
math::Swap(m.ElementAt(i, j), m.ElementAt(j, i)); math::Swap(m.ElementAt(i, j), m.ElementAt(j, i));
} }
return m; return m;
} }
template <class T> Matrix44<T> Inverse(const Matrix44<T> &m) { template <class T> Matrix44<T> Inverse(const Matrix44<T> &m) {
Eigen::Matrix4d mm,mmi; Eigen::Matrix4d mm,mmi;
m.ToEigenMatrix(mm); m.ToEigenMatrix(mm);
mmi=mm.inverse(); mmi=mm.inverse();
Matrix44<T> res; Matrix44<T> res;
res.FromEigenMatrix(mmi); res.FromEigenMatrix(mmi);
return res; return res;
} }
} //namespace } //namespace