From 07f2e976ea7b2396675b4b1f2160884e21dd2010 Mon Sep 17 00:00:00 2001 From: cnr-isti-vclab Date: Wed, 29 Oct 2008 00:05:44 +0000 Subject: [PATCH] * change all remaining Transpose to transpose, * update the gl/math wrappers to make them more Eigen friendly (and remove the useless, and not used, and somehow dangerous *Direct and *E functions) * add automatic reinterpret_casting from Eigen::Matrix to vcg specialized types --- apps/ptx2ply/ptx2ply.cpp | 36 ++-- apps/test/quaternion/quat_test.cpp | 29 ++- vcg/math/deprecated_matrix44.h | 287 +++++++++++++++-------------- vcg/math/eigen.h | 57 ++++++ vcg/math/eigen_matrix_addons.h | 54 +++--- vcg/math/matrix.h | 5 +- vcg/math/matrix33.h | 20 +- vcg/math/matrix44.h | 65 +------ vcg/math/shot.h | 40 ++-- vcg/space/point2.h | 2 + vcg/space/point3.h | 9 +- vcg/space/point4.h | 4 +- wrap/gl/deprecated_math.h | 144 +++++++++++++++ wrap/gl/math.h | 90 ++++----- wrap/gui/trackball.cpp | 68 +++---- 15 files changed, 533 insertions(+), 377 deletions(-) create mode 100644 wrap/gl/deprecated_math.h diff --git a/apps/ptx2ply/ptx2ply.cpp b/apps/ptx2ply/ptx2ply.cpp index d4dc0d02..9bffdc9b 100644 --- a/apps/ptx2ply/ptx2ply.cpp +++ b/apps/ptx2ply/ptx2ply.cpp @@ -43,7 +43,7 @@ using namespace tri; class MyFace; class MyEdgeC; class MyFaceC; - + class MyVertexC:public VertexVCVN{}; class MyFaceC :public FaceFN{}; class MyMeshC: public tri::TriMesh< std::vector, std::vector >{}; @@ -148,11 +148,11 @@ int readmesh(FILE* fp) if(hascolor && savecolor) { - viC = Allocator::AddVertices(currentmeshC,(rownum*colnum)); + viC = Allocator::AddVertices(currentmeshC,(rownum*colnum)); } else { - vi = Allocator::AddVertices(currentmesh,(rownum*colnum)); + vi = Allocator::AddVertices(currentmesh,(rownum*colnum)); } // parse the first line.... @@ -360,7 +360,7 @@ int readmesh(FILE* fp) if(hascolor && savecolor) printf("\nV: %8i F: %8i \n", currentmeshC.vn, currentmeshC.fn); - else + else printf("\nV: %8i F: %8i \n", currentmesh.vn, currentmesh.fn); } @@ -419,7 +419,7 @@ int readmesh(FILE* fp) if(hascolor && savecolor) printf("V: %8i F: %8i \n", currentmeshC.vn, currentmeshC.fn); - else + else printf("V: %8i F: %8i \n", currentmesh.vn, currentmesh.fn); @@ -464,7 +464,7 @@ int readmesh(FILE* fp) } } - Transpose(currtrasf); + currtrasf = currtrasf.transpose().eval(); // apply tranformation if(hascolor && savecolor) @@ -496,7 +496,7 @@ int readmesh(FILE* fp) if(! onlypoints) int unref = tri::Clean::RemoveUnreferencedVertex(currentmesh); } - + return 0; } @@ -568,7 +568,7 @@ void dounpack(FILE* fp) } -// skip a mesh +// skip a mesh int skipmesh(FILE* fp) { int colnum; @@ -595,7 +595,7 @@ int skipmesh(FILE* fp) fread(&linebuf,1,1,fp); while(linebuf != '\n') fread(&linebuf,1,1,fp); - } + } return 0; } @@ -643,17 +643,17 @@ void parseparams(int argn, char** argvect) dumpit = true; printf("dumping # %i chars\n",todump); } - if(argvect[pit][1] == 'u') // unpack the file in different + if(argvect[pit][1] == 'u') // unpack the file in different { unpack = true; printf("UNPACKING \n"); } - if(argvect[pit][1] == 'c') // save color if present + if(argvect[pit][1] == 'c') // save color if present { savecolor = true; printf("saving color \n"); } - if(argvect[pit][1] == 'k') // keep all + if(argvect[pit][1] == 'k') // keep all { saveall = true; printf("keeping all elements \n"); @@ -673,7 +673,7 @@ void parseparams(int argn, char** argvect) switchside = true; printf("swapped triangulation \n"); } - + } } @@ -750,19 +750,19 @@ int main(int argc, char *argv[]) dumpit = false; unpack = false; savecolor = false; - saveall = false; + saveall = false; flipfaces = false; onlypoints = false; - switchside = false; + switchside = false; if(argc < 2) - printhelp(); + printhelp(); //-- parseparams(argc, argv); - + strcpy(modelname,argv[1]); modelname[strlen(argv[1])-4] = '\0'; @@ -794,7 +794,7 @@ int main(int argc, char *argv[]) fclose(fp); exit(0); } - + printf("reading "); readmesh(fp); diff --git a/apps/test/quaternion/quat_test.cpp b/apps/test/quaternion/quat_test.cpp index a0666295..6fe2e1e7 100644 --- a/apps/test/quaternion/quat_test.cpp +++ b/apps/test/quaternion/quat_test.cpp @@ -47,13 +47,13 @@ ostream &operator<<(ostream &o, Matrix44f &m) { return o; } -bool verify(Quaternionf &q){ +bool verify(Quaternionf &q){ cout << "Quaternion: " << q << endl; Matrix33f m; q.ToMatrix(m); cout << "To Matrix: " << m << endl; - cout << "Row norms: " << m.GetRow(0).Norm() << " " - << m.GetRow(1).Norm() << " " + cout << "Row norms: " << m.GetRow(0).Norm() << " " + << m.GetRow(1).Norm() << " " << m.GetRow(2).Norm() << endl; Point3f p(3, 4, 5); Point3f qp = q.Rotate(p); @@ -68,14 +68,13 @@ bool verify(Quaternionf &q){ bool verify(Matrix33f &m) { cout << "Matrix: " << m << endl; cout << "Det: " << m.Determinant() << endl; - cout << "Row norms: " << m.GetRow(0).Norm() << " " - << m.GetRow(1).Norm() << " " + cout << "Row norms: " << m.GetRow(0).Norm() << " " + << m.GetRow(1).Norm() << " " << m.GetRow(2).Norm() << endl; - cout << "Column norms: " << m.GetColumn(0).Norm() << " " - << m.GetColumn(1).Norm() << " " + cout << "Column norms: " << m.GetColumn(0).Norm() << " " + << m.GetColumn(1).Norm() << " " << m.GetColumn(2).Norm() << endl; - Matrix33f im = m; - im.Transpose(); + Matrix33f im = m.transpose() * m; im = im*m; cout << "Check ortonormality: " << im << endl; Quaternionf q; @@ -87,7 +86,7 @@ bool verify(Matrix33f &m) { cout << "Norm: " << axis.SquaredNorm() + q[0]*q[0] << endl; axis.Normalize(); cout << "angle: " << 2*acos(q[0]) << " Axis: " << axis << endl; - + Point3f p(3, 4, 5); Point3f qp = q.Rotate(p); Point3f mp = m*p; @@ -112,8 +111,8 @@ int main() { verify(q); Matrix33f m; - m[0][0] = 0.70145; m[0][1] = 0.372035; m[0][2] = 0.607913; - m[1][0] = -0.628023; m[1][1] = 0.725922; m[1][2] = 0.2804; + m[0][0] = 0.70145; m[0][1] = 0.372035; m[0][2] = 0.607913; + m[1][0] = -0.628023; m[1][1] = 0.725922; m[1][2] = 0.2804; m[2][0] = -0.336978; m[2][1] = -0.57847; m[2][2] = 0.742845; cout << "verify matrix: " << endl; @@ -129,7 +128,7 @@ int main() { //Matrix33f m; cout << "matrix: " << m << endl; - + Point3f p(3, 4, 5); //test this matrix: cout << "norms: " << m.GetRow(0).Norm() << " " << m.GetRow(1).Norm() << " " << m.GetRow(2).Norm() << endl; @@ -143,14 +142,14 @@ Point3f p(3, 4, 5); } cout << endl; } - cout << endl; + cout << endl; cout << "Point: " << p[0] << " " << p[1] << " " << p[2] << endl; Point3f r = q.Rotate(p); cout << "Rotated by q: " << r[0] << " " << r[1] << " " << r[2] << endl; r = m*p; - cout << "Rotated by m: " << r[0] << " " << r[1] << " " << r[2] << endl; + cout << "Rotated by m: " << r[0] << " " << r[1] << " " << r[2] << endl; q.FromMatrix(m); cout << "quaternion: " << q[0] << " " << q[1] << " " << q[2] << " " << q[3] < #include #include +#include namespace vcg { @@ -120,13 +121,13 @@ Opengl stores matrix in column-major order. That is, the matrix is stored as: a2 a6 a10 a14 a3 a7 a11 a15 - Usually in opengl (see opengl specs) vectors are 'column' vectors + Usually in opengl (see opengl specs) vectors are 'column' vectors so usually matrix are PRE-multiplied for a vector. - So the command glTranslate generate a matrix that + So the command glTranslate generate a matrix that is ready to be premultipled for a vector: 1 0 0 tx - 0 1 0 ty + 0 1 0 ty 0 0 1 tz 0 0 0 1 @@ -138,7 +139,7 @@ Matrix44 stores matrix in row-major order i.e. a12 a13 a14 a15 So for the use of that matrix in opengl with their supposed meaning you have to transpose them before feeding to glMultMatrix. -This mechanism is hidden by the templated function defined in wrap/gl/math.h; +This mechanism is hidden by the templated function defined in wrap/gl/math.h; If your machine has the ARB_transpose_matrix extension it will use the appropriate; The various gl-like command SetRotate, SetTranslate assume that you are making matrix for 'column' vectors. @@ -150,7 +151,7 @@ class Matrix44Diag:public Point4{ public: /** @name Matrix33 Class Matrix33Diag. - This is the class for definition of a diagonal matrix 4x4. + 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(p0,p1,p2,p3){}; @@ -160,11 +161,11 @@ public: /** This class represent a 4x4 matrix. T is the kind of element in the matrix. */ -template class Matrix44 { +template class Matrix44 { protected: T _a[16]; -public: +public: typedef T ScalarType; ///@{ @@ -172,7 +173,7 @@ public: /** $name Constructors * No automatic casting and default constructor is empty */ - Matrix44() {}; + Matrix44() {}; ~Matrix44() {}; Matrix44(const Matrix44 &m); Matrix44(const T v[]); @@ -195,7 +196,7 @@ public: //const T &operator[](const int i) const; T *V(); const T *V() const ; - + T *operator[](const int i); const T *operator[](const int i) const; @@ -216,7 +217,7 @@ public: return Point4(ElementAt(i,0),ElementAt(i,1),ElementAt(i,2),ElementAt(i,3)); // return *((Point4*)(&_a[i<<2])); alternativa forse + efficiente } - + Point3 GetRow3(const int& i)const{ assert(i>=0 && i<4); return Point3(ElementAt(i,0),ElementAt(i,1),ElementAt(i,2)); @@ -227,18 +228,18 @@ public: Matrix44 operator-(const Matrix44 &m) const; Matrix44 operator*(const Matrix44 &m) const; Matrix44 operator*(const Matrix44Diag &m) const; - Point4 operator*(const Point4 &v) const; + Point4 operator*(const Point4 &v) 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; 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 ); - + template void ToMatrix(Matrix44Type & m) const {for(int i = 0; i < 16; i++) m.V()[i]=V()[i];} @@ -250,31 +251,31 @@ public: void SetZero(); void SetIdentity(); 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); Matrix44 &SetTranslate(const Point3 &t); Matrix44 &SetTranslate(const T sx, const T sy, const T sz); - Matrix44 &SetShearXY(const T sz); - Matrix44 &SetShearXZ(const T sy); - Matrix44 &SetShearYZ(const T sx); + Matrix44 &SetShearXY(const T sz); + Matrix44 &SetShearXZ(const T sy); + Matrix44 &SetShearYZ(const T sx); ///use radiants for angle. - Matrix44 &SetRotateDeg(T AngleDeg, const Point3 & axis); - Matrix44 &SetRotateRad(T AngleRad, const Point3 & axis); + Matrix44 &SetRotateDeg(T AngleDeg, const Point3 & axis); + Matrix44 &SetRotateRad(T AngleRad, const Point3 & axis); T Determinant() const; template void Import(const Matrix44 &m) { - for(int i = 0; i < 16; i++) + for(int i = 0; i < 16; i++) _a[i] = (T)(m.V()[i]); } - template + template static inline Matrix44 Construct( const Matrix44 & b ) { Matrix44 tmp; tmp.FromMatrix(b); return tmp; } - + static inline const Matrix44 &Identity( ) { static Matrix44 tmp; tmp.SetIdentity(); @@ -291,6 +292,18 @@ public: // for the transistion to eigen const Matrix44& eval() { return *this; } + void print() { + unsigned int i, j, p; + for (i=0, p=0; i<4; i++, p+=4) + { + std::cout << "[\t"; + for (j=0; j<4; j++) + std::cout << _a[p+j] << "\t"; + + std::cout << "]\n"; + } + std::cout << "\n"; + } }; @@ -299,7 +312,7 @@ public: template class LinearSolve: public Matrix44 { public: LinearSolve(const Matrix44 &m); - Point4 Solve(const Point4 &b); // solve A � x = b + Point4 Solve(const Point4 &b); // solve A � x = b ///If you need to solve some equation you can use this function instead of Matrix44 one for speed. T Determinant() const; protected: @@ -313,13 +326,13 @@ protected: /*** Postmultiply */ //template Point3 operator*(const Point3 &p, const Matrix44 &m); -///Premultiply +///Premultiply template Point3 operator*(const Matrix44 &m, const Point3 &p); template Matrix44 &Transpose(Matrix44 &m); //return NULL matrix if not invertible -template Matrix44 &Invert(Matrix44 &m); -template Matrix44 Inverse(const Matrix44 &m); +template Matrix44 &Invert(Matrix44 &m); +template Matrix44 Inverse(const Matrix44 &m); typedef Matrix44 Matrix44s; typedef Matrix44 Matrix44i; @@ -329,11 +342,11 @@ typedef Matrix44 Matrix44d; template Matrix44::Matrix44(const Matrix44 &m) { - memcpy((T *)_a, (T *)m._a, 16 * sizeof(T)); + memcpy((T *)_a, (T *)m._a, 16 * sizeof(T)); } template Matrix44::Matrix44(const T v[]) { - memcpy((T *)_a, v, 16 * sizeof(T)); + memcpy((T *)_a, v, 16 * sizeof(T)); } template T &Matrix44::ElementAt(const int row, const int col) { @@ -356,7 +369,7 @@ template T Matrix44::ElementAt(const int row, const int col) const //template const T &Matrix44::operator[](const int i) const { // assert(i >= 0 && i < 16); // return ((T *)_a)[i]; -//} +//} template T *Matrix44::operator[](const int i) { assert(i >= 0 && i < 4); return _a+i*4; @@ -366,26 +379,26 @@ template const T *Matrix44::operator[](const int i) const { assert(i >= 0 && i < 4); return _a+i*4; } -template T *Matrix44::V() { return _a;} -template const T *Matrix44::V() const { return _a;} +template T *Matrix44::V() { return _a;} +template const T *Matrix44::V() const { return _a;} template Matrix44 Matrix44::operator+(const Matrix44 &m) const { Matrix44 ret; - for(int i = 0; i < 16; i++) - ret.V()[i] = V()[i] + m.V()[i]; + for(int i = 0; i < 16; i++) + ret.V()[i] = V()[i] + m.V()[i]; return ret; } template Matrix44 Matrix44::operator-(const Matrix44 &m) const { Matrix44 ret; - for(int i = 0; i < 16; i++) - ret.V()[i] = V()[i] - m.V()[i]; + for(int i = 0; i < 16; i++) + ret.V()[i] = V()[i] - m.V()[i]; return ret; } template Matrix44 Matrix44::operator*(const Matrix44 &m) const { - Matrix44 ret; + Matrix44 ret; for(int i = 0; i < 4; i++) for(int j = 0; j < 4; j++) { T t = 0.0; @@ -397,15 +410,15 @@ template Matrix44 Matrix44::operator*(const Matrix44 &m) const { } template Matrix44 Matrix44::operator*(const Matrix44Diag &m) const { - Matrix44 ret = (*this); + Matrix44 ret = (*this); for(int i = 0; i < 4; ++i) for(int j = 0; j < 4; ++j) - ret[i][j]*=m[i]; + ret[i][j]*=m[i]; return ret; } template Point4 Matrix44::operator*(const Point4 &v) const { - Point4 ret; + Point4 ret; for(int i = 0; i < 4; i++){ T t = 0.0; for(int k = 0; k < 4; k++) @@ -418,14 +431,14 @@ template Point4 Matrix44::operator*(const Point4 &v) const { template bool Matrix44::operator==(const Matrix44 &m) const { 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)) return false; return true; } template bool Matrix44::operator!=(const Matrix44 &m) const { 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)) return true; return false; @@ -447,7 +460,7 @@ template Matrix44 Matrix44::operator*(const T k) const { template void Matrix44::operator+=(const Matrix44 &m) { for(int i = 0; i < 16; i++) - V()[i] += m.V()[i]; + V()[i] += m.V()[i]; } template void Matrix44::operator-=(const Matrix44 &m) { for(int i = 0; i < 16; i++) @@ -455,7 +468,7 @@ template void Matrix44::operator-=(const Matrix44 &m) { } template void Matrix44::operator*=( const Matrix44 & m ) { *this = *this *m; - + /*for(int i = 0; i < 4; i++) { //sbagliato Point4 t(0, 0, 0, 0); for(int k = 0; k < 4; k++) { @@ -487,7 +500,7 @@ void Matrix44::ToEulerAngles(T &alpha, T &beta, T &gamma) gamma = atan2(ElementAt(0,1), ElementAt(1,1)); } -template +template void Matrix44::FromEulerAngles(T alpha, T beta, T gamma) { this->SetZero(); @@ -499,18 +512,18 @@ void Matrix44::FromEulerAngles(T alpha, T beta, T gamma) T sinbeta = sin(beta); T singamma = sin(gamma); - ElementAt(0,0) = cosbeta * cosgamma; - ElementAt(1,0) = -cosalpha * singamma + sinalpha * sinbeta * cosgamma; + ElementAt(0,0) = cosbeta * cosgamma; + ElementAt(1,0) = -cosalpha * singamma + sinalpha * sinbeta * cosgamma; ElementAt(2,0) = sinalpha * singamma + cosalpha * sinbeta * cosgamma; - + ElementAt(0,1) = cosbeta * singamma; - ElementAt(1,1) = cosalpha * cosgamma + sinalpha * sinbeta * singamma; + ElementAt(1,1) = cosalpha * cosgamma + sinalpha * sinbeta * singamma; ElementAt(2,1) = -sinalpha * cosgamma + cosalpha * sinbeta * singamma; - - ElementAt(0,2) = -sinbeta; - ElementAt(1,2) = sinalpha * cosbeta; + + ElementAt(0,2) = -sinbeta; + ElementAt(1,2) = sinalpha * cosbeta; ElementAt(2,2) = cosalpha * cosbeta; - + ElementAt(3,3) = 1; } @@ -518,8 +531,8 @@ template void Matrix44::SetZero() { memset((T *)_a, 0, 16 * sizeof(T)); } -template void Matrix44::SetIdentity() { - SetDiagonal(1); +template void Matrix44::SetIdentity() { + SetDiagonal(1); } template void Matrix44::SetDiagonal(const T k) { @@ -527,7 +540,7 @@ template void Matrix44::SetDiagonal(const T k) { ElementAt(0, 0) = k; ElementAt(1, 1) = k; ElementAt(2, 2) = k; - ElementAt(3, 3) = 1; + ElementAt(3, 3) = 1; } template Matrix44 &Matrix44::SetScale(const Point3 &t) { @@ -555,15 +568,15 @@ template Matrix44 &Matrix44::SetTranslate(const T tx, const T ty return *this; } -template Matrix44 &Matrix44::SetRotateDeg(T AngleDeg, const Point3 & axis) { +template Matrix44 &Matrix44::SetRotateDeg(T AngleDeg, const Point3 & axis) { return SetRotateRad(math::ToRad(AngleDeg),axis); } -template Matrix44 &Matrix44::SetRotateRad(T AngleRad, const Point3 & axis) { +template Matrix44 &Matrix44::SetRotateRad(T AngleRad, const Point3 & axis) { //angle = angle*(T)3.14159265358979323846/180; e' in radianti! T c = math::Cos(AngleRad); T s = math::Sin(AngleRad); - T q = 1-c; + T q = 1-c; Point3 t = axis; t.Normalize(); ElementAt(0,0) = t[0]*t[0]*q + c; @@ -579,14 +592,14 @@ template Matrix44 &Matrix44::SetRotateRad(T AngleRad, const Poin ElementAt(2,2) = t[2]*t[2]*q +c; ElementAt(2,3) = 0; ElementAt(3,0) = 0; - ElementAt(3,1) = 0; + ElementAt(3,1) = 0; ElementAt(3,2) = 0; - ElementAt(3,3) = 1; + ElementAt(3,3) = 1; return *this; } /* Shear Matrixes - XY + XY 1 k 0 0 x x+ky 0 1 0 0 y y 0 0 1 0 z z @@ -604,19 +617,19 @@ template Matrix44 &Matrix44::SetRotateRad(T AngleRad, const Poin */ - template Matrix44 & Matrix44:: SetShearXY( const T sh) {// shear the X coordinate as the Y coordinate change + template Matrix44 & Matrix44:: SetShearXY( const T sh) {// shear the X coordinate as the Y coordinate change SetIdentity(); ElementAt(0,1) = sh; return *this; } - - template Matrix44 & Matrix44:: SetShearXZ( const T sh) {// shear the X coordinate as the Z coordinate change + + template Matrix44 & Matrix44:: SetShearXZ( const T sh) {// shear the X coordinate as the Z coordinate change SetIdentity(); ElementAt(0,2) = sh; return *this; } - - template Matrix44 &Matrix44:: SetShearYZ( const T sh) {// shear the Y coordinate as the Z coordinate change + + template Matrix44 &Matrix44:: SetShearYZ( const T sh) {// shear the Y coordinate as the Z coordinate change SetIdentity(); ElementAt(1,2) = sh; return *this; @@ -625,7 +638,7 @@ template Matrix44 &Matrix44::SetRotateRad(T AngleRad, const Poin /* 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 - ScaleV and Tranv are obiviously scaling and translation. @@ -635,7 +648,7 @@ This procedure decompose it in a sequence of The input matrix is modified leaving inside it a simple roto translation. To obtain the original matrix the above transformation have to be applied in the strict following way: - + OriginalMatrix = Trn * Rtx*Rty*Rtz * ShearYZ*ShearXZ*ShearXY * Scl Example Code: @@ -655,7 +668,7 @@ double srv() { return (double(rand()%40)-20)/2.0; } // small random value 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 Trn; Trn.SetTranslate(TrV); - + Matrix44d StartM = Trn * Rtx*Rty*Rtz * Syz*Sxz*Sxy *Scl; Matrix44d ResultM=StartM; Decompose(ResultM,ScVOut,ShVOut,RtVOut,TrVOut); @@ -670,75 +683,75 @@ double srv() { return (double(rand()%40)-20)/2.0; } // small random value Trn.SetTranslate(TrVOut); // 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 -bool Decompose(Matrix44 &M, Point3 &ScaleV, Point3 &ShearV, Point3 &RotV,Point3 &TranV) +bool Decompose(Matrix44 &M, Point3 &ScaleV, Point3 &ShearV, Point3 &RotV,Point3 &TranV) { if(!(M[3][0]==0 && M[3][1]==0 && M[3][2]==0 && M[3][3]==1) ) // the matrix is projective return false; 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); - - // Second Step Recover Scale and Shearing interleaved + + // Second Step Recover Scale and Shearing interleaved ScaleV[0]=Norm(M.GetColumn3(0)); Point3 R[3]; R[0]=M.GetColumn3(0); R[0].Normalize(); - - 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]; 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]; - 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]; assert(math::Abs(R[2]*R[0])<1e-10); - + R[2] = R[2]-R[1]*(R[2]*R[1]); assert(math::Abs(R[2]*R[1])<1e-10); assert(math::Abs(R[2]*R[0])<1e-10); - + ScaleV[2]=Norm(R[2]); - ShearV[1]=ShearV[1]/ScaleV[2]; + ShearV[1]=ShearV[1]/ScaleV[2]; R[2]=R[2]/ScaleV[2]; assert(math::Abs(R[2]*R[1])<1e-10); assert(math::Abs(R[2]*R[0])<1e-10); - + ShearV[2]=R[1]*M.GetColumn3(2); // yz shearing ShearV[2]=ShearV[2]/ScaleV[2]; - int i,j; + int i,j; for(i=0;i<3;++i) for(j=0;j<3;++j) M[i][j]=R[j][i]; // 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(); 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... - if(det<0) { + if(det<0) { ScaleV *= -1; M *= -1; } 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); if(math::Abs(cosbeta) > 1e-5) - { + { alpha=asin(-M[1][2]/cosbeta); if((M[2][2]/cosbeta) < 0 ) alpha=M_PI-alpha; gamma=asin(-M[0][1]/cosbeta); if((M[0][0]/cosbeta)<0) gamma = M_PI-gamma; } else - { + { alpha=asin(-M[1][0]); if(M[1][1]<0) alpha=M_PI-alpha; gamma=0; @@ -754,7 +767,7 @@ bool Decompose(Matrix44 &M, Point3 &ScaleV, Point3 &ShearV, Point3 & -template T Matrix44::Determinant() const { +template T Matrix44::Determinant() const { LinearSolve solve(*this); return solve.Determinant(); } @@ -785,45 +798,45 @@ template Point3 operator*(const Matrix44 &m, const Point3 &p) template Matrix44 &Transpose(Matrix44 &m) { for(int i = 1; i < 4; i++) 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; } /* - To invert a matrix you can + 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 Matrix44 & Invert(Matrix44 &m) { +template Matrix44 & Invert(Matrix44 &m) { LinearSolve solve(m); for(int j = 0; j < 4; j++) { //Find inverse by columns. Point4 col(0, 0, 0, 0); col[j] = 1.0; col = solve.Solve(col); - for(int i = 0; i < 4; i++) + for(int i = 0; i < 4; i++) m.ElementAt(i, j) = col[i]; - } + } return m; } -template Matrix44 Inverse(const Matrix44 &m) { +template Matrix44 Inverse(const Matrix44 &m) { LinearSolve solve(m); Matrix44 res; for(int j = 0; j < 4; j++) { //Find inverse by columns. Point4 col(0, 0, 0, 0); col[j] = 1.0; col = solve.Solve(col); - for(int i = 0; i < 4; i++) + for(int i = 0; i < 4; i++) res.ElementAt(i, j) = col[i]; - } + } return res; } @@ -842,8 +855,8 @@ template LinearSolve::LinearSolve(const Matrix44 &m): Matrix44 T LinearSolve::Determinant() const { T det = d; - for(int j = 0; j < 4; j++) - det *= this-> ElementAt(j, j); + for(int j = 0; j < 4; j++) + det *= this-> ElementAt(j, j); return det; } @@ -852,18 +865,18 @@ template T LinearSolve::Determinant() const { d is +or -1 depeneing of row permutation even or odd.*/ #define TINY 1e-100 -template bool LinearSolve::Decompose() { - +template bool LinearSolve::Decompose() { + /* Matrix44 A; for(int i = 0; i < 16; i++) - A[i] = operator[](i); - SetIdentity(); + A[i] = operator[](i); + SetIdentity(); Point4 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++) + 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; } @@ -895,23 +908,23 @@ template bool LinearSolve::Decompose() { 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++ ) + //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; + 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++) { + for(i = 0; i < 4; i++) { T largest = 0.0; for(j = 0; j < 4; j++) { T t = math::Abs(this->ElementAt(i, j)); @@ -920,70 +933,70 @@ template bool LinearSolve::Decompose() { if (largest == 0.0) { //oooppps there is a zero row! return false; - } - scaling[i] = (T)1.0 / largest; + } + scaling[i] = (T)1.0 / largest; } int imax = 0; - for(j = 0; j < 4; j++) { + 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++) + 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 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) { + if(t >= largest) { largest = t; imax = i; } } - if (j != imax) { - for (int k = 0; k < 4; k++) { + 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]; + scaling[imax] = scaling[j]; } index[j]=imax; if (this->ElementAt(j,j) == 0.0) this->ElementAt(j,j) = (T)TINY; - if (j != 3) { + if (j != 3) { T dum = (T)1.0 / (this->ElementAt(j,j)); - for (i = j+1; i < 4; i++) + for (i = j+1; i < 4; i++) this->ElementAt(i,j) *= dum; } } - return true; + return true; } -template Point4 LinearSolve::Solve(const Point4 &b) { +template Point4 LinearSolve::Solve(const Point4 &b) { Point4 x(b); - int first = -1, ip; - for(int i = 0; i < 4; i++) { + 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++) + for(int j = first; j <= i-1; j++) sum -= this->ElementAt(i,j) * x[j]; - else - if(sum) first = i; + else + if(sum) first = i; x[i] = sum; } - for (int i = 3; i >= 0; i--) { + for (int i = 3; i >= 0; i--) { T sum = x[i]; - for (int j = i+1; j < 4; j++) + for (int j = i+1; j < 4; j++) sum -= this->ElementAt(i, j) * x[j]; - x[i] = sum / this->ElementAt(i, i); + x[i] = sum / this->ElementAt(i, i); } return x; } diff --git a/vcg/math/eigen.h b/vcg/math/eigen.h index 4e15336a..d2eb1235 100644 --- a/vcg/math/eigen.h +++ b/vcg/math/eigen.h @@ -31,7 +31,24 @@ // forward declarations namespace Eigen { + +#include "../Eigen/src/Core/util/Meta.h" + template struct ei_lexi_comparison; + +template::ret, + bool SameSize = Derived1::SizeAtCompileTime==Derived2::SizeAtCompileTime> +struct ei_import_selector; + +template +struct ei_to_vcgtype; + } #include "base.h" @@ -156,6 +173,46 @@ template struct ei_lexi_comparison +struct ei_import_selector +{ + static void run(Derived1& a, const Derived2& b) { a = b; } +}; + +template +struct ei_import_selector +{ + static void run(Derived1& a, const Derived2& b) + { a = b.template cast(); } +}; + +template +struct ei_import_selector +{ + static void run(Derived1& a, const Derived2& b) + { + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived1); + EIGEN_STATIC_ASSERT_FIXED_SIZE(Derived1); + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived2); + EIGEN_STATIC_ASSERT_FIXED_SIZE(Derived2); + enum { + Size1 = Derived1::SizeAtCompileTime, + Size2 = Derived2::SizeAtCompileTime + }; + assert(Size1<=4 && Size2<=4); + a.coeffRef(0) = Scalar(b.coeff(0)); + if (Size1>1) { if (Size2>1) a.coeffRef(1) = Scalar(b.coeff(1)); else a.coeffRef(1) = 0; } + if (Size1>2) { if (Size2>2) a.coeffRef(2) = Scalar(b.coeff(2)); else a.coeffRef(2) = 0; } + if (Size1>3) { if (Size2>3) a.coeffRef(3) = Scalar(b.coeff(3)); else a.coeffRef(3) = 0; } + } +}; + +// default implementation of ei_to_vcgtype +// the specialization are with +template +struct ei_to_vcgtype { typedef Matrix type; }; + } #define VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATOR(Derived, Op) \ diff --git a/vcg/math/eigen_matrix_addons.h b/vcg/math/eigen_matrix_addons.h index 09da01b0..7d38e860 100644 --- a/vcg/math/eigen_matrix_addons.h +++ b/vcg/math/eigen_matrix_addons.h @@ -24,28 +24,31 @@ #warning You are including deprecated math stuff enum {Dimension = SizeAtCompileTime}; +typedef typename ei_to_vcgtype::type EquivVcgType; typedef vcg::VoidType ParamType; typedef Matrix PointType; using Base::V; +// automatic conversion to similar vcg types +// the otherway round is implicit because they inherits this Matrix tyoe +operator EquivVcgType& () { return *reinterpret_cast(this); } +operator const EquivVcgType& () const { return *reinterpret_cast(this); } + +/** \deprecated use m.cast() */ /// importer for points with different scalar type and-or dimensionality // FIXME the Point3/Point4 specialization were only for same sizes ?? // while the Point version was generic like this one template inline void Import(const MatrixBase& b) { - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Matrix); - EIGEN_STATIC_ASSERT_FIXED_SIZE(Matrix); - EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived); - EIGEN_STATIC_ASSERT_FIXED_SIZE(OtherDerived); - enum { OtherSize = OtherDerived::SizeAtCompileTime }; - assert(SizeAtCompileTime<=4 && OtherSize<=4); - data()[0] = Scalar(b[0]); - if (SizeAtCompileTime>1) { if (OtherSize>1) data()[1] = Scalar(b[1]); else data()[1] = 0; } - if (SizeAtCompileTime>2) { if (OtherSize>2) data()[2] = Scalar(b[2]); else data()[2] = 0; } - if (SizeAtCompileTime>3) { if (OtherSize>3) data()[3] = Scalar(b[3]); else data()[3] = 0; } + ei_import_selector::run(*this,b.derived()); } +/// constructor for points with different scalar type and-or dimensionality +template +static inline Matrix Construct(const MatrixBase& b) +{ Matrix p; p.Import(b); return p; } + /// importer for homogeneous points template inline void ImportHomo(const MatrixBase& b) @@ -53,26 +56,15 @@ inline void ImportHomo(const MatrixBase& b) EIGEN_STATIC_ASSERT_VECTOR_ONLY(Matrix); EIGEN_STATIC_ASSERT_FIXED_SIZE(Matrix); EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(OtherDerived,SizeAtCompileTime-1); - + this->template start = b; data()[SizeAtCompileTime-1] = Scalar(1.0); } -/// constructor for points with different scalar type and-or dimensionality -template -static inline Matrix Construct(const MatrixBase& b) -{ - Matrix p; p.Import(b); - return p; -} - - /// constructor for homogeneus point. +/// constructor for homogeneus point. template static inline Matrix ConstructHomo(const MatrixBase& b) -{ - Matrix p; p.ImportHomo(b); - return p; -} +{ Matrix p; p.ImportHomo(b); return p; } inline const Scalar &X() const { return data()[0]; } inline const Scalar &Y() const { return data()[1]; } @@ -81,15 +73,19 @@ inline Scalar &X() { return data()[0]; } inline Scalar &Y() { return data()[1]; } inline Scalar &Z() { assert(SizeAtCompileTime>2); return data()[2]; } -// note, W always returns the last entry +/** note, W always returns the last entry */ inline Scalar& W() { return data()[SizeAtCompileTime-1]; } +/** note, W always returns the last entry */ inline const Scalar& W() const { return data()[SizeAtCompileTime-1]; } -Scalar* V() { return data(); } -const Scalar* V() const { return data(); } +/** \deprecated use .data() */ +EIGEN_DEPRECATED Scalar* V() { return data(); } +/** \deprecated use .data() */ +EIGEN_DEPRECATED const Scalar* V() const { return data(); } +/** \deprecated use m.coeff(i) or m[i] or m(i) */ // overloaded to return a const reference -inline const Scalar& V( const int i ) const +EIGEN_DEPRECATED inline const Scalar& V( const int i ) const { assert(i>=0 && i class Matrix; namespace Eigen{ template struct ei_traits > : ei_traits > {}; +template struct ei_to_vcgtype +{ typedef vcg::ndim::Matrix type; }; } namespace vcg{ @@ -50,7 +52,7 @@ namespace ndim{ /* @{ */ /*! - * \deprecated use Matrix or Matrix + * \deprecated use Matrix or Matrix or any typedef * This class represent a generic m�n matrix. The class is templated over the scalar type field. * @param Scalar (Templete Parameter) Specifies the ScalarType field. */ @@ -63,7 +65,6 @@ public: _EIGEN_GENERIC_PUBLIC_INTERFACE(Matrix,_Base); typedef _Scalar ScalarType; - VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Matrix) /*! diff --git a/vcg/math/matrix33.h b/vcg/math/matrix33.h index 726274bc..0912d94b 100644 --- a/vcg/math/matrix33.h +++ b/vcg/math/matrix33.h @@ -38,12 +38,14 @@ template class Matrix33; namespace Eigen{ template struct ei_traits > : ei_traits > {}; +template struct ei_to_vcgtype +{ typedef vcg::Matrix33 type; }; } namespace vcg { - -/** @name Matrix33 +/** \deprecated use Matrix + @name Matrix33 Class Matrix33. This is the class for definition of a matrix 3x3. @param S (Templete Parameter) Specifies the ScalarType field. @@ -98,7 +100,6 @@ public: return SetRotateRad(math::ToRad(angle),axis); } - // Warning, this Inversion code can be HIGHLY NUMERICALLY UNSTABLE! // In most case you are advised to use the Invert() method based on SVD decomposition. /** \deprecated */ @@ -110,18 +111,13 @@ public: printf("| %g \t%g \t%g |\n",coeff(i,0),coeff(i,1),coeff(i,2)); } - /* + /** \deprecated use a * b.transpose() compute the matrix generated by the product of a * b^T */ - void ExternalProduct(const Point3 &a, const Point3 &b) - { - for(int i=0;i<3;++i) - for(int j=0;j<3;++j) - (*this)[i][j] = a[i]*b[j]; - } + // hm.... this is the outer product + void ExternalProduct(const Point3 &a, const Point3 &b) { *this = a * b.transpose(); } - /* Compute the Frobenius Norm of the Matrix - */ + /** Compute the Frobenius Norm of the Matrix */ Scalar Norm() { return Base::cwise().abs2().sum(); } // { // // FIXME looks like there was a bug: j is not used !!! diff --git a/vcg/math/matrix44.h b/vcg/math/matrix44.h index 89c64b50..5fc6f397 100644 --- a/vcg/math/matrix44.h +++ b/vcg/math/matrix44.h @@ -41,6 +41,8 @@ template class Matrix44; namespace Eigen{ template struct ei_traits > : ei_traits > {}; +template struct ei_to_vcgtype +{ typedef vcg::Matrix44 type; }; } namespace vcg { @@ -140,17 +142,6 @@ public: Matrix44 &SetRotateDeg(Scalar AngleDeg, const Point3 & axis); Matrix44 &SetRotateRad(Scalar AngleRad, const Point3 & axis); - template void Import(const Matrix44 &m) { - for(int i = 0; i < 16; i++) - Base::data()[i] = (Scalar)(m.data()[i]); - } - template - static inline Matrix44 Construct( const Matrix44 & b ) - { - Matrix44 tmp; tmp.FromMatrix(b); - return tmp; - } - // template Point3 operator*(const Point3 &p) { // T w; // Point3 s; @@ -173,30 +164,10 @@ public: return s; } + void print() {std::cout << *this << "\n\n";} + }; - -/** Class for solving A * x = b. */ -template class LinearSolve: public Matrix44 { -public: - LinearSolve(const Matrix44 &m); - Point4 Solve(const Point4 &b); // solve A � x = b - ///If you need to solve some equation you can use this function instead of Matrix44 one for speed. - T Determinant() const; -protected: - ///Holds row permutation. - int index[4]; //hold permutation - ///Hold sign of row permutation (used for determinant sign) - T d; - bool Decompose(); -}; - -/*** Postmultiply */ -//template Point3 operator*(const Point3 &p, const Matrix44 &m); - -///Premultiply -// template Point3 operator*(const Matrix44 &m, const Point3 &p); - //return NULL matrix if not invertible template Matrix44 &Invert(Matrix44 &m); template Matrix44 Inverse(const Matrix44 &m); @@ -206,18 +177,6 @@ typedef Matrix44 Matrix44i; typedef Matrix44 Matrix44f; typedef Matrix44 Matrix44d; - -//template T &Matrix44::operator[](const int i) { -// assert(i >= 0 && i < 16); -// return ((T *)_a)[i]; -//} -// -//template const T &Matrix44::operator[](const int i) const { -// assert(i >= 0 && i < 16); -// return ((T *)_a)[i]; -//} - - template < class PointType , class T > void operator*=( std::vector &vert, const Matrix44 & m ) { typename std::vector::iterator ii; for(ii=vert.begin();ii!=vert.end();++ii) @@ -487,22 +446,6 @@ bool Decompose(Matrix44 &M, Point3 &ScaleV, Point3 &ShearV, Point3 & return true; } - - - - -// template Point3 operator*(const Point3 &p, const Matrix44 &m) { -// T w; -// Point3 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; -// } - - /* To invert a matrix you can either invert the matrix inplace calling diff --git a/vcg/math/shot.h b/vcg/math/shot.h index 4cbd125c..e992661d 100644 --- a/vcg/math/shot.h +++ b/vcg/math/shot.h @@ -8,7 +8,7 @@ * \ * * All rights reserved. * * * -* This program is free software; you can redistribute it and/or modify * +* 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. * @@ -93,9 +93,9 @@ creation ****************************************************************************/ -/** class Shot +/** class Shot The shot is made of two things: -* the Instrinsics paramaters, which are stored as a Camera type (check vcg/math/camera) and that +* the Instrinsics paramaters, which are stored as a Camera type (check vcg/math/camera) and that determines how a point in the frame of the camera is projected in the 2D projection plane * the Extrinsics parameters, which are stored in the class Shot (che the type ReferenceFrame) @@ -103,7 +103,7 @@ and that tell viewpoint and view direction. The Extrinsics parameters are kept as a rotation matrix "rot" and a translation vector "tra" NOTE: the translation matrix "tra" corresponds to -viewpoint while the rotation matrix -"rot" corresponds to the axis of the reference frame by row, i.e. +"rot" corresponds to the axis of the reference frame by row, i.e. rot[0][0|1|2] == X axis rot[1][0|1|2] == Y axis rot[2][0|1|2] == Z axis @@ -136,7 +136,7 @@ public: class ReferenceFrame { friend class Shot; RotoType rot; // rotation - Point3 tra; // viewpoint + Point3 tra; // viewpoint public: void SetIdentity(){ rot.SetIdentity(); tra = Point3(0.0,0.0,0.0);} void SetTra(const Point3 & tr) {tra = tr;} @@ -146,7 +146,7 @@ public: }; Camera Intrinsics; // the camera that made the shot - ReferenceFrame Extrinsics; // the position and orientation of the camera + ReferenceFrame Extrinsics; // the position and orientation of the camera Shot(Camera c) @@ -186,13 +186,13 @@ public: /// convert a 3d point from camera to world coordinates vcg::Point3 ConvertCameraToWorldCoordinates(const vcg::Point3 & p) const; - /// project a 3d point from world coordinates to 2d camera viewport (pixel) + /// project a 3d point from world coordinates to 2d camera viewport (pixel) vcg::Point2 Project(const vcg::Point3 & p) const; /// inverse projection from 2d camera viewport (pixel) + Zdepth to 3d world coordinates vcg::Point3 UnProject(const vcg::Point2 & p, const S & d) const; - /// returns distance of point p from camera plane (z depth), used for unprojection + /// returns distance of point p from camera plane (z depth), used for unprojection S Depth(const vcg::Point3 & p)const; @@ -205,8 +205,7 @@ public: Matrix44 GetExtrinsicsToWorldMatrix() const { Matrix44 rotM ; Extrinsics.rot.ToMatrix(rotM); - Transpose(rotM); - return Matrix44().SetTranslate(Extrinsics.tra) * rotM; + return Matrix44().SetTranslate(Extrinsics.tra) * rotM.transpose(); } /* Returns the matrix M such that @@ -240,7 +239,7 @@ public: //--- /// GET the viewpoint template -const vcg::Point3 Shot::GetViewPoint() const +const vcg::Point3 Shot::GetViewPoint() const { return Extrinsics.tra; } @@ -254,10 +253,10 @@ void Shot::SetViewPoint(const vcg::Point3 & viewpoint) /// GET the i-th axis of the coordinate system of the camera template -vcg::Point3 Shot::Axis(const int & i) const -{ - vcg::Matrix44 m; - Extrinsics.rot.ToMatrix(m); +vcg::Point3 Shot::Axis(const int & i) const +{ + vcg::Matrix44 m; + Extrinsics.rot.ToMatrix(m); vcg::Point3 aa = m.GetRow3(i); return aa; } @@ -285,7 +284,7 @@ void Shot::LookTowards(const vcg::Point3 & z_dir,const vcg::P { vcg::Point3 x_dir = up ^-z_dir; vcg::Point3 y_dir = -z_dir ^x_dir; - + Matrix44 m; m.SetIdentity(); *(vcg::Point3 *)&m[0][0] = x_dir/x_dir.Norm(); @@ -316,12 +315,11 @@ vcg::Point3 Shot::ConvertCameraToWorldCoordinates(const vcg:: vcg::Point3 cp = p; cp[2]=-cp[2]; // note: World reference system is left handed Extrinsics.rot.ToMatrix(rotM); - Transpose(rotM); - cp = rotM * cp+ GetViewPoint(); + cp = rotM.transpose() * cp + GetViewPoint(); return cp; } -/// project a 3d point from world coordinates to 2d camera viewport (pixel) +/// project a 3d point from world coordinates to 2d camera viewport (pixel) template vcg::Point2 Shot::Project(const vcg::Point3 & p) const { @@ -341,9 +339,9 @@ vcg::Point3 Shot::UnProject(const vcg::Point2 & p, const S return wp; } -/// returns distance of point p from camera plane (z depth), used for unprojection +/// returns distance of point p from camera plane (z depth), used for unprojection template -S Shot::Depth(const vcg::Point3 & p)const +S Shot::Depth(const vcg::Point3 & p)const { return ConvertWorldToCameraCoordinates(p).Z(); } diff --git a/vcg/space/point2.h b/vcg/space/point2.h index 749e751f..19129d00 100644 --- a/vcg/space/point2.h +++ b/vcg/space/point2.h @@ -37,6 +37,8 @@ template class Point2; namespace Eigen { template struct ei_traits > : ei_traits > {}; +template struct ei_to_vcgtype +{ typedef vcg::Point2 type; }; } namespace vcg { diff --git a/vcg/space/point3.h b/vcg/space/point3.h index d1b8e28a..f1b81077 100644 --- a/vcg/space/point3.h +++ b/vcg/space/point3.h @@ -8,7 +8,7 @@ * \ * * All rights reserved. * * * -* This program is free software; you can redistribute it and/or modify * +* 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. * @@ -38,6 +38,9 @@ namespace Eigen{ template struct ei_traits > : ei_traits > {}; +template struct ei_to_vcgtype +{ typedef vcg::Point3 type; }; + template struct NumTraits > : NumTraits { @@ -79,7 +82,7 @@ public: _EIGEN_GENERIC_PUBLIC_INTERFACE(Point3,_Base); VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Point3) - /** @name Standard Constructors and Initializers + /** @name Standard Constructors and Initializers No casting operators have been introduced to avoid automatic unattended (and costly) conversion between different point types **/ @@ -92,7 +95,7 @@ public: // this one is very useless - template + template static inline Point3 Construct( const Q & P0, const Q & P1, const Q & P2) { return Point3(Scalar(P0),Scalar(P1),Scalar(P2)); diff --git a/vcg/space/point4.h b/vcg/space/point4.h index 71219c37..3c753811 100644 --- a/vcg/space/point4.h +++ b/vcg/space/point4.h @@ -8,7 +8,7 @@ * \ * * All rights reserved. * * * -* This program is free software; you can redistribute it and/or modify * +* 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. * @@ -36,6 +36,8 @@ template class Point4; namespace Eigen { template struct ei_traits > : ei_traits > {}; +template struct ei_to_vcgtype +{ typedef vcg::Point4 type; }; } namespace vcg { diff --git a/wrap/gl/deprecated_math.h b/wrap/gl/deprecated_math.h new file mode 100644 index 00000000..b90998ab --- /dev/null +++ b/wrap/gl/deprecated_math.h @@ -0,0 +1,144 @@ +/**************************************************************************** +* 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. * +* * +****************************************************************************/ +/**************************************************************************** + History + +$Log: not supported by cvs2svn $ +Revision 1.10 2006/02/13 13:05:05 cignoni +Removed glew inclusion + +Revision 1.9 2004/09/30 00:48:07 ponchio + -> + +Revision 1.8 2004/09/28 14:04:36 ganovelli +glGet added + +Revision 1.7 2004/07/13 15:55:57 cignoni +Added test on presence of glTranspose extension (for old hw support) + +Revision 1.6 2004/05/26 15:12:39 cignoni +Removed inclusion of gl extension stuff + +Revision 1.5 2004/05/12 20:54:55 ponchio +*** empty log message *** + +Revision 1.4 2004/05/12 13:07:47 ponchio +Added #include + +Revision 1.3 2004/05/04 23:36:23 cignoni +remove include of gl and added glextgension exploiting, + +Revision 1.2 2004/04/07 10:47:03 cignoni +inlined functions for avoid multiple linking errors + +Revision 1.1 2004/03/31 15:27:17 ponchio +*** empty log message *** + + +****************************************************************************/ + +#ifndef VCG_GL_MATH_H +#define VCG_GL_MATH_H + +// Please note that this file assume that you have already included your +// gl-extension wrapping utility, and that therefore all the extension symbol are already defined. + +#include +#include +//#include // please do not include it! + + +namespace vcg { + +inline void glMultMatrixE(const Matrix44f &matrix) { + //glMultMatrixf((const GLfloat *)(matrix[0])); + if(glMultTransposeMatrixf) glMultTransposeMatrixf((const GLfloat *)(matrix.V())); + else { + glMultMatrixf((const GLfloat *)(matrix.transpose().eval().V())); + } +} + +inline void glMultMatrix(const Matrix44f &matrix) { + glMultMatrixf((const GLfloat *)(matrix.transpose().eval().V())); +} + +inline void glMultMatrixE(const Matrix44d &matrix) { + if(glMultTransposeMatrixd) glMultTransposeMatrixd((const GLdouble *)(matrix.V())); + else { + glMultMatrixd((const GLdouble *)(matrix.transpose().eval().V())); + } +} +inline void glMultMatrix(const Matrix44d &matrix) { + glMultMatrixd((const GLdouble *)(matrix.transpose().eval().V())); +} + +inline void glMultMatrixDirect(const Matrix44f &matrix) { + glMultMatrixf((const GLfloat *)(matrix.V())); +} + +inline void glMultMatrixDirect(const Matrix44d &matrix) { + glMultMatrixd((const GLdouble *)(matrix.V())); +} + +inline void glMultMatrix(const Similarityf &s) { + glTranslatef(s.tra[0], s.tra[1], s.tra[2]); + glScalef(s.sca, s.sca, s.sca); + float alpha; + Point3f axis; + s.rot.ToAxis(alpha, axis); + glRotatef(math::ToDeg(alpha), axis[0], axis[1], axis[2]); + +} + +inline void glMultMatrix(const Similarityd &s) { + glTranslated(s.tra[0], s.tra[1], s.tra[2]); + double alpha; + Point3d axis; + s.rot.ToAxis(alpha, axis); + glRotated(math::ToDeg(alpha), axis[0], axis[1], axis[2]); + glScaled(s.sca, s.sca, s.sca); +} + +inline void glGetv(const GLenum pname, Matrix44f & m){ + Matrix44f tmp; + glGetFloatv(pname,tmp.V()); + m = tmp.transpose(); +} + +inline void glGetv(const GLenum pname, Matrix44d & m){ + Matrix44d tmp; + glGetDoublev(pname,tmp.V()); + m = tmp.transpose(); +} + +inline void glGetDirectv(const GLenum pname, Matrix44f & m){ + glGetFloatv(pname,m.V()); +} + +inline void glGetDirecv(const GLenum pname, Matrix44d & m){ + glGetDoublev(pname,m.V()); +} + + +}//namespace +#endif diff --git a/wrap/gl/math.h b/wrap/gl/math.h index b90998ab..5fef8b58 100644 --- a/wrap/gl/math.h +++ b/wrap/gl/math.h @@ -8,7 +8,7 @@ * \ * * All rights reserved. * * * -* This program is free software; you can redistribute it and/or modify * +* 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. * @@ -57,10 +57,14 @@ Revision 1.1 2004/03/31 15:27:17 ponchio ****************************************************************************/ +#ifndef VCG_USE_EIGEN +#include "deprecated_math.h" +#else + #ifndef VCG_GL_MATH_H #define VCG_GL_MATH_H -// Please note that this file assume that you have already included your +// Please note that this file assume that you have already included your // gl-extension wrapping utility, and that therefore all the extension symbol are already defined. #include @@ -70,44 +74,41 @@ Revision 1.1 2004/03/31 15:27:17 ponchio namespace vcg { -inline void glMultMatrixE(const Matrix44f &matrix) { - //glMultMatrixf((const GLfloat *)(matrix[0])); - if(glMultTransposeMatrixf) glMultTransposeMatrixf((const GLfloat *)(matrix.V())); - else { - glMultMatrixf((const GLfloat *)(matrix.transpose().eval().V())); - } +inline void glLoadMatrix(const Eigen::Matrix& matrix) { glLoadMatrixf(matrix.data()); } +inline void glLoadMatrix(const Eigen::Matrix& matrix) { + Eigen::Matrix4f tmp(matrix); + glLoadMatrixf(tmp.data()); } +inline void glLoadMatrix(const Eigen::Matrix& matrix) { glLoadMatrixd(matrix.data()); } +inline void glLoadMatrix(const Eigen::Matrix& matrix) { + Eigen::Matrix4d tmp(matrix); + glLoadMatrixd(tmp.data()); +} +template +inline void glLoadMatrix(const Eigen::Transform& t) { glLoadMatrix(t.matrix()); } -inline void glMultMatrix(const Matrix44f &matrix) { - glMultMatrixf((const GLfloat *)(matrix.transpose().eval().V())); -} -inline void glMultMatrixE(const Matrix44d &matrix) { - if(glMultTransposeMatrixd) glMultTransposeMatrixd((const GLdouble *)(matrix.V())); - else { - glMultMatrixd((const GLdouble *)(matrix.transpose().eval().V())); - } +inline void glMultMatrix(const Eigen::Matrix& matrix) { glMultMatrixf(matrix.data()); } +inline void glMultMatrix(const Eigen::Matrix& matrix) { + Eigen::Matrix tmp(matrix); + glMultMatrixf(tmp.data()); } -inline void glMultMatrix(const Matrix44d &matrix) { - glMultMatrixd((const GLdouble *)(matrix.transpose().eval().V())); +inline void glMultMatrix(const Eigen::Matrix& matrix) { glMultMatrixd(matrix.data()); } +inline void glMultMatrix(const Eigen::Matrix& matrix) { + Eigen::Matrix tmp(matrix); + glMultMatrixd(tmp.data()); } +template +inline void glMultMatrix(const Eigen::Transform& t) { glMultMatrix(t.matrix()); } -inline void glMultMatrixDirect(const Matrix44f &matrix) { - glMultMatrixf((const GLfloat *)(matrix.V())); -} - -inline void glMultMatrixDirect(const Matrix44d &matrix) { - glMultMatrixd((const GLdouble *)(matrix.V())); -} inline void glMultMatrix(const Similarityf &s) { glTranslatef(s.tra[0], s.tra[1], s.tra[2]); glScalef(s.sca, s.sca, s.sca); float alpha; Point3f axis; - s.rot.ToAxis(alpha, axis); - glRotatef(math::ToDeg(alpha), axis[0], axis[1], axis[2]); - + s.rot.ToAxis(alpha, axis); + glRotatef(math::ToDeg(alpha), axis[0], axis[1], axis[2]); } inline void glMultMatrix(const Similarityd &s) { @@ -119,26 +120,27 @@ inline void glMultMatrix(const Similarityd &s) { glScaled(s.sca, s.sca, s.sca); } -inline void glGetv(const GLenum pname, Matrix44f & m){ - Matrix44f tmp; - glGetFloatv(pname,tmp.V()); - m = tmp.transpose(); +inline void glGetv(const GLenum pname, Eigen::Matrix& matrix){ + glGetFloatv(pname,matrix.data()); } - -inline void glGetv(const GLenum pname, Matrix44d & m){ - Matrix44d tmp; - glGetDoublev(pname,tmp.V()); - m = tmp.transpose(); +inline void glGetv(const GLenum pname, Eigen::Matrix& matrix){ + glGetDoublev(pname,matrix.data()); } - -inline void glGetDirectv(const GLenum pname, Matrix44f & m){ - glGetFloatv(pname,m.V()); +inline void glGetv(const GLenum pname, Eigen::Matrix& matrix){ + Eigen::Matrix4f tmp; + glGetFloatv(pname,tmp.data()); + matrix = tmp; } - -inline void glGetDirecv(const GLenum pname, Matrix44d & m){ - glGetDoublev(pname,m.V()); +inline void glGetv(const GLenum pname, Eigen::Matrix& matrix){ + Eigen::Matrix4d tmp; + glGetDoublev(pname,tmp.data()); + matrix = tmp; } - +template +inline void glGetv(const GLenum pname, const Eigen::Transform& t) +{ glGetv(pname,t.matrix()); } }//namespace #endif + +#endif diff --git a/wrap/gui/trackball.cpp b/wrap/gui/trackball.cpp index ce56ec7a..8f131315 100644 --- a/wrap/gui/trackball.cpp +++ b/wrap/gui/trackball.cpp @@ -8,7 +8,7 @@ * \ * * All rights reserved. * * * -* This program is free software; you can redistribute it and/or modify * +* 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. * @@ -124,12 +124,12 @@ Transform::Transform() { } Trackball::Trackball(): current_button(0), current_mode(NULL), inactive_mode(NULL), - dragging(false), spinnable(true), spinning(false), + dragging(false), spinnable(true), spinning(false), history_size(10){ setDefaultMapping (); } -Trackball::~Trackball() +Trackball::~Trackball() { std::map::iterator it; for(it = modes.begin(); it != modes.end(); it++) @@ -167,7 +167,7 @@ void Trackball::GetView() { } // the drawing code has been moved to the trackmodes -void Trackball::DrawPostApply() { +void Trackball::DrawPostApply() { if(current_mode !=NULL){ current_mode->Draw(this); }else{ @@ -177,12 +177,12 @@ void Trackball::DrawPostApply() { } void Trackball::Apply () { - glTranslate (center); + glTranslate (center);track.Matrix().print(); glMultMatrix (track.Matrix()); glTranslate (-center); } -void Trackball::Apply(bool ToDraw) { +void Trackball::Apply(bool ToDraw) { Apply(); if(ToDraw){ DrawPostApply(); @@ -190,7 +190,7 @@ void Trackball::Apply(bool ToDraw) { } -void Trackball::ApplyInverse() { +void Trackball::ApplyInverse() { glTranslate(center); glMultMatrix(track.InverseMatrix()); glTranslate(-center); @@ -242,15 +242,15 @@ void Trackball::DrawPlane() { void Trackball::ToAscii(char* result){ float * f = (float*) &track; - sprintf(result, "trackball(%f,%f,%f,%f,%f,%f,%f,%f,%f)", + sprintf(result, "trackball(%f,%f,%f,%f,%f,%f,%f,%f,%f)", f[0],f[1],f[2],f[3],f[4],f[5],f[6],f[7],f[8] ); } bool Trackball::SetFromAscii(const char * st){ float * f = (float*) &track; - int res= sscanf(st, "trackball(%f,%f,%f,%f,%f,%f,%f,%f,%f)", + int res= sscanf(st, "trackball(%f,%f,%f,%f,%f,%f,%f,%f,%f)", f+0,f+1,f+2,f+3,f+4,f+5,f+6,f+7,f+8 ); - + return (res==9); } @@ -279,7 +279,7 @@ void Trackball::DrawPlaneHandle() { void Trackball::DrawIcon() { glPushMatrix(); - + glScale(radius); /// Here start the real drawing stuff float amb[4] ={.3f,.3f,.3f,1.0f}; @@ -287,31 +287,31 @@ void Trackball::DrawIcon() { //float col2[4]={.9f,.9f,1.0f,1.0f}; glPushAttrib(GL_ENABLE_BIT | GL_LINE_BIT | GL_CURRENT_BIT | GL_LIGHTING_BIT); - + if(current_mode == NULL ) glLineWidth(DH.LineWidthStill); else glLineWidth(DH.LineWidthMoving); - + glEnable(GL_LIGHTING); glEnable(GL_LIGHT0); glEnable(GL_LINE_SMOOTH); glEnable(GL_BLEND); glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); glColor(DH.color); - + glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,amb); glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,col); glPushMatrix(); DrawCircle(); glPushMatrix(); - + glRotatef(90,1,0,0); DrawCircle(); glRotatef(90,0,1,0); DrawCircle(); - + glPopMatrix(); glPopMatrix(); - + //glColor4f(1.0,.8f,.8f,1.0f); glPopAttrib(); @@ -336,36 +336,36 @@ void Trackball::Reset() { //interface void Trackball::MouseDown(int button) { undo_track = track; - current_button |= button; + current_button |= button; SetCurrentAction(); Hits.clear(); } void Trackball::MouseDown(int x, int y, int button) { undo_track = track; - current_button |= button; + current_button |= button; SetCurrentAction(); last_point = Point3f((float)x, (float)y, 0); Hits.clear(); } -void Trackball::MouseMove(int x, int y) { - if(current_mode == NULL) return; +void Trackball::MouseMove(int x, int y) { + if(current_mode == NULL) return; if(last_point[2] == -1) { //changed mode in the middle of moving last_point = Point3f((float)x, (float)y, 0); return; } undo_track = track; current_mode->Apply(this, Point3f(float(x), float(y), 0)); -} +} -void Trackball::MouseUp(int /* x */, int /* y */, int button) { +void Trackball::MouseUp(int /* x */, int /* y */, int button) { undo_track = track; current_button &= (~button); SetCurrentAction(); -} +} // it assumes that a notch of 1.0 is a single step of the wheel -void Trackball::MouseWheel(float notch) +void Trackball::MouseWheel(float notch) { undo_track = track; int buttons = current_button; @@ -375,7 +375,7 @@ void Trackball::MouseWheel(float notch) { ScaleMode scalemode; scalemode.Apply (this, notch); - } + } else { current_mode->Apply(this, notch); @@ -405,7 +405,7 @@ void Trackball::ButtonDown(Trackball::Button button) { if ( ( modes.count (current_button) ) && ( modes[current_button] != NULL ) ) { old_sticky = modes[current_button]->isSticky(); } - current_button |= button; + current_button |= button; if ( ( modes.count (current_button) ) && ( modes[current_button] != NULL ) ) { new_sticky = modes[current_button]->isSticky(); } @@ -414,13 +414,13 @@ void Trackball::ButtonDown(Trackball::Button button) { SetCurrentAction(); } -void Trackball::ButtonUp(Trackball::Button button) { +void Trackball::ButtonUp(Trackball::Button button) { bool old_sticky=false, new_sticky=false; assert ( modes.count (0) ); if ( ( modes.count (current_button) ) && ( modes[current_button] != NULL ) ) { old_sticky = modes[current_button]->isSticky(); } - current_button &= (~button); + current_button &= (~button); if ( ( modes.count (current_button) ) && ( modes[current_button] != NULL ) ) { new_sticky = modes[current_button]->isSticky(); } @@ -440,12 +440,12 @@ void Trackball::Undo(){ void Trackball::SetSpinnable(bool /* on*/ ){} bool Trackball::IsSpinnable() { return spinnable; -} +} void Trackball::SetSpinning(Quaternionf &/* spin*/){} void Trackball::StopSpinning(){} bool Trackball::IsSpinning() { return spinning; -} +} //navigation interface: void Trackball::Back(){} @@ -469,14 +469,14 @@ void Trackball::SetCurrentAction () } ////return center of trackball in Window coordinates. -//Point3f Trackball::ScreenOrigin() { -// return camera.Project(ModelOrigin()); +//Point3f Trackball::ScreenOrigin() { +// return camera.Project(ModelOrigin()); //} //return center of trackball in Model coordinates //Point3f Trackball::ModelOrigin() { -// return center; +// return center; //} //Matrix44f Trackball::ScreenToModel() {