From 0361427bc092cbd26401575af2a2ac7ed7d0eb43 Mon Sep 17 00:00:00 2001 From: cnr-isti-vclab Date: Tue, 28 Oct 2008 20:06:17 +0000 Subject: [PATCH] big cleaning in Point* and Matrix*, now they are very closed to be simple typedef of Eigen's Matrix. Now the dilema is how to mimic those typedefs, using inheritence ? or using the classic workaround: typename Point3::Type; with Point3::Type defined to Eigen::Matrix. Anyway currently I support both (and the inheritence scheme has to be preserved for compatibility). The advantage of the second approach is that when eigen has to evaluate an expression it uses an Eigen::Matrix<>, so it is probably better to only use Eigen::Matrix but I'm not 100% sure that makes a big difference especially if we add some automatic reinterpret_cast between Eigen::Matrix and vcg::Point*.... --- vcg/math/eigen.h | 17 +- vcg/math/eigen_matrix_addons.h | 143 +++++++++ ..._vcgaddons.h => eigen_matrixbase_addons.h} | 15 +- vcg/math/matrix33.h | 14 +- vcg/math/matrix44.h | 272 +++++++++--------- vcg/space/point.h | 242 +++------------- vcg/space/point2.h | 78 ++--- vcg/space/point3.h | 197 +++---------- vcg/space/point4.h | 85 +++--- 9 files changed, 461 insertions(+), 602 deletions(-) create mode 100644 vcg/math/eigen_matrix_addons.h rename vcg/math/{eigen_vcgaddons.h => eigen_matrixbase_addons.h} (98%) 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