diff --git a/vcg/math/eigen_vcgaddons.h b/vcg/math/eigen_vcgaddons.h index 9be8b7e6..fa520472 100644 --- a/vcg/math/eigen_vcgaddons.h +++ b/vcg/math/eigen_vcgaddons.h @@ -24,46 +24,14 @@ #warning You are including deprecated math stuff /*! * \deprecated use cols() -* Number of columns */ -EIGEN_DEPRECATED inline unsigned int ColumnsNumber() const -{ - return cols(); -}; +EIGEN_DEPRECATED inline unsigned int ColumnsNumber() const { return cols(); }; /*! * \deprecated use rows() -* Number of rows */ -EIGEN_DEPRECATED inline unsigned int RowsNumber() const -{ - return rows(); -}; - -/* -* \deprecated use *this.isApprox(m) or *this.cwise() == m -* Equality operator. -* \param m -* \return true iff the matrices have same size and its elements have same values. -*/ -// template -// EIGEN_DEPRECATED bool operator==(const MatrixBase &m) const -// { -// return (this->cwise() == m); -// } - -/* -* \deprecated use !*this.isApprox(m) or *this.cwise() != m -* Inequality operator -* \param m -* \return true iff the matrices have different size or if their elements have different values. -*/ -// template -// EIGEN_DEPRECATED bool operator!=(const MatrixBase &m) const -// { -// return (this->cwise() != m); -// }; +EIGEN_DEPRECATED inline unsigned int RowsNumber() const { return rows(); }; /*! * \deprecated use *this(i,j) (or *this.coeff(i,j)) @@ -72,25 +40,15 @@ EIGEN_DEPRECATED inline unsigned int RowsNumber() const * \param j the column index * \return the element */ -EIGEN_DEPRECATED inline Scalar ElementAt(unsigned int i, unsigned int j) const -{ - return (*this)(i,j); -} - -EIGEN_DEPRECATED inline Scalar& ElementAt(unsigned int i, unsigned int j) -{ - return (*this)(i,j); -} +EIGEN_DEPRECATED inline Scalar ElementAt(unsigned int i, unsigned int j) const { return (*this)(i,j); } +EIGEN_DEPRECATED inline Scalar& ElementAt(unsigned int i, unsigned int j) { return (*this)(i,j); } /*! * \deprecated use *this.determinant() (or *this.lu().determinant() for large matrices) * Calculate and return the matrix determinant (Laplace) * \return the matrix determinant */ -EIGEN_DEPRECATED Scalar Determinant() const -{ - return determinant(); -}; +EIGEN_DEPRECATED Scalar Determinant() const { return determinant(); }; /*! * Return the cofactor Ai,j of the ai,j element @@ -103,53 +61,23 @@ EIGEN_DEPRECATED Scalar Cofactor(unsigned int i, unsigned int j) const return (((i+j)%2==0) ? 1. : -1.) * minor(i,j).determinant(); }; -/*! -* \deprecated use *this.col(j) -* Get the j-th column on the matrix. -* \param j the column index. -* \return the reference to the column elements. This pointer must be deallocated by the caller. -*/ -EIGEN_DEPRECATED ColXpr GetColumn(const unsigned int j) -{ - return col(j); -}; +/*! \deprecated use *this.col(j) */ +EIGEN_DEPRECATED ColXpr GetColumn(const unsigned int j) { return col(j); }; -/*! -* \deprecated use *this.row(i) -* Get the i-th row on the matrix. -* \param i the column index. -* \return the reference to the row elements. This pointer must be deallocated by the caller. -*/ -EIGEN_DEPRECATED RowXpr GetRow(const unsigned int i) -{ - return row(i); -}; +/*! \deprecated use *this.row(i) */ +EIGEN_DEPRECATED RowXpr GetRow(const unsigned int i) { return row(i); }; -/*! -* \deprecated use m1.col(i).swap(m1.col(j)); -* Swaps the values of the elements between the i-th and the j-th column. -* \param i the index of the first column -* \param j the index of the second column -*/ +/*! \deprecated use m1.col(i).swap(m1.col(j)); */ EIGEN_DEPRECATED void SwapColumns(const unsigned int i, const unsigned int j) { - if (i==j) - return; - + if (i==j) return; col(i).swap(col(j)); }; -/*! -* \deprecated use m1.col(i).swap(m1.col(j)) -* Swaps the values of the elements between the i-th and the j-th row. -* \param i the index of the first row -* \param j the index of the second row -*/ +/*! \deprecated use m1.col(i).swap(m1.col(j)) */ EIGEN_DEPRECATED void SwapRows(const unsigned int i, const unsigned int j) { - if (i==j) - return; - + if (i==j) return; row(i).swap(row(j)); }; @@ -198,94 +126,44 @@ EIGEN_DEPRECATED Derived& operator-=(const Scalar k) // }; -/*! -* \deprecated use *this = a * b.transpose() -* Matrix from outer product. -*/ +/*! \deprecated use *this = a * b.transpose() (or *this = a * b.adjoint() for complexes) */ template EIGEN_DEPRECATED void OuterProduct(const MatrixBase& a, const MatrixBase& b) -{ - *this = a * b.transpose(); -} +{ *this = a * b.adjoint(); } typedef CwiseUnaryOp, Derived> ScalarAddReturnType; -/*! -* \deprecated use *this.cwise() + k -* Scalar sum. -* \param k -* \return the resultant matrix -*/ +/*! \deprecated use *this.cwise() + k */ EIGEN_DEPRECATED const ScalarAddReturnType operator+(const Scalar k) { return cwise() + k; } -/*! -* \deprecated use *this.cwise() - k -* Scalar difference. -* \param k -* \return the resultant matrix -*/ +/*! \deprecated use *this.cwise() - k */ EIGEN_DEPRECATED const ScalarAddReturnType operator-(const Scalar k) { return cwise() - k; } +/*! \deprecated use *this.setZero() or *this = MatrixType::Zero(rows,cols), etc. */ +EIGEN_DEPRECATED void SetZero() { setZero(); }; -/*! -* \deprecated use *this.setZero() or *this = MatrixType::Zero(rows,cols), etc. -* Set all the matrix elements to zero. -*/ -EIGEN_DEPRECATED void SetZero() -{ - setZero(); -}; +/*! \deprecated use *this.setIdentity() or *this = MatrixType::Identity(rows,cols), etc. */ +EIGEN_DEPRECATED void SetIdentity() { setIdentity(); }; -/*! -* \deprecated use *this.setIdentity() or *this = MatrixType::Identity(rows,cols), etc. -* Set the matrix to identity. -*/ -EIGEN_DEPRECATED void SetIdentity() -{ - setIdentity(); -}; - -/*! -* \deprecated use *this.col(j) = expression -* Set the values of j-th column to v[j] -* \param j the column index -* \param v ... -*/ +/*! \deprecated use *this.col(j) = expression */ EIGEN_DEPRECATED void SetColumn(unsigned int j, Scalar* v) -{ - col(j) = Map >(v,cols(),1); -}; +{ col(j) = Map >(v,cols(),1); }; /** \deprecated use *this.col(i) = other */ template EIGEN_DEPRECATED void SetColumn(unsigned int j, const MatrixBase& other) -{ - col(j) = other; -}; +{ col(j) = other; }; -/*! -* \deprecated use *this.row(i) = expression -* Set the elements of the i-th row to v[j] -* \param i the row index -* \param v ... -*/ +/*! \deprecated use *this.row(i) = expression */ EIGEN_DEPRECATED void SetRow(unsigned int i, Scalar* v) -{ - row(i) = Map >(v,1,rows()); -}; +{ row(i) = Map >(v,1,rows()); }; /** \deprecated use *this.row(i) = other */ template EIGEN_DEPRECATED void SetRow(unsigned int j, const MatrixBase& other) -{ - row(j) = other; -}; +{ row(j) = other; }; -/*! -* \deprecated use *this.diagonal() = expression -* Set the diagonal elements vi,i to v[i] -* \param v -*/ +/*! \deprecated use *this.diagonal() = expression */ EIGEN_DEPRECATED void SetDiagonal(Scalar *v) { assert(rows() == cols()); @@ -295,20 +173,7 @@ EIGEN_DEPRECATED void SetDiagonal(Scalar *v) /** \deprecated use trace() */ EIGEN_DEPRECATED Scalar Trace() const { return trace(); } -/*! -* \deprecated use *this = *this.transpose() -*/ -// Transpose already exist -// EIGEN_DEPRECATED void Transpose() -// { -// assert(0 && "dangerous use of deprecated Transpose function, please use: m = m.transpose();"); -// }; - - -/*! -* \deprecated use ostream << *this or ostream << *this.withFormat(...) -* Print all matrix elements -*/ +/*! \deprecated use ostream << *this or even ostream << *this.withFormat(...) */ EIGEN_DEPRECATED void Dump() { unsigned int i, j; diff --git a/vcg/math/matrix44.h b/vcg/math/matrix44.h index 5f8fe529..3a43117e 100644 --- a/vcg/math/matrix44.h +++ b/vcg/math/matrix44.h @@ -54,13 +54,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 @@ -72,7 +72,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. @@ -99,9 +99,6 @@ public: VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Matrix44) - /** \name Constructors - * No automatic casting and default constructor is empty - */ Matrix44() : Base() {} ~Matrix44() {} Matrix44(const Matrix44 &m) : Base(m) {} @@ -113,13 +110,12 @@ public: const typename Base::RowXpr operator[](int i) const { return Base::row(i); } typename Base::RowXpr operator[](int i) { return Base::row(i); } - // return a copy of the i-th column typename Base::ColXpr GetColumn4(const int& i) const { return Base::col(i); } const Eigen::Block GetColumn3(const int& i) const { return this->template block<3,1>(0,i); } typename Base::RowXpr GetRow4(const int& i) const { return Base::row(i); } Eigen::Block GetRow3(const int& i) const { return this->template block<1,3>(i,0); } - + template void ToMatrix(Matrix44Type & m) const { m = (*this).template cast(); } @@ -127,7 +123,7 @@ public: template void FromMatrix(const Matrix44Type & m) { for(int i = 0; i < 16; i++) Base::data()[i] = m.data()[i]; } - + void FromEulerAngles(Scalar alpha, Scalar beta, Scalar gamma); void SetDiagonal(const Scalar k); Matrix44 &SetScale(const Scalar sx, const Scalar sy, const Scalar sz); @@ -143,10 +139,10 @@ public: Matrix44 &SetRotateRad(Scalar AngleRad, const Point3 & axis); template void Import(const Matrix44 &m) { - for(int i = 0; i < 16; i++) + for(int i = 0; i < 16; i++) Base::data()[i] = (Scalar)(m.data()[i]); } - template + template static inline Matrix44 Construct( const Matrix44 & b ) { Matrix44 tmp; tmp.FromMatrix(b); @@ -161,7 +157,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: @@ -175,12 +171,12 @@ protected: /*** Postmultiply */ //template Point3 operator*(const Point3 &p, const Matrix44 &m); -///Premultiply +///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); +template Matrix44 &Invert(Matrix44 &m); +template Matrix44 Inverse(const Matrix44 &m); typedef Matrix44 Matrix44s; typedef Matrix44 Matrix44i; @@ -196,7 +192,7 @@ typedef Matrix44 Matrix44d; //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 ) { @@ -213,7 +209,7 @@ void Matrix44::ToEulerAngles(Scalar &alpha, Scalar &beta, Scalar &gamma) gamma = atan2(coeff(0,1), coeff(1,1)); } -template +template void Matrix44::FromEulerAngles(Scalar alpha, Scalar beta, Scalar gamma) { this->SetZero(); @@ -225,18 +221,18 @@ void Matrix44::FromEulerAngles(Scalar alpha, Scalar beta, Scalar 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; } @@ -281,7 +277,7 @@ template Matrix44 &Matrix44::SetRotateRad(Scalar AngleRad, const //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; @@ -304,7 +300,7 @@ template Matrix44 &Matrix44::SetRotateRad(Scalar AngleRad, const } /* Shear Matrixes - XY + XY 1 k 0 0 x x+ky 0 1 0 0 y y 0 0 1 0 z z @@ -327,13 +323,13 @@ template Matrix44 &Matrix44::SetRotateRad(Scalar AngleRad, const ElementAt(0,1) = sh; return *this; } - + template Matrix44 & Matrix44::SetShearXZ( const Scalar sh) {// shear the X coordinate as the Z coordinate change Base::setIdentity(); ElementAt(0,2) = sh; return *this; } - + template Matrix44 &Matrix44::SetShearYZ( const Scalar sh) {// shear the Y coordinate as the Z coordinate change Base::setIdentity(); ElementAt(1,2) = sh; @@ -343,7 +339,7 @@ template Matrix44 &Matrix44::SetRotateRad(Scalar AngleRad, const /* 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. @@ -353,7 +349,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: @@ -373,7 +369,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); @@ -388,77 +384,76 @@ 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 + 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].dot(M.GetColumn3(1)); // xy shearing R[1]= M.GetColumn3(1)-R[0]*ShearV[0]; assert(math::Abs(R[1].dot(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].dot(M.GetColumn3(2)); // xz shearing R[2]= M.GetColumn3(2)-R[0]*ShearV[1]; assert(math::Abs(R[2].dot(R[0]))<1e-10); - + R[2] = R[2]-R[1]*(R[2].dot(R[1])); assert(math::Abs(R[2].dot(R[1]))<1e-10); assert(math::Abs(R[2].dot(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].dot(R[1]))<1e-10); assert(math::Abs(R[2].dot(R[0]))<1e-10); - + ShearV[2]=R[1].dot(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]; + 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; + { + 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; + { + alpha=asin(-M(1,0)); + if(M(1,1)<0) alpha=M_PI-alpha; gamma=0; } @@ -495,15 +490,15 @@ template Point3 operator*(const Matrix44 &m, const Point3 &p) /* - 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) { return m = m.lu().inverse();