diff --git a/vcg/complex/trimesh/crease_cut.h b/vcg/complex/trimesh/crease_cut.h index e81b238f..65823559 100644 --- a/vcg/complex/trimesh/crease_cut.h +++ b/vcg/complex/trimesh/crease_cut.h @@ -91,7 +91,7 @@ void CreaseCut(MESH_TYPE &m, float angleRad) if(!isBorderVertex) // for internal vertex we search the first crease and start from it { do { - ScalarType dotProd = iPos.FFlip()->cN() * iPos.f->N(); + ScalarType dotProd = iPos.FFlip()->cN().dot(iPos.f->N()); iPos.NextFE(); if(dotProdcN() * iPos.f->N(); // test normal with the next face (fflip) + ScalarType dotProd=iPos.FFlip()->cN().dot(iPos.f->N()); // test normal with the next face (fflip) size_t faceInd = Index(m,iPos.f); indVec[faceInd*3+ iPos.VInd()] = curVertexCounter; diff --git a/vcg/complex/trimesh/create/ball_pivoting.h b/vcg/complex/trimesh/create/ball_pivoting.h index c0ca9765..6991b20b 100644 --- a/vcg/complex/trimesh/create/ball_pivoting.h +++ b/vcg/complex/trimesh/create/ball_pivoting.h @@ -147,7 +147,7 @@ template class BallPivoting: public AdvancingFront { if(d0 < min_edge || d0 > max_edge) continue; Point3x normal = (p1 - p0)^(p2 - p0); - if(normal * (p0 - baricenter) < 0) continue; + if(normal.dot(p0 - baricenter) < 0) continue; /* if(use_normals) { if(normal * vv0->N() < 0) continue; if(normal * vv1->N() < 0) continue; @@ -169,7 +169,7 @@ template class BallPivoting: public AdvancingFront { } //check on the other side there is not a surface - Point3x opposite = center + normal*(((center - p0)*normal)*2/normal.SquaredNorm()); + Point3x opposite = center + normal*(((center - p0).dot(normal))*2/normal.SquaredNorm()); for(t = 0; t < n; t++) { VertexType &v = *(targets[t]); if((v.IsV()) && (opposite - v.P()).Norm() <= radius) @@ -291,7 +291,7 @@ template class BallPivoting: public AdvancingFront { if(min_angle >= M_PI - 0.1) { return -1; } - + if(candidate == NULL) { return -1; } @@ -304,7 +304,7 @@ template class BallPivoting: public AdvancingFront { assert(id != edge.v0 && id != edge.v1); Point3x newnormal = ((candidate->P() - v0)^(v1 - v0)).Normalize(); - if(normal * newnormal < max_angle || this->nb[id] >= 2) { + if(normal.dot(newnormal) < max_angle || this->nb[id] >= 2) { return -1; } @@ -315,7 +315,7 @@ template class BallPivoting: public AdvancingFront { if((*k).v0 == id) touch = k; //mark vertices close to candidate - Mark(candidate); + Mark(candidate); return id; } @@ -360,9 +360,9 @@ template class BallPivoting: public AdvancingFront { up /= uplen; - ScalarType a11 = q1*q1; - ScalarType a12 = q1*q2; - ScalarType a22 = q2*q2; + ScalarType a11 = q1.dot(q1); + ScalarType a12 = q1.dot(q2); + ScalarType a22 = q2.dot(q2); ScalarType m = 4*(a11*a22 - a12*a12); ScalarType l1 = 2*(a11*a22 - a22*a12)/m; @@ -385,8 +385,8 @@ template class BallPivoting: public AdvancingFront { p.Normalize(); q.Normalize(); Point3x vec = p^q; - ScalarType angle = acos(p*q); - if(vec*axis < 0) angle = -angle; + ScalarType angle = acos(p.dot(q)); + if(vec.dot(axis) < 0) angle = -angle; if(angle < 0) angle += 2*M_PI; return angle; } diff --git a/vcg/complex/trimesh/inertia.h b/vcg/complex/trimesh/inertia.h index 4907cee7..97183b3d 100644 --- a/vcg/complex/trimesh/inertia.h +++ b/vcg/complex/trimesh/inertia.h @@ -334,7 +334,7 @@ static void Covariance(const MeshType & m, vcg::Point3 & bary, vcg:: // integral of (x,y,0) in the same triangle CoordType X(1/6.0,1/6.0,0); vcg::Matrix33 A, // matrix that bring the vertices to (v1-v0,v2-v0,n) - At,DC; + DC; for(fi = m.face.begin(); fi != m.face.end(); ++fi) if(!(*fi).IsD()) { @@ -350,18 +350,15 @@ static void Covariance(const MeshType & m, vcg::Point3 & bary, vcg:: A.SetColumn(2,n); CoordType delta = P0 - bary; - At= A; - At.Transpose(); /* DC is calculated as integral of (A*x+delta) * (A*x+delta)^T over the triangle, where delta = v0-bary */ DC.SetZero(); - DC+= A*C0*At; + DC+= A*C0*A.transpose(); vcg::Matrix33 tmp; tmp.OuterProduct(A*X,delta); - DC+= tmp; - tmp.Transpose(); + DC += tmp + tmp.transpose(); DC+= tmp; tmp.OuterProduct(delta,delta); DC+=tmp*0.5; diff --git a/vcg/complex/trimesh/smooth.h b/vcg/complex/trimesh/smooth.h index 2df36f6c..73158c5c 100644 --- a/vcg/complex/trimesh/smooth.h +++ b/vcg/complex/trimesh/smooth.h @@ -488,7 +488,7 @@ static void VertexCoordLaplacianQuality(MeshType &m, int step, bool SmoothSelect /* Improved Laplacian Smoothing of Noisy Surface Meshes - J. Vollmer, R. Mencl, and H. Müller + J. Vollmer, R. Mencl, and H. M�ller EUROGRAPHICS Volume 18 (1999), Number 3 */ @@ -1099,7 +1099,7 @@ static void FaceNormalAngleThreshold(MeshType &m, { if(! (*ep.f).IsV() ) { - ScalarType cosang=ep.f->N()*(*fi).N(); + ScalarType cosang=ep.f->N().dot((*fi).N()); if(cosang >= sigma) { ScalarType w = cosang-sigma; @@ -1223,7 +1223,7 @@ static void FastFitMesh(MeshType &m, for (;!ep.End();++ep) { CoordType bc=Barycenter(*ep.F()); - Sum += ep.F()->N()*(ep.F()->N()*(bc - (*vi).P())); + Sum += ep.F()->N()*(ep.F()->N().dot(bc - (*vi).P())); ++cnt; } TDV[*vi].np=(*vi).P()+ Sum*(1.0/cnt); diff --git a/vcg/complex/trimesh/update/curvature.h b/vcg/complex/trimesh/update/curvature.h index d1ceb397..ac0778a2 100644 --- a/vcg/complex/trimesh/update/curvature.h +++ b/vcg/complex/trimesh/update/curvature.h @@ -211,11 +211,8 @@ public: tempMatrix.ExternalProduct(W,W); Q -= tempMatrix * 2.0f; - Matrix33 Qt(Q); - Qt.Transpose(); - // compute matrix Q^t M Q - Matrix33 QtMQ = (Qt * M * Q); + Matrix33 QtMQ = (Q.transpose() * M * Q); CoordType bl = Q.GetColumn(0); CoordType T1 = Q.GetColumn(1); @@ -644,7 +641,7 @@ public: int n_rot; Jacobi(m33,lambda, vect,n_rot); - vect.Transpose(); + vect = vect.transpose().eval(); ScalarType normal = std::numeric_limits::min(); int normI = 0; for(int i = 0; i < 3; ++i) diff --git a/vcg/complex/trimesh/update/edges.h b/vcg/complex/trimesh/update/edges.h index d0c6a93c..7c772ba4 100644 --- a/vcg/complex/trimesh/update/edges.h +++ b/vcg/complex/trimesh/update/edges.h @@ -74,7 +74,7 @@ namespace tri { f.Edge(2) = f.V(0)->P(); f.Edge(2) -= f.V(2)->P(); // Calcolo di plane f.Plane().SetDirection(f.Edge(0)^f.Edge(1)); - f.Plane().SetOffset(f.Plane().Direction() * f.V(0)->P()); + f.Plane().SetOffset(f.Plane().Direction().dot(f.V(0)->P())); f.Plane().Normalize(); // Calcolo migliore proiezione ScalarType nx = math::Abs(f.Plane().Direction()[0]); diff --git a/vcg/math/deprecated_matrix.h b/vcg/math/deprecated_matrix.h index 2113cdd9..b799fba5 100644 --- a/vcg/math/deprecated_matrix.h +++ b/vcg/math/deprecated_matrix.h @@ -721,12 +721,15 @@ namespace vcg{ delete []temp; }; + // for the transistion to eigen Matrix transpose() { Matrix res = *this; res.Transpose(); return res; } + // for the transistion to eigen + const Matrix& eval() { return *this; } /*! diff --git a/vcg/math/deprecated_matrix33.h b/vcg/math/deprecated_matrix33.h index 9a8d3b1a..77e6576f 100644 --- a/vcg/math/deprecated_matrix33.h +++ b/vcg/math/deprecated_matrix33.h @@ -403,6 +403,8 @@ public: res.Transpose(); return res; } + // for the transistion to eigen + const Matrix33& eval() { return *this; } /// Funzione per costruire una matrice diagonale dati i tre elem. Matrix33 & SetDiagonal(S *v) diff --git a/vcg/math/deprecated_matrix44.h b/vcg/math/deprecated_matrix44.h index 037846f8..e0c55fda 100644 --- a/vcg/math/deprecated_matrix44.h +++ b/vcg/math/deprecated_matrix44.h @@ -281,6 +281,16 @@ public: return tmp; } + // for the transistion to eigen + Matrix44 transpose() const + { + Matrix44 res = *this; + Transpose(res); + return res; + } + // for the transistion to eigen + const Matrix44& eval() { return *this; } + }; @@ -289,7 +299,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: diff --git a/vcg/math/eigen.h b/vcg/math/eigen.h index 23d47fd2..34a81b9b 100644 --- a/vcg/math/eigen.h +++ b/vcg/math/eigen.h @@ -26,8 +26,40 @@ #define EIGEN_MATRIXBASE_PLUGIN +#include "../Eigen/LU" +#include "../Eigen/Geometry" #include "../Eigen/Array" #include "../Eigen/Core" +#include "base.h" + +namespace Eigen { +template<> struct NumTraits +{ + typedef unsigned char Real; + typedef float FloatingPoint; + enum { + IsComplex = 0, + HasFloatingPoint = 0, + ReadCost = 1, + AddCost = 1, + MulCost = 1 + }; +}; + +template<> struct NumTraits +{ + typedef short int Real; + typedef float FloatingPoint; + enum { + IsComplex = 0, + HasFloatingPoint = 0, + ReadCost = 1, + AddCost = 1, + MulCost = 1 + }; +}; + +} #define VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATOR(Derived, Op) \ template \ @@ -54,4 +86,46 @@ VCG_EIGEN_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, *=) \ VCG_EIGEN_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, /=) + +namespace vcg { + +template +typename Eigen::ei_traits::Scalar +Angle(const Eigen::MatrixBase& p1, const Eigen::MatrixBase & p2) +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived1) + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived2) + EIGEN_STATIC_ASSERT_FIXED_SIZE(Derived1) + EIGEN_STATIC_ASSERT_FIXED_SIZE(Derived2) + EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(Derived1,Derived2) + typedef typename Eigen::ei_traits::Scalar Scalar; + + Scalar w = p1.norm()*p2.norm(); + if(w==0) return Scalar(-1); + Scalar t = (p1.dot(p2))/w; + if(t>1) t = 1; + else if(t<-1) t = -1; + return vcg::math::Acos(t); +} + +template +inline typename Eigen::ei_traits::Scalar Norm( const Eigen::MatrixBase& p) +{ return p.norm(); } + +template +inline typename Eigen::ei_traits::Scalar SquaredNorm( const Eigen::MatrixBase& p) +{ return p.norm2(); } + +template +inline typename Eigen::ei_traits::Scalar +Distance(const Eigen::MatrixBase& p1, const Eigen::MatrixBase & p2) +{ return (p1-p2).norm(); } + +template +inline typename Eigen::ei_traits::Scalar +SquaredDistance(const Eigen::MatrixBase& p1, const Eigen::MatrixBase & p2) +{ return (p1-p2).norm2(); } + +} + #endif diff --git a/vcg/math/eigen_vcgaddons.h b/vcg/math/eigen_vcgaddons.h index 3aa2cffe..f9211b25 100644 --- a/vcg/math/eigen_vcgaddons.h +++ b/vcg/math/eigen_vcgaddons.h @@ -72,10 +72,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) +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) @@ -148,6 +153,8 @@ EIGEN_DEPRECATED void SwapRows(const unsigned int i, const unsigned int j) row(i).swap(row(j)); }; +Scalar* V() { return derived().data(); } +const Scalar* V() const { return derived().data(); } /*! * \deprecated use *this.cwise() += k @@ -288,6 +295,8 @@ EIGEN_DEPRECATED void SetDiagonal(Scalar *v) diagonal() = Map >(v,cols(),1); } +/** \deprecated use trace() */ +EIGEN_DEPRECATED Scalar Trace() const { return trace(); } /*! * \deprecated use *this = *this.transpose() @@ -316,4 +325,13 @@ EIGEN_DEPRECATED void Dump() printf("\n"); } +/** \deprecated use norm() */ +EIGEN_DEPRECATED inline Scalar Norm() const { return norm(); } +/** \deprecated use squaredNorm() */ +EIGEN_DEPRECATED inline Scalar SquaredNorm() const { return norm2(); } +/** \deprecated use normalize() or normalized() */ +EIGEN_DEPRECATED inline Derived& Normalize() { normalize(); return derived(); } + +/** \deprecated use .cross(p) */ +inline EvalType operator ^ (const Derived& p ) const { return this->cross(p); } diff --git a/vcg/math/gen_normal.h b/vcg/math/gen_normal.h index aefcbf31..2fed9807 100644 --- a/vcg/math/gen_normal.h +++ b/vcg/math/gen_normal.h @@ -73,7 +73,7 @@ static void UniformCone(int vn, std::vector > &NN, ScalarTyp ScalarType DotProd = cos(AngleRad); for(vi=NNT.begin();vi!=NNT.end();++vi) { - if(dir*(*vi) >= DotProd) NN.push_back(*vi); + if(dir.dot(*vi) >= DotProd) NN.push_back(*vi); } } @@ -124,7 +124,7 @@ static int BestMatchingNormal(const Point3x &n, std::vector &nv) typename std::vector::iterator ni; for(ni=nv.begin();ni!=nv.end();++ni) { - cosang=(*ni)*n; + cosang=(*ni).dot(n); if(cosang>bestang) { bestang=cosang; ret=ni-nv.begin(); diff --git a/vcg/math/matrix.h b/vcg/math/matrix.h index a16d188e..58e2973d 100644 --- a/vcg/math/matrix.h +++ b/vcg/math/matrix.h @@ -23,8 +23,8 @@ #ifndef VCG_USE_EIGEN #include "deprecated_matrix.h" -#endif +#else #ifndef MATRIX_VCGLIB #define MATRIX_VCGLIB @@ -69,10 +69,9 @@ public: * @param Scalar (Templete Parameter) Specifies the ScalarType field. */ template -class Matrix : public Eigen::Matrix<_Scalar,Eigen::Dynamic,Eigen::Dynamic> // FIXME col or row major ? +class Matrix : public Eigen::Matrix<_Scalar,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> // FIXME col or row major ? { - - typedef Eigen::Matrix<_Scalar,Eigen::Dynamic,Eigen::Dynamic> _Base; + typedef Eigen::Matrix<_Scalar,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> _Base; public: @@ -87,7 +86,7 @@ public: * \param m the number of matrix rows * \param n the number of matrix columns */ - Matrix(unsigned int m, unsigned int n) + Matrix(int m, int n) : Base(m,n) { memset(Base::data(), 0, m*n*sizeof(Scalar)); @@ -100,10 +99,10 @@ public: * \param n the number of matrix columns * \param values the values of the matrix elements */ - Matrix(unsigned int m, unsigned int n, Scalar *values) - : Base(m.n) + Matrix(int m, int n, Scalar *values) + : Base(m,n) { - *this = Eigen::Map >(values, m , n); + *this = Eigen::Map >(values, m , n); } /*! @@ -197,3 +196,6 @@ void Invert(MatrixType & m) } // end of namespace #endif + +#endif + diff --git a/vcg/math/matrix33.h b/vcg/math/matrix33.h index fc319fad..013fdd4e 100644 --- a/vcg/math/matrix33.h +++ b/vcg/math/matrix33.h @@ -21,14 +21,15 @@ * * ****************************************************************************/ -// #ifndef VCG_USE_EIGEN +#ifndef VCG_USE_EIGEN #include "deprecated_matrix33.h" -// #endif +#else #ifndef __VCGLIB_MATRIX33_H #define __VCGLIB_MATRIX33_H #include "eigen.h" +#include "matrix44.h" namespace vcg{ template class Matrix33; @@ -48,17 +49,20 @@ namespace vcg { @param S (Templete Parameter) Specifies the ScalarType field. */ template -class Matrix33 : public Eigen::Matrix<_Scalar,3,3,RowMajor> // FIXME col or row major ? +class Matrix33 : public Eigen::Matrix<_Scalar,3,3,Eigen::RowMajor> // FIXME col or row major ? { - typedef Eigen::Matrix<_Scalar,Eigen::Dynamic,Eigen::Dynamic> _Base; + typedef Eigen::Matrix<_Scalar,3,3,Eigen::RowMajor> _Base; + using _Base::coeff; + using _Base::coeffRef; + using _Base::setZero; public: - _EIGEN_GENERIC_PUBLIC_INTERFACE(Matrix,_Base); + _EIGEN_GENERIC_PUBLIC_INTERFACE(Matrix33,_Base); typedef _Scalar ScalarType; - VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Matrix) + VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Matrix33) /// Default constructor inline Matrix33() : Base() {} @@ -67,10 +71,13 @@ public: Matrix33(const Matrix33& m ) : Base(m) {} /// create from a \b row-major array - Matrix33(const Scalar * v ) : Base(Map(v)) {} + Matrix33(const Scalar * v ) : Base(Eigen::Map >(v)) {} /// create from Matrix44 excluding row and column k - Matrix33(const Matrix44 & m, const int & k) : Base(m.minor(k,k)) {} + Matrix33(const Matrix44 & m, const int & k) : Base(m.minor(k,k)) {} + + template + Matrix33(const Eigen::MatrixBase& other) : Base(other) {} /*! * \deprecated use *this.row(i) @@ -86,13 +93,13 @@ public: /** \deprecated */ - Matrix33 & SetRotateRad(S angle, const Point3 & axis ) + Matrix33 & SetRotateRad(Scalar angle, const Point3 & axis ) { *this = Eigen::AngleAxis(angle,axis).toRotationMatrix(); return (*this); } /** \deprecated */ - Matrix33 & SetRotateDeg(S angle, const Point3 & axis ){ + Matrix33 & SetRotateDeg(Scalar angle, const Point3 & axis ){ return SetRotateRad(math::ToRad(angle),axis); } @@ -100,7 +107,7 @@ public: // 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 */ - Matrix33 & FastInvert() { return *this = inverse(); } + Matrix33 & FastInvert() { return *this = Base::inverse(); } void show(FILE * fp) { @@ -111,7 +118,7 @@ public: /* compute the matrix generated by the product of a * b^T */ - void ExternalProduct(const Point3 &a, const Point3 &b) + void ExternalProduct(const Point3 &a, const Point3 &b) { for(int i=0;i<3;++i) for(int j=0;j<3;++j) @@ -120,15 +127,15 @@ public: /* Compute the Frobenius Norm of the Matrix */ - ScalarType Norm() - { - // FIXME looks like there is a bug: j is not used !!! - ScalarType SQsum=0; - for(int i=0;i<3;++i) - for(int j=0;j<3;++j) - SQsum += a[i]*a[i]; - return (math::Sqrt(SQsum)); - } + Scalar Norm() { return Base::cwise().abs2().sum(); } +// { +// // FIXME looks like there is a bug: j is not used !!! +// Scalar SQsum=0; +// for(int i=0;i<3;++i) +// for(int j=0;j<3;++j) +// SQsum += a[i]*a[i]; +// return (math::Sqrt(SQsum)); +// } /** @@ -136,7 +143,7 @@ public: */ // FIXME should be outside Matrix template - void Covariance(const STLPOINTCONTAINER &points, Point3 &bp) { + void Covariance(const STLPOINTCONTAINER &points, Point3 &bp) { assert(!points.empty()); typedef typename STLPOINTCONTAINER::const_iterator PointIte; // first cycle: compute the barycenter @@ -147,14 +154,12 @@ public: this->setZero(); vcg::Matrix33 A; for( PointIte pi = points.begin(); pi != points.end(); ++pi) { - Point3 p = (*pi)-bp; + Point3 p = (*pi)-bp; A.OuterProduct(p,p); (*this)+= A; } } - - /** It computes the cross covariance matrix of two set of 3d points P and X; it returns also the barycenters of P and X. @@ -168,14 +173,14 @@ public: // FIXME should be outside Matrix template void CrossCovariance(const STLPOINTCONTAINER &P, const STLPOINTCONTAINER &X, - Point3 &bp, Point3 &bx) + Point3 &bp, Point3 &bx) { setZero(); assert(P.size()==X.size()); bx.setZero(); bp.setZero(); - Matrix33 tmp; - typename std::vector >::const_iterator pi,xi; + Matrix33 tmp; + typename std::vector >::const_iterator pi,xi; for(pi=P.begin(),xi=X.begin();pi!=P.end();++pi,++xi){ bp+=*pi; bx+=*xi; @@ -193,15 +198,15 @@ public: void WeightedCrossCovariance(const STLREALCONTAINER & weights, const STLPOINTCONTAINER &P, const STLPOINTCONTAINER &X, - Point3 &bp, - Point3 &bx) + Point3 &bp, + Point3 &bx) { - SetZero(); + setZero(); assert(P.size()==X.size()); bx.SetZero(); bp.SetZero(); - Matrix33 tmp; - typename std::vector >::const_iterator pi,xi; + Matrix33 tmp; + typename std::vector >::const_iterator pi,xi; typename STLREALCONTAINER::const_iterator pw; for(pi=P.begin(),xi=X.begin();pi!=P.end();++pi,++xi){ @@ -317,10 +322,12 @@ template /// typedef Matrix33 Matrix33s; -typedef Matrix33 Matrix33i; +typedef Matrix33 Matrix33i; typedef Matrix33 Matrix33f; typedef Matrix33 Matrix33d; } // end of namespace #endif + +#endif diff --git a/vcg/math/matrix44.h b/vcg/math/matrix44.h index 0b6b0e26..5f8fe529 100644 --- a/vcg/math/matrix44.h +++ b/vcg/math/matrix44.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. * @@ -21,9 +21,9 @@ * * ****************************************************************************/ -// #ifndef VCG_USE_EIGEN +#ifndef VCG_USE_EIGEN #include "deprecated_matrix44.h" -// #endif +#else #ifndef __VCGLIB_MATRIX44 #define __VCGLIB_MATRIX44 @@ -83,17 +83,21 @@ for 'column' vectors. /** This class represents a 4x4 matrix. T is the kind of element in the matrix. */ template -class Matrix44 : public Eigen::Matrix<_Scalar,4,4,RowMajor> // FIXME col or row major ! +class Matrix44 : public Eigen::Matrix<_Scalar,4,4,Eigen::RowMajor> // FIXME col or row major ! { - typedef Eigen::Matrix<_Scalar,Eigen::Dynamic,Eigen::Dynamic> _Base; + typedef Eigen::Matrix<_Scalar,4,4,Eigen::RowMajor> _Base; + using _Base::coeff; + using _Base::coeffRef; + using _Base::ElementAt; + using _Base::setZero; public: - _EIGEN_GENERIC_PUBLIC_INTERFACE(Matrix,_Base); + _EIGEN_GENERIC_PUBLIC_INTERFACE(Matrix44,_Base); typedef _Scalar ScalarType; - VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Matrix) + VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Matrix44) /** \name Constructors * No automatic casting and default constructor is empty @@ -101,50 +105,51 @@ public: Matrix44() : Base() {} ~Matrix44() {} Matrix44(const Matrix44 &m) : Base(m) {} - Matrix33(const Scalar * v ) : Base(Map(v)) {} + Matrix44(const Scalar * v ) : Base(Eigen::Map >(v)) {} - Scalar *V() { return Base::data(); } - const Scalar *V() const { return Base::data(); } + template + Matrix44(const Eigen::MatrixBase& other) : Base(other) {} + + 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 col(i); } - Block GetColumn3(const int& i) const { return block<3,1>(0,i); } + 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 col(i); } - Block GetRow3(const int& i) const { return block<1,3>(i,0); } + 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).cast(); } + void ToMatrix(Matrix44Type & m) const { m = (*this).template cast(); } void ToEulerAngles(Scalar &alpha, Scalar &beta, Scalar &gamma); template - void FromMatrix(const Matrix44Type & m) { for(int i = 0; i < 16; i++) data()[i] = m.data()[i]; } + void FromMatrix(const Matrix44Type & m) { for(int i = 0; i < 16; i++) Base::data()[i] = m.data()[i]; } - void FromEulerAngles(T alpha, T beta, T gamma); - void SetDiagonal(const T k); - 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); + void FromEulerAngles(Scalar alpha, Scalar beta, Scalar gamma); + void SetDiagonal(const Scalar k); + Matrix44 &SetScale(const Scalar sx, const Scalar sy, const Scalar sz); + Matrix44 &SetScale(const Point3 &t); + Matrix44 &SetTranslate(const Point3 &t); + Matrix44 &SetTranslate(const Scalar sx, const Scalar sy, const Scalar sz); + Matrix44 &SetShearXY(const Scalar sz); + Matrix44 &SetShearXZ(const Scalar sy); + Matrix44 &SetShearYZ(const Scalar sx); ///use radiants for angle. - Matrix44 &SetRotateDeg(T AngleDeg, const Point3 & axis); - Matrix44 &SetRotateRad(T AngleRad, const Point3 & axis); - - T Determinant() const; + 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++) - _a[i] = (T)(m.V()[i]); + Base::data()[i] = (Scalar)(m.data()[i]); } template static inline Matrix44 Construct( const Matrix44 & b ) { - Matrix44 tmp; tmp.FromMatrix(b); + Matrix44 tmp; tmp.FromMatrix(b); return tmp; } @@ -173,7 +178,6 @@ protected: ///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); @@ -184,27 +188,6 @@ typedef Matrix44 Matrix44f; typedef Matrix44 Matrix44d; - -template Matrix44::Matrix44(const Matrix44 &m) { - memcpy((T *)_a, (T *)m._a, 16 * sizeof(T)); -} - -template Matrix44::Matrix44(const T v[]) { - memcpy((T *)_a, v, 16 * sizeof(T)); -} - -template T &Matrix44::ElementAt(const int row, const int col) { - assert(row >= 0 && row < 4); - assert(col >= 0 && col < 4); - return _a[(row<<2) + col]; -} - -template T Matrix44::ElementAt(const int row, const int col) const { - assert(row >= 0 && row < 4); - assert(col >= 0 && col < 4); - return _a[(row<<2) + col]; -} - //template T &Matrix44::operator[](const int i) { // assert(i >= 0 && i < 16); // return ((T *)_a)[i]; @@ -214,116 +197,7 @@ template T Matrix44::ElementAt(const int row, const int col) 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; -} -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 Matrix44 Matrix44::operator+(const Matrix44 &m) const { - Matrix44 ret; - 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]; - return ret; -} - -template Matrix44 Matrix44::operator*(const Matrix44 &m) const { - Matrix44 ret; - for(int i = 0; i < 4; i++) - for(int j = 0; j < 4; j++) { - T t = 0.0; - for(int k = 0; k < 4; k++) - t += ElementAt(i, k) * m.ElementAt(k, j); - ret.ElementAt(i, j) = t; - } - return ret; -} - -template Matrix44 Matrix44::operator*(const Matrix44Diag &m) const { - Matrix44 ret = (*this); - for(int i = 0; i < 4; ++i) - for(int j = 0; j < 4; ++j) - ret[i][j]*=m[i]; - return ret; -} - -template Point4 Matrix44::operator*(const Point4 &v) const { - Point4 ret; - for(int i = 0; i < 4; i++){ - T t = 0.0; - for(int k = 0; k < 4; k++) - t += ElementAt(i,k) * v[k]; - ret[i] = t; - } - return ret; -} - - -template bool Matrix44::operator==(const Matrix44 &m) const { - for(int i = 0; i < 4; ++i) - 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) - if(ElementAt(i,j) != m.ElementAt(i,j)) - return true; - return false; -} - -template Matrix44 Matrix44::operator-() const { - Matrix44 res; - for(int i = 0; i < 16; i++) - res.V()[i] = -V()[i]; - return res; -} - -template Matrix44 Matrix44::operator*(const T k) const { - Matrix44 res; - for(int i = 0; i < 16; i++) - res.V()[i] =V()[i] * k; - return res; -} - -template void Matrix44::operator+=(const Matrix44 &m) { - for(int i = 0; i < 16; i++) - V()[i] += m.V()[i]; -} -template void Matrix44::operator-=(const Matrix44 &m) { - for(int i = 0; i < 16; i++) - V()[i] -= m.V()[i]; -} -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++) { - for(int j = 0; j < 4; j++) { - t[k] += ElementAt(i, k) * m.ElementAt(k, j); - } - } - for(int l = 0; l < 4; l++) - ElementAt(i, l) = t[l]; - } */ -} template < class PointType , class T > void operator*=( std::vector &vert, const Matrix44 & m ) { typename std::vector::iterator ii; @@ -331,21 +205,16 @@ template < class PointType , class T > void operator*=( std::vector & (*ii).P()=m * (*ii).P(); } -template void Matrix44::operator*=( const T k ) { - for(int i = 0; i < 16; i++) - _a[i] *= k; -} - template -void Matrix44::ToEulerAngles(T &alpha, T &beta, T &gamma) +void Matrix44::ToEulerAngles(Scalar &alpha, Scalar &beta, Scalar &gamma) { - alpha = atan2(ElementAt(1,2), ElementAt(2,2)); - beta = asin(-ElementAt(0,2)); - gamma = atan2(ElementAt(0,1), ElementAt(1,1)); + alpha = atan2(coeff(1,2), coeff(2,2)); + beta = asin(-coeff(0,2)); + gamma = atan2(coeff(0,1), coeff(1,1)); } template -void Matrix44::FromEulerAngles(T alpha, T beta, T gamma) +void Matrix44::FromEulerAngles(Scalar alpha, Scalar beta, Scalar gamma) { this->SetZero(); @@ -371,28 +240,20 @@ void Matrix44::FromEulerAngles(T alpha, T beta, T gamma) ElementAt(3,3) = 1; } -template void Matrix44::SetZero() { - memset((T *)_a, 0, 16 * sizeof(T)); -} - -template void Matrix44::SetIdentity() { - SetDiagonal(1); -} - -template void Matrix44::SetDiagonal(const T k) { - SetZero(); +template void Matrix44::SetDiagonal(const Scalar k) { + setZero(); 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) { +template Matrix44 &Matrix44::SetScale(const Point3 &t) { SetScale(t[0], t[1], t[2]); return *this; } -template Matrix44 &Matrix44::SetScale(const T sx, const T sy, const T sz) { - SetZero(); +template Matrix44 &Matrix44::SetScale(const Scalar sx, const Scalar sy, const Scalar sz) { + setZero(); ElementAt(0, 0) = sx; ElementAt(1, 1) = sy; ElementAt(2, 2) = sz; @@ -400,23 +261,23 @@ template Matrix44 &Matrix44::SetScale(const T sx, const T sy, co return *this; } -template Matrix44 &Matrix44::SetTranslate(const Point3 &t) { +template Matrix44 &Matrix44::SetTranslate(const Point3 &t) { SetTranslate(t[0], t[1], t[2]); return *this; } -template Matrix44 &Matrix44::SetTranslate(const T tx, const T ty, const T tz) { - SetIdentity(); +template Matrix44 &Matrix44::SetTranslate(const Scalar tx, const Scalar ty, const Scalar tz) { + Base::setIdentity(); ElementAt(0, 3) = tx; ElementAt(1, 3) = ty; ElementAt(2, 3) = tz; return *this; } -template Matrix44 &Matrix44::SetRotateDeg(T AngleDeg, const Point3 & axis) { +template Matrix44 &Matrix44::SetRotateDeg(Scalar AngleDeg, const Point3 & axis) { return SetRotateRad(math::ToRad(AngleDeg),axis); } -template Matrix44 &Matrix44::SetRotateRad(T AngleRad, const Point3 & axis) { +template Matrix44 &Matrix44::SetRotateRad(Scalar AngleRad, const Point3 & axis) { //angle = angle*(T)3.14159265358979323846/180; e' in radianti! T c = math::Cos(AngleRad); T s = math::Sin(AngleRad); @@ -436,9 +297,9 @@ 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; } @@ -461,20 +322,20 @@ template Matrix44 &Matrix44::SetRotateRad(T AngleRad, const Poin */ - template Matrix44 & Matrix44:: SetShearXY( const T sh) {// shear the X coordinate as the Y coordinate change - SetIdentity(); + template Matrix44 & Matrix44::SetShearXY( const Scalar sh) {// shear the X coordinate as the Y coordinate change + Base::setIdentity(); ElementAt(0,1) = sh; return *this; } - template Matrix44 & Matrix44:: SetShearXZ( const T sh) {// shear the X coordinate as the Z coordinate change - SetIdentity(); + 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 T sh) {// shear the Y coordinate as the Z coordinate change - SetIdentity(); + template Matrix44 &Matrix44::SetShearYZ( const Scalar sh) {// shear the Y coordinate as the Z coordinate change + Base::setIdentity(); ElementAt(1,2) = sh; return *this; } @@ -546,28 +407,28 @@ bool Decompose(Matrix44 &M, Point3 &ScaleV, Point3 &ShearV, Point3 & R[0]=M.GetColumn3(0); R[0].Normalize(); - ShearV[0]=R[0]*M.GetColumn3(1); // xy shearing + ShearV[0]=R[0].dot(M.GetColumn3(1)); // xy shearing R[1]= M.GetColumn3(1)-R[0]*ShearV[0]; - assert(math::Abs(R[1]*R[0])<1e-10); + assert(math::Abs(R[1].dot(R[0]))<1e-10); ScaleV[1]=Norm(R[1]); // y scaling R[1]=R[1]/ScaleV[1]; ShearV[0]=ShearV[0]/ScaleV[1]; - ShearV[1]=R[0]*M.GetColumn3(2); // xz shearing + ShearV[1]=R[0].dot(M.GetColumn3(2)); // xz shearing R[2]= M.GetColumn3(2)-R[0]*ShearV[1]; - assert(math::Abs(R[2]*R[0])<1e-10); + assert(math::Abs(R[2].dot(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); + 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]; R[2]=R[2]/ScaleV[2]; - assert(math::Abs(R[2]*R[1])<1e-10); - assert(math::Abs(R[2]*R[0])<1e-10); + assert(math::Abs(R[2].dot(R[1]))<1e-10); + assert(math::Abs(R[2].dot(R[0]))<1e-10); - ShearV[2]=R[1]*M.GetColumn3(2); // yz shearing + ShearV[2]=R[1].dot(M.GetColumn3(2)); // yz shearing ShearV[2]=ShearV[2]/ScaleV[2]; int i,j; for(i=0;i<3;++i) @@ -608,15 +469,6 @@ bool Decompose(Matrix44 &M, Point3 &ScaleV, Point3 &ShearV, Point3 & return true; } - - - -template T Matrix44::Determinant() const { - LinearSolve solve(*this); - return solve.Determinant(); -} - - template Point3 operator*(const Matrix44 &m, const Point3 &p) { T w; Point3 s; @@ -628,24 +480,19 @@ template Point3 operator*(const Matrix44 &m, const Point3 &p) return s; } -//template Point3 operator*(const Point3 &p, const Matrix44 &m) { + + +// 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; +// if(w != 0) s /= w; // return s; -//} +// } -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)); - } - return m; -} /* To invert a matrix you can @@ -658,194 +505,16 @@ template Matrix44 &Transpose(Matrix44 &m) { invertedMatrix = vcg::Inverse(untouchedMatrix); */ -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++) - m.ElementAt(i, j) = col[i]; - } - return m; +template Matrix44 & Invert(Matrix44 &m) { + return m = m.lu().inverse(); } -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++) - res.ElementAt(i, j) = col[i]; - } - return res; +template Matrix44 Inverse(const Matrix44 &m) { + return m.lu().inverse(); } - -/* LINEAR SOLVE IMPLEMENTATION */ - -template LinearSolve::LinearSolve(const Matrix44 &m): Matrix44(m) { - if(!Decompose()) { - for(int i = 0; i < 4; i++) - index[i] = i; - Matrix44::SetZero(); - } -} - - -template T LinearSolve::Determinant() const { - T det = d; - for(int j = 0; j < 4; j++) - det *= this-> ElementAt(j, j); - return det; -} - - -/*replaces a matrix by its LU decomposition of a rowwise permutation. -d is +or -1 depeneing of row permutation even or odd.*/ -#define TINY 1e-100 - -template bool LinearSolve::Decompose() { - - /* Matrix44 A; - for(int i = 0; i < 16; i++) - 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++) - scalemax = (scalemax > math::Abs(A.ElementAt(i,j))) ? scalemax : math::Abs(A.ElementAt(i,j)); - scale[i] = scalemax; - } - - // Loop over rows k = 1, ..., (N-1) - int signDet = 1; - for(int k = 0; k < 3; k++ ) { - // Select pivot row from max( |a(j,k)/s(j)| ) - T ratiomax = (T)0.0; - int jPivot = k; - for(int i = k; i < 4; i++ ) { - T ratio = math::Abs(A.ElementAt(index[i], k))/scale[index[i]]; - if(ratio > ratiomax) { - jPivot = i; - ratiomax = ratio; - } - } - // Perform pivoting using row index list - int indexJ = index[k]; - if( jPivot != k ) { // Pivot - indexJ = index[jPivot]; - index[jPivot] = index[k]; // Swap index jPivot and k - index[k] = indexJ; - signDet *= -1; // Flip sign of determinant - } - // Perform forward elimination - for(int i=k+1; i < 4; i++ ) { - T coeff = A.ElementAt(index[i],k)/A.ElementAt(indexJ,k); - for(int j=k+1; j < 4; j++ ) - A.ElementAt(index[i],j) -= coeff*A.ElementAt(indexJ,j); - A.ElementAt(index[i],k) = coeff; - //for( j=1; j< 4; j++ ) - // ElementAt(index[i],j) -= A.ElementAt(index[i], k)*ElementAt(indexJ, j); - } - } - for(int i = 0; i < 16; i++) - operator[](i) = A[i]; - - d = signDet; - // this = A; - return true; */ - - d = 1; //no permutation still - - T scaling[4]; - int i,j,k; - //Saving the scvaling information per row - for(i = 0; i < 4; i++) { - T largest = 0.0; - for(j = 0; j < 4; j++) { - T t = math::Abs(this->ElementAt(i, j)); - if (t > largest) largest = t; - } - - if (largest == 0.0) { //oooppps there is a zero row! - return false; - } - scaling[i] = (T)1.0 / largest; - } - - int imax = 0; - for(j = 0; j < 4; j++) { - for(i = 0; i < j; i++) { - T sum = this->ElementAt(i,j); - for(int k = 0; k < i; k++) - sum -= this->ElementAt(i,k)*this->ElementAt(k,j); - this->ElementAt(i,j) = sum; - } - T largest = 0.0; - for(i = j; i < 4; i++) { - T sum = this->ElementAt(i,j); - for(k = 0; k < j; k++) - sum -= this->ElementAt(i,k)*this->ElementAt(k,j); - this->ElementAt(i,j) = sum; - T t = scaling[i] * math::Abs(sum); - if(t >= largest) { - largest = t; - imax = i; - } - } - if (j != imax) { - for (int k = 0; k < 4; k++) { - T dum = this->ElementAt(imax,k); - this->ElementAt(imax,k) = this->ElementAt(j,k); - this->ElementAt(j,k) = dum; - } - d = -(d); - scaling[imax] = scaling[j]; - } - index[j]=imax; - if (this->ElementAt(j,j) == 0.0) this->ElementAt(j,j) = (T)TINY; - if (j != 3) { - T dum = (T)1.0 / (this->ElementAt(j,j)); - for (i = j+1; i < 4; i++) - this->ElementAt(i,j) *= dum; - } - } - return true; -} - - -template Point4 LinearSolve::Solve(const Point4 &b) { - Point4 x(b); - int first = -1, ip; - for(int i = 0; i < 4; i++) { - ip = index[i]; - T sum = x[ip]; - x[ip] = x[i]; - if(first!= -1) - for(int j = first; j <= i-1; j++) - sum -= this->ElementAt(i,j) * x[j]; - else - if(sum) first = i; - x[i] = sum; - } - for (int i = 3; i >= 0; i--) { - T sum = x[i]; - for (int j = i+1; j < 4; j++) - sum -= this->ElementAt(i, j) * x[j]; - x[i] = sum / this->ElementAt(i, i); - } - return x; -} - } //namespace #endif - +#endif diff --git a/vcg/math/point_matching.h b/vcg/math/point_matching.h index ec7b70fd..572afa91 100644 --- a/vcg/math/point_matching.h +++ b/vcg/math/point_matching.h @@ -137,16 +137,12 @@ bool ComputeWeightedRigidMatchMatrix(Matrix44x &res, ) { - Matrix33x tmp; Matrix33x ccm; Point3x bfix,bmov; // baricenter of src e trg ccm.WeightedCrossCovariance(weights,Pmov,Pfix,bmov,bfix); Matrix33x cyc; // the cyclic components of the cross covariance matrix. - cyc=ccm; - tmp=ccm; - tmp.Transpose(); - cyc-=tmp; + cyc=ccm - ccm.transpose(); Matrix44x QQ; QQ.SetZero(); @@ -157,9 +153,7 @@ bool ComputeWeightedRigidMatchMatrix(Matrix44x &res, RM[0][0]=-ccm.Trace(); RM[1][1]=-ccm.Trace(); RM[2][2]=-ccm.Trace(); - RM+=ccm; - ccm.Transpose(); - RM+=ccm; + RM += ccm + ccm.transpose(); QQ[0][0] = ccm.Trace(); QQ[0][1] = D[0]; QQ[0][2] = D[1]; QQ[0][3] = D[2]; @@ -208,16 +202,12 @@ bool ComputeRigidMatchMatrix(Matrix44x &res, Point3x &tr) { - Matrix33x tmp; Matrix33x ccm; Point3x bfix,bmov; // baricenter of src e trg ccm.CrossCovariance(Pmov,Pfix,bmov,bfix); Matrix33x cyc; // the cyclic components of the cross covariance matrix. - cyc=ccm; - tmp=ccm; - tmp.Transpose(); - cyc-=tmp; + cyc=ccm-ccm.transpose(); Matrix44x QQ; QQ.SetZero(); @@ -228,9 +218,7 @@ bool ComputeRigidMatchMatrix(Matrix44x &res, RM[0][0]=-ccm.Trace(); RM[1][1]=-ccm.Trace(); RM[2][2]=-ccm.Trace(); - RM+=ccm; - ccm.Transpose(); - RM+=ccm; + RM += ccm + ccm.transpose(); QQ[0][0] = ccm.Trace(); QQ[0][1] = D[0]; QQ[0][2] = D[1]; QQ[0][3] = D[2]; @@ -260,7 +248,7 @@ bool ComputeRigidMatchMatrix(Matrix44x &res, maxv=d[i]; } // The corresponding eigenvector define the searched rotation, - Matrix44x Rot; + Matrix44x Rot; q.ToMatrix(Rot); // the translation (last row) is simply the difference between the transformed src barycenter and the trg baricenter tr= (bfix - Rot*bmov); diff --git a/vcg/math/polar_decomposition.h b/vcg/math/polar_decomposition.h index 97b51ab8..beea1453 100644 --- a/vcg/math/polar_decomposition.h +++ b/vcg/math/polar_decomposition.h @@ -49,35 +49,31 @@ void RotationalPartByPolarDecomposition( const vcg::Matrix33 & m, vcg::Matrix r.SetZero(); s.SetZero(); - tmp= m; - tmp.Transpose(); - tmp = m*tmp; + tmp = m*m.transpose(); Matrix33 res; Point3 e; bool ss = SingularValueDecomposition >(tmp,&e[0],res); - res.Transpose(); e[0]=math::Sqrt(e[0]); e[1]=math::Sqrt(e[1]); e[2]=math::Sqrt(e[2]); #ifdef VCG_USE_EIGEN - tmp = tmp*e.asDiagonal()*res; + tmp = tmp*e.asDiagonal()*res.transpose(); #else - tmp = tmp*Matrix33Diag(e)*res; + tmp = tmp*Matrix33Diag(e)*res.transpose(); #endif - bool s1 = SingularValueDecomposition >(tmp,&e[0],res); - tmp.Transpose(); + bool s1 = SingularValueDecomposition >(tmp,&e[0],res.transpose()); e[0]=1/e[0]; e[1]=1/e[1]; e[2]=1/e[2]; #ifdef VCG_USE_EIGEN - tmp = res*e.asDiagonal()*tmp; + tmp = res*e.asDiagonal()*tmp.transpose(); #else - tmp = res*Matrix33Diag(e)*tmp; + tmp = res*Matrix33Diag(e)*tmp.transpose(); #endif r = m*tmp; diff --git a/vcg/math/quadric.h b/vcg/math/quadric.h index 97d2ce6b..117db161 100644 --- a/vcg/math/quadric.h +++ b/vcg/math/quadric.h @@ -88,7 +88,7 @@ template< class PlaneType > c = (ScalarType)p.Offset()*p.Offset(); } - void Zero() // Azzera la quadrica + void SetZero() // Azzera la quadrica { a[0] = 0; a[1] = 0; diff --git a/vcg/math/quaternion.h b/vcg/math/quaternion.h index 24eeca64..0c895b82 100644 --- a/vcg/math/quaternion.h +++ b/vcg/math/quaternion.h @@ -184,9 +184,9 @@ template Quaternion Quaternion::operator*(const S &s) const { template Quaternion Quaternion::operator*(const Quaternion &q) const { Point3 t1(V(1), V(2), V(3)); - Point3 t2(q.V(1), q.V(2), q.V(3)); + Point3 t2(q.V(1), q.V(2), q.V(3)); - S d = t2 * t1; + S d = t2.dot(t1); Point3 t3 = t1 ^ t2; t1 *= q.V(0); diff --git a/vcg/math/similarity.h b/vcg/math/similarity.h index 212148ff..831abbb5 100644 --- a/vcg/math/similarity.h +++ b/vcg/math/similarity.h @@ -191,7 +191,7 @@ template Matrix44 Similarity::Ma rot.ToMatrix(r); Matrix44 s = Matrix44().SetScale(sca, sca, sca); Matrix44 t = Matrix44().SetTranslate(tra[0], tra[1], tra[2]); - return s*r*t; // trans * scale * rot; + return Matrix44(s*r*t); // trans * scale * rot; } template Matrix44 Similarity::InverseMatrix() const { diff --git a/vcg/space/color4.h b/vcg/space/color4.h index 138b111b..7db71034 100644 --- a/vcg/space/color4.h +++ b/vcg/space/color4.h @@ -89,6 +89,7 @@ namespace vcg { template class Color4 : public Point4 { + typedef Point4 Base; public: /// Constant for storing standard colors. /// Each color is stored in a simple in so that the bit pattern match with the one of Color4b. @@ -121,22 +122,30 @@ public: inline Color4 ( const Point4 &c) :Point4(c) {}; inline Color4 (){}; inline Color4 (ColorConstant cc); + #ifdef VCG_USE_EIGEN + template + inline Color4 (const Eigen::MatrixBase& other) : Base(other) + { + // TODO make sure the types are the same + } + #endif + template inline void Import(const Color4 & b ) { - Point4::_v[0] = T(b[0]); - Point4::_v[1] = T(b[1]); - Point4::_v[2] = T(b[2]); - Point4::_v[3] = T(b[3]); + Point4::V()[0] = T(b[0]); + Point4::V()[1] = T(b[1]); + Point4::V()[2] = T(b[2]); + Point4::V()[3] = T(b[3]); } template inline void Import(const Point4 & b ) { - Point4::_v[0] = T(b[0]); - Point4::_v[1] = T(b[1]); - Point4::_v[2] = T(b[2]); - Point4::_v[3] = T(b[3]); + Point4::V()[0] = T(b[0]); + Point4::V()[1] = T(b[1]); + Point4::V()[2] = T(b[2]); + Point4::V()[3] = T(b[3]); } template @@ -150,7 +159,7 @@ public: inline Color4 operator + ( const Color4 & p) const { - return Color4( this->_v[0]+p._v[0], this->_v[1]+p._v[1], this->_v[2]+p._v[2], this->_v[3]+p._v[3] ); + return Color4( (*this)[0]+p.V()[0], (*this)[1]+p.V()[1], (*this)[2]+p.V()[2], (*this)[3]+p.V()[3] ); } @@ -161,20 +170,20 @@ public: inline void SetRGB( unsigned char r, unsigned char g, unsigned char b ) { - Point4::_v[0] = r; - Point4::_v[1] = g; - Point4::_v[2] = b; - Point4::_v[3] = 0; + Point4::V()[0] = r; + Point4::V()[1] = g; + Point4::V()[2] = b; + Point4::V()[3] = 0; } void SetHSVColor( float h, float s, float v){ float r,g,b; if(s==0.0){ // gray color r = g = b = v; - Point4::_v[0]=(unsigned char)(255*r); - Point4::_v[1]=(unsigned char)(255*g); - Point4::_v[2]=(unsigned char)(255*b); - Point4::_v[3]=255; + Point4::V()[0]=(unsigned char)(255*r); + Point4::V()[1]=(unsigned char)(255*g); + Point4::V()[2]=(unsigned char)(255*b); + Point4::V()[3]=255; return; } if(h==1.0) h = 0.0; @@ -194,11 +203,11 @@ public: case 4: r=t; g=p; b=v; break; case 5: r=v; g=p; b=q; break; } - Point4::_v[0]=(unsigned char)(255*r); - Point4::_v[1]=(unsigned char)(255*g); - Point4::_v[2]=(unsigned char)(255*b); - Point4::_v[3]=255; -// _v[0]=r*256;_v[1]=g*256;_v[2]=b*256; + Point4::V()[0]=(unsigned char)(255*r); + Point4::V()[1]=(unsigned char)(255*g); + Point4::V()[2]=(unsigned char)(255*b); + Point4::V()[3]=255; +// V()[0]=r*256;V()[1]=g*256;V()[2]=b*256; } inline static Color4 GrayShade(float f) @@ -246,10 +255,10 @@ inline void Color4::lerp(const Color4 &c0, const Color4 &c1, const floa assert(x>=0); assert(x<=1); - Point4::_v[0]=(T)(c1._v[0]*x + c0._v[0]*(1.0f-x)); - Point4::_v[1]=(T)(c1._v[1]*x + c0._v[1]*(1.0f-x)); - Point4::_v[2]=(T)(c1._v[2]*x + c0._v[2]*(1.0f-x)); - Point4::_v[3]=(T)(c1._v[3]*x + c0._v[3]*(1.0f-x)); + Point4::V()[0]=(T)(c1.V()[0]*x + c0.V()[0]*(1.0f-x)); + Point4::V()[1]=(T)(c1.V()[1]*x + c0.V()[1]*(1.0f-x)); + Point4::V()[2]=(T)(c1.V()[2]*x + c0.V()[2]*(1.0f-x)); + Point4::V()[3]=(T)(c1.V()[3]*x + c0.V()[3]*(1.0f-x)); } template @@ -257,10 +266,10 @@ inline void Color4::lerp(const Color4 &c0, const Color4 &c1, const Colo { assert(fabs(ip[0]+ip[1]+ip[2]-1)<0.00001); - Point4::_v[0]=(T)(c0[0]*ip[0] + c1[0]*ip[1]+ c2[0]*ip[2]); - Point4::_v[1]=(T)(c0[1]*ip[0] + c1[1]*ip[1]+ c2[1]*ip[2]); - Point4::_v[2]=(T)(c0[2]*ip[0] + c1[2]*ip[1]+ c2[2]*ip[2]); - Point4::_v[3]=(T)(c0[3]*ip[0] + c1[3]*ip[1]+ c2[3]*ip[2]); + Point4::V()[0]=(T)(c0[0]*ip[0] + c1[0]*ip[1]+ c2[0]*ip[2]); + Point4::V()[1]=(T)(c0[1]*ip[0] + c1[1]*ip[1]+ c2[1]*ip[2]); + Point4::V()[2]=(T)(c0[2]*ip[0] + c1[2]*ip[1]+ c2[2]*ip[2]); + Point4::V()[3]=(T)(c0[3]*ip[0] + c1[3]*ip[1]+ c2[3]*ip[2]); } @@ -292,10 +301,10 @@ template <> template <> inline void Color4::Import(const Color4 &b) { - this->_v[0]=b[0]/255.0f; - this->_v[1]=b[1]/255.0f; - this->_v[2]=b[2]/255.0f; - this->_v[3]=b[3]/255.0f; + (*this)[0]=b[0]/255.0f; + (*this)[1]=b[1]/255.0f; + (*this)[2]=b[2]/255.0f; + (*this)[3]=b[3]/255.0f; } #if !defined(__GNUC__) || (__GNUC__ > 3) @@ -304,10 +313,10 @@ template <> // [Bug c++/14479] enum definition in template class with template m template <> inline void Color4::Import(const Color4 &b) { - this->_v[0]=(unsigned char)(b[0]*255.0f); - this->_v[1]=(unsigned char)(b[1]*255.0f); - this->_v[2]=(unsigned char)(b[2]*255.0f); - this->_v[3]=(unsigned char)(b[3]*255.0f); + (*this)[0]=(unsigned char)(b[0]*255.0f); + (*this)[1]=(unsigned char)(b[1]*255.0f); + (*this)[2]=(unsigned char)(b[2]*255.0f); + (*this)[3]=(unsigned char)(b[3]*255.0f); } #if !defined(__GNUC__) || (__GNUC__ > 3) @@ -316,10 +325,10 @@ template <> // [Bug c++/14479] enum definition in template class with template m template <> inline void Color4::Import(const Point4 &b) { - this->_v[0]=(unsigned char)(b[0]*255.0f); - this->_v[1]=(unsigned char)(b[1]*255.0f); - this->_v[2]=(unsigned char)(b[2]*255.0f); - this->_v[3]=(unsigned char)(b[3]*255.0f); + (*this)[0]=(unsigned char)(b[0]*255.0f); + (*this)[1]=(unsigned char)(b[1]*255.0f); + (*this)[2]=(unsigned char)(b[2]*255.0f); + (*this)[3]=(unsigned char)(b[3]*255.0f); } #if !defined(__GNUC__) || (__GNUC__ > 3) @@ -351,10 +360,10 @@ inline Color4 Color4::Construct( const Color4 & b ) //template //inline void Color4::Import(const Color4 &b) //{ -// _v[0] = T(b[0]); -// _v[1] = T(b[1]); -// _v[2] = T(b[2]); -// _v[3] = T(b[3]); +// V()[0] = T(b[0]); +// V()[1] = T(b[1]); +// V()[2] = T(b[2]); +// V()[3] = T(b[3]); //} // template<> @@ -382,10 +391,10 @@ template<> inline Color4 Color4::operator + ( const Color4 & p) const { return Color4( - math::Clamp(int(this->_v[0])+int(p._v[0]),0,255), - math::Clamp(int(this->_v[1])+int(p._v[1]),0,255), - math::Clamp(int(this->_v[2])+int(p._v[2]),0,255), - math::Clamp(int(this->_v[3])+int(p._v[3]),0,255) + math::Clamp(int((*this)[0])+int(p[0]),0,255), + math::Clamp(int((*this)[1])+int(p[1]),0,255), + math::Clamp(int((*this)[2])+int(p[2]),0,255), + math::Clamp(int((*this)[3])+int(p[3]),0,255) ); } diff --git a/vcg/space/deprecated_point3.h b/vcg/space/deprecated_point3.h new file mode 100644 index 00000000..1e05d37c --- /dev/null +++ b/vcg/space/deprecated_point3.h @@ -0,0 +1,572 @@ +/**************************************************************************** +* 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.26 2006/11/13 13:03:45 ponchio +Added GetBBox in Point3 (declaration) the body of the function is in box3.h + +Revision 1.25 2006/10/13 12:59:24 cignoni +Added **explicit** constructor from three coords of a different scalartype + +Revision 1.24 2006/09/28 13:37:35 m_di_benedetto +added non const * V() + +Revision 1.23 2005/11/09 16:11:55 cignoni +Added Abs and LowClampToZero + +Revision 1.22 2005/09/14 14:09:21 m_di_benedetto +Added specialized Convert() for the same scalar type. + +Revision 1.21 2005/05/06 14:45:33 spinelli +cambiato parentesi nel costruttore di GetUV per rendere compatibile tale costruttore con MVC e borland + +Revision 1.20 2005/04/27 16:05:19 callieri +line 466, added parentesis on default value creator getUV [borland] + +Revision 1.19 2004/11/09 15:49:07 ganovelli +added GetUV + +Revision 1.18 2004/10/13 12:45:51 cignoni +Better Doxygen documentation + +Revision 1.17 2004/09/10 14:01:40 cignoni +Added polar to cartesian + +Revision 1.16 2004/03/21 17:14:56 ponchio +Added a math:: + +Revision 1.15 2004/03/07 22:45:32 cignoni +Moved quality and normal functions to the triangle class. + +Revision 1.14 2004/03/05 17:55:01 tarini +errorino: upper case in Zero() + +Revision 1.13 2004/03/03 14:22:48 cignoni +Yet against cr lf mismatch + +Revision 1.12 2004/02/23 23:42:26 cignoni +Translated comments, removed unusued stuff. corrected linefeed/cr + +Revision 1.11 2004/02/19 16:12:28 cignoni +cr lf mismatch 2 + +Revision 1.10 2004/02/19 16:06:24 cignoni +cr lf mismatch + +Revision 1.8 2004/02/19 15:13:40 cignoni +corrected sqrt and added doxygen groups + +Revision 1.7 2004/02/17 02:08:47 cignoni +Di prova... + +Revision 1.6 2004/02/15 23:35:47 cignoni +Cambiato nome type template in accordo alla styleguide + +Revision 1.5 2004/02/10 01:07:15 cignoni +Edited Comments and GPL license + +Revision 1.4 2004/02/09 13:48:02 cignoni +Edited doxygen comments +****************************************************************************/ + +#ifndef __VCGLIB_POINT3 +#define __VCGLIB_POINT3 + +#include +#include + +namespace vcg { + +/** \addtogroup space */ +/*@{*/ + /** + The templated class for representing a point in 3D space. + The class is templated over the ScalarType class that is used to represent coordinates. All the usual + operator overloading (* + - ...) is present. + */ +template class Box3; + +template class Point3 +{ +protected: + /// The only data member. Hidden to user. + P3ScalarType _v[3]; + +public: + typedef P3ScalarType ScalarType; + enum {Dimension = 3}; + + +//@{ + + /** @name Standard Constructors and Initializers + No casting operators have been introduced to avoid automatic unattended (and costly) conversion between different point types + **/ + + inline Point3 () { } + inline Point3 ( const P3ScalarType nx, const P3ScalarType ny, const P3ScalarType nz ) + { + _v[0] = nx; + _v[1] = ny; + _v[2] = nz; + } + inline Point3 ( Point3 const & p ) + { + _v[0]= p._v[0]; + _v[1]= p._v[1]; + _v[2]= p._v[2]; + } + inline Point3 ( const P3ScalarType nv[3] ) + { + _v[0] = nv[0]; + _v[1] = nv[1]; + _v[2] = nv[2]; + } + inline Point3 & operator =( Point3 const & p ) + { + _v[0]= p._v[0]; _v[1]= p._v[1]; _v[2]= p._v[2]; + return *this; + } + inline void SetZero() + { + _v[0] = 0; + _v[1] = 0; + _v[2] = 0; + } + + /// Padding function: give a default 0 value to all the elements that are not in the [0..2] range. + /// Useful for managing in a consistent way object that could have point2 / point3 / point4 + inline P3ScalarType Ext( const int i ) const + { + if(i>=0 && i<=2) return _v[i]; + else return 0; + } + + template + inline void Import( const Point3 & b ) + { + _v[0] = P3ScalarType(b[0]); + _v[1] = P3ScalarType(b[1]); + _v[2] = P3ScalarType(b[2]); + } + + template + static inline Point3 Construct( const Point3 & b ) + { + return Point3(P3ScalarType(b[0]),P3ScalarType(b[1]),P3ScalarType(b[2])); + } + + template + static inline Point3 Construct( const Q & P0, const Q & P1, const Q & P2) + { + return Point3(P3ScalarType(P0),P3ScalarType(P1),P3ScalarType(P2)); + } + + static inline Point3 Construct( const Point3 & b ) + { + return b; + } + +//@} + +//@{ + + /** @name Data Access. + access to data is done by overloading of [] or explicit naming of coords (x,y,z)**/ + + inline P3ScalarType & operator [] ( const int i ) + { + assert(i>=0 && i<3); + return _v[i]; + } + inline const P3ScalarType & operator [] ( const int i ) const + { + assert(i>=0 && i<3); + return _v[i]; + } + inline const P3ScalarType &X() const { return _v[0]; } + inline const P3ScalarType &Y() const { return _v[1]; } + inline const P3ScalarType &Z() const { return _v[2]; } + inline P3ScalarType &X() { return _v[0]; } + inline P3ScalarType &Y() { return _v[1]; } + inline P3ScalarType &Z() { return _v[2]; } + inline const P3ScalarType * V() const + { + return _v; + } + inline P3ScalarType * V() + { + return _v; + } + inline P3ScalarType & V( const int i ) + { + assert(i>=0 && i<3); + return _v[i]; + } + inline const P3ScalarType & V( const int i ) const + { + assert(i>=0 && i<3); + return _v[i]; + } +//@} +//@{ + + /** @name Classical overloading of operators + Note + **/ + + inline Point3 operator + ( Point3 const & p) const + { + return Point3( _v[0]+p._v[0], _v[1]+p._v[1], _v[2]+p._v[2] ); + } + inline Point3 operator - ( Point3 const & p) const + { + return Point3( _v[0]-p._v[0], _v[1]-p._v[1], _v[2]-p._v[2] ); + } + inline Point3 operator * ( const P3ScalarType s ) const + { + return Point3( _v[0]*s, _v[1]*s, _v[2]*s ); + } + inline Point3 operator / ( const P3ScalarType s ) const + { + return Point3( _v[0]/s, _v[1]/s, _v[2]/s ); + } + /// Dot product + inline P3ScalarType operator * ( Point3 const & p ) const + { + return ( _v[0]*p._v[0] + _v[1]*p._v[1] + _v[2]*p._v[2] ); + } + inline P3ScalarType dot( const Point3 & p ) const { return (*this) * p; } + /// Cross product + inline Point3 operator ^ ( Point3 const & p ) const + { + return Point3 + ( + _v[1]*p._v[2] - _v[2]*p._v[1], + _v[2]*p._v[0] - _v[0]*p._v[2], + _v[0]*p._v[1] - _v[1]*p._v[0] + ); + } + + inline Point3 & operator += ( Point3 const & p) + { + _v[0] += p._v[0]; + _v[1] += p._v[1]; + _v[2] += p._v[2]; + return *this; + } + inline Point3 & operator -= ( Point3 const & p) + { + _v[0] -= p._v[0]; + _v[1] -= p._v[1]; + _v[2] -= p._v[2]; + return *this; + } + inline Point3 & operator *= ( const P3ScalarType s ) + { + _v[0] *= s; + _v[1] *= s; + _v[2] *= s; + return *this; + } + inline Point3 & operator /= ( const P3ScalarType s ) + { + _v[0] /= s; + _v[1] /= s; + _v[2] /= s; + return *this; + } + // Norme + inline P3ScalarType Norm() const + { + return math::Sqrt( _v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] ); + } + inline P3ScalarType SquaredNorm() const + { + return ( _v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] ); + } + // Scalatura differenziata + inline Point3 & Scale( const P3ScalarType sx, const P3ScalarType sy, const P3ScalarType sz ) + { + _v[0] *= sx; + _v[1] *= sy; + _v[2] *= sz; + return *this; + } + inline Point3 & Scale( const Point3 & p ) + { + _v[0] *= p._v[0]; + _v[1] *= p._v[1]; + _v[2] *= p._v[2]; + return *this; + } + + // Normalizzazione + inline Point3 & Normalize() + { + P3ScalarType n = math::Sqrt(_v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2]); + if(n>0.0) { _v[0] /= n; _v[1] /= n; _v[2] /= n; } + return *this; + } + + // for compatibility with eigen port + inline Point3 & normalized() { return Normalize(); } + + /** + * Convert to polar coordinates from cartesian coordinates. + * + * Theta is the azimuth angle and ranges between [0, 360) degrees. + * Phi is the elevation angle (not the polar angle) and ranges between [-90, 90] degrees. + * + * /note Note that instead of the classical polar angle, which ranges between + * 0 and 180 degrees we opt for the elevation angle to obtain a more + * intuitive spherical coordinate system. + */ + void ToPolar(P3ScalarType &ro, P3ScalarType &theta, P3ScalarType &phi) const + { + ro = Norm(); + theta = (P3ScalarType)atan2(_v[2], _v[0]); + phi = (P3ScalarType)asin(_v[1]/ro); + } + + /** + * Convert from polar coordinates to cartesian coordinates. + * + * Theta is the azimuth angle and ranges between [0, 360) degrees. + * Phi is the elevation angle (not the polar angle) and ranges between [-90, 90] degrees. + * + * \note Note that instead of the classical polar angle, which ranges between + * 0 and 180 degrees, we opt for the elevation angle to obtain a more + * intuitive spherical coordinate system. + */ + void FromPolar(const P3ScalarType &ro, const P3ScalarType &theta, const P3ScalarType &phi) + { + _v[0]= ro*cos(theta)*cos(phi); + _v[1]= ro*sin(phi); + _v[2]= ro*sin(theta)*cos(phi); + } + + Box3 GetBBox(Box3 &bb) const; +//@} +//@{ + + /** @name Comparison Operators. + Note that the reverse z prioritized ordering, useful in many situations. + **/ + +inline bool operator == ( Point3 const & p ) const + { + return _v[0]==p._v[0] && _v[1]==p._v[1] && _v[2]==p._v[2]; + } + inline bool operator != ( Point3 const & p ) const + { + return _v[0]!=p._v[0] || _v[1]!=p._v[1] || _v[2]!=p._v[2]; + } + inline bool operator < ( Point3 const & p ) const + { + return (_v[2]!=p._v[2])?(_v[2] ( Point3 const & p ) const + { + return (_v[2]!=p._v[2])?(_v[2]>p._v[2]): + (_v[1]!=p._v[1])?(_v[1]>p._v[1]): + (_v[0]>p._v[0]); + } + inline bool operator <= ( Point3 const & p ) const + { + return (_v[2]!=p._v[2])?(_v[2]< p._v[2]): + (_v[1]!=p._v[1])?(_v[1]< p._v[1]): + (_v[0]<=p._v[0]); + } + inline bool operator >= ( Point3 const & p ) const + { + return (_v[2]!=p._v[2])?(_v[2]> p._v[2]): + (_v[1]!=p._v[1])?(_v[1]> p._v[1]): + (_v[0]>=p._v[0]); + } + + + inline Point3 operator - () const + { + return Point3 ( -_v[0], -_v[1], -_v[2] ); + } + //@} + +}; // end class definition + + +template +inline P3ScalarType Angle( Point3 const & p1, Point3 const & p2 ) +{ + P3ScalarType w = p1.Norm()*p2.Norm(); + if(w==0) return -1; + P3ScalarType t = (p1*p2)/w; + if(t>1) t = 1; + else if(t<-1) t = -1; + return (P3ScalarType) acos(t); +} + +// versione uguale alla precedente ma che assume che i due vettori sono unitari +template +inline P3ScalarType AngleN( Point3 const & p1, Point3 const & p2 ) +{ + P3ScalarType w = p1*p2; + if(w>1) + w = 1; + else if(w<-1) + w=-1; + return (P3ScalarType) acos(w); +} + + +template +inline P3ScalarType Norm( Point3 const & p ) +{ + return p.Norm(); +} + +template +inline P3ScalarType SquaredNorm( Point3 const & p ) +{ + return p.SquaredNorm(); +} + +template +inline Point3 & Normalize( Point3 & p ) +{ + p.Normalize(); + return p; +} + +template +inline P3ScalarType Distance( Point3 const & p1,Point3 const & p2 ) +{ + return (p1-p2).Norm(); +} + +template +inline P3ScalarType SquaredDistance( Point3 const & p1,Point3 const & p2 ) +{ + return (p1-p2).SquaredNorm(); +} + + // Dot product preciso numericamente (solo double!!) + // Implementazione: si sommano i prodotti per ordine di esponente + // (prima le piu' grandi) +template +double stable_dot ( Point3 const & p0, Point3 const & p1 ) +{ + P3ScalarType k0 = p0._v[0]*p1._v[0]; + P3ScalarType k1 = p0._v[1]*p1._v[1]; + P3ScalarType k2 = p0._v[2]*p1._v[2]; + + int exp0,exp1,exp2; + + frexp( double(k0), &exp0 ); + frexp( double(k1), &exp1 ); + frexp( double(k2), &exp2 ); + + if( exp0 +P3ScalarType PSDist( const Point3 & p, + const Point3 & v1, + const Point3 & v2, + Point3 & q ) +{ + Point3 e = v2-v1; + P3ScalarType t = ((p-v1)*e)/e.SquaredNorm(); + if(t<0) t = 0; + else if(t>1) t = 1; + q = v1+e*t; + return Distance(p,q); +} + + +template +void GetUV( Point3 &n,Point3 &u, Point3 &v, Point3 up=(Point3(0,1,0)) ) +{ + n.Normalize(); + const double LocEps=double(1e-7); + u=n^up; + double len = u.Norm(); + if(len < LocEps) + { + if(fabs(n[0])(1,0,0); // x is the min + else up=Point3(0,0,1); // z is the min + }else { + if(fabs(n[1])(0,1,0); // y is the min + else up=Point3(0,0,1); // z is the min + } + u=n^up; + } + u.Normalize(); + v=n^u; + v.Normalize(); + Point3 uv=u^v; +} + + +template +inline Point3 Abs(const Point3 & p) { + return (Point3(math::Abs(p[0]), math::Abs(p[1]), math::Abs(p[2]))); +} +// probably a more uniform naming should be defined... +template +inline Point3 LowClampToZero(const Point3 & p) { + return (Point3(math::Max(p[0], (SCALARTYPE)0), math::Max(p[1], (SCALARTYPE)0), math::Max(p[2], (SCALARTYPE)0))); +} + +typedef Point3 Point3s; +typedef Point3 Point3i; +typedef Point3 Point3f; +typedef Point3 Point3d; + +/*@}*/ + +} // end namespace + +#endif + diff --git a/vcg/space/deprecated_point4.h b/vcg/space/deprecated_point4.h new file mode 100644 index 00000000..73375ade --- /dev/null +++ b/vcg/space/deprecated_point4.h @@ -0,0 +1,383 @@ +/**************************************************************************** +* 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.12 2006/06/21 11:06:16 ganovelli +changed return type of Zero() (to void) + +Revision 1.11 2005/04/13 09:40:30 ponchio +Including math/bash.h + +Revision 1.10 2005/03/18 16:34:42 fiorin +minor changes to comply gcc compiler + +Revision 1.9 2005/01/21 18:02:11 ponchio +Removed dependence from matrix44 and changed VectProd + +Revision 1.8 2005/01/12 11:25:02 ganovelli +added Dimension + +Revision 1.7 2004/10/11 17:46:11 ganovelli +added definition of vector product (not implemented) + +Revision 1.6 2004/05/10 11:16:19 ganovelli +include assert.h added + +Revision 1.5 2004/03/31 10:09:58 cignoni +missing return value in zero() + +Revision 1.4 2004/03/11 17:17:49 tarini +added commets (doxy), uniformed with new style, now using math::, ... +added HomoNormalize(), Zero()... remade StableDot() (hand made sort). + +Revision 1.1 2004/02/10 01:11:28 cignoni +Edited Comments and GPL license + +****************************************************************************/ + +#ifndef __VCGLIB_POINT4 +#define __VCGLIB_POINT4 +#include + +#include + +namespace vcg { +/** \addtogroup space */ +/*@{*/ + /** + The templated class for representing a point in 4D space. + The class is templated over the ScalarType class that is used to represent coordinates. + All the usual operator (* + - ...) are defined. + */ + +template class Point4 +{ +public: + /// The only data member. Hidden to user. + T _v[4]; + +public: + typedef T ScalarType; + enum {Dimension = 4}; + +//@{ + + /** @name Standard Constructors and Initializers + No casting operators have been introduced to avoid automatic unattended (and costly) conversion between different point types + **/ + + inline Point4 () { } + inline Point4 ( const T nx, const T ny, const T nz , const T nw ) + { + _v[0] = nx; _v[1] = ny; _v[2] = nz; _v[3] = nw; + } + inline Point4 ( const T p[4] ) + { + _v[0] = p[0]; _v[1]= p[1]; _v[2] = p[2]; _v[3]= p[3]; + } + inline Point4 ( const Point4 & p ) + { + _v[0]= p._v[0]; _v[1]= p._v[1]; _v[2]= p._v[2]; _v[3]= p._v[3]; + } + inline void Zero() + { + _v[0] = _v[1] = _v[2] = _v[3]= 0; + } + template + /// importer from different Point4 types + inline void Import( const Point4 & b ) + { + _v[0] = T(b[0]); _v[1] = T(b[1]); + _v[2] = T(b[2]); + _v[3] = T(b[3]); + } + /// constuctor that imports from different Point4 types + template + static inline Point4 Construct( const Point4 & b ) + { + return Point4(T(b[0]),T(b[1]),T(b[2]),T(b[3])); + } + +//@} + +//@{ + + /** @name Data Access. + access to data is done by overloading of [] or explicit naming of coords (x,y,z,w) + **/ + inline const T & operator [] ( const int i ) const + { + assert(i>=0 && i<4); + return _v[i]; + } + inline T & operator [] ( const int i ) + { + assert(i>=0 && i<4); + return _v[i]; + } + inline T &X() {return _v[0];} + inline T &Y() {return _v[1];} + inline T &Z() {return _v[2];} + inline T &W() {return _v[3];} + inline T const * V() const + { + return _v; + } + inline const T & V ( const int i ) const + { + assert(i>=0 && i<4); + return _v[i]; + } + inline T & V ( const int i ) + { + assert(i>=0 && i<4); + return _v[i]; + } + /// Padding function: give a default 0 value to all the elements that are not in the [0..2] range. + /// Useful for managing in a consistent way object that could have point2 / point3 / point4 + inline T Ext( const int i ) const + { + if(i>=0 && i<=3) return _v[i]; + else return 0; + } +//@} + +//@{ + /** @name Linear operators and the likes + **/ + inline Point4 operator + ( const Point4 & p) const + { + return Point4( _v[0]+p._v[0], _v[1]+p._v[1], _v[2]+p._v[2], _v[3]+p._v[3] ); + } + inline Point4 operator - ( const Point4 & p) const + { + return Point4( _v[0]-p._v[0], _v[1]-p._v[1], _v[2]-p._v[2], _v[3]-p._v[3] ); + } + inline Point4 operator * ( const T s ) const + { + return Point4( _v[0]*s, _v[1]*s, _v[2]*s, _v[3]*s ); + } + inline Point4 operator / ( const T s ) const + { + return Point4( _v[0]/s, _v[1]/s, _v[2]/s, _v[3]/s ); + } + inline Point4 & operator += ( const Point4 & p) + { + _v[0] += p._v[0]; _v[1] += p._v[1]; _v[2] += p._v[2]; _v[3] += p._v[3]; + return *this; + } + inline Point4 & operator -= ( const Point4 & p ) + { + _v[0] -= p._v[0]; _v[1] -= p._v[1]; _v[2] -= p._v[2]; _v[3] -= p._v[3]; + return *this; + } + inline Point4 & operator *= ( const T s ) + { + _v[0] *= s; _v[1] *= s; _v[2] *= s; _v[3] *= s; + return *this; + } + inline Point4 & operator /= ( const T s ) + { + _v[0] /= s; _v[1] /= s; _v[2] /= s; _v[3] /= s; + return *this; + } + inline Point4 operator - () const + { + return Point4( -_v[0], -_v[1], -_v[2], -_v[3] ); + } + inline Point4 VectProd ( const Point4 &x, const Point4 &z ) const + { + Point4 res; + const Point4 &y = *this; + + res[0] = y[1]*x[2]*z[3]-y[1]*x[3]*z[2]-x[1]*y[2]*z[3]+ + x[1]*y[3]*z[2]+z[1]*y[2]*x[3]-z[1]*y[3]*x[2]; + res[1] = y[0]*x[3]*z[2]-z[0]*y[2]*x[3]-y[0]*x[2]* + z[3]+z[0]*y[3]*x[2]+x[0]*y[2]*z[3]-x[0]*y[3]*z[2]; + res[2] = -y[0]*z[1]*x[3]+x[0]*z[1]*y[3]+y[0]*x[1]* + z[3]-x[0]*y[1]*z[3]-z[0]*x[1]*y[3]+z[0]*y[1]*x[3]; + res[3] = -z[0]*y[1]*x[2]-y[0]*x[1]*z[2]+x[0]*y[1]* + z[2]+y[0]*z[1]*x[2]-x[0]*z[1]*y[2]+z[0]*x[1]*y[2]; + return res; + } +//@} + +//@{ + /** @name Norms and normalizations + **/ + /// Euclidian normal + inline T Norm() const + { + return math::Sqrt( _v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] + _v[3]*_v[3] ); + } + /// Squared euclidian normal + inline T SquaredNorm() const + { + return _v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] + _v[3]*_v[3]; + } + /// Euclidian normalization + inline Point4 & Normalize() + { + T n = sqrt(_v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] + _v[3]*_v[3] ); + if(n>0.0) { _v[0] /= n; _v[1] /= n; _v[2] /= n; _v[3] /= n; } + return *this; + } + /// Homogeneous normalization (division by W) + inline Point4 & HomoNormalize(){ + if (_v[3]!=0.0) { _v[0] /= _v[3]; _v[1] /= _v[3]; _v[2] /= _v[3]; _v[3]=1.0; } + return *this; + }; + +//@} + +//@{ + /** @name Comparison operators (lexicographical order) + **/ + inline bool operator == ( const Point4& p ) const + { + return _v[0]==p._v[0] && _v[1]==p._v[1] && _v[2]==p._v[2] && _v[3]==p._v[3]; + } + inline bool operator != ( const Point4 & p ) const + { + return _v[0]!=p._v[0] || _v[1]!=p._v[1] || _v[2]!=p._v[2] || _v[3]!=p._v[3]; + } + inline bool operator < ( Point4 const & p ) const + { + return (_v[3]!=p._v[3])?(_v[3] ( const Point4 & p ) const + { + return (_v[3]!=p._v[3])?(_v[3]>p._v[3]): + (_v[2]!=p._v[2])?(_v[2]>p._v[2]): + (_v[1]!=p._v[1])?(_v[1]>p._v[1]): + (_v[0]>p._v[0]); + } + inline bool operator <= ( const Point4 & p ) const + { + return (_v[3]!=p._v[3])?(_v[3]< p._v[3]): + (_v[2]!=p._v[2])?(_v[2]< p._v[2]): + (_v[1]!=p._v[1])?(_v[1]< p._v[1]): + (_v[0]<=p._v[0]); + } + inline bool operator >= ( const Point4 & p ) const + { + return (_v[3]!=p._v[3])?(_v[3]> p._v[3]): + (_v[2]!=p._v[2])?(_v[2]> p._v[2]): + (_v[1]!=p._v[1])?(_v[1]> p._v[1]): + (_v[0]>=p._v[0]); + } +//@} + +//@{ + /** @name Dot products + **/ + + // dot product + inline T operator * ( const Point4 & p ) const + { + return _v[0]*p._v[0] + _v[1]*p._v[1] + _v[2]*p._v[2] + _v[3]*p._v[3]; + } + inline T dot( const Point4 & p ) const { return (*this) * p; } + inline Point4 operator ^ ( const Point4& p ) const + { + assert(0);// not defined by two vectors (only put for metaprogramming) + return Point4(); + } + + /// slower version, more stable (double precision only) + T StableDot ( const Point4 & p ) const + { + + T k0=_v[0]*p._v[0], k1=_v[1]*p._v[1], k2=_v[2]*p._v[2], k3=_v[3]*p._v[3]; + int exp0,exp1,exp2,exp3; + + frexp( double(k0), &exp0 );frexp( double(k1), &exp1 ); + frexp( double(k2), &exp2 );frexp( double(k3), &exp3 ); + + if (exp0>exp1) { math::Swap(k0,k1); math::Swap(exp0,exp1); } + if (exp2>exp3) { math::Swap(k2,k3); math::Swap(exp2,exp3); } + if (exp0>exp2) { math::Swap(k0,k2); math::Swap(exp0,exp2); } + if (exp1>exp3) { math::Swap(k1,k3); math::Swap(exp1,exp3); } + if (exp2>exp3) { math::Swap(k2,k3); math::Swap(exp2,exp3); } + + return ( (k0 + k1) + k2 ) +k3; + } +//@} + + +}; // end class definition + +template +T Angle( const Point4& p1, const Point4 & p2 ) +{ + T w = p1.Norm()*p2.Norm(); + if(w==0) return -1; + T t = (p1*p2)/w; + if(t>1) t=1; + return T( math::Acos(t) ); +} + +template +inline T Norm( const Point4 & p ) +{ + return p.Norm(); +} + +template +inline T SquaredNorm( const Point4 & p ) +{ + return p.SquaredNorm(); +} + +template +inline T Distance( const Point4 & p1, const Point4 & p2 ) +{ + return Norm(p1-p2); +} + +template +inline T SquaredDistance( const Point4 & p1, const Point4 & p2 ) +{ + return SquaredNorm(p1-p2); +} + + /// slower version of dot product, more stable (double precision only) +template +double StableDot ( Point4 const & p0, Point4 const & p1 ) +{ + return p0.StableDot(p1); +} + +typedef Point4 Point4s; +typedef Point4 Point4i; +typedef Point4 Point4f; +typedef Point4 Point4d; + +/*@}*/ +} // end namespace +#endif diff --git a/vcg/space/intersection3.h b/vcg/space/intersection3.h index 0033ecd6..6e1ef6f4 100644 --- a/vcg/space/intersection3.h +++ b/vcg/space/intersection3.h @@ -339,10 +339,10 @@ namespace vcg { inline bool IntersectionLinePlane( const Plane3 & pl, const Line3 & li, Point3 & po){ const T epsilon = T(1e-8); - T k = pl.Direction() * li.Direction(); // Compute 'k' factor + T k = pl.Direction().dot(li.Direction()); // Compute 'k' factor if( (k > -epsilon) && (k < epsilon)) return false; - T r = (pl.Offset() - pl.Direction()*li.Origin())/k; // Compute ray distance + T r = (pl.Offset() - pl.Direction().dot(li.Origin()))/k; // Compute ray distance po = li.Origin() + li.Direction()*r; return true; } diff --git a/vcg/space/line3.h b/vcg/space/line3.h index e1d8a1c5..b1f0e6f2 100644 --- a/vcg/space/line3.h +++ b/vcg/space/line3.h @@ -121,8 +121,8 @@ public: { return _ori!=p._ori || _dir!=p._dir; } /// Projects a point on the line inline ScalarType Projection( const PointType &p ) const - { if (NORM) return ScalarType((p-_ori)*_dir); - else return ScalarType((p-_ori)*_dir/_dir.SquaredNorm()); + { if (NORM) return ScalarType((p-_ori).dot(_dir)); + else return ScalarType((p-_ori).dot(_dir)/_dir.SquaredNorm()); } /// returns wheter this type is normalized or not static bool IsNormalized() {return NORM;}; diff --git a/vcg/space/plane3.h b/vcg/space/plane3.h index 924c0a26..5620338d 100644 --- a/vcg/space/plane3.h +++ b/vcg/space/plane3.h @@ -141,7 +141,7 @@ public: ///Project a point on the plane PointType Projection(const PointType &p) const { - ScalarType k = p * _dir - _offset; + ScalarType k = p.dot(_dir) - _offset; return p - _dir * k; } @@ -154,13 +154,13 @@ public: void Init(const PointType &p0, const PointType &p1, const PointType &p2) { _dir = (p2 - p0) ^ (p1 - p0); if(NORM) Normalize(); - _offset = p0 * _dir; + _offset = p0.dot(_dir); } /// Calculates the plane passing through a point and the normal (Rename this method inline void Init(const PointType &p0, const PointType &norm) { _dir = norm; - _offset = p0 * _dir; + _offset = p0.dot(_dir); } }; // end class Plane3 @@ -169,12 +169,12 @@ typedef Plane3 Plane3d; ///Distance plane - point (Move this function to somewhere else) template T Distance(const Plane3 &plane, const Point3 &point) { - return plane.Direction() * point - plane.Offset(); + return plane.Direction().dot(point) - plane.Offset(); } ///Distance point-plane (Move this function to somewhere else) template T Distance(const Point3 &point, const Plane3 &plane) { - return plane.Direction() * point - plane.Offset(); + return plane.Direction().dot(point) - plane.Offset(); } } // end namespace diff --git a/vcg/space/point3.h b/vcg/space/point3.h index 189e0abd..1e338daf 100644 --- a/vcg/space/point3.h +++ b/vcg/space/point3.h @@ -20,83 +20,26 @@ * for more details. * * * ****************************************************************************/ -/**************************************************************************** - History -$Log: not supported by cvs2svn $ -Revision 1.26 2006/11/13 13:03:45 ponchio -Added GetBBox in Point3 (declaration) the body of the function is in box3.h - -Revision 1.25 2006/10/13 12:59:24 cignoni -Added **explicit** constructor from three coords of a different scalartype - -Revision 1.24 2006/09/28 13:37:35 m_di_benedetto -added non const * V() - -Revision 1.23 2005/11/09 16:11:55 cignoni -Added Abs and LowClampToZero - -Revision 1.22 2005/09/14 14:09:21 m_di_benedetto -Added specialized Convert() for the same scalar type. - -Revision 1.21 2005/05/06 14:45:33 spinelli -cambiato parentesi nel costruttore di GetUV per rendere compatibile tale costruttore con MVC e borland - -Revision 1.20 2005/04/27 16:05:19 callieri -line 466, added parentesis on default value creator getUV [borland] - -Revision 1.19 2004/11/09 15:49:07 ganovelli -added GetUV - -Revision 1.18 2004/10/13 12:45:51 cignoni -Better Doxygen documentation - -Revision 1.17 2004/09/10 14:01:40 cignoni -Added polar to cartesian - -Revision 1.16 2004/03/21 17:14:56 ponchio -Added a math:: - -Revision 1.15 2004/03/07 22:45:32 cignoni -Moved quality and normal functions to the triangle class. - -Revision 1.14 2004/03/05 17:55:01 tarini -errorino: upper case in Zero() - -Revision 1.13 2004/03/03 14:22:48 cignoni -Yet against cr lf mismatch - -Revision 1.12 2004/02/23 23:42:26 cignoni -Translated comments, removed unusued stuff. corrected linefeed/cr - -Revision 1.11 2004/02/19 16:12:28 cignoni -cr lf mismatch 2 - -Revision 1.10 2004/02/19 16:06:24 cignoni -cr lf mismatch - -Revision 1.8 2004/02/19 15:13:40 cignoni -corrected sqrt and added doxygen groups - -Revision 1.7 2004/02/17 02:08:47 cignoni -Di prova... - -Revision 1.6 2004/02/15 23:35:47 cignoni -Cambiato nome type template in accordo alla styleguide - -Revision 1.5 2004/02/10 01:07:15 cignoni -Edited Comments and GPL license - -Revision 1.4 2004/02/09 13:48:02 cignoni -Edited doxygen comments -****************************************************************************/ +#ifndef VCG_USE_EIGEN +#include "deprecated_point3.h" +#else #ifndef __VCGLIB_POINT3 #define __VCGLIB_POINT3 -#include +#include "../math/eigen.h" #include +namespace vcg{ +template class Point3; +} + +namespace Eigen{ +template +struct ei_traits > : ei_traits > {}; +} + namespace vcg { /** \addtogroup space */ @@ -108,14 +51,22 @@ namespace vcg { */ template class Box3; -template class Point3 +template class Point3 : public Eigen::Matrix<_Scalar,3,1> { -protected: - /// The only data member. Hidden to user. - P3ScalarType _v[3]; + typedef Eigen::Matrix<_Scalar,3,1> _Base; + using _Base::coeff; + using _Base::coeffRef; + using _Base::setZero; + using _Base::data; + using _Base::V; public: - typedef P3ScalarType ScalarType; + + _EIGEN_GENERIC_PUBLIC_INTERFACE(Point3,_Base); + typedef Scalar ScalarType; + + VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Point3) + enum {Dimension = 3}; @@ -125,63 +76,40 @@ public: No casting operators have been introduced to avoid automatic unattended (and costly) conversion between different point types **/ - inline Point3 () { } - inline Point3 ( const P3ScalarType nx, const P3ScalarType ny, const P3ScalarType nz ) - { - _v[0] = nx; - _v[1] = ny; - _v[2] = nz; - } - inline Point3 ( Point3 const & p ) - { - _v[0]= p._v[0]; - _v[1]= p._v[1]; - _v[2]= p._v[2]; - } - inline Point3 ( const P3ScalarType nv[3] ) - { - _v[0] = nv[0]; - _v[1] = nv[1]; - _v[2] = nv[2]; - } - inline Point3 & operator =( Point3 const & p ) - { - _v[0]= p._v[0]; _v[1]= p._v[1]; _v[2]= p._v[2]; - return *this; - } - inline void SetZero() - { - _v[0] = 0; - _v[1] = 0; - _v[2] = 0; - } - + inline Point3 () {} + inline Point3 ( const Scalar nx, const Scalar ny, const Scalar nz ) : Base(nx,ny,nz) {} + inline Point3 ( Point3 const & p ) : Base(p) {} + inline Point3 ( const Scalar nv[3] ) : Base(nv) {} + template + inline Point3(const Eigen::MatrixBase& other) : Base(other) {} + /// Padding function: give a default 0 value to all the elements that are not in the [0..2] range. /// Useful for managing in a consistent way object that could have point2 / point3 / point4 - inline P3ScalarType Ext( const int i ) const + inline Scalar Ext( const int i ) const { - if(i>=0 && i<=2) return _v[i]; + if(i>=0 && i<=2) return data()[i]; else return 0; } - template - inline void Import( const Point3 & b ) + template + inline void Import( const Eigen::MatrixBase& b ) { - _v[0] = P3ScalarType(b[0]); - _v[1] = P3ScalarType(b[1]); - _v[2] = P3ScalarType(b[2]); + EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(OtherDerived,3); + data()[0] = Scalar(b[0]); + data()[1] = Scalar(b[1]); + data()[2] = Scalar(b[2]); } template static inline Point3 Construct( const Point3 & b ) { - return Point3(P3ScalarType(b[0]),P3ScalarType(b[1]),P3ScalarType(b[2])); + return Point3(Scalar(b[0]),Scalar(b[1]),Scalar(b[2])); } template static inline Point3 Construct( const Q & P0, const Q & P1, const Q & P2) { - return Point3(P3ScalarType(P0),P3ScalarType(P1),P3ScalarType(P2)); + return Point3(Scalar(P0),Scalar(P1),Scalar(P2)); } static inline Point3 Construct( const Point3 & b ) @@ -196,39 +124,21 @@ public: /** @name Data Access. access to data is done by overloading of [] or explicit naming of coords (x,y,z)**/ - inline P3ScalarType & operator [] ( const int i ) + inline const Scalar &X() const { return data()[0]; } + inline const Scalar &Y() const { return data()[1]; } + inline const Scalar &Z() const { return data()[2]; } + inline Scalar &X() { return data()[0]; } + inline Scalar &Y() { return data()[1]; } + inline Scalar &Z() { return data()[2]; } + inline Scalar & V( const int i ) { assert(i>=0 && i<3); - return _v[i]; + return data()[i]; } - inline const P3ScalarType & operator [] ( const int i ) const + inline const Scalar & V( const int i ) const { assert(i>=0 && i<3); - return _v[i]; - } - inline const P3ScalarType &X() const { return _v[0]; } - inline const P3ScalarType &Y() const { return _v[1]; } - inline const P3ScalarType &Z() const { return _v[2]; } - inline P3ScalarType &X() { return _v[0]; } - inline P3ScalarType &Y() { return _v[1]; } - inline P3ScalarType &Z() { return _v[2]; } - inline const P3ScalarType * V() const - { - return _v; - } - inline P3ScalarType * V() - { - return _v; - } - inline P3ScalarType & V( const int i ) - { - assert(i>=0 && i<3); - return _v[i]; - } - inline const P3ScalarType & V( const int i ) const - { - assert(i>=0 && i<3); - return _v[i]; + return data()[i]; } //@} //@{ @@ -237,96 +147,19 @@ public: Note **/ - inline Point3 operator + ( Point3 const & p) const + // Scalatura differenziata + inline Point3 & Scale( const Scalar sx, const Scalar sy, const Scalar sz ) { - return Point3( _v[0]+p._v[0], _v[1]+p._v[1], _v[2]+p._v[2] ); - } - inline Point3 operator - ( Point3 const & p) const - { - return Point3( _v[0]-p._v[0], _v[1]-p._v[1], _v[2]-p._v[2] ); - } - inline Point3 operator * ( const P3ScalarType s ) const - { - return Point3( _v[0]*s, _v[1]*s, _v[2]*s ); - } - inline Point3 operator / ( const P3ScalarType s ) const - { - return Point3( _v[0]/s, _v[1]/s, _v[2]/s ); - } - /// Dot product - inline P3ScalarType operator * ( Point3 const & p ) const - { - return ( _v[0]*p._v[0] + _v[1]*p._v[1] + _v[2]*p._v[2] ); - } - /// Cross product - inline Point3 operator ^ ( Point3 const & p ) const - { - return Point3 - ( - _v[1]*p._v[2] - _v[2]*p._v[1], - _v[2]*p._v[0] - _v[0]*p._v[2], - _v[0]*p._v[1] - _v[1]*p._v[0] - ); - } - - inline Point3 & operator += ( Point3 const & p) - { - _v[0] += p._v[0]; - _v[1] += p._v[1]; - _v[2] += p._v[2]; - return *this; - } - inline Point3 & operator -= ( Point3 const & p) - { - _v[0] -= p._v[0]; - _v[1] -= p._v[1]; - _v[2] -= p._v[2]; - return *this; - } - inline Point3 & operator *= ( const P3ScalarType s ) - { - _v[0] *= s; - _v[1] *= s; - _v[2] *= s; - return *this; - } - inline Point3 & operator /= ( const P3ScalarType s ) - { - _v[0] /= s; - _v[1] /= s; - _v[2] /= s; - return *this; - } - // Norme - inline P3ScalarType Norm() const - { - return math::Sqrt( _v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] ); - } - inline P3ScalarType SquaredNorm() const - { - return ( _v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] ); - } - // Scalatura differenziata - inline Point3 & Scale( const P3ScalarType sx, const P3ScalarType sy, const P3ScalarType sz ) - { - _v[0] *= sx; - _v[1] *= sy; - _v[2] *= sz; + data()[0] *= sx; + data()[1] *= sy; + data()[2] *= sz; return *this; } inline Point3 & Scale( const Point3 & p ) { - _v[0] *= p._v[0]; - _v[1] *= p._v[1]; - _v[2] *= p._v[2]; - return *this; - } - - // Normalizzazione - inline Point3 & Normalize() - { - P3ScalarType n = math::Sqrt(_v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2]); - if(n>0.0) { _v[0] /= n; _v[1] /= n; _v[2] /= n; } + data()[0] *= p.data()[0]; + data()[1] *= p.data()[1]; + data()[2] *= p.data()[2]; return *this; } @@ -340,11 +173,11 @@ public: * 0 and 180 degrees we opt for the elevation angle to obtain a more * intuitive spherical coordinate system. */ - void ToPolar(P3ScalarType &ro, P3ScalarType &theta, P3ScalarType &phi) const + void ToPolar(Scalar &ro, Scalar &theta, Scalar &phi) const { - ro = Norm(); - theta = (P3ScalarType)atan2(_v[2], _v[0]); - phi = (P3ScalarType)asin(_v[1]/ro); + ro = this->norm(); + theta = (Scalar)atan2(data()[2], data()[0]); + phi = (Scalar)asin(data()[1]/ro); } /** @@ -357,14 +190,14 @@ public: * 0 and 180 degrees, we opt for the elevation angle to obtain a more * intuitive spherical coordinate system. */ - void FromPolar(const P3ScalarType &ro, const P3ScalarType &theta, const P3ScalarType &phi) + void FromPolar(const Scalar &ro, const Scalar &theta, const Scalar &phi) { - _v[0]= ro*cos(theta)*cos(phi); - _v[1]= ro*sin(phi); - _v[2]= ro*sin(theta)*cos(phi); + data()[0]= ro*cos(theta)*cos(phi); + data()[1]= ro*sin(phi); + data()[2]= ro*sin(theta)*cos(phi); } - Box3 GetBBox(Box3 &bb) const; + Box3<_Scalar> GetBBox(Box3<_Scalar> &bb) const; //@} //@{ @@ -372,113 +205,64 @@ public: Note that the reverse z prioritized ordering, useful in many situations. **/ -inline bool operator == ( Point3 const & p ) const - { - return _v[0]==p._v[0] && _v[1]==p._v[1] && _v[2]==p._v[2]; - } - inline bool operator != ( Point3 const & p ) const - { - return _v[0]!=p._v[0] || _v[1]!=p._v[1] || _v[2]!=p._v[2]; - } inline bool operator < ( Point3 const & p ) const { - return (_v[2]!=p._v[2])?(_v[2] ( Point3 const & p ) const { - return (_v[2]!=p._v[2])?(_v[2]>p._v[2]): - (_v[1]!=p._v[1])?(_v[1]>p._v[1]): - (_v[0]>p._v[0]); + return (data()[2]!=p.data()[2])?(data()[2]>p.data()[2]): + (data()[1]!=p.data()[1])?(data()[1]>p.data()[1]): + (data()[0]>p.data()[0]); } inline bool operator <= ( Point3 const & p ) const { - return (_v[2]!=p._v[2])?(_v[2]< p._v[2]): - (_v[1]!=p._v[1])?(_v[1]< p._v[1]): - (_v[0]<=p._v[0]); + return (data()[2]!=p.data()[2])?(data()[2]< p.data()[2]): + (data()[1]!=p.data()[1])?(data()[1]< p.data()[1]): + (data()[0]<=p.data()[0]); } inline bool operator >= ( Point3 const & p ) const { - return (_v[2]!=p._v[2])?(_v[2]> p._v[2]): - (_v[1]!=p._v[1])?(_v[1]> p._v[1]): - (_v[0]>=p._v[0]); - } - - - inline Point3 operator - () const - { - return Point3 ( -_v[0], -_v[1], -_v[2] ); + return (data()[2]!=p.data()[2])?(data()[2]> p.data()[2]): + (data()[1]!=p.data()[1])?(data()[1]> p.data()[1]): + (data()[0]>=p.data()[0]); } //@} }; // end class definition -template -inline P3ScalarType Angle( Point3 const & p1, Point3 const & p2 ) -{ - P3ScalarType w = p1.Norm()*p2.Norm(); - if(w==0) return -1; - P3ScalarType t = (p1*p2)/w; - if(t>1) t = 1; - else if(t<-1) t = -1; - return (P3ScalarType) acos(t); -} - // versione uguale alla precedente ma che assume che i due vettori sono unitari -template -inline P3ScalarType AngleN( Point3 const & p1, Point3 const & p2 ) +template +inline Scalar AngleN( Point3 const & p1, Point3 const & p2 ) { - P3ScalarType w = p1*p2; + Scalar w = p1*p2; if(w>1) w = 1; else if(w<-1) w=-1; - return (P3ScalarType) acos(w); + return (Scalar) acos(w); } -template -inline P3ScalarType Norm( Point3 const & p ) -{ - return p.Norm(); -} - -template -inline P3ScalarType SquaredNorm( Point3 const & p ) -{ - return p.SquaredNorm(); -} - -template -inline Point3 & Normalize( Point3 & p ) +template +inline Point3 & Normalize( Point3 & p ) { p.Normalize(); return p; } -template -inline P3ScalarType Distance( Point3 const & p1,Point3 const & p2 ) -{ - return (p1-p2).Norm(); -} - -template -inline P3ScalarType SquaredDistance( Point3 const & p1,Point3 const & p2 ) -{ - return (p1-p2).SquaredNorm(); -} - // Dot product preciso numericamente (solo double!!) // Implementazione: si sommano i prodotti per ordine di esponente // (prima le piu' grandi) -template -double stable_dot ( Point3 const & p0, Point3 const & p1 ) +template +double stable_dot ( Point3 const & p0, Point3 const & p1 ) { - P3ScalarType k0 = p0._v[0]*p1._v[0]; - P3ScalarType k1 = p0._v[1]*p1._v[1]; - P3ScalarType k2 = p0._v[2]*p1._v[2]; + Scalar k0 = p0.data()[0]*p1.data()[0]; + Scalar k1 = p0.data()[1]*p1.data()[1]; + Scalar k2 = p0.data()[2]*p1.data()[2]; int exp0,exp1,exp2; @@ -505,14 +289,14 @@ double stable_dot ( Point3 const & p0, Point3 const /// Point(p) Edge(v1-v2) dist, q is the point in v1-v2 with min dist -template -P3ScalarType PSDist( const Point3 & p, - const Point3 & v1, - const Point3 & v2, - Point3 & q ) +template +Scalar PSDist( const Point3 & p, + const Point3 & v1, + const Point3 & v2, + Point3 & q ) { - Point3 e = v2-v1; - P3ScalarType t = ((p-v1)*e)/e.SquaredNorm(); + Point3 e = v2-v1; + Scalar t = ((p-v1).dot(e))/e.SquaredNorm(); if(t<0) t = 0; else if(t>1) t = 1; q = v1+e*t; @@ -520,8 +304,8 @@ P3ScalarType PSDist( const Point3 & p, } -template -void GetUV( Point3 &n,Point3 &u, Point3 &v, Point3 up=(Point3(0,1,0)) ) +template +void GetUV( Point3 &n,Point3 &u, Point3 &v, Point3 up=(Point3(0,1,0)) ) { n.Normalize(); const double LocEps=double(1e-7); @@ -530,18 +314,18 @@ void GetUV( Point3 &n,Point3 &u, Point3(1,0,0); // x is the min - else up=Point3(0,0,1); // z is the min + if(fabs(n[0])(1,0,0); // x is the min + else up=Point3(0,0,1); // z is the min }else { - if(fabs(n[1])(0,1,0); // y is the min - else up=Point3(0,0,1); // z is the min + if(fabs(n[1])(0,1,0); // y is the min + else up=Point3(0,0,1); // z is the min } u=n^up; } u.Normalize(); v=n^u; v.Normalize(); - Point3 uv=u^v; + Point3 uv=u^v; } @@ -566,3 +350,4 @@ typedef Point3 Point3d; #endif +#endif diff --git a/vcg/space/point4.h b/vcg/space/point4.h index f0892e9e..5e8eda55 100644 --- a/vcg/space/point4.h +++ b/vcg/space/point4.h @@ -20,48 +20,24 @@ * for more details. * * * ****************************************************************************/ -/**************************************************************************** - History -$Log: not supported by cvs2svn $ -Revision 1.12 2006/06/21 11:06:16 ganovelli -changed return type of Zero() (to void) - -Revision 1.11 2005/04/13 09:40:30 ponchio -Including math/bash.h - -Revision 1.10 2005/03/18 16:34:42 fiorin -minor changes to comply gcc compiler - -Revision 1.9 2005/01/21 18:02:11 ponchio -Removed dependence from matrix44 and changed VectProd - -Revision 1.8 2005/01/12 11:25:02 ganovelli -added Dimension - -Revision 1.7 2004/10/11 17:46:11 ganovelli -added definition of vector product (not implemented) - -Revision 1.6 2004/05/10 11:16:19 ganovelli -include assert.h added - -Revision 1.5 2004/03/31 10:09:58 cignoni -missing return value in zero() - -Revision 1.4 2004/03/11 17:17:49 tarini -added commets (doxy), uniformed with new style, now using math::, ... -added HomoNormalize(), Zero()... remade StableDot() (hand made sort). - -Revision 1.1 2004/02/10 01:11:28 cignoni -Edited Comments and GPL license - -****************************************************************************/ +#ifndef VCG_USE_EIGEN +#include "deprecated_point4.h" +#else #ifndef __VCGLIB_POINT4 #define __VCGLIB_POINT4 -#include -#include +#include "../math/eigen.h" + +namespace vcg{ +template class Point4; +} + +namespace Eigen{ +template +struct ei_traits > : ei_traits > {}; +} namespace vcg { /** \addtogroup space */ @@ -72,224 +48,119 @@ namespace vcg { All the usual operator (* + - ...) are defined. */ -template class Point4 +template class Point4 : public Eigen::Matrix { -public: - /// The only data member. Hidden to user. - T _v[4]; + typedef Eigen::Matrix _Base; + using _Base::coeff; + using _Base::coeffRef; + using _Base::setZero; + using _Base::data; + using _Base::V; public: - typedef T ScalarType; + + _EIGEN_GENERIC_PUBLIC_INTERFACE(Point4,_Base); + typedef Scalar ScalarType; + + VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Point4) + enum {Dimension = 4}; -//@{ + inline Point4() : Base() {} + inline Point4( const T nx, const T ny, const T nz , const T nw ) : Base(nx,ny,nz,nw) {} + inline Point4(const T p[4]) : Base(p) {} + inline Point4(const Point4& p) : Base(p) {} + template + inline Point4(const Eigen::MatrixBase& other) : Base(other) {} + - /** @name Standard Constructors and Initializers - No casting operators have been introduced to avoid automatic unattended (and costly) conversion between different point types - **/ - - inline Point4 () { } - inline Point4 ( const T nx, const T ny, const T nz , const T nw ) - { - _v[0] = nx; _v[1] = ny; _v[2] = nz; _v[3] = nw; - } - inline Point4 ( const T p[4] ) - { - _v[0] = p[0]; _v[1]= p[1]; _v[2] = p[2]; _v[3]= p[3]; - } - inline Point4 ( const Point4 & p ) - { - _v[0]= p._v[0]; _v[1]= p._v[1]; _v[2]= p._v[2]; _v[3]= p._v[3]; - } - inline void Zero() - { - _v[0] = _v[1] = _v[2] = _v[3]= 0; - } - template /// importer from different Point4 types - inline void Import( const Point4 & b ) - { - _v[0] = T(b[0]); _v[1] = T(b[1]); - _v[2] = T(b[2]); - _v[3] = T(b[3]); - } + template inline void Import( const Point4 & b ) { *this = b.template cast(); } + /// constuctor that imports from different Point4 types - template - static inline Point4 Construct( const Point4 & b ) - { - return Point4(T(b[0]),T(b[1]),T(b[2]),T(b[3])); - } + template + static inline Point4 Construct( const Point4 & b ) { return b.template cast(); } -//@} //@{ - /** @name Data Access. - access to data is done by overloading of [] or explicit naming of coords (x,y,z,w) - **/ - inline const T & operator [] ( const int i ) const - { - assert(i>=0 && i<4); - return _v[i]; - } - inline T & operator [] ( const int i ) - { - assert(i>=0 && i<4); - return _v[i]; - } - inline T &X() {return _v[0];} - inline T &Y() {return _v[1];} - inline T &Z() {return _v[2];} - inline T &W() {return _v[3];} - inline T const * V() const - { - return _v; - } + inline T &X() {return Base::x();} + inline T &Y() {return Base::y();} + inline T &Z() {return Base::z();} + inline T &W() {return Base::w();} + inline const T & V ( const int i ) const { assert(i>=0 && i<4); - return _v[i]; + return data()[i]; } inline T & V ( const int i ) { assert(i>=0 && i<4); - return _v[i]; + return data()[i]; } - /// Padding function: give a default 0 value to all the elements that are not in the [0..2] range. - /// Useful for managing in a consistent way object that could have point2 / point3 / point4 + + /// Padding function: give a default 0 value to all the elements that are not in the [0..2] range. + /// Useful for managing in a consistent way object that could have point2 / point3 / point4 inline T Ext( const int i ) const { - if(i>=0 && i<=3) return _v[i]; + if(i>=0 && i<=3) return data()[i]; else return 0; } //@} -//@{ - /** @name Linear operators and the likes - **/ - inline Point4 operator + ( const Point4 & p) const - { - return Point4( _v[0]+p._v[0], _v[1]+p._v[1], _v[2]+p._v[2], _v[3]+p._v[3] ); - } - inline Point4 operator - ( const Point4 & p) const - { - return Point4( _v[0]-p._v[0], _v[1]-p._v[1], _v[2]-p._v[2], _v[3]-p._v[3] ); - } - inline Point4 operator * ( const T s ) const - { - return Point4( _v[0]*s, _v[1]*s, _v[2]*s, _v[3]*s ); - } - inline Point4 operator / ( const T s ) const - { - return Point4( _v[0]/s, _v[1]/s, _v[2]/s, _v[3]/s ); - } - inline Point4 & operator += ( const Point4 & p) - { - _v[0] += p._v[0]; _v[1] += p._v[1]; _v[2] += p._v[2]; _v[3] += p._v[3]; - return *this; - } - inline Point4 & operator -= ( const Point4 & p ) - { - _v[0] -= p._v[0]; _v[1] -= p._v[1]; _v[2] -= p._v[2]; _v[3] -= p._v[3]; - return *this; - } - inline Point4 & operator *= ( const T s ) - { - _v[0] *= s; _v[1] *= s; _v[2] *= s; _v[3] *= s; - return *this; - } - inline Point4 & operator /= ( const T s ) - { - _v[0] /= s; _v[1] /= s; _v[2] /= s; _v[3] /= s; - return *this; - } - inline Point4 operator - () const - { - return Point4( -_v[0], -_v[1], -_v[2], -_v[3] ); - } inline Point4 VectProd ( const Point4 &x, const Point4 &z ) const { Point4 res; const Point4 &y = *this; - res[0] = y[1]*x[2]*z[3]-y[1]*x[3]*z[2]-x[1]*y[2]*z[3]+ - x[1]*y[3]*z[2]+z[1]*y[2]*x[3]-z[1]*y[3]*x[2]; - res[1] = y[0]*x[3]*z[2]-z[0]*y[2]*x[3]-y[0]*x[2]* - z[3]+z[0]*y[3]*x[2]+x[0]*y[2]*z[3]-x[0]*y[3]*z[2]; + res[0] = y[1]*x[2]*z[3]-y[1]*x[3]*z[2]-x[1]*y[2]*z[3]+ + x[1]*y[3]*z[2]+z[1]*y[2]*x[3]-z[1]*y[3]*x[2]; + res[1] = y[0]*x[3]*z[2]-z[0]*y[2]*x[3]-y[0]*x[2]* + z[3]+z[0]*y[3]*x[2]+x[0]*y[2]*z[3]-x[0]*y[3]*z[2]; res[2] = -y[0]*z[1]*x[3]+x[0]*z[1]*y[3]+y[0]*x[1]* - z[3]-x[0]*y[1]*z[3]-z[0]*x[1]*y[3]+z[0]*y[1]*x[3]; + z[3]-x[0]*y[1]*z[3]-z[0]*x[1]*y[3]+z[0]*y[1]*x[3]; res[3] = -z[0]*y[1]*x[2]-y[0]*x[1]*z[2]+x[0]*y[1]* - z[2]+y[0]*z[1]*x[2]-x[0]*z[1]*y[2]+z[0]*x[1]*y[2]; + z[2]+y[0]*z[1]*x[2]-x[0]*z[1]*y[2]+z[0]*x[1]*y[2]; return res; } -//@} -//@{ - /** @name Norms and normalizations - **/ - /// Euclidian normal - inline T Norm() const - { - return math::Sqrt( _v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] + _v[3]*_v[3] ); - } - /// Squared euclidian normal - inline T SquaredNorm() const - { - return _v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] + _v[3]*_v[3]; - } - /// Euclidian normalization - inline Point4 & Normalize() - { - T n = sqrt(_v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] + _v[3]*_v[3] ); - if(n>0.0) { _v[0] /= n; _v[1] /= n; _v[2] /= n; _v[3] /= n; } - return *this; - } /// Homogeneous normalization (division by W) - inline Point4 & HomoNormalize(){ - if (_v[3]!=0.0) { _v[0] /= _v[3]; _v[1] /= _v[3]; _v[2] /= _v[3]; _v[3]=1.0; } + inline Point4 & HomoNormalize() { + if (data()[3]!=0.0) { Base::template start<3>() /= coeff(3); coeffRef(3) = 1.0; } return *this; - }; - -//@} + } //@{ /** @name Comparison operators (lexicographical order) **/ - inline bool operator == ( const Point4& p ) const - { - return _v[0]==p._v[0] && _v[1]==p._v[1] && _v[2]==p._v[2] && _v[3]==p._v[3]; - } - inline bool operator != ( const Point4 & p ) const - { - return _v[0]!=p._v[0] || _v[1]!=p._v[1] || _v[2]!=p._v[2] || _v[3]!=p._v[3]; - } inline bool operator < ( Point4 const & p ) const { - return (_v[3]!=p._v[3])?(_v[3] ( const Point4 & p ) const { - return (_v[3]!=p._v[3])?(_v[3]>p._v[3]): - (_v[2]!=p._v[2])?(_v[2]>p._v[2]): - (_v[1]!=p._v[1])?(_v[1]>p._v[1]): - (_v[0]>p._v[0]); + return (data()[3]!=p.data()[3])?(data()[3]>p.data()[3]): + (data()[2]!=p.data()[2])?(data()[2]>p.data()[2]): + (data()[1]!=p.data()[1])?(data()[1]>p.data()[1]): + (data()[0]>p.data()[0]); } inline bool operator <= ( const Point4 & p ) const { - return (_v[3]!=p._v[3])?(_v[3]< p._v[3]): - (_v[2]!=p._v[2])?(_v[2]< p._v[2]): - (_v[1]!=p._v[1])?(_v[1]< p._v[1]): - (_v[0]<=p._v[0]); + return (data()[3]!=p.data()[3])?(data()[3]< p.data()[3]): + (data()[2]!=p.data()[2])?(data()[2]< p.data()[2]): + (data()[1]!=p.data()[1])?(data()[1]< p.data()[1]): + (data()[0]<=p.data()[0]); } inline bool operator >= ( const Point4 & p ) const { - return (_v[3]!=p._v[3])?(_v[3]> p._v[3]): - (_v[2]!=p._v[2])?(_v[2]> p._v[2]): - (_v[1]!=p._v[1])?(_v[1]> p._v[1]): - (_v[0]>=p._v[0]); + return (data()[3]!=p.data()[3])?(data()[3]> p.data()[3]): + (data()[2]!=p.data()[2])?(data()[2]> p.data()[2]): + (data()[1]!=p.data()[1])?(data()[1]> p.data()[1]): + (data()[0]>=p.data()[0]); } //@} @@ -297,22 +168,16 @@ public: /** @name Dot products **/ - // dot product - inline T operator * ( const Point4 & p ) const - { - return _v[0]*p._v[0] + _v[1]*p._v[1] + _v[2]*p._v[2] + _v[3]*p._v[3]; - } inline Point4 operator ^ ( const Point4& p ) const { - assert(0);// not defined by two vectors (only put for metaprogramming) + assert(0 && "not defined by two vectors (only put for metaprogramming)"); return Point4(); } /// slower version, more stable (double precision only) T StableDot ( const Point4 & p ) const { - - T k0=_v[0]*p._v[0], k1=_v[1]*p._v[1], k2=_v[2]*p._v[2], k3=_v[3]*p._v[3]; + T k0=data()[0]*p.data()[0], k1=data()[1]*p.data()[1], k2=data()[2]*p.data()[2], k3=data()[3]*p.data()[3]; int exp0,exp1,exp2,exp3; frexp( double(k0), &exp0 );frexp( double(k1), &exp1 ); @@ -331,39 +196,6 @@ public: }; // end class definition -template -T Angle( const Point4& p1, const Point4 & p2 ) -{ - T w = p1.Norm()*p2.Norm(); - if(w==0) return -1; - T t = (p1*p2)/w; - if(t>1) t=1; - return T( math::Acos(t) ); -} - -template -inline T Norm( const Point4 & p ) -{ - return p.Norm(); -} - -template -inline T SquaredNorm( const Point4 & p ) -{ - return p.SquaredNorm(); -} - -template -inline T Distance( const Point4 & p1, const Point4 & p2 ) -{ - return Norm(p1-p2); -} - -template -inline T SquaredDistance( const Point4 & p1, const Point4 & p2 ) -{ - return SquaredNorm(p1-p2); -} /// slower version of dot product, more stable (double precision only) template @@ -380,3 +212,5 @@ typedef Point4 Point4d; /*@}*/ } // end namespace #endif + +#endif diff --git a/vcg/space/ray3.h b/vcg/space/ray3.h index 829e40b1..2a6fb8ae 100644 --- a/vcg/space/ray3.h +++ b/vcg/space/ray3.h @@ -112,8 +112,8 @@ public: { return _ori!=p._ori || _dir!=p._dir; } /// Projects a point on the ray inline ScalarType Projection( const PointType &p ) const - { if (NORM) return ScalarType((p-_ori)*_dir); - else return ScalarType((p-_ori)*_dir/_dir.SquaredNorm()); + { if (NORM) return ScalarType((p-_ori).dot(_dir)); + else return ScalarType((p-_ori).dot(_dir)/_dir.SquaredNorm()); } /// returns wheter this type is normalized or not static bool IsNormalized() {return NORM;}; diff --git a/wrap/gl/math.h b/wrap/gl/math.h index c9f343f5..b90998ab 100644 --- a/wrap/gl/math.h +++ b/wrap/gl/math.h @@ -72,40 +72,32 @@ namespace vcg { inline void glMultMatrixE(const Matrix44f &matrix) { //glMultMatrixf((const GLfloat *)(matrix[0])); - if(glMultTransposeMatrixf) glMultTransposeMatrixf((const GLfloat *)(matrix[0])); + if(glMultTransposeMatrixf) glMultTransposeMatrixf((const GLfloat *)(matrix.V())); else { - Matrix44f tmp(matrix); - Transpose(tmp); - glMultMatrixf((const GLfloat *)(tmp[0])); + glMultMatrixf((const GLfloat *)(matrix.transpose().eval().V())); } } inline void glMultMatrix(const Matrix44f &matrix) { - Matrix44f tmp(matrix); - Transpose(tmp); - glMultMatrixf((const GLfloat *)(tmp[0])); + glMultMatrixf((const GLfloat *)(matrix.transpose().eval().V())); } inline void glMultMatrixE(const Matrix44d &matrix) { - if(glMultTransposeMatrixd) glMultTransposeMatrixd((const GLdouble *)(matrix[0])); + if(glMultTransposeMatrixd) glMultTransposeMatrixd((const GLdouble *)(matrix.V())); else { - Matrix44d tmp(matrix); - Transpose(tmp); - glMultMatrixd((const GLdouble *)(tmp[0])); + glMultMatrixd((const GLdouble *)(matrix.transpose().eval().V())); } } inline void glMultMatrix(const Matrix44d &matrix) { - Matrix44d tmp(matrix); - Transpose(tmp); - glMultMatrixd((const GLdouble *)(tmp[0])); + glMultMatrixd((const GLdouble *)(matrix.transpose().eval().V())); } inline void glMultMatrixDirect(const Matrix44f &matrix) { - glMultMatrixf((const GLfloat *)(matrix[0])); + glMultMatrixf((const GLfloat *)(matrix.V())); } inline void glMultMatrixDirect(const Matrix44d &matrix) { - glMultMatrixd((const GLdouble *)(matrix[0])); + glMultMatrixd((const GLdouble *)(matrix.V())); } inline void glMultMatrix(const Similarityf &s) { @@ -128,21 +120,23 @@ inline void glMultMatrix(const Similarityd &s) { } inline void glGetv(const GLenum pname, Matrix44f & m){ - glGetFloatv(pname,&m[0][0]); - Transpose(m); + Matrix44f tmp; + glGetFloatv(pname,tmp.V()); + m = tmp.transpose(); } inline void glGetv(const GLenum pname, Matrix44d & m){ - glGetDoublev(pname,&m[0][0]); - Transpose(m); + Matrix44d tmp; + glGetDoublev(pname,tmp.V()); + m = tmp.transpose(); } inline void glGetDirectv(const GLenum pname, Matrix44f & m){ - glGetFloatv(pname,&m[0][0]); + glGetFloatv(pname,m.V()); } inline void glGetDirecv(const GLenum pname, Matrix44d & m){ - glGetDoublev(pname,&m[0][0]); + glGetDoublev(pname,m.V()); } diff --git a/wrap/gl/space.h b/wrap/gl/space.h index 0ec4ff2f..dea392f2 100644 --- a/wrap/gl/space.h +++ b/wrap/gl/space.h @@ -243,9 +243,8 @@ inline void glBoxClip(const Box3 & b) *(Point3*)&m[2][0] = *(Point3*)&v[0];m[2][3]=0; *(Point3*)&m[3][0] = *(Point3*)&c1[0];m[3][3]=1; - vcg::Transpose(m); glPushMatrix(); - glMultMatrix(m); + glMultMatrix(m.transpose()); glBegin(GL_QUADS); glNormal(Point3(0,1,0)); @@ -279,5 +278,18 @@ template glTriangle3(Triangle3(c.P(1),c.P(0),c.P(3))); } +#ifdef VCG_USE_EIGEN + +#define _WRAP_EIGEN_XPR(FUNC) template \ + inline void FUNC(const Eigen::MatrixBase& p) { FUNC(p.eval()); } + +_WRAP_EIGEN_XPR(glVertex) +_WRAP_EIGEN_XPR(glNormal) +_WRAP_EIGEN_XPR(glTexCoord) +_WRAP_EIGEN_XPR(glTranslate) +_WRAP_EIGEN_XPR(glScale) + +#endif + }//namespace #endif diff --git a/wrap/gui/trackmode.cpp b/wrap/gui/trackmode.cpp index c5c7da63..91a8eb6c 100644 --- a/wrap/gui/trackmode.cpp +++ b/wrap/gui/trackmode.cpp @@ -486,8 +486,8 @@ int PathMode::Verse(Point3f reference_point,Point3f current_point,Point3f prev_p prev_dir.Normalize(); next_dir.Normalize(); float prev_coeff,next_coeff; - prev_coeff = prev_dir * reference_dir; - next_coeff = next_dir * reference_dir; + prev_coeff = prev_dir.dot(reference_dir); + next_coeff = next_dir.dot(reference_dir); if (prev_coeff < 0.0f) prev_coeff = 0.0f; if (next_coeff < 0.0f) @@ -573,8 +573,8 @@ void AreaMode::Init(const std::vector < Point3f > &pts) bool pts_not_in_line=false; Point3f a,b; for(unsigned int i=0;i EPSILON; if(pts_not_in_line){ plane.Init( pts[i%npts], diff --git a/wrap/gui/trackutils.h b/wrap/gui/trackutils.h index 79dd0b1e..c7462c25 100644 --- a/wrap/gui/trackutils.h +++ b/wrap/gui/trackutils.h @@ -103,8 +103,8 @@ Plane3f GetViewPlane (const View < float >&camera, const Point3f & center) Point3f vp = camera.ViewPoint (); Plane3f pl; Point3f plnorm = vp - center; - plnorm.Normalize (); - pl.Set (plnorm, plnorm * center); + plnorm.Normalize(); + pl.Set(plnorm, plnorm.dot(center)); return pl; } @@ -331,16 +331,16 @@ Point3f HitSphere (Trackball * tb, const Point3f & p) std::pair< float, bool > LineLineDistance(const Line3f & P,const Line3f & Q,Point3f & P_s, Point3f & Q_t){ Point3f p0 = P.Origin (), Vp = P.Direction (); Point3f q0 = Q.Origin (), Vq = Q.Direction (); - float VPVP = Vp * Vp; - float VQVQ = Vq * Vq; - float VPVQ = Vp * Vq; + float VPVP = Vp.dot(Vp); + float VQVQ = Vq.dot(Vq); + float VPVQ = Vp.dot(Vq); const float det = ( VPVP * VQVQ ) - ( VPVQ * VPVQ ); const float EPSILON = 0.00001f; if ( fabs(det) < EPSILON ) { return std::make_pair(Distance(P,q0), true); } - float b1= (q0 - p0) * Vp; - float b2= (p0 - q0) * Vq; + float b1= (q0 - p0).dot(Vp); + float b2= (p0 - q0).dot(Vq); float s = ( (VQVQ * b1) + (VPVQ * b2) ) / det; float t = ( (VPVQ * b1) + (VPVP * b2) ) / det; P_s = p0 + (Vp * s); @@ -367,16 +367,16 @@ std::pair< float, bool > LineLineDistance(const Line3f & P,const Line3f & Q,Poin std::pair< float, bool > RayLineDistance(const Ray3f & R,const Line3f & Q,Point3f & R_s, Point3f & Q_t){ Point3f r0 = R.Origin (), Vr = R.Direction (); Point3f q0 = Q.Origin (), Vq = Q.Direction (); - float VRVR = Vr * Vr; - float VQVQ = Vq * Vq; - float VRVQ = Vr * Vq; + float VRVR = Vr.dot(Vr); + float VQVQ = Vq.dot(Vq); + float VRVQ = Vr.dot(Vq); const float det = ( VRVR * VQVQ ) - ( VRVQ * VRVQ ); const float EPSILON = 0.00001f; if ( ( det >= 0.0f ? det : -det) < EPSILON ) { return std::make_pair(Distance(Q,r0), true); } - float b1= (q0 - r0) * Vr; - float b2= (r0 - q0) * Vq; + float b1= (q0 - r0).dot(Vr); + float b2= (r0 - q0).dot(Vq); float s = ( (VQVQ * b1) + (VRVQ * b2) ) / det; float t = ( (VRVQ * b1) + (VRVR * b2) ) / det; if(s<0){ @@ -418,17 +418,17 @@ std::pair< float, bool > SegmentSegmentDistance(const Segment3f & R, const Segme R_s=ClosestPoint(R,Q_t); return std::make_pair(Distance(R_s,Q_t),true); } - Point3f r0 = R.P0(), Vr = (R.P1()-R.P0()).Normalize(); - Point3f q0 = Q.P0(), Vq = (Q.P1()-Q.P0()).Normalize(); - float VRVR = Vr * Vr; - float VQVQ = Vq * Vq; - float VRVQ = Vr * Vq; + Point3f r0 = R.P0(), Vr = (R.P1()-R.P0()).normalized(); + Point3f q0 = Q.P0(), Vq = (Q.P1()-Q.P0()).normalized(); + float VRVR = Vr.dot(Vr); + float VQVQ = Vq.dot(Vq); + float VRVQ = Vr.dot(Vq); const float det = ( VRVR * VQVQ ) - ( VRVQ * VRVQ ); const float EPSILON = 0.00001f; if ( ( det >= 0.0f ? det : -det) < EPSILON ) { Line3f lR(R.P0(),R.P1()); float qa=lR.Projection(Q.P0()); - float qb=lR.Projection(Q.P1()); + float qb=lR.Projection(Q.P1()); if( (qa<=0.0f) && qb<=(0.0f)){ R_s=R.P0(); Q_t=ClosestPoint(Q,R_s); @@ -453,8 +453,8 @@ std::pair< float, bool > SegmentSegmentDistance(const Segment3f & R, const Segme } return std::make_pair(Distance(R_s,Q_t),true); } - float b1= (q0 - r0) * Vr; - float b2= (r0 - q0) * Vq; + float b1= (q0 - r0).dot(Vr); + float b2= (r0 - q0).dot(Vq); float s = ( (VQVQ * b1) + (VRVQ * b2) ) / det; float t = ( (VRVQ * b1) + (VRVR * b2) ) / det; if( s < 0 ){ @@ -529,7 +529,7 @@ Line3f ProjectLineOnPlane(const Line3f & ln, const Plane3f & pl) */ float signedDistance(Line3f line,Point3f pt,Point3f positive_dir) { - return Distance(line,pt) * ((((pt-ClosestPoint(line,pt)) * positive_dir) >= 0.0f )? 1.0f: -1.0f); + return Distance(line,pt) * ((((pt-ClosestPoint(line,pt)).dot(positive_dir)) >= 0.0f )? 1.0f: -1.0f); } /*! @@ -557,10 +557,10 @@ template inline bool IntersectionRayPlane( const Plane3 & pl, const Ray3 & ray, Point3 &po){ const T epsilon = T(1e-8); - T k = pl.Direction() * ray.Direction(); // Compute 'k' factor + T k = pl.Direction().dot(ray.Direction()); // Compute 'k' factor if( (k > -epsilon) && (k < epsilon)) return false; - T r = (pl.Offset() - pl.Direction()*ray.Origin())/k; // Compute ray distance + T r = (pl.Offset() - pl.Direction().dot(ray.Origin()))/k; // Compute ray distance if (r < 0) return false; po = ray.Origin() + ray.Direction()*r; @@ -887,8 +887,8 @@ void DrawUglyPlaneMode(Trackball * tb,Plane3f plane) if(norm == d1 || norm == -d1) d1 = Point3f(1,0,0); d2=plane.Projection(d1); - d1=(d2 - p0).Normalize(); - d2=(d1 ^ norm).Normalize(); + d1=(d2 - p0).normalized(); + d2=(d1 ^ norm).normalized(); glLineWidth(3.0); glColor3f(0.2f, 0.2f, 0.9f); glBegin(GL_LINES); @@ -945,8 +945,8 @@ void DrawUglyCylinderMode(Trackball * tb,Line3f axis) if(norm == d1 || norm == -d1) d1 = Point3f(1,0,0); d2=plane.Projection(d1); - d1=(d2 - p0).Normalize(); - d2=(d1 ^ norm).Normalize(); + d1=(d2 - p0).normalized(); + d2=(d1 ^ norm).normalized(); glLineWidth(1.0); glColor3f(0.2f, 0.2f, 0.9f); for(int i=-100;i<100;i++){ @@ -1098,8 +1098,8 @@ void DrawUglyAreaMode(Trackball * tb,const std::vector < Point3f > &points, if(norm == d1 || norm == -d1) d1 = Point3f(1,0,0); d2=plane.Projection(d1); - d1=(d2 - p0).Normalize(); - d2=(d1 ^ norm).Normalize(); + d1=(d2 - p0).normalized(); + d2=(d1 ^ norm).normalized(); glLineWidth(3.0); glColor3f(0.2f, 0.2f, 0.9f); glBegin(GL_LINES); diff --git a/wrap/gui/view.h b/wrap/gui/view.h index 2ed115da..75794387 100644 --- a/wrap/gui/view.h +++ b/wrap/gui/view.h @@ -151,10 +151,10 @@ template void View::GetView() { glGetDoublev(GL_PROJECTION_MATRIX, m); proj.Import(Matrix44d(m)); - Transpose(proj); + proj = proj.transpose().eval(); glGetDoublev(GL_MODELVIEW_MATRIX, m); model.Import(Matrix44d(m)); - Transpose(model); + model = model.transpose().eval(); glGetIntegerv(GL_VIEWPORT, (GLint*)viewport); @@ -205,14 +205,14 @@ template Line3 View::ViewLineFromWindow(const Point3 &p) template Point3 View::Project(const Point3 &p) const { Point3 r; - r = matrix * p; + r = matrix * p; return NormDevCoordToWindowCoord(r); } template Point3 View::UnProject(const Point3 &p) const { Point3 s = WindowCoordToNormDevCoord(p); s = inverse * s ; - return s; + return s; } // Come spiegato nelle glspec diff --git a/wrap/io_trimesh/import_ptx.h b/wrap/io_trimesh/import_ptx.h index 7d1d89a3..e955ba36 100644 --- a/wrap/io_trimesh/import_ptx.h +++ b/wrap/io_trimesh/import_ptx.h @@ -230,7 +230,7 @@ namespace io { else return false; // PTX transformation matrix is transposed - Transpose(currtrasf); + currtrasf = currtrasf.transpose().eval(); // allocating vertex space int vn = rownum*colnum;