diff --git a/vcg/math/eigen.h b/vcg/math/eigen.h index 17359802..4e15336a 100644 --- a/vcg/math/eigen.h +++ b/vcg/math/eigen.h @@ -26,18 +26,19 @@ // TODO enable the vectorization #define EIGEN_DONT_VECTORIZE -#define EIGEN_MATRIXBASE_PLUGIN +#define EIGEN_MATRIXBASE_PLUGIN +#define EIGEN_MATRIX_PLUGIN // forward declarations namespace Eigen { template struct ei_lexi_comparison; } +#include "base.h" #include "../Eigen/LU" #include "../Eigen/Geometry" #include "../Eigen/Array" #include "../Eigen/Core" -#include "base.h" // add support for unsigned char and short int namespace Eigen { @@ -239,6 +240,18 @@ inline typename Eigen::ei_traits::Scalar SquaredDistance(const Eigen::MatrixBase& p1, const Eigen::MatrixBase & p2) { return (p1-p2).norm2(); } +template +inline const Eigen::CwiseUnaryOp::Scalar>, Derived> +Abs(const Eigen::MatrixBase& p) +{ return p.cwise().abs(); } + +template +inline const Eigen::CwiseBinaryOp::Scalar>, + Derived, + Eigen::NestByValue > +LowClampToZero(const Eigen::MatrixBase& p) +{ return p.cwise().max(Derived::Zero().nestByValue()); } + } #endif diff --git a/vcg/math/eigen_matrix_addons.h b/vcg/math/eigen_matrix_addons.h new file mode 100644 index 00000000..09da01b0 --- /dev/null +++ b/vcg/math/eigen_matrix_addons.h @@ -0,0 +1,143 @@ +/**************************************************************************** +* 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. * +* * +****************************************************************************/ + +#warning You are including deprecated math stuff + +enum {Dimension = SizeAtCompileTime}; +typedef vcg::VoidType ParamType; +typedef Matrix PointType; +using Base::V; + +/// importer for points with different scalar type and-or dimensionality +// FIXME the Point3/Point4 specialization were only for same sizes ?? +// while the Point version was generic like this one +template +inline void Import(const MatrixBase& b) +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Matrix); + EIGEN_STATIC_ASSERT_FIXED_SIZE(Matrix); + EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived); + EIGEN_STATIC_ASSERT_FIXED_SIZE(OtherDerived); + enum { OtherSize = OtherDerived::SizeAtCompileTime }; + assert(SizeAtCompileTime<=4 && OtherSize<=4); + data()[0] = Scalar(b[0]); + if (SizeAtCompileTime>1) { if (OtherSize>1) data()[1] = Scalar(b[1]); else data()[1] = 0; } + if (SizeAtCompileTime>2) { if (OtherSize>2) data()[2] = Scalar(b[2]); else data()[2] = 0; } + if (SizeAtCompileTime>3) { if (OtherSize>3) data()[3] = Scalar(b[3]); else data()[3] = 0; } +} + +/// importer for homogeneous points +template +inline void ImportHomo(const MatrixBase& b) +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Matrix); + EIGEN_STATIC_ASSERT_FIXED_SIZE(Matrix); + EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(OtherDerived,SizeAtCompileTime-1); + + this->template start = b; + data()[SizeAtCompileTime-1] = Scalar(1.0); +} + +/// constructor for points with different scalar type and-or dimensionality +template +static inline Matrix Construct(const MatrixBase& b) +{ + Matrix p; p.Import(b); + return p; +} + + /// constructor for homogeneus point. +template +static inline Matrix ConstructHomo(const MatrixBase& b) +{ + Matrix p; p.ImportHomo(b); + return p; +} + +inline const Scalar &X() const { return data()[0]; } +inline const Scalar &Y() const { return data()[1]; } +inline const Scalar &Z() const { assert(SizeAtCompileTime>2); return data()[2]; } +inline Scalar &X() { return data()[0]; } +inline Scalar &Y() { return data()[1]; } +inline Scalar &Z() { assert(SizeAtCompileTime>2); return data()[2]; } + +// note, W always returns the last entry +inline Scalar& W() { return data()[SizeAtCompileTime-1]; } +inline const Scalar& W() const { return data()[SizeAtCompileTime-1]; } + +Scalar* V() { return data(); } +const Scalar* V() const { return data(); } + +// overloaded to return a const reference +inline const Scalar& V( const int i ) const +{ + assert(i>=0 && inorm(); + theta = (Scalar)atan2(data()[2], data()[0]); + phi = (Scalar)asin(data()[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 Scalar &ro, const Scalar &theta, const Scalar &phi) +{ + EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Matrix,3); + data()[0]= ro*ei_cos(theta)*ei_cos(phi); + data()[1]= ro*ei_sin(phi); + data()[2]= ro*ei_sin(theta)*ei_cos(phi); +} diff --git a/vcg/math/eigen_vcgaddons.h b/vcg/math/eigen_matrixbase_addons.h similarity index 98% rename from vcg/math/eigen_vcgaddons.h rename to vcg/math/eigen_matrixbase_addons.h index aba3c936..b2e4dbe5 100644 --- a/vcg/math/eigen_vcgaddons.h +++ b/vcg/math/eigen_matrixbase_addons.h @@ -22,15 +22,15 @@ ****************************************************************************/ #warning You are including deprecated math stuff -/*! -* \deprecated use cols() -*/ + + +typedef Scalar ScalarType; + +/*! \deprecated use cols() */ EIGEN_DEPRECATED inline unsigned int ColumnsNumber() const { return cols(); }; -/*! -* \deprecated use rows() -*/ +/*! \deprecated use rows() */ EIGEN_DEPRECATED inline unsigned int RowsNumber() const { return rows(); }; /*! @@ -83,9 +83,6 @@ 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 * (Modifier) Add to each element of this matrix the scalar constant k. diff --git a/vcg/math/matrix33.h b/vcg/math/matrix33.h index cbbf5b05..726274bc 100644 --- a/vcg/math/matrix33.h +++ b/vcg/math/matrix33.h @@ -79,19 +79,14 @@ public: template Matrix33(const Eigen::MatrixBase& other) : Base(other) {} - /*! - * \deprecated use *this.row(i) - */ + /*! \deprecated use *this.row(i) */ inline typename Base::RowXpr operator[](const unsigned int i) { return Base::row(i); } - /*! - * \deprecated use *this.row(i) - */ + /*! \deprecated use *this.row(i) */ inline const typename Base::RowXpr operator[](const unsigned int i) const { return Base::row(i); } - /** \deprecated */ Matrix33 & SetRotateRad(Scalar angle, const Point3 & axis ) { @@ -129,7 +124,7 @@ public: */ Scalar Norm() { return Base::cwise().abs2().sum(); } // { -// // FIXME looks like there is a bug: j is not used !!! +// // FIXME looks like there was a bug: j is not used !!! // Scalar SQsum=0; // for(int i=0;i<3;++i) // for(int j=0;j<3;++j) @@ -138,8 +133,7 @@ public: // } - /** - It computes the covariance matrix of a set of 3d points. Returns the barycenter + /** Computes the covariance matrix of a set of 3d points. Returns the barycenter. */ // FIXME should be outside Matrix template diff --git a/vcg/math/matrix44.h b/vcg/math/matrix44.h index 3a43117e..89c64b50 100644 --- a/vcg/math/matrix44.h +++ b/vcg/math/matrix44.h @@ -45,7 +45,7 @@ struct ei_traits > : ei_traits (or the typedef) you want a real 4x4 matrix, or use Eigen::Transform if you want a transformation matrix for a 3D space (a Eigen::Transform is internally a 4x4 col-major matrix) + * + * This class represents a 4x4 matrix. T is the kind of element in the matrix. + */ template class Matrix44 : public Eigen::Matrix<_Scalar,4,4,Eigen::RowMajor> // FIXME col or row major ! { @@ -91,6 +93,7 @@ class Matrix44 : public Eigen::Matrix<_Scalar,4,4,Eigen::RowMajor> // FIXME col using _Base::coeffRef; using _Base::ElementAt; using _Base::setZero; + using _Base::operator*; public: @@ -103,7 +106,6 @@ public: ~Matrix44() {} Matrix44(const Matrix44 &m) : Base(m) {} Matrix44(const Scalar * v ) : Base(Eigen::Map >(v)) {} - template Matrix44(const Eigen::MatrixBase& other) : Base(other) {} @@ -116,7 +118,7 @@ public: 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 + template void ToMatrix(Matrix44Type & m) const { m = (*this).template cast(); } void ToEulerAngles(Scalar &alpha, Scalar &beta, Scalar &gamma); @@ -125,30 +127,51 @@ public: void FromMatrix(const Matrix44Type & m) { for(int i = 0; i < 16; i++) Base::data()[i] = m.data()[i]; } void FromEulerAngles(Scalar alpha, Scalar beta, Scalar gamma); - void SetDiagonal(const Scalar k); + 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 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); + Matrix44 &SetShearXY(const Scalar sz); + Matrix44 &SetShearXZ(const Scalar sy); + Matrix44 &SetShearYZ(const Scalar sx); - ///use radiants for angle. - Matrix44 &SetRotateDeg(Scalar AngleDeg, const Point3 & axis); - Matrix44 &SetRotateRad(Scalar AngleRad, const Point3 & axis); + ///use radiants for angle. + Matrix44 &SetRotateDeg(Scalar AngleDeg, const Point3 & axis); + Matrix44 &SetRotateRad(Scalar AngleRad, const Point3 & axis); - template void Import(const Matrix44 &m) { - for(int i = 0; i < 16; i++) - Base::data()[i] = (Scalar)(m.data()[i]); - } - template - static inline Matrix44 Construct( const Matrix44 & b ) - { - Matrix44 tmp; tmp.FromMatrix(b); - return tmp; - } + template void Import(const Matrix44 &m) { + for(int i = 0; i < 16; i++) + Base::data()[i] = (Scalar)(m.data()[i]); + } + template + static inline Matrix44 Construct( const Matrix44 & b ) + { + Matrix44 tmp; tmp.FromMatrix(b); + return tmp; + } +// template Point3 operator*(const Point3 &p) { +// T w; +// Point3 s; +// s[0] = m.ElementAt(0, 0)*p[0] + m.ElementAt(0, 1)*p[1] + m.ElementAt(0, 2)*p[2] + m.ElementAt(0, 3); +// s[1] = m.ElementAt(1, 0)*p[0] + m.ElementAt(1, 1)*p[1] + m.ElementAt(1, 2)*p[2] + m.ElementAt(1, 3); +// s[2] = m.ElementAt(2, 0)*p[0] + m.ElementAt(2, 1)*p[1] + m.ElementAt(2, 2)*p[2] + m.ElementAt(2, 3); +// w = m.ElementAt(3, 0)*p[0] + m.ElementAt(3, 1)*p[1] + m.ElementAt(3, 2)*p[2] + m.ElementAt(3, 3); +// if(w!= 0) s /= w; +// return s; +// } + + Eigen::Matrix operator * (const Eigen::Matrix& p) const { + Scalar w; + Eigen::Matrix s; + s[0] = ElementAt(0, 0)*p[0] + ElementAt(0, 1)*p[1] + ElementAt(0, 2)*p[2] + ElementAt(0, 3); + s[1] = ElementAt(1, 0)*p[0] + ElementAt(1, 1)*p[1] + ElementAt(1, 2)*p[2] + ElementAt(1, 3); + s[2] = ElementAt(2, 0)*p[0] + ElementAt(2, 1)*p[1] + ElementAt(2, 2)*p[2] + ElementAt(2, 3); + w = ElementAt(3, 0)*p[0] + ElementAt(3, 1)*p[1] + ElementAt(3, 2)*p[2] + ElementAt(3, 3); + if(w!= 0) s /= w; + return s; + } }; @@ -156,23 +179,23 @@ public: /** Class for solving A * x = b. */ template class LinearSolve: public Matrix44 { public: - LinearSolve(const Matrix44 &m); - Point4 Solve(const Point4 &b); // solve A � x = b - ///If you need to solve some equation you can use this function instead of Matrix44 one for speed. - T Determinant() const; + LinearSolve(const Matrix44 &m); + Point4 Solve(const Point4 &b); // solve A � x = b + ///If you need to solve some equation you can use this function instead of Matrix44 one for speed. + T Determinant() const; protected: - ///Holds row permutation. - int index[4]; //hold permutation - ///Hold sign of row permutation (used for determinant sign) - T d; - bool Decompose(); + ///Holds row permutation. + int index[4]; //hold permutation + ///Hold sign of row permutation (used for determinant sign) + T d; + bool Decompose(); }; /*** Postmultiply */ //template Point3 operator*(const Point3 &p, const Matrix44 &m); ///Premultiply -template Point3 operator*(const Matrix44 &m, const Point3 &p); +// template Point3 operator*(const Matrix44 &m, const Point3 &p); //return NULL matrix if not invertible template Matrix44 &Invert(Matrix44 &m); @@ -196,9 +219,9 @@ typedef Matrix44 Matrix44d; template < class PointType , class T > void operator*=( std::vector &vert, const Matrix44 & m ) { - typename std::vector::iterator ii; - for(ii=vert.begin();ii!=vert.end();++ii) - (*ii).P()=m * (*ii).P(); + typename std::vector::iterator ii; + for(ii=vert.begin();ii!=vert.end();++ii) + (*ii).P()=m * (*ii).P(); } template @@ -237,36 +260,36 @@ void Matrix44::FromEulerAngles(Scalar alpha, Scalar beta, Scalar gamma) } template void Matrix44::SetDiagonal(const Scalar k) { - setZero(); - ElementAt(0, 0) = k; - ElementAt(1, 1) = k; - ElementAt(2, 2) = k; - ElementAt(3, 3) = 1; + setZero(); + ElementAt(0, 0) = k; + ElementAt(1, 1) = k; + ElementAt(2, 2) = k; + ElementAt(3, 3) = 1; } template Matrix44 &Matrix44::SetScale(const Point3 &t) { - SetScale(t[0], t[1], t[2]); - return *this; + SetScale(t[0], t[1], t[2]); + return *this; } 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; - ElementAt(3, 3) = 1; - return *this; + setZero(); + ElementAt(0, 0) = sx; + ElementAt(1, 1) = sy; + ElementAt(2, 2) = sz; + ElementAt(3, 3) = 1; + return *this; } template Matrix44 &Matrix44::SetTranslate(const Point3 &t) { - SetTranslate(t[0], t[1], t[2]); - return *this; + SetTranslate(t[0], t[1], t[2]); + return *this; } template Matrix44 &Matrix44::SetTranslate(const Scalar tx, const Scalar ty, const Scalar tz) { - Base::setIdentity(); + Base::setIdentity(); ElementAt(0, 3) = tx; - ElementAt(1, 3) = ty; - ElementAt(2, 3) = tz; - return *this; + ElementAt(1, 3) = ty; + ElementAt(2, 3) = tz; + return *this; } template Matrix44 &Matrix44::SetRotateDeg(Scalar AngleDeg, const Point3 & axis) { @@ -274,9 +297,9 @@ template Matrix44 &Matrix44::SetRotateDeg(Scalar AngleDeg, const } 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); + //angle = angle*(T)3.14159265358979323846/180; e' in radianti! + T c = math::Cos(AngleRad); + T s = math::Sin(AngleRad); T q = 1-c; Point3 t = axis; t.Normalize(); @@ -296,94 +319,94 @@ template Matrix44 &Matrix44::SetRotateRad(Scalar AngleRad, const ElementAt(3,1) = 0; ElementAt(3,2) = 0; ElementAt(3,3) = 1; - return *this; + return *this; } - /* Shear Matrixes - XY - 1 k 0 0 x x+ky - 0 1 0 0 y y - 0 0 1 0 z z - 0 0 0 1 1 1 +/* Shear Matrixes +XY +1 k 0 0 x x+ky +0 1 0 0 y y +0 0 1 0 z z +0 0 0 1 1 1 - 1 0 k 0 x x+kz - 0 1 0 0 y y - 0 0 1 0 z z - 0 0 0 1 1 1 +1 0 k 0 x x+kz +0 1 0 0 y y +0 0 1 0 z z +0 0 0 1 1 1 - 1 1 0 0 x x - 0 1 k 0 y y+kz - 0 0 1 0 z z - 0 0 0 1 1 1 +1 1 0 0 x x +0 1 k 0 y y+kz +0 0 1 0 z z +0 0 0 1 1 1 - */ +*/ template Matrix44 & Matrix44::SetShearXY( const Scalar sh) {// shear the X coordinate as the Y coordinate change Base::setIdentity(); ElementAt(0,1) = sh; - return *this; + return *this; } template Matrix44 & Matrix44::SetShearXZ( const Scalar sh) {// shear the X coordinate as the Z coordinate change Base::setIdentity(); ElementAt(0,2) = sh; - return *this; + return *this; } template Matrix44 &Matrix44::SetShearYZ( const Scalar sh) {// shear the Y coordinate as the Z coordinate change Base::setIdentity(); ElementAt(1,2) = sh; - return *this; + return *this; } /* Given a non singular, non projective matrix (e.g. with the last row equal to [0,0,0,1] ) This procedure decompose it in a sequence of - Scale,Shear,Rotation e Translation + Scale,Shear,Rotation e Translation - ScaleV and Tranv are obiviously scaling and translation. - ShearV contains three scalars with, respectively - ShearXY, ShearXZ e ShearYZ + ShearXY, ShearXZ e ShearYZ - RotateV contains the rotations (in degree!) around the x,y,z axis - The input matrix is modified leaving inside it a simple roto translation. + The input matrix is modified leaving inside it a simple roto translation. - To obtain the original matrix the above transformation have to be applied in the strict following way: + To obtain the original matrix the above transformation have to be applied in the strict following way: - OriginalMatrix = Trn * Rtx*Rty*Rtz * ShearYZ*ShearXZ*ShearXY * Scl + OriginalMatrix = Trn * Rtx*Rty*Rtz * ShearYZ*ShearXZ*ShearXY * Scl Example Code: double srv() { return (double(rand()%40)-20)/2.0; } // small random value - srand(time(0)); - Point3d ScV(10+srv(),10+srv(),10+srv()),ScVOut(-1,-1,-1); - Point3d ShV(srv(),srv(),srv()),ShVOut(-1,-1,-1); - Point3d RtV(10+srv(),srv(),srv()),RtVOut(-1,-1,-1); - Point3d TrV(srv(),srv(),srv()),TrVOut(-1,-1,-1); + srand(time(0)); + Point3d ScV(10+srv(),10+srv(),10+srv()),ScVOut(-1,-1,-1); + Point3d ShV(srv(),srv(),srv()),ShVOut(-1,-1,-1); + Point3d RtV(10+srv(),srv(),srv()),RtVOut(-1,-1,-1); + Point3d TrV(srv(),srv(),srv()),TrVOut(-1,-1,-1); - Matrix44d Scl; Scl.SetScale(ScV); - Matrix44d Sxy; Sxy.SetShearXY(ShV[0]); + Matrix44d Scl; Scl.SetScale(ScV); + Matrix44d Sxy; Sxy.SetShearXY(ShV[0]); Matrix44d Sxz; Sxz.SetShearXZ(ShV[1]); Matrix44d Syz; Syz.SetShearYZ(ShV[2]); - Matrix44d Rtx; Rtx.SetRotate(math::ToRad(RtV[0]),Point3d(1,0,0)); + Matrix44d Rtx; Rtx.SetRotate(math::ToRad(RtV[0]),Point3d(1,0,0)); Matrix44d Rty; Rty.SetRotate(math::ToRad(RtV[1]),Point3d(0,1,0)); Matrix44d Rtz; Rtz.SetRotate(math::ToRad(RtV[2]),Point3d(0,0,1)); Matrix44d Trn; Trn.SetTranslate(TrV); Matrix44d StartM = Trn * Rtx*Rty*Rtz * Syz*Sxz*Sxy *Scl; - Matrix44d ResultM=StartM; - Decompose(ResultM,ScVOut,ShVOut,RtVOut,TrVOut); + Matrix44d ResultM=StartM; + Decompose(ResultM,ScVOut,ShVOut,RtVOut,TrVOut); - Scl.SetScale(ScVOut); - Sxy.SetShearXY(ShVOut[0]); - Sxz.SetShearXZ(ShVOut[1]); - Syz.SetShearYZ(ShVOut[2]); - Rtx.SetRotate(math::ToRad(RtVOut[0]),Point3d(1,0,0)); - Rty.SetRotate(math::ToRad(RtVOut[1]),Point3d(0,1,0)); - Rtz.SetRotate(math::ToRad(RtVOut[2]),Point3d(0,0,1)); - Trn.SetTranslate(TrVOut); + Scl.SetScale(ScVOut); + Sxy.SetShearXY(ShVOut[0]); + Sxz.SetShearXZ(ShVOut[1]); + Syz.SetShearYZ(ShVOut[2]); + Rtx.SetRotate(math::ToRad(RtVOut[0]),Point3d(1,0,0)); + Rty.SetRotate(math::ToRad(RtVOut[1]),Point3d(0,1,0)); + Rtz.SetRotate(math::ToRad(RtVOut[2]),Point3d(0,0,1)); + Trn.SetTranslate(TrVOut); - // Now Rebuild is equal to StartM + // Now Rebuild is equal to StartM Matrix44d RebuildM = Trn * Rtx*Rty*Rtz * Syz*Sxz*Sxy * Scl ; */ template @@ -393,7 +416,7 @@ bool Decompose(Matrix44 &M, Point3 &ScaleV, Point3 &ShearV, Point3 & return false; if(math::Abs(M.Determinant())<1e-10) return false; // matrix should be at least invertible... - // First Step recover the traslation + // First Step recover the traslation TranV=M.GetColumn3(3); // Second Step Recover Scale and Shearing interleaved @@ -404,7 +427,7 @@ bool Decompose(Matrix44 &M, Point3 &ScaleV, Point3 &ShearV, Point3 & ShearV[0]=R[0].dot(M.GetColumn3(1)); // xy shearing R[1]= M.GetColumn3(1)-R[0]*ShearV[0]; - assert(math::Abs(R[1].dot(R[0]))<1e-10); + 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]; @@ -425,16 +448,16 @@ bool Decompose(Matrix44 &M, Point3 &ScaleV, Point3 &ShearV, Point3 & ShearV[2]=R[1].dot(M.GetColumn3(2)); // yz shearing ShearV[2]=ShearV[2]/ScaleV[2]; - int i,j; + int i,j; for(i=0;i<3;++i) for(j=0;j<3;++j) M(i,j)=R[j][i]; // Third and last step: Recover the rotation //now the matrix should be a pure rotation matrix so its determinant is +-1 - double det=M.Determinant(); - if(math::Abs(det)<1e-10) return false; // matrix should be at least invertible... - assert(math::Abs(math::Abs(det)-1.0)<1e-10); // it should be +-1... + double det=M.Determinant(); + if(math::Abs(det)<1e-10) return false; // matrix should be at least invertible... + assert(math::Abs(math::Abs(det)-1.0)<1e-10); // it should be +-1... if(det<0) { ScaleV *= -1; M *= -1; @@ -443,37 +466,28 @@ bool Decompose(Matrix44 &M, Point3 &ScaleV, Point3 &ShearV, Point3 & double alpha,beta,gamma; // rotations around the x,y and z axis beta=asin( M(0,2)); double cosbeta=cos(beta); - if(math::Abs(cosbeta) > 1e-5) + if(math::Abs(cosbeta) > 1e-5) { alpha=asin(-M(1,2)/cosbeta); if((M(2,2)/cosbeta) < 0 ) alpha=M_PI-alpha; gamma=asin(-M(0,1)/cosbeta); if((M(0,0)/cosbeta)<0) gamma = M_PI-gamma; } - else + else { alpha=asin(-M(1,0)); if(M(1,1)<0) alpha=M_PI-alpha; gamma=0; } - RotV[0]=math::ToDeg(alpha); + RotV[0]=math::ToDeg(alpha); RotV[1]=math::ToDeg(beta); RotV[2]=math::ToDeg(gamma); return true; } -template Point3 operator*(const Matrix44 &m, const Point3 &p) { - T w; - Point3 s; - s[0] = m.ElementAt(0, 0)*p[0] + m.ElementAt(0, 1)*p[1] + m.ElementAt(0, 2)*p[2] + m.ElementAt(0, 3); - s[1] = m.ElementAt(1, 0)*p[0] + m.ElementAt(1, 1)*p[1] + m.ElementAt(1, 2)*p[2] + m.ElementAt(1, 3); - s[2] = m.ElementAt(2, 0)*p[0] + m.ElementAt(2, 1)*p[1] + m.ElementAt(2, 2)*p[2] + m.ElementAt(2, 3); - w = m.ElementAt(3, 0)*p[0] + m.ElementAt(3, 1)*p[1] + m.ElementAt(3, 2)*p[2] + m.ElementAt(3, 3); - if(w!= 0) s /= w; - return s; -} + @@ -490,14 +504,14 @@ template Point3 operator*(const Matrix44 &m, const Point3 &p) /* - To invert a matrix you can - either invert the matrix inplace calling +To invert a matrix you can +either invert the matrix inplace calling - vcg::Invert(yourMatrix); +vcg::Invert(yourMatrix); - or get the inverse matrix of a given matrix without touching it: +or get the inverse matrix of a given matrix without touching it: - invertedMatrix = vcg::Inverse(untouchedMatrix); +invertedMatrix = vcg::Inverse(untouchedMatrix); */ template Matrix44 & Invert(Matrix44 &m) { diff --git a/vcg/space/point.h b/vcg/space/point.h index 13f0912c..b62d2045 100644 --- a/vcg/space/point.h +++ b/vcg/space/point.h @@ -33,18 +33,18 @@ #include namespace vcg{ -template class Point4; -} +template class Point2; +template class Point3; +template class Point4; -namespace Eigen{ -template -struct ei_traits > : ei_traits > {}; +namespace ndim{ +template class Point; +} } namespace vcg { namespace ndim{ - /** \addtogroup space */ /*@{*/ /** @@ -55,7 +55,17 @@ namespace ndim{ */ template class Point : public Eigen::Matrix { - typedef Eigen::Matrix _Base; +//---------------------------------------- +// template typedef part +// use it as follow: typename Point::Type instead of simply Point +//---------------------------------------- +public: + typedef Eigen::Matrix Type; +//---------------------------------------- +// inheritence part +//---------------------------------------- +private: + typedef Eigen::Matrix _Base; using _Base::coeff; using _Base::coeffRef; using _Base::setZero; @@ -65,212 +75,50 @@ template class Point : public Eigen::Matrix public: _EIGEN_GENERIC_PUBLIC_INTERFACE(Point,_Base); - - typedef S ScalarType; - typedef VoidType ParamType; - typedef Point PointType; - enum {Dimension = N}; - VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Point) -//@{ - - /** @name Standard Constructors and Initializers - No casting operators have been introduced to avoid automatic unattended (and costly) conversion between different PointType types - **/ inline Point() : Base() {} template - inline Point(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 S Ext( const int i ) const - { - if(i>=0 && i<=N) return data()[i]; - else return 0; - } - - /// importer for points with different scalar type and-or dimensionality - template - inline void Import( const Point & b ) - { - data()[0] = ScalarType(b[0]); - data()[1] = ScalarType(b[1]); - if (N>2) { if (N2>2) data()[2] = ScalarType(b[2]); else data()[2] = 0;}; - if (N>3) { if (N2>3) data()[3] = ScalarType(b[3]); else data()[3] = 0;}; - } - - /// constructor for points with different scalar type and-or dimensionality - template - static inline PointType Construct( const Point & b ) - { - PointType p; p.Import(b); - return p; - } - - /// importer for homogeneous points - template - inline void ImportHomo( const Point & b ) - { - data()[0] = ScalarType(b[0]); - data()[1] = ScalarType(b[1]); - if (N>2) { data()[2] = ScalarType(data()[2]); }; - data()[N-1] = 1.0; - } - - /// constructor for homogeneus point. - template - static inline PointType Construct( const Point & b ) - { - PointType p; p.ImportHomo(b); - return p; - } - -//@} - -//@{ - - /** @name Data Access. - access to data is done by overloading of [] or explicit naming of coords (x,y,z)**/ - - inline const S &X() const { return data()[0]; } - inline const S &Y() const { return data()[1]; } - inline const S &Z() const { static_assert(N>2); return data()[2]; } - /// W is in any case the last coordinate. - /// (in a 2D point, W() == Y(). In a 3D point, W()==Z() - /// in a 4D point, W() is a separate component) - inline const S &W() const { return data()[N-1]; } - inline S &X() { return data()[0]; } - inline S &Y() { return data()[1]; } - inline S &Z() { static_assert(N>2); return data()[2]; } - inline S &W() { return data()[N-1]; } -//@} - - -//@{ - - /** @name Dot products (cross product "%" is defined olny for 3D points) - **/ + inline Point(const Eigen::MatrixBase& other) : Base(other) {} /// slower version, more stable (double precision only) - inline S StableDot ( const PointType & p ) const; - -//@} - -//@{ - - /** @name Norms - **/ - - /// Euclidean norm, static version - template static S Norm(const PT &p ); - /// Squared Euclidean norm, static version - template static S SquaredNorm(const PT &p ); - /// Normalization (division by norm), static version - template static PointType & Normalize(const PT &p); - -//@} + inline S StableDot (const Point& p) const; /// Signed area operator - /// a % b returns the signed area of the parallelogram inside a and b - // inline S operator % ( PointType const & p ) const; - - /// Convert to polar coordinates - void ToPolar( S & ro, S & tetha, S & fi ) const - { - ro = Norm(); - tetha = (S)atan2( data()[1], data()[0] ); - fi = (S)acos( data()[2]/ro ); - } - -//@{ - - /** @name Comparison Operators. - Lexicographic order. - **/ - - inline bool operator == ( PointType const & p ) const; - inline bool operator != ( PointType const & p ) const; - inline bool operator < ( PointType const & p ) const; - inline bool operator > ( PointType const & p ) const; - inline bool operator <= ( PointType const & p ) const; - inline bool operator >= ( PointType const & p ) const; - //@} - -//@{ - - /** @name - Glocal to Local and viceversa - (provided for uniformity with other spatial classes. trivial for points) - **/ - - inline PointType LocalToGlobal(ParamType p) const { return *this; } - inline ParamType GlobalToLocal(PointType p) const { ParamType p(); return p; } -//@} + /// a % b returns the signed area of the parallelogram inside a and b + // inline S operator % ( PointType const & p ) const; }; // end class definition -// workaround the lack of template typedef (the next c++ standard will support them :) ) +typedef Eigen::Matrix Point2s; +typedef Eigen::Matrix Point2i; +typedef Eigen::Matrix Point2f; +typedef Eigen::Matrix Point2d; +typedef Eigen::Matrix Vector2s; +typedef Eigen::Matrix Vector2i; +typedef Eigen::Matrix Vector2f; +typedef Eigen::Matrix Vector2d; -template -struct Point2:public Point<2,S>{ - typedef Point<3,S> Base; - inline Point2() : Base() {}; - inline Point2(const Point2& p):Base(p){}; - inline Point2(S a, S b):Base(a,b){}; - template - inline Point2(const Eigen::MatrixBase& other) : Base(other) {} -}; - -template -struct Point3:public Point<3,S> { - typedef Point<3,S> Base; - inline Point3() : Base() {}; - inline Point3(const Point3& p):Base(p){} - inline Point3(S a, S b, S c):Base(a,b,c){}; - template - inline Point3(const Eigen::MatrixBase& other) : Base(other) {} -}; +typedef Eigen::Matrix Point3s; +typedef Eigen::Matrix Point3i; +typedef Eigen::Matrix Point3f; +typedef Eigen::Matrix Point3d; +typedef Eigen::Matrix Vector3s; +typedef Eigen::Matrix Vector3i; +typedef Eigen::Matrix Vector3f; +typedef Eigen::Matrix Vector3d; -template -struct Point4:public Point<4,S>{ - typedef Point<4,S> Base; - inline Point4() : Base() {}; - inline Point4(const Point4& p):Base(p) {} - inline Point4(S a, S b, S c, S d):Base(a,b,c,d){}; - template - inline Point4(const Eigen::MatrixBase& other) : Base(other) {} -}; +typedef Eigen::Matrix Point4s; +typedef Eigen::Matrix Point4i; +typedef Eigen::Matrix Point4f; +typedef Eigen::Matrix Point4d; +typedef Eigen::Matrix Vector4s; +typedef Eigen::Matrix Vector4i; +typedef Eigen::Matrix Vector4f; +typedef Eigen::Matrix Vector4d; -typedef Point<2,short> Point2s; -typedef Point<2,int> Point2i; -typedef Point<2,float> Point2f; -typedef Point<2,double> Point2d; -typedef Point<2,short> Vector2s; -typedef Point<2,int> Vector2i; -typedef Point<2,float> Vector2f; -typedef Point<2,double> Vector2d; - -typedef Point<3,short> Point3s; -typedef Point<3,int> Point3i; -typedef Point<3,float> Point3f; -typedef Point<3,double> Point3d; -typedef Point<3,short> Vector3s; -typedef Point<3,int> Vector3i; -typedef Point<3,float> Vector3f; -typedef Point<3,double> Vector3d; - - -typedef Point<4,short> Point4s; -typedef Point<4,int> Point4i; -typedef Point<4,float> Point4f; -typedef Point<4,double> Point4d; -typedef Point<4,short> Vector4s; -typedef Point<4,int> Vector4i; -typedef Point<4,float> Vector4f; -typedef Point<4,double> Vector4d; /*@}*/ diff --git a/vcg/space/point2.h b/vcg/space/point2.h index 69be34db..749e751f 100644 --- a/vcg/space/point2.h +++ b/vcg/space/point2.h @@ -29,28 +29,37 @@ #define __VCGLIB_POINT2 #include "../math/eigen.h" -#include +// #include "point.h" namespace vcg{ -template class Point2; +template class Point2; } -namespace Eigen{ -template -struct ei_traits > : ei_traits > {}; +namespace Eigen { +template struct ei_traits > : ei_traits > {}; } namespace vcg { /** \addtogroup space */ /*@{*/ - /** - The templated class for representing a point in 2D space. - The class is templated over the Scalar class that is used to represent coordinates. - All the usual operator overloading (* + - ...) is present. - */ +/** + The templated class for representing a point in 2D space. + The class is templated over the Scalar class that is used to represent coordinates. + All the usual operator overloading (* + - ...) is present. + */ template class Point2 : public Eigen::Matrix<_Scalar,2,1> { +//---------------------------------------- +// template typedef part +// use it as follow: typename Point2::Type instead of simply Point2 +//---------------------------------------- +public: + typedef Eigen::Matrix<_Scalar,2,1> Type; +//---------------------------------------- +// inheritence part +//---------------------------------------- +private: typedef Eigen::Matrix<_Scalar,2,1> _Base; using _Base::coeff; using _Base::coeffRef; @@ -61,29 +70,8 @@ template class Point2 : public Eigen::Matrix<_Scalar,2,1> public: _EIGEN_GENERIC_PUBLIC_INTERFACE(Point2,_Base); - typedef Scalar ScalarType; - VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Point2) - enum {Dimension = 2}; - -//@{ - /** @name Access to Coords. - access to coords is done by overloading of [] or explicit naming of coords (X,Y,) - ("p[0]" or "p.X()" are equivalent) **/ - inline const Scalar &X() const {return data()[0];} - inline const Scalar &Y() const {return data()[1];} - inline Scalar &X() {return data()[0];} - inline Scalar &Y() {return data()[1];} - - // overloaded to return a const reference - inline const Scalar & V( const int i ) const - { - assert(i>=0 && i<2); - return data()[i]; - } -//@} - /// empty constructor (does nothing) inline Point2 () { } /// x,y constructor @@ -105,7 +93,7 @@ public: { return math::Atan2(data()[1],data()[0]); } - /// transform the point in cartesian coords into polar coords + /// transform the point in cartesian coords into polar coords inline Point2 & Cartesian2Polar() { Scalar t = Angle(); @@ -113,7 +101,7 @@ public: data()[1] = t; return *this; } - /// transform the point in polar coords into cartesian coords + /// transform the point in polar coords into cartesian coords inline Point2 & Polar2Cartesian() { Scalar l = data()[0]; @@ -121,7 +109,7 @@ public: data()[1] = (Scalar)(l*math::Sin(data()[1])); return *this; } - /// rotates the point of an angle (radiants, counterclockwise) + /// rotates the point of an angle (radiants, counterclockwise) inline Point2 & Rotate( const Scalar rad ) { Scalar t = data()[0]; @@ -133,19 +121,6 @@ public: return *this; } - - /// imports from 2D points of different types - template - inline void Import( const Point2 & b ) - { - data()[0] = b.X(); data()[1] = b.Y(); - } - /// constructs a 2D points from an existing one of different type - template - static Point2 Construct( const Point2 & b ) - { - return Point2(b.X(),b.Y()); - } }; // end class definition typedef Point2 Point2s; @@ -153,6 +128,15 @@ typedef Point2 Point2i; typedef Point2 Point2f; typedef Point2 Point2d; +// typedef Eigen::Matrix Point2s; +// typedef Eigen::Matrix Point2i; +// typedef Eigen::Matrix Point2f; +// typedef Eigen::Matrix Point2d; +// typedef Eigen::Matrix Vector2s; +// typedef Eigen::Matrix Vector2i; +// typedef Eigen::Matrix Vector2f; +// typedef Eigen::Matrix Vector2d; + /*@}*/ } // end namespace #endif diff --git a/vcg/space/point3.h b/vcg/space/point3.h index adfce8ea..d1b8e28a 100644 --- a/vcg/space/point3.h +++ b/vcg/space/point3.h @@ -29,15 +29,14 @@ #define __VCGLIB_POINT3 #include "../math/eigen.h" -#include namespace vcg{ -template class Point3; +template class Point3; } namespace Eigen{ -template -struct ei_traits > : ei_traits > {}; + +template struct ei_traits > : ei_traits > {}; template struct NumTraits > : NumTraits @@ -53,35 +52,32 @@ struct NumTraits > : NumTraits namespace vcg { +template class Box3; + /** \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; - +/** + 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 Point3 : public Eigen::Matrix<_Scalar,3,1> { - typedef Eigen::Matrix<_Scalar,3,1> _Base; - using _Base::coeff; - using _Base::coeffRef; - using _Base::setZero; - using _Base::data; - using _Base::V; - +//---------------------------------------- +// template typedef part +// use it as follow: typename Point3::Type instead of simply Point3 +//---------------------------------------- +public: + typedef Eigen::Matrix<_Scalar,3,1> Type; +//---------------------------------------- +// inheritence part +//---------------------------------------- +private: + typedef Eigen::Matrix<_Scalar,3,1> _Base; + using _Base::Construct; public: - _EIGEN_GENERIC_PUBLIC_INTERFACE(Point3,_Base); - typedef Scalar ScalarType; - VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Point3) - - 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 @@ -95,133 +91,15 @@ public: inline Point3(const Eigen::MatrixBase& other) : Base(other) {} - template - inline void Import( const Eigen::MatrixBase& b ) - { - 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(Scalar(b[0]),Scalar(b[1]),Scalar(b[2])); - } - + // this one is very useless template static inline Point3 Construct( const Q & P0, const Q & P1, const Q & P2) { return Point3(Scalar(P0),Scalar(P1),Scalar(P2)); } + vcg::Box3<_Scalar> GetBBox(vcg::Box3<_Scalar> &bb) const; - 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 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]; } - // overloaded to return a const reference - inline const Scalar & V( const int i ) const - { - assert(i>=0 && i<3); - return data()[i]; - } -//@} -//@{ - - /** @name Classical overloading of operators - Note - **/ - - // Scalatura differenziata - inline Point3 & Scale( const Scalar sx, const Scalar sy, const Scalar sz ) - { - data()[0] *= sx; - data()[1] *= sy; - data()[2] *= sz; - return *this; - } - inline Point3 & Scale( const Point3 & p ) - { - data()[0] *= p.data()[0]; - data()[1] *= p.data()[1]; - data()[2] *= p.data()[2]; - return *this; - } - - /** - * 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(Scalar &ro, Scalar &theta, Scalar &phi) const - { - ro = this->norm(); - theta = (Scalar)atan2(data()[2], data()[0]); - phi = (Scalar)asin(data()[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 Scalar &ro, const Scalar &theta, const Scalar &phi) - { - data()[0]= ro*cos(theta)*cos(phi); - data()[1]= ro*sin(phi); - data()[2]= ro*sin(theta)*cos(phi); - } - - Box3<_Scalar> GetBBox(Box3<_Scalar> &bb) const; -//@} - -}; // end class definition - - -// versione uguale alla precedente ma che assume che i due vettori sono unitari -template -inline Scalar AngleN( Point3 const & p1, Point3 const & p2 ) -{ - Scalar w = p1*p2; - if(w>1) - w = 1; - else if(w<-1) - w=-1; - return (Scalar) acos(w); -} - - -template -inline Point3 & Normalize( Point3 & p ) -{ - p.Normalize(); - return p; -} +}; // end class definition (Point3) // Dot product preciso numericamente (solo double!!) // Implementazione: si sommano i prodotti per ordine di esponente @@ -253,9 +131,7 @@ double stable_dot ( Point3 const & p0, Point3 const & p1 ) else return (k0+k1)+k2; } -} - - +} /// Point(p) Edge(v1-v2) dist, q is the point in v1-v2 with min dist template @@ -272,7 +148,6 @@ Scalar PSDist( const Point3 & p, return Distance(p,q); } - template void GetUV( Point3 &n,Point3 &u, Point3 &v, Point3 up=(Point3(0,1,0)) ) { @@ -297,23 +172,21 @@ void GetUV( Point3 &n,Point3 &u, Point3 &v, 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; -/*@}*/ +// typedef Eigen::Matrix Point3s; +// typedef Eigen::Matrix Point3i; +// typedef Eigen::Matrix Point3f; +// typedef Eigen::Matrix Point3d; +// typedef Eigen::Matrix Vector3s; +// typedef Eigen::Matrix Vector3i; +// typedef Eigen::Matrix Vector3f; +// typedef Eigen::Matrix Vector3d; } // end namespace diff --git a/vcg/space/point4.h b/vcg/space/point4.h index 17a4038a..71219c37 100644 --- a/vcg/space/point4.h +++ b/vcg/space/point4.h @@ -31,25 +31,35 @@ #include "../math/eigen.h" namespace vcg{ -template class Point4; +template class Point4; } -namespace Eigen{ -template -struct ei_traits > : ei_traits > {}; +namespace Eigen { +template struct ei_traits > : ei_traits > {}; } 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. - */ - +/** + 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 Eigen::Matrix { +//---------------------------------------- +// template typedef part +// use it as follow: typename Point4::Type instead of simply Point4 +//---------------------------------------- +public: + typedef Eigen::Matrix Type; +//---------------------------------------- +// inheritence part +//---------------------------------------- +private: typedef Eigen::Matrix _Base; using _Base::coeff; using _Base::coeffRef; @@ -63,8 +73,6 @@ public: 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) {} @@ -72,34 +80,10 @@ public: inline Point4(const Point4& p) : Base(p) {} template inline Point4(const Eigen::MatrixBase& other) : Base(other) {} - - - /// importer from different Point4 types - 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 b.template cast(); } -//@{ - - inline T &X() {return Base::x();} - inline T &Y() {return Base::y();} - inline T &Z() {return Base::z();} - inline T &W() {return Base::w();} - - // overloaded to return a const reference - inline const T & V (int i) const - { - assert(i>=0 && i<4); - return data()[i]; - } - -//@} - inline Point4 VectProd ( const Point4 &x, const Point4 &z ) const - { + { Point4 res; const Point4 &y = *this; @@ -113,7 +97,7 @@ public: z[2]+y[0]*z[1]*x[2]-x[0]*z[1]*y[2]+z[0]*x[1]*y[2]; return res; } - + //@{ /** @name Dot products **/ @@ -140,24 +124,33 @@ public: if (exp2>exp3) { math::Swap(k2,k3); math::Swap(exp2,exp3); } return ( (k0 + k1) + k2 ) +k3; - } + } //@} }; // end class definition - /// slower version of dot product, more stable (double precision only) +typedef Point4 Point4s; +typedef Point4 Point4i; +typedef Point4 Point4f; +typedef Point4 Point4d; + +// typedef Eigen::Matrix Point4s; +// typedef Eigen::Matrix Point4i; +// typedef Eigen::Matrix Point4f; +// typedef Eigen::Matrix Point4d; +// typedef Eigen::Matrix Vector4s; +// typedef Eigen::Matrix Vector4i; +// typedef Eigen::Matrix Vector4f; +// typedef Eigen::Matrix Vector4d; + +/// 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