manolas
/
vcglib
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

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<float>::Type; with Point3<T>::Type defined
to Eigen::Matrix<T,3,1>. 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*....
This commit is contained in:
Paolo Cignoni 2008-10-28 20:06:17 +00:00
parent 632f4842f1
commit 0361427bc0
9 changed files with 461 additions and 602 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
vcg/math/eigen.h
View File

@ -26,18 +26,19 @@
// TODO enable the vectorization // TODO enable the vectorization
#define EIGEN_DONT_VECTORIZE #define EIGEN_DONT_VECTORIZE
#define EIGEN_MATRIXBASE_PLUGIN <vcg/math/eigen_vcgaddons.h> #define EIGEN_MATRIXBASE_PLUGIN <vcg/math/eigen_matrixbase_addons.h>
#define EIGEN_MATRIX_PLUGIN <vcg/math/eigen_matrix_addons.h>
// forward declarations // forward declarations
namespace Eigen { namespace Eigen {
template<typename Derived1, typename Derived2, int Size> struct ei_lexi_comparison; template<typename Derived1, typename Derived2, int Size> struct ei_lexi_comparison;
} }
#include "base.h"
#include "../Eigen/LU" #include "../Eigen/LU"
#include "../Eigen/Geometry" #include "../Eigen/Geometry"
#include "../Eigen/Array" #include "../Eigen/Array"
#include "../Eigen/Core" #include "../Eigen/Core"
#include "base.h"
// add support for unsigned char and short int // add support for unsigned char and short int
namespace Eigen { namespace Eigen {
@ -239,6 +240,18 @@ inline typename Eigen::ei_traits<Derived1>::Scalar
SquaredDistance(const Eigen::MatrixBase<Derived1>& p1, const Eigen::MatrixBase<Derived2> & p2) SquaredDistance(const Eigen::MatrixBase<Derived1>& p1, const Eigen::MatrixBase<Derived2> & p2)
{ return (p1-p2).norm2(); } { return (p1-p2).norm2(); }
template<typename Derived>
inline const Eigen::CwiseUnaryOp<Eigen::ei_scalar_abs_op<typename Eigen::ei_traits<Derived>::Scalar>, Derived>
Abs(const Eigen::MatrixBase<Derived>& p)
{ return p.cwise().abs(); }
template<typename Derived>
inline const Eigen::CwiseBinaryOp<Eigen::ei_scalar_max_op<typename Eigen::ei_traits<Derived>::Scalar>,
Derived,
Eigen::NestByValue<typename Derived::ConstantReturnType> >
LowClampToZero(const Eigen::MatrixBase<Derived>& p)
{ return p.cwise().max(Derived::Zero().nestByValue()); }
} }
#endif #endif

143
vcg/math/eigen_matrix_addons.h Normal file
View File

@ -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<typename OtherDerived>
inline void Import(const MatrixBase<OtherDerived>& 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<typename OtherDerived>
inline void ImportHomo(const MatrixBase<OtherDerived>& 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<SizeAtCompileTime-1> = b;
data()[SizeAtCompileTime-1] = Scalar(1.0);
}
/// constructor for points with different scalar type and-or dimensionality
template<typename OtherDerived>
static inline Matrix Construct(const MatrixBase<OtherDerived>& b)
{
Matrix p; p.Import(b);
return p;
}
/// constructor for homogeneus point.
template<typename OtherDerived>
static inline Matrix ConstructHomo(const MatrixBase<OtherDerived>& 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 && i<SizeAtCompileTime);
return data()[i];
}
//--------------------------------------------------------------------------------
// SPACE
//--------------------------------------------------------------------------------
/** Local to Glocal
* (provided for uniformity with other spatial classes. trivial for points) */
inline Matrix LocalToGlobal(ParamType p) const { return *this; }
/** Glocal to Local
* (provided for uniformity with other spatial classes. trivial for points) */
inline ParamType GlobalToLocal(PointType /*p*/) const { return ParamType(); }
/**
* 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
{
EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Matrix,3);
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)
{
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);
}

15
vcg/math/eigen_vcgaddons.h → vcg/math/eigen_matrixbase_addons.h
View File

@ -22,15 +22,15 @@
****************************************************************************/ ****************************************************************************/
#warning You are including deprecated math stuff #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(); }; EIGEN_DEPRECATED inline unsigned int ColumnsNumber() const { return cols(); };
/*! /*! \deprecated use rows() */
* \deprecated use rows()
*/
EIGEN_DEPRECATED inline unsigned int RowsNumber() const { return 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)); row(i).swap(row(j));
}; };
Scalar* V() { return derived().data(); }
const Scalar* V() const { return derived().data(); }
/*! /*!
* \deprecated use *this.cwise() += k * \deprecated use *this.cwise() += k
* (Modifier) Add to each element of this matrix the scalar constant <I>k</I>. * (Modifier) Add to each element of this matrix the scalar constant <I>k</I>.

14
vcg/math/matrix33.h
View File

@ -79,19 +79,14 @@ public:
template<typename OtherDerived> template<typename OtherDerived>
Matrix33(const Eigen::MatrixBase<OtherDerived>& other) : Base(other) {} Matrix33(const Eigen::MatrixBase<OtherDerived>& other) : Base(other) {}
/*! /*! \deprecated use *this.row(i) */
* \deprecated use *this.row(i)
*/
inline typename Base::RowXpr operator[](const unsigned int i) inline typename Base::RowXpr operator[](const unsigned int i)
{ return Base::row(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 inline const typename Base::RowXpr operator[](const unsigned int i) const
{ return Base::row(i); } { return Base::row(i); }
/** \deprecated */ /** \deprecated */
Matrix33 & SetRotateRad(Scalar angle, const Point3<Scalar> & axis ) Matrix33 & SetRotateRad(Scalar angle, const Point3<Scalar> & axis )
{ {
@ -129,7 +124,7 @@ public:
*/ */
Scalar Norm() { return Base::cwise().abs2().sum(); } 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; // Scalar SQsum=0;
// for(int i=0;i<3;++i) // for(int i=0;i<3;++i)
// for(int j=0;j<3;++j) // for(int j=0;j<3;++j)
@ -138,8 +133,7 @@ public:
// } // }
/** /** Computes the covariance matrix of a set of 3d points. Returns the barycenter.
It computes the covariance matrix of a set of 3d points. Returns the barycenter
*/ */
// FIXME should be outside Matrix // FIXME should be outside Matrix
template <class STLPOINTCONTAINER > template <class STLPOINTCONTAINER >

272
vcg/math/matrix44.h
View File

@ -45,7 +45,7 @@ struct ei_traits<vcg::Matrix44<Scalar> > : ei_traits<Eigen::Matrix<Scalar,4,4,Ro
namespace vcg { namespace vcg {
/* /*
Annotations: Annotations:
Opengl stores matrix in column-major order. That is, the matrix is stored as: Opengl stores matrix in column-major order. That is, the matrix is stored as:
@ -54,10 +54,10 @@ Opengl stores matrix in column-major order. That is, the matrix is stored as:
a2 a6 a10 a14 a2 a6 a10 a14
a3 a7 a11 a15 a3 a7 a11 a15
Usually in opengl (see opengl specs) vectors are 'column' vectors Usually in opengl (see opengl specs) vectors are 'column' vectors
so usually matrix are PRE-multiplied for a vector. so usually matrix are PRE-multiplied for a vector.
So the command glTranslate generate a matrix that So the command glTranslate generate a matrix that
is ready to be premultipled for a vector: is ready to be premultipled for a vector:
1 0 0 tx 1 0 0 tx
0 1 0 ty 0 1 0 ty
@ -80,8 +80,10 @@ for 'column' vectors.
*/ */
/** This class represents a 4x4 matrix. T is the kind of element in the matrix. /** \deprecated use Eigen::Matrix<Scalar,4,4> (or the typedef) you want a real 4x4 matrix, or use Eigen::Transform<Scalar,3> if you want a transformation matrix for a 3D space (a Eigen::Transform<Scalar,3> is internally a 4x4 col-major matrix)
*/ *
* This class represents a 4x4 matrix. T is the kind of element in the matrix.
*/
template<typename _Scalar> template<typename _Scalar>
class Matrix44 : public Eigen::Matrix<_Scalar,4,4,Eigen::RowMajor> // FIXME col or row major ! 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::coeffRef;
using _Base::ElementAt; using _Base::ElementAt;
using _Base::setZero; using _Base::setZero;
using _Base::operator*;
public: public:
@ -103,7 +106,6 @@ public:
~Matrix44() {} ~Matrix44() {}
Matrix44(const Matrix44 &m) : Base(m) {} Matrix44(const Matrix44 &m) : Base(m) {}
Matrix44(const Scalar * v ) : Base(Eigen::Map<Eigen::Matrix<Scalar,4,4,Eigen::RowMajor> >(v)) {} Matrix44(const Scalar * v ) : Base(Eigen::Map<Eigen::Matrix<Scalar,4,4,Eigen::RowMajor> >(v)) {}
template<typename OtherDerived> template<typename OtherDerived>
Matrix44(const Eigen::MatrixBase<OtherDerived>& other) : Base(other) {} Matrix44(const Eigen::MatrixBase<OtherDerived>& other) : Base(other) {}
@ -116,7 +118,7 @@ public:
typename Base::RowXpr GetRow4(const int& i) const { return Base::row(i); } typename Base::RowXpr GetRow4(const int& i) const { return Base::row(i); }
Eigen::Block<Base,1,3> GetRow3(const int& i) const { return this->template block<1,3>(i,0); } Eigen::Block<Base,1,3> GetRow3(const int& i) const { return this->template block<1,3>(i,0); }
template <class Matrix44Type> template <class Matrix44Type>
void ToMatrix(Matrix44Type & m) const { m = (*this).template cast<typename Matrix44Type::Scalar>(); } void ToMatrix(Matrix44Type & m) const { m = (*this).template cast<typename Matrix44Type::Scalar>(); }
void ToEulerAngles(Scalar &alpha, Scalar &beta, Scalar &gamma); 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 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 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 Scalar sx, const Scalar sy, const Scalar sz);
Matrix44 &SetScale(const Point3<Scalar> &t); Matrix44 &SetScale(const Point3<Scalar> &t);
Matrix44 &SetTranslate(const Point3<Scalar> &t); Matrix44 &SetTranslate(const Point3<Scalar> &t);
Matrix44 &SetTranslate(const Scalar sx, const Scalar sy, const Scalar sz); Matrix44 &SetTranslate(const Scalar sx, const Scalar sy, const Scalar sz);
Matrix44 &SetShearXY(const Scalar sz); Matrix44 &SetShearXY(const Scalar sz);
Matrix44 &SetShearXZ(const Scalar sy); Matrix44 &SetShearXZ(const Scalar sy);
Matrix44 &SetShearYZ(const Scalar sx); Matrix44 &SetShearYZ(const Scalar sx);
///use radiants for angle. ///use radiants for angle.
Matrix44 &SetRotateDeg(Scalar AngleDeg, const Point3<Scalar> & axis); Matrix44 &SetRotateDeg(Scalar AngleDeg, const Point3<Scalar> & axis);
Matrix44 &SetRotateRad(Scalar AngleRad, const Point3<Scalar> & axis); Matrix44 &SetRotateRad(Scalar AngleRad, const Point3<Scalar> & axis);
template <class Q> void Import(const Matrix44<Q> &m) { template <class Q> void Import(const Matrix44<Q> &m) {
for(int i = 0; i < 16; i++) for(int i = 0; i < 16; i++)
Base::data()[i] = (Scalar)(m.data()[i]); Base::data()[i] = (Scalar)(m.data()[i]);
} }
template <class Q> template <class Q>
static inline Matrix44 Construct( const Matrix44<Q> & b ) static inline Matrix44 Construct( const Matrix44<Q> & b )
{ {
Matrix44 tmp; tmp.FromMatrix(b); Matrix44 tmp; tmp.FromMatrix(b);
return tmp; return tmp;
} }
// template <class T> Point3<T> operator*(const Point3<T> &p) {
// T w;
// Point3<T> 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<Scalar,3,1> operator * (const Eigen::Matrix<Scalar,3,1>& p) const {
Scalar w;
Eigen::Matrix<Scalar,3,1> 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. */ /** Class for solving A * x = b. */
template <class T> class LinearSolve: public Matrix44<T> { template <class T> class LinearSolve: public Matrix44<T> {
public: public:
LinearSolve(const Matrix44<T> &m); LinearSolve(const Matrix44<T> &m);
Point4<T> Solve(const Point4<T> &b); // solve A <20> x = b Point4<T> Solve(const Point4<T> &b); // solve A <20> x = b
///If you need to solve some equation you can use this function instead of Matrix44 one for speed. ///If you need to solve some equation you can use this function instead of Matrix44 one for speed.
T Determinant() const; T Determinant() const;
protected: protected:
///Holds row permutation. ///Holds row permutation.
int index[4]; //hold permutation int index[4]; //hold permutation
///Hold sign of row permutation (used for determinant sign) ///Hold sign of row permutation (used for determinant sign)
T d; T d;
bool Decompose(); bool Decompose();
}; };
/*** Postmultiply */ /*** Postmultiply */
//template <class T> Point3<T> operator*(const Point3<T> &p, const Matrix44<T> &m); //template <class T> Point3<T> operator*(const Point3<T> &p, const Matrix44<T> &m);
///Premultiply ///Premultiply
template <class T> Point3<T> operator*(const Matrix44<T> &m, const Point3<T> &p); // template <class T> Point3<T> operator*(const Matrix44<T> &m, const Point3<T> &p);
//return NULL matrix if not invertible //return NULL matrix if not invertible
template <class T> Matrix44<T> &Invert(Matrix44<T> &m); template <class T> Matrix44<T> &Invert(Matrix44<T> &m);
@ -196,9 +219,9 @@ typedef Matrix44<double> Matrix44d;
template < class PointType , class T > void operator*=( std::vector<PointType> &vert, const Matrix44<T> & m ) { template < class PointType , class T > void operator*=( std::vector<PointType> &vert, const Matrix44<T> & m ) {
typename std::vector<PointType>::iterator ii; typename std::vector<PointType>::iterator ii;
for(ii=vert.begin();ii!=vert.end();++ii) for(ii=vert.begin();ii!=vert.end();++ii)
(*ii).P()=m * (*ii).P(); (*ii).P()=m * (*ii).P();
} }
template <class T> template <class T>
@ -237,36 +260,36 @@ void Matrix44<T>::FromEulerAngles(Scalar alpha, Scalar beta, Scalar gamma)
} }
template <class T> void Matrix44<T>::SetDiagonal(const Scalar k) { template <class T> void Matrix44<T>::SetDiagonal(const Scalar k) {
setZero(); setZero();
ElementAt(0, 0) = k; ElementAt(0, 0) = k;
ElementAt(1, 1) = k; ElementAt(1, 1) = k;
ElementAt(2, 2) = k; ElementAt(2, 2) = k;
ElementAt(3, 3) = 1; ElementAt(3, 3) = 1;
} }
template <class T> Matrix44<T> &Matrix44<T>::SetScale(const Point3<Scalar> &t) { template <class T> Matrix44<T> &Matrix44<T>::SetScale(const Point3<Scalar> &t) {
SetScale(t[0], t[1], t[2]); SetScale(t[0], t[1], t[2]);
return *this; return *this;
} }
template <class T> Matrix44<T> &Matrix44<T>::SetScale(const Scalar sx, const Scalar sy, const Scalar sz) { template <class T> Matrix44<T> &Matrix44<T>::SetScale(const Scalar sx, const Scalar sy, const Scalar sz) {
setZero(); setZero();
ElementAt(0, 0) = sx; ElementAt(0, 0) = sx;
ElementAt(1, 1) = sy; ElementAt(1, 1) = sy;
ElementAt(2, 2) = sz; ElementAt(2, 2) = sz;
ElementAt(3, 3) = 1; ElementAt(3, 3) = 1;
return *this; return *this;
} }
template <class T> Matrix44<T> &Matrix44<T>::SetTranslate(const Point3<Scalar> &t) { template <class T> Matrix44<T> &Matrix44<T>::SetTranslate(const Point3<Scalar> &t) {
SetTranslate(t[0], t[1], t[2]); SetTranslate(t[0], t[1], t[2]);
return *this; return *this;
} }
template <class T> Matrix44<T> &Matrix44<T>::SetTranslate(const Scalar tx, const Scalar ty, const Scalar tz) { template <class T> Matrix44<T> &Matrix44<T>::SetTranslate(const Scalar tx, const Scalar ty, const Scalar tz) {
Base::setIdentity(); Base::setIdentity();
ElementAt(0, 3) = tx; ElementAt(0, 3) = tx;
ElementAt(1, 3) = ty; ElementAt(1, 3) = ty;
ElementAt(2, 3) = tz; ElementAt(2, 3) = tz;
return *this; return *this;
} }
template <class T> Matrix44<T> &Matrix44<T>::SetRotateDeg(Scalar AngleDeg, const Point3<Scalar> & axis) { template <class T> Matrix44<T> &Matrix44<T>::SetRotateDeg(Scalar AngleDeg, const Point3<Scalar> & axis) {
@ -274,9 +297,9 @@ template <class T> Matrix44<T> &Matrix44<T>::SetRotateDeg(Scalar AngleDeg, const
} }
template <class T> Matrix44<T> &Matrix44<T>::SetRotateRad(Scalar AngleRad, const Point3<Scalar> & axis) { template <class T> Matrix44<T> &Matrix44<T>::SetRotateRad(Scalar AngleRad, const Point3<Scalar> & axis) {
//angle = angle*(T)3.14159265358979323846/180; e' in radianti! //angle = angle*(T)3.14159265358979323846/180; e' in radianti!
T c = math::Cos(AngleRad); T c = math::Cos(AngleRad);
T s = math::Sin(AngleRad); T s = math::Sin(AngleRad);
T q = 1-c; T q = 1-c;
Point3<T> t = axis; Point3<T> t = axis;
t.Normalize(); t.Normalize();
@ -296,94 +319,94 @@ template <class T> Matrix44<T> &Matrix44<T>::SetRotateRad(Scalar AngleRad, const
ElementAt(3,1) = 0; ElementAt(3,1) = 0;
ElementAt(3,2) = 0; ElementAt(3,2) = 0;
ElementAt(3,3) = 1; ElementAt(3,3) = 1;
return *this; return *this;
} }
/* Shear Matrixes /* Shear Matrixes
XY XY
1 k 0 0 x x+ky 1 k 0 0 x x+ky
0 1 0 0 y y 0 1 0 0 y y
0 0 1 0 z z 0 0 1 0 z z
0 0 0 1 1 1 0 0 0 1 1 1
1 0 k 0 x x+kz 1 0 k 0 x x+kz
0 1 0 0 y y 0 1 0 0 y y
0 0 1 0 z z 0 0 1 0 z z
0 0 0 1 1 1 0 0 0 1 1 1
1 1 0 0 x x 1 1 0 0 x x
0 1 k 0 y y+kz 0 1 k 0 y y+kz
0 0 1 0 z z 0 0 1 0 z z
0 0 0 1 1 1 0 0 0 1 1 1
*/ */
template <class T> Matrix44<T> & Matrix44<T>::SetShearXY( const Scalar sh) {// shear the X coordinate as the Y coordinate change template <class T> Matrix44<T> & Matrix44<T>::SetShearXY( const Scalar sh) {// shear the X coordinate as the Y coordinate change
Base::setIdentity(); Base::setIdentity();
ElementAt(0,1) = sh; ElementAt(0,1) = sh;
return *this; return *this;
} }
template <class T> Matrix44<T> & Matrix44<T>::SetShearXZ( const Scalar sh) {// shear the X coordinate as the Z coordinate change template <class T> Matrix44<T> & Matrix44<T>::SetShearXZ( const Scalar sh) {// shear the X coordinate as the Z coordinate change
Base::setIdentity(); Base::setIdentity();
ElementAt(0,2) = sh; ElementAt(0,2) = sh;
return *this; return *this;
} }
template <class T> Matrix44<T> &Matrix44<T>::SetShearYZ( const Scalar sh) {// shear the Y coordinate as the Z coordinate change template <class T> Matrix44<T> &Matrix44<T>::SetShearYZ( const Scalar sh) {// shear the Y coordinate as the Z coordinate change
Base::setIdentity(); Base::setIdentity();
ElementAt(1,2) = sh; 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] ) Given a non singular, non projective matrix (e.g. with the last row equal to [0,0,0,1] )
This procedure decompose it in a sequence of This procedure decompose it in a sequence of
Scale,Shear,Rotation e Translation Scale,Shear,Rotation e Translation
- ScaleV and Tranv are obiviously scaling and translation. - ScaleV and Tranv are obiviously scaling and translation.
- ShearV contains three scalars with, respectively - 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 - 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: Example Code:
double srv() { return (double(rand()%40)-20)/2.0; } // small random value double srv() { return (double(rand()%40)-20)/2.0; } // small random value
srand(time(0)); srand(time(0));
Point3d ScV(10+srv(),10+srv(),10+srv()),ScVOut(-1,-1,-1); Point3d ScV(10+srv(),10+srv(),10+srv()),ScVOut(-1,-1,-1);
Point3d ShV(srv(),srv(),srv()),ShVOut(-1,-1,-1); Point3d ShV(srv(),srv(),srv()),ShVOut(-1,-1,-1);
Point3d RtV(10+srv(),srv(),srv()),RtVOut(-1,-1,-1); Point3d RtV(10+srv(),srv(),srv()),RtVOut(-1,-1,-1);
Point3d TrV(srv(),srv(),srv()),TrVOut(-1,-1,-1); Point3d TrV(srv(),srv(),srv()),TrVOut(-1,-1,-1);
Matrix44d Scl; Scl.SetScale(ScV); Matrix44d Scl; Scl.SetScale(ScV);
Matrix44d Sxy; Sxy.SetShearXY(ShV[0]); Matrix44d Sxy; Sxy.SetShearXY(ShV[0]);
Matrix44d Sxz; Sxz.SetShearXZ(ShV[1]); Matrix44d Sxz; Sxz.SetShearXZ(ShV[1]);
Matrix44d Syz; Syz.SetShearYZ(ShV[2]); 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 Rty; Rty.SetRotate(math::ToRad(RtV[1]),Point3d(0,1,0));
Matrix44d Rtz; Rtz.SetRotate(math::ToRad(RtV[2]),Point3d(0,0,1)); Matrix44d Rtz; Rtz.SetRotate(math::ToRad(RtV[2]),Point3d(0,0,1));
Matrix44d Trn; Trn.SetTranslate(TrV); Matrix44d Trn; Trn.SetTranslate(TrV);
Matrix44d StartM = Trn * Rtx*Rty*Rtz * Syz*Sxz*Sxy *Scl; Matrix44d StartM = Trn * Rtx*Rty*Rtz * Syz*Sxz*Sxy *Scl;
Matrix44d ResultM=StartM; Matrix44d ResultM=StartM;
Decompose(ResultM,ScVOut,ShVOut,RtVOut,TrVOut); Decompose(ResultM,ScVOut,ShVOut,RtVOut,TrVOut);
Scl.SetScale(ScVOut); Scl.SetScale(ScVOut);
Sxy.SetShearXY(ShVOut[0]); Sxy.SetShearXY(ShVOut[0]);
Sxz.SetShearXZ(ShVOut[1]); Sxz.SetShearXZ(ShVOut[1]);
Syz.SetShearYZ(ShVOut[2]); Syz.SetShearYZ(ShVOut[2]);
Rtx.SetRotate(math::ToRad(RtVOut[0]),Point3d(1,0,0)); Rtx.SetRotate(math::ToRad(RtVOut[0]),Point3d(1,0,0));
Rty.SetRotate(math::ToRad(RtVOut[1]),Point3d(0,1,0)); Rty.SetRotate(math::ToRad(RtVOut[1]),Point3d(0,1,0));
Rtz.SetRotate(math::ToRad(RtVOut[2]),Point3d(0,0,1)); Rtz.SetRotate(math::ToRad(RtVOut[2]),Point3d(0,0,1));
Trn.SetTranslate(TrVOut); 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 ; Matrix44d RebuildM = Trn * Rtx*Rty*Rtz * Syz*Sxz*Sxy * Scl ;
*/ */
template <class T> template <class T>
@ -393,7 +416,7 @@ bool Decompose(Matrix44<T> &M, Point3<T> &ScaleV, Point3<T> &ShearV, Point3<T> &
return false; return false;
if(math::Abs(M.Determinant())<1e-10) return false; // matrix should be at least invertible... 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); TranV=M.GetColumn3(3);
// Second Step Recover Scale and Shearing interleaved // Second Step Recover Scale and Shearing interleaved
@ -404,7 +427,7 @@ bool Decompose(Matrix44<T> &M, Point3<T> &ScaleV, Point3<T> &ShearV, Point3<T> &
ShearV[0]=R[0].dot(M.GetColumn3(1)); // xy shearing ShearV[0]=R[0].dot(M.GetColumn3(1)); // xy shearing
R[1]= M.GetColumn3(1)-R[0]*ShearV[0]; 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 ScaleV[1]=Norm(R[1]); // y scaling
R[1]=R[1]/ScaleV[1]; R[1]=R[1]/ScaleV[1];
ShearV[0]=ShearV[0]/ScaleV[1]; ShearV[0]=ShearV[0]/ScaleV[1];
@ -425,16 +448,16 @@ bool Decompose(Matrix44<T> &M, Point3<T> &ScaleV, Point3<T> &ShearV, Point3<T> &
ShearV[2]=R[1].dot(M.GetColumn3(2)); // yz shearing ShearV[2]=R[1].dot(M.GetColumn3(2)); // yz shearing
ShearV[2]=ShearV[2]/ScaleV[2]; ShearV[2]=ShearV[2]/ScaleV[2];
int i,j; int i,j;
for(i=0;i<3;++i) for(i=0;i<3;++i)
for(j=0;j<3;++j) for(j=0;j<3;++j)
M(i,j)=R[j][i]; M(i,j)=R[j][i];
// Third and last step: Recover the rotation // Third and last step: Recover the rotation
//now the matrix should be a pure rotation matrix so its determinant is +-1 //now the matrix should be a pure rotation matrix so its determinant is +-1
double det=M.Determinant(); double det=M.Determinant();
if(math::Abs(det)<1e-10) return false; // matrix should be at least invertible... 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... assert(math::Abs(math::Abs(det)-1.0)<1e-10); // it should be +-1...
if(det<0) { if(det<0) {
ScaleV *= -1; ScaleV *= -1;
M *= -1; M *= -1;
@ -443,37 +466,28 @@ bool Decompose(Matrix44<T> &M, Point3<T> &ScaleV, Point3<T> &ShearV, Point3<T> &
double alpha,beta,gamma; // rotations around the x,y and z axis double alpha,beta,gamma; // rotations around the x,y and z axis
beta=asin( M(0,2)); beta=asin( M(0,2));
double cosbeta=cos(beta); double cosbeta=cos(beta);
if(math::Abs(cosbeta) > 1e-5) if(math::Abs(cosbeta) > 1e-5)
{ {
alpha=asin(-M(1,2)/cosbeta); alpha=asin(-M(1,2)/cosbeta);
if((M(2,2)/cosbeta) < 0 ) alpha=M_PI-alpha; if((M(2,2)/cosbeta) < 0 ) alpha=M_PI-alpha;
gamma=asin(-M(0,1)/cosbeta); gamma=asin(-M(0,1)/cosbeta);
if((M(0,0)/cosbeta)<0) gamma = M_PI-gamma; if((M(0,0)/cosbeta)<0) gamma = M_PI-gamma;
} }
else else
{ {
alpha=asin(-M(1,0)); alpha=asin(-M(1,0));
if(M(1,1)<0) alpha=M_PI-alpha; if(M(1,1)<0) alpha=M_PI-alpha;
gamma=0; gamma=0;
} }
RotV[0]=math::ToDeg(alpha); RotV[0]=math::ToDeg(alpha);
RotV[1]=math::ToDeg(beta); RotV[1]=math::ToDeg(beta);
RotV[2]=math::ToDeg(gamma); RotV[2]=math::ToDeg(gamma);
return true; return true;
} }
template <class T> Point3<T> operator*(const Matrix44<T> &m, const Point3<T> &p) {
T w;
Point3<T> 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 <class T> Point3<T> operator*(const Matrix44<T> &m, const Point3<T> &p)
/* /*
To invert a matrix you can To invert a matrix you can
either invert the matrix inplace calling 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 <class T> Matrix44<T> & Invert(Matrix44<T> &m) { template <class T> Matrix44<T> & Invert(Matrix44<T> &m) {

242
vcg/space/point.h
View File

@ -33,18 +33,18 @@
#include <vcg/space/space.h> #include <vcg/space/space.h>
namespace vcg{ namespace vcg{
template<class Scalar> class Point4; template<typename Scalar> class Point2;
} template<typename Scalar> class Point3;
template<typename Scalar> class Point4;
namespace Eigen{ namespace ndim{
template<typename Scalar,int Size> template <int Size, typename Scalar> class Point;
struct ei_traits<vcg::Point<Scalar,Size> > : ei_traits<Eigen::Matrix<Scalar,Size,1> > {}; }
} }
namespace vcg { namespace vcg {
namespace ndim{ namespace ndim{
/** \addtogroup space */ /** \addtogroup space */
/*@{*/ /*@{*/
/** /**
@ -55,7 +55,17 @@ namespace ndim{
*/ */
template <int N, class S> class Point : public Eigen::Matrix<S,N,1> template <int N, class S> class Point : public Eigen::Matrix<S,N,1>
{ {
typedef Eigen::Matrix<T,N,1> _Base; //----------------------------------------
// template typedef part
// use it as follow: typename Point<N,S>::Type instead of simply Point<N,S>
//----------------------------------------
public:
typedef Eigen::Matrix<S,N,1> Type;
//----------------------------------------
// inheritence part
//----------------------------------------
private:
typedef Eigen::Matrix<S,N,1> _Base;
using _Base::coeff; using _Base::coeff;
using _Base::coeffRef; using _Base::coeffRef;
using _Base::setZero; using _Base::setZero;
@ -65,212 +75,50 @@ template <int N, class S> class Point : public Eigen::Matrix<S,N,1>
public: public:
_EIGEN_GENERIC_PUBLIC_INTERFACE(Point,_Base); _EIGEN_GENERIC_PUBLIC_INTERFACE(Point,_Base);
typedef S ScalarType;
typedef VoidType ParamType;
typedef Point<N,S> PointType;
enum {Dimension = N};
VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Point) 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() {} inline Point() : Base() {}
template<typename OtherDerived> template<typename OtherDerived>
inline Point(const Eigen::MatrixBase<OtherDerived>& other) : Base(other) {} inline Point(const Eigen::MatrixBase<OtherDerived>& 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 <int N2, class S2>
inline void Import( const Point<N2,S2> & 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 <int N2, class S2>
static inline PointType Construct( const Point<N2,S2> & b )
{
PointType p; p.Import(b);
return p;
}
/// importer for homogeneous points
template <class S2>
inline void ImportHomo( const Point<N-1,S2> & 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 <int N2, class S2>
static inline PointType Construct( const Point<N-1,S2> & 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)
**/
/// slower version, more stable (double precision only) /// slower version, more stable (double precision only)
inline S StableDot ( const PointType & p ) const; inline S StableDot (const Point& p) const;
//@}
//@{
/** @name Norms
**/
/// Euclidean norm, static version
template <class PT> static S Norm(const PT &p );
/// Squared Euclidean norm, static version
template <class PT> static S SquaredNorm(const PT &p );
/// Normalization (division by norm), static version
template <class PT> static PointType & Normalize(const PT &p);
//@}
/// Signed area operator /// Signed area operator
/// a % b returns the signed area of the parallelogram inside a and b /// a % b returns the signed area of the parallelogram inside a and b
// inline S operator % ( PointType const & p ) const; // 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; }
//@}
}; // end class definition }; // end class definition
// workaround the lack of template typedef (the next c++ standard will support them :) ) typedef Eigen::Matrix<short ,2,1> Point2s;
typedef Eigen::Matrix<int ,2,1> Point2i;
typedef Eigen::Matrix<float ,2,1> Point2f;
typedef Eigen::Matrix<double,2,1> Point2d;
typedef Eigen::Matrix<short ,2,1> Vector2s;
typedef Eigen::Matrix<int ,2,1> Vector2i;
typedef Eigen::Matrix<float ,2,1> Vector2f;
typedef Eigen::Matrix<double,2,1> Vector2d;
template <typename S> typedef Eigen::Matrix<short ,3,1> Point3s;
struct Point2:public Point<2,S>{ typedef Eigen::Matrix<int ,3,1> Point3i;
typedef Point<3,S> Base; typedef Eigen::Matrix<float ,3,1> Point3f;
inline Point2() : Base() {}; typedef Eigen::Matrix<double,3,1> Point3d;
inline Point2(const Point2& p):Base(p){}; typedef Eigen::Matrix<short ,3,1> Vector3s;
inline Point2(S a, S b):Base(a,b){}; typedef Eigen::Matrix<int ,3,1> Vector3i;
template<typename OtherDerived> typedef Eigen::Matrix<float ,3,1> Vector3f;
inline Point2(const Eigen::MatrixBase<OtherDerived>& other) : Base(other) {} typedef Eigen::Matrix<double,3,1> Vector3d;
};
template <typename S>
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<typename OtherDerived>
inline Point3(const Eigen::MatrixBase<OtherDerived>& other) : Base(other) {}
};
template <typename S> typedef Eigen::Matrix<short ,4,1> Point4s;
struct Point4:public Point<4,S>{ typedef Eigen::Matrix<int ,4,1> Point4i;
typedef Point<4,S> Base; typedef Eigen::Matrix<float ,4,1> Point4f;
inline Point4() : Base() {}; typedef Eigen::Matrix<double,4,1> Point4d;
inline Point4(const Point4& p):Base(p) {} typedef Eigen::Matrix<short ,4,1> Vector4s;
inline Point4(S a, S b, S c, S d):Base(a,b,c,d){}; typedef Eigen::Matrix<int ,4,1> Vector4i;
template<typename OtherDerived> typedef Eigen::Matrix<float ,4,1> Vector4f;
inline Point4(const Eigen::MatrixBase<OtherDerived>& other) : Base(other) {} typedef Eigen::Matrix<double,4,1> 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;
/*@}*/ /*@}*/

78
vcg/space/point2.h
View File

@ -29,28 +29,37 @@
#define __VCGLIB_POINT2 #define __VCGLIB_POINT2
#include "../math/eigen.h" #include "../math/eigen.h"
#include <vcg/math/base.h> // #include "point.h"
namespace vcg{ namespace vcg{
template<class Scalar> class Point2; template<typename Scalar> class Point2;
} }
namespace Eigen{ namespace Eigen {
template<typename Scalar> template<typename Scalar> struct ei_traits<vcg::Point2<Scalar> > : ei_traits<Eigen::Matrix<Scalar,2,1> > {};
struct ei_traits<vcg::Point2<Scalar> > : ei_traits<Eigen::Matrix<Scalar,2,1> > {};
} }
namespace vcg { namespace vcg {
/** \addtogroup space */ /** \addtogroup space */
/*@{*/ /*@{*/
/** /**
The templated class for representing a point in 2D 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. The class is templated over the Scalar class that is used to represent coordinates.
All the usual operator overloading (* + - ...) is present. All the usual operator overloading (* + - ...) is present.
*/ */
template <class _Scalar> class Point2 : public Eigen::Matrix<_Scalar,2,1> template <class _Scalar> class Point2 : public Eigen::Matrix<_Scalar,2,1>
{ {
//----------------------------------------
// template typedef part
// use it as follow: typename Point2<S>::Type instead of simply Point2<S>
//----------------------------------------
public:
typedef Eigen::Matrix<_Scalar,2,1> Type;
//----------------------------------------
// inheritence part
//----------------------------------------
private:
typedef Eigen::Matrix<_Scalar,2,1> _Base; typedef Eigen::Matrix<_Scalar,2,1> _Base;
using _Base::coeff; using _Base::coeff;
using _Base::coeffRef; using _Base::coeffRef;
@ -61,29 +70,8 @@ template <class _Scalar> class Point2 : public Eigen::Matrix<_Scalar,2,1>
public: public:
_EIGEN_GENERIC_PUBLIC_INTERFACE(Point2,_Base); _EIGEN_GENERIC_PUBLIC_INTERFACE(Point2,_Base);
typedef Scalar ScalarType;
VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Point2) 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) /// empty constructor (does nothing)
inline Point2 () { } inline Point2 () { }
/// x,y constructor /// x,y constructor
@ -105,7 +93,7 @@ public:
{ {
return math::Atan2(data()[1],data()[0]); 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() inline Point2 & Cartesian2Polar()
{ {
Scalar t = Angle(); Scalar t = Angle();
@ -113,7 +101,7 @@ public:
data()[1] = t; data()[1] = t;
return *this; return *this;
} }
/// transform the point in polar coords into cartesian coords /// transform the point in polar coords into cartesian coords
inline Point2 & Polar2Cartesian() inline Point2 & Polar2Cartesian()
{ {
Scalar l = data()[0]; Scalar l = data()[0];
@ -121,7 +109,7 @@ public:
data()[1] = (Scalar)(l*math::Sin(data()[1])); data()[1] = (Scalar)(l*math::Sin(data()[1]));
return *this; return *this;
} }
/// rotates the point of an angle (radiants, counterclockwise) /// rotates the point of an angle (radiants, counterclockwise)
inline Point2 & Rotate( const Scalar rad ) inline Point2 & Rotate( const Scalar rad )
{ {
Scalar t = data()[0]; Scalar t = data()[0];
@ -133,19 +121,6 @@ public:
return *this; return *this;
} }
/// imports from 2D points of different types
template <class T>
inline void Import( const Point2<T> & b )
{
data()[0] = b.X(); data()[1] = b.Y();
}
/// constructs a 2D points from an existing one of different type
template <class T>
static Point2 Construct( const Point2<T> & b )
{
return Point2(b.X(),b.Y());
}
}; // end class definition }; // end class definition
typedef Point2<short> Point2s; typedef Point2<short> Point2s;
@ -153,6 +128,15 @@ typedef Point2<int> Point2i;
typedef Point2<float> Point2f; typedef Point2<float> Point2f;
typedef Point2<double> Point2d; typedef Point2<double> Point2d;
// typedef Eigen::Matrix<short ,2,1> Point2s;
// typedef Eigen::Matrix<int ,2,1> Point2i;
// typedef Eigen::Matrix<float ,2,1> Point2f;
// typedef Eigen::Matrix<double,2,1> Point2d;
// typedef Eigen::Matrix<short ,2,1> Vector2s;
// typedef Eigen::Matrix<int ,2,1> Vector2i;
// typedef Eigen::Matrix<float ,2,1> Vector2f;
// typedef Eigen::Matrix<double,2,1> Vector2d;
/*@}*/ /*@}*/
} // end namespace } // end namespace
#endif #endif

195
vcg/space/point3.h
View File

@ -29,15 +29,14 @@
#define __VCGLIB_POINT3 #define __VCGLIB_POINT3
#include "../math/eigen.h" #include "../math/eigen.h"
#include <vcg/math/base.h>
namespace vcg{ namespace vcg{
template<class Scalar> class Point3; template<typename Scalar> class Point3;
} }
namespace Eigen{ namespace Eigen{
template<typename Scalar>
struct ei_traits<vcg::Point3<Scalar> > : ei_traits<Eigen::Matrix<Scalar,3,1> > {}; template<typename Scalar> struct ei_traits<vcg::Point3<Scalar> > : ei_traits<Eigen::Matrix<Scalar,3,1> > {};
template<typename Scalar> template<typename Scalar>
struct NumTraits<vcg::Point3<Scalar> > : NumTraits<Scalar> struct NumTraits<vcg::Point3<Scalar> > : NumTraits<Scalar>
@ -53,36 +52,33 @@ struct NumTraits<vcg::Point3<Scalar> > : NumTraits<Scalar>
namespace vcg { namespace vcg {
template<typename Scalar> class Box3;
/** \addtogroup space */ /** \addtogroup space */
/*@{*/ /*@{*/
/** /**
The templated class for representing a point in 3D 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 The class is templated over the ScalarType class that is used to represent coordinates. All the usual
operator overloading (* + - ...) is present. operator overloading (* + - ...) is present.
*/ */
template <class T> class Box3;
template <class _Scalar> class Point3 : public Eigen::Matrix<_Scalar,3,1> template <class _Scalar> class Point3 : public Eigen::Matrix<_Scalar,3,1>
{ {
typedef Eigen::Matrix<_Scalar,3,1> _Base; //----------------------------------------
using _Base::coeff; // template typedef part
using _Base::coeffRef; // use it as follow: typename Point3<S>::Type instead of simply Point3<S>
using _Base::setZero; //----------------------------------------
using _Base::data; public:
using _Base::V; typedef Eigen::Matrix<_Scalar,3,1> Type;
//----------------------------------------
// inheritence part
//----------------------------------------
private:
typedef Eigen::Matrix<_Scalar,3,1> _Base;
using _Base::Construct;
public: public:
_EIGEN_GENERIC_PUBLIC_INTERFACE(Point3,_Base); _EIGEN_GENERIC_PUBLIC_INTERFACE(Point3,_Base);
typedef Scalar ScalarType;
VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Point3) VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Point3)
enum {Dimension = 3};
//@{
/** @name Standard Constructors and Initializers /** @name Standard Constructors and Initializers
No casting operators have been introduced to avoid automatic unattended (and costly) conversion between different point types 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<OtherDerived>& other) : Base(other) {} inline Point3(const Eigen::MatrixBase<OtherDerived>& other) : Base(other) {}
template<class OtherDerived> // this one is very useless
inline void Import( const Eigen::MatrixBase<OtherDerived>& 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 <class Q>
static inline Point3 Construct( const Point3<Q> & b )
{
return Point3(Scalar(b[0]),Scalar(b[1]),Scalar(b[2]));
}
template <class Q> template <class Q>
static inline Point3 Construct( const Q & P0, const Q & P1, const Q & P2) static inline Point3 Construct( const Q & P0, const Q & P1, const Q & P2)
{ {
return Point3(Scalar(P0),Scalar(P1),Scalar(P2)); return Point3(Scalar(P0),Scalar(P1),Scalar(P2));
} }
vcg::Box3<_Scalar> GetBBox(vcg::Box3<_Scalar> &bb) const;
static inline Point3 Construct( const Point3<ScalarType> & b ) }; // end class definition (Point3)
{
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 <class Scalar>
inline Scalar AngleN( Point3<Scalar> const & p1, Point3<Scalar> const & p2 )
{
Scalar w = p1*p2;
if(w>1)
w = 1;
else if(w<-1)
w=-1;
return (Scalar) acos(w);
}
template <class Scalar>
inline Point3<Scalar> & Normalize( Point3<Scalar> & p )
{
p.Normalize();
return p;
}
// Dot product preciso numericamente (solo double!!) // Dot product preciso numericamente (solo double!!)
// Implementazione: si sommano i prodotti per ordine di esponente // Implementazione: si sommano i prodotti per ordine di esponente
@ -255,8 +133,6 @@ double stable_dot ( Point3<Scalar> const & p0, Point3<Scalar> const & p1 )
} }
} }
/// Point(p) Edge(v1-v2) dist, q is the point in v1-v2 with min dist /// Point(p) Edge(v1-v2) dist, q is the point in v1-v2 with min dist
template<class Scalar> template<class Scalar>
Scalar PSDist( const Point3<Scalar> & p, Scalar PSDist( const Point3<Scalar> & p,
@ -272,7 +148,6 @@ Scalar PSDist( const Point3<Scalar> & p,
return Distance(p,q); return Distance(p,q);
} }
template <class Scalar> template <class Scalar>
void GetUV( Point3<Scalar> &n,Point3<Scalar> &u, Point3<Scalar> &v, Point3<Scalar> up=(Point3<Scalar>(0,1,0)) ) void GetUV( Point3<Scalar> &n,Point3<Scalar> &u, Point3<Scalar> &v, Point3<Scalar> up=(Point3<Scalar>(0,1,0)) )
{ {
@ -297,23 +172,21 @@ void GetUV( Point3<Scalar> &n,Point3<Scalar> &u, Point3<Scalar> &v, Point3<Scala
Point3<Scalar> uv=u^v; Point3<Scalar> uv=u^v;
} }
/*@}*/
template <class SCALARTYPE>
inline Point3<SCALARTYPE> Abs(const Point3<SCALARTYPE> & p) {
return (Point3<SCALARTYPE>(math::Abs(p[0]), math::Abs(p[1]), math::Abs(p[2])));
}
// probably a more uniform naming should be defined...
template <class SCALARTYPE>
inline Point3<SCALARTYPE> LowClampToZero(const Point3<SCALARTYPE> & p) {
return (Point3<SCALARTYPE>(math::Max(p[0], (SCALARTYPE)0), math::Max(p[1], (SCALARTYPE)0), math::Max(p[2], (SCALARTYPE)0)));
}
typedef Point3<short> Point3s; typedef Point3<short> Point3s;
typedef Point3<int> Point3i; typedef Point3<int> Point3i;
typedef Point3<float> Point3f; typedef Point3<float> Point3f;
typedef Point3<double> Point3d; typedef Point3<double> Point3d;
/*@}*/ // typedef Eigen::Matrix<short ,3,1> Point3s;
// typedef Eigen::Matrix<int ,3,1> Point3i;
// typedef Eigen::Matrix<float ,3,1> Point3f;
// typedef Eigen::Matrix<double,3,1> Point3d;
// typedef Eigen::Matrix<short ,3,1> Vector3s;
// typedef Eigen::Matrix<int ,3,1> Vector3i;
// typedef Eigen::Matrix<float ,3,1> Vector3f;
// typedef Eigen::Matrix<double,3,1> Vector3d;
} // end namespace } // end namespace

77
vcg/space/point4.h
View File

@ -31,25 +31,35 @@
#include "../math/eigen.h" #include "../math/eigen.h"
namespace vcg{ namespace vcg{
template<class Scalar> class Point4; template<typename Scalar> class Point4;
} }
namespace Eigen{ namespace Eigen {
template<typename Scalar> template<typename Scalar> struct ei_traits<vcg::Point4<Scalar> > : ei_traits<Eigen::Matrix<Scalar,4,1> > {};
struct ei_traits<vcg::Point4<Scalar> > : ei_traits<Eigen::Matrix<Scalar,4,1> > {};
} }
namespace vcg { namespace vcg {
/** \addtogroup space */ /** \addtogroup space */
/*@{*/ /*@{*/
/** /**
The templated class for representing a point in 4D 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. The class is templated over the ScalarType class that is used to represent coordinates.
All the usual operator (* + - ...) are defined. All the usual operator (* + - ...) are defined.
*/ */
template <class T> class Point4 : public Eigen::Matrix<T,4,1> template <class T> class Point4 : public Eigen::Matrix<T,4,1>
{ {
//----------------------------------------
// template typedef part
// use it as follow: typename Point4<S>::Type instead of simply Point4<S>
//----------------------------------------
public:
typedef Eigen::Matrix<T,4,1> Type;
//----------------------------------------
// inheritence part
//----------------------------------------
private:
typedef Eigen::Matrix<T,4,1> _Base; typedef Eigen::Matrix<T,4,1> _Base;
using _Base::coeff; using _Base::coeff;
using _Base::coeffRef; using _Base::coeffRef;
@ -64,8 +74,6 @@ public:
VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Point4) VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Point4)
enum {Dimension = 4};
inline Point4() : Base() {} 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 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 T p[4]) : Base(p) {}
@ -74,30 +82,6 @@ public:
inline Point4(const Eigen::MatrixBase<OtherDerived>& other) : Base(other) {} inline Point4(const Eigen::MatrixBase<OtherDerived>& other) : Base(other) {}
/// importer from different Point4 types
template <class Q> inline void Import( const Point4<Q> & b ) { *this = b.template cast<T>(); }
/// constuctor that imports from different Point4 types
template <class Q>
static inline Point4 Construct( const Point4<Q> & b ) { return b.template cast<T>(); }
//@{
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 inline Point4 VectProd ( const Point4 &x, const Point4 &z ) const
{ {
Point4 res; Point4 res;
@ -147,18 +131,27 @@ public:
}; // end class definition }; // end class definition
/// slower version of dot product, more stable (double precision only) typedef Point4<short> Point4s;
typedef Point4<int> Point4i;
typedef Point4<float> Point4f;
typedef Point4<double> Point4d;
// typedef Eigen::Matrix<short ,4,1> Point4s;
// typedef Eigen::Matrix<int ,4,1> Point4i;
// typedef Eigen::Matrix<float ,4,1> Point4f;
// typedef Eigen::Matrix<double,4,1> Point4d;
// typedef Eigen::Matrix<short ,4,1> Vector4s;
// typedef Eigen::Matrix<int ,4,1> Vector4i;
// typedef Eigen::Matrix<float ,4,1> Vector4f;
// typedef Eigen::Matrix<double,4,1> Vector4d;
/// slower version of dot product, more stable (double precision only)
template<class T> template<class T>
double StableDot ( Point4<T> const & p0, Point4<T> const & p1 ) double StableDot ( Point4<T> const & p0, Point4<T> const & p1 )
{ {
return p0.StableDot(p1); return p0.StableDot(p1);
} }
typedef Point4<short> Point4s;
typedef Point4<int> Point4i;
typedef Point4<float> Point4f;
typedef Point4<double> Point4d;
/*@}*/ /*@}*/
} // end namespace } // end namespace
#endif #endif