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

Port to eigen2: state of the mess: 

* curently nothing change if you don't define VCG_USE_EIGEN
* make Matrix*, Point3 and Point4 derive Eigen::Matrix (still ugly)
* now catching all the dot products to replace them by .dot()
  note that most of meshlab already compile
This commit is contained in:
Paolo Cignoni 2008-10-27 19:35:17 +00:00
parent 69ebba04b9
commit ab200fc950
35 changed files with 1557 additions and 1205 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
vcg/complex/trimesh/crease_cut.h
View File

@ -91,7 +91,7 @@ void CreaseCut(MESH_TYPE &m, float angleRad)
if(!isBorderVertex) // for internal vertex we search the first crease and start from it
{
do {
ScalarType dotProd = iPos.FFlip()->cN() * iPos.f->N();
ScalarType dotProd = iPos.FFlip()->cN().dot(iPos.f->N());
iPos.NextFE();
if(dotProd<cosangle) break;
} while (startPos!=iPos);
@ -102,7 +102,7 @@ void CreaseCut(MESH_TYPE &m, float angleRad)
int curVertexCounter =vertInd;
do { // The real Loop
ScalarType dotProd=iPos.FFlip()->cN() * iPos.f->N(); // test normal with the next face (fflip)
ScalarType dotProd=iPos.FFlip()->cN().dot(iPos.f->N()); // test normal with the next face (fflip)
size_t faceInd = Index(m,iPos.f);
indVec[faceInd*3+ iPos.VInd()] = curVertexCounter;

20
vcg/complex/trimesh/create/ball_pivoting.h
View File

@ -147,7 +147,7 @@ template <class MESH> class BallPivoting: public AdvancingFront<MESH> {
if(d0 < min_edge || d0 > max_edge) continue;
Point3x normal = (p1 - p0)^(p2 - p0);
if(normal * (p0 - baricenter) < 0) continue;
if(normal.dot(p0 - baricenter) < 0) continue;
/* if(use_normals) {
if(normal * vv0->N() < 0) continue;
if(normal * vv1->N() < 0) continue;
@ -169,7 +169,7 @@ template <class MESH> class BallPivoting: public AdvancingFront<MESH> {
}
//check on the other side there is not a surface
Point3x opposite = center + normal*(((center - p0)*normal)*2/normal.SquaredNorm());
Point3x opposite = center + normal*(((center - p0).dot(normal))*2/normal.SquaredNorm());
for(t = 0; t < n; t++) {
VertexType &v = *(targets[t]);
if((v.IsV()) && (opposite - v.P()).Norm() <= radius)
@ -291,7 +291,7 @@ template <class MESH> class BallPivoting: public AdvancingFront<MESH> {
if(min_angle >= M_PI - 0.1) {
return -1;
}
if(candidate == NULL) {
return -1;
}
@ -304,7 +304,7 @@ template <class MESH> class BallPivoting: public AdvancingFront<MESH> {
assert(id != edge.v0 && id != edge.v1);
Point3x newnormal = ((candidate->P() - v0)^(v1 - v0)).Normalize();
if(normal * newnormal < max_angle || this->nb[id] >= 2) {
if(normal.dot(newnormal) < max_angle || this->nb[id] >= 2) {
return -1;
}
@ -315,7 +315,7 @@ template <class MESH> class BallPivoting: public AdvancingFront<MESH> {
if((*k).v0 == id) touch = k;
//mark vertices close to candidate
Mark(candidate);
Mark(candidate);
return id;
}
@ -360,9 +360,9 @@ template <class MESH> class BallPivoting: public AdvancingFront<MESH> {
up /= uplen;
ScalarType a11 = q1*q1;
ScalarType a12 = q1*q2;
ScalarType a22 = q2*q2;
ScalarType a11 = q1.dot(q1);
ScalarType a12 = q1.dot(q2);
ScalarType a22 = q2.dot(q2);
ScalarType m = 4*(a11*a22 - a12*a12);
ScalarType l1 = 2*(a11*a22 - a22*a12)/m;
@ -385,8 +385,8 @@ template <class MESH> class BallPivoting: public AdvancingFront<MESH> {
p.Normalize();
q.Normalize();
Point3x vec = p^q;
ScalarType angle = acos(p*q);
if(vec*axis < 0) angle = -angle;
ScalarType angle = acos(p.dot(q));
if(vec.dot(axis) < 0) angle = -angle;
if(angle < 0) angle += 2*M_PI;
return angle;
}

9
vcg/complex/trimesh/inertia.h
View File

@ -334,7 +334,7 @@ static void Covariance(const MeshType & m, vcg::Point3<ScalarType> & bary, vcg::
// integral of (x,y,0) in the same triangle
CoordType X(1/6.0,1/6.0,0);
vcg::Matrix33<ScalarType> A, // matrix that bring the vertices to (v1-v0,v2-v0,n)
At,DC;
DC;
for(fi = m.face.begin(); fi != m.face.end(); ++fi)
if(!(*fi).IsD())
{
@ -350,18 +350,15 @@ static void Covariance(const MeshType & m, vcg::Point3<ScalarType> & bary, vcg::
A.SetColumn(2,n);
CoordType delta = P0 - bary;
At= A;
At.Transpose();
/* DC is calculated as integral of (A*x+delta) * (A*x+delta)^T over the triangle,
where delta = v0-bary
*/
DC.SetZero();
DC+= A*C0*At;
DC+= A*C0*A.transpose();
vcg::Matrix33<ScalarType> tmp;
tmp.OuterProduct(A*X,delta);
DC+= tmp;
tmp.Transpose();
DC += tmp + tmp.transpose();
DC+= tmp;
tmp.OuterProduct(delta,delta);
DC+=tmp*0.5;

6
vcg/complex/trimesh/smooth.h
View File

@ -488,7 +488,7 @@ static void VertexCoordLaplacianQuality(MeshType &m, int step, bool SmoothSelect
/*
Improved Laplacian Smoothing of Noisy Surface Meshes
J. Vollmer, R. Mencl, and H. Müller
J. Vollmer, R. Mencl, and H. M<EFBFBD>ller
EUROGRAPHICS Volume 18 (1999), Number 3
*/
@ -1099,7 +1099,7 @@ static void FaceNormalAngleThreshold(MeshType &m,
{
if(! (*ep.f).IsV() )
{
ScalarType cosang=ep.f->N()*(*fi).N();
ScalarType cosang=ep.f->N().dot((*fi).N());
if(cosang >= sigma)
{
ScalarType w = cosang-sigma;
@ -1223,7 +1223,7 @@ static void FastFitMesh(MeshType &m,
for (;!ep.End();++ep)
{
CoordType bc=Barycenter<FaceType>(*ep.F());
Sum += ep.F()->N()*(ep.F()->N()*(bc - (*vi).P()));
Sum += ep.F()->N()*(ep.F()->N().dot(bc - (*vi).P()));
++cnt;
}
TDV[*vi].np=(*vi).P()+ Sum*(1.0/cnt);

7
vcg/complex/trimesh/update/curvature.h
View File

@ -211,11 +211,8 @@ public:
tempMatrix.ExternalProduct(W,W);
Q -= tempMatrix * 2.0f;
Matrix33<ScalarType> Qt(Q);
Qt.Transpose();
// compute matrix Q^t M Q
Matrix33<ScalarType> QtMQ = (Qt * M * Q);
Matrix33<ScalarType> QtMQ = (Q.transpose() * M * Q);
CoordType bl = Q.GetColumn(0);
CoordType T1 = Q.GetColumn(1);
@ -644,7 +641,7 @@ public:
int n_rot;
Jacobi(m33,lambda, vect,n_rot);
vect.Transpose();
vect = vect.transpose().eval();
ScalarType normal = std::numeric_limits<ScalarType>::min();
int normI = 0;
for(int i = 0; i < 3; ++i)

2
vcg/complex/trimesh/update/edges.h
View File

@ -74,7 +74,7 @@ namespace tri {
f.Edge(2) = f.V(0)->P(); f.Edge(2) -= f.V(2)->P();
// Calcolo di plane
f.Plane().SetDirection(f.Edge(0)^f.Edge(1));
f.Plane().SetOffset(f.Plane().Direction() * f.V(0)->P());
f.Plane().SetOffset(f.Plane().Direction().dot(f.V(0)->P()));
f.Plane().Normalize();
// Calcolo migliore proiezione
ScalarType nx = math::Abs(f.Plane().Direction()[0]);

3
vcg/math/deprecated_matrix.h
View File

@ -721,12 +721,15 @@ namespace vcg{
delete []temp;
};
// for the transistion to eigen
Matrix transpose()
{
Matrix res = *this;
res.Transpose();
return res;
}
// for the transistion to eigen
const Matrix& eval() { return *this; }
/*!

2
vcg/math/deprecated_matrix33.h
View File

@ -403,6 +403,8 @@ public:
res.Transpose();
return res;
}
// for the transistion to eigen
const Matrix33& eval() { return *this; }
/// Funzione per costruire una matrice diagonale dati i tre elem.
Matrix33 & SetDiagonal(S *v)

12
vcg/math/deprecated_matrix44.h
View File

@ -281,6 +281,16 @@ public:
return tmp;
}
// for the transistion to eigen
Matrix44 transpose() const
{
Matrix44 res = *this;
Transpose(res);
return res;
}
// for the transistion to eigen
const Matrix44& eval() { return *this; }
};
@ -289,7 +299,7 @@ public:
template <class T> class LinearSolve: public Matrix44<T> {
public:
LinearSolve(const Matrix44<T> &m);
Point4<T> Solve(const Point4<T> &b); // solve A · x = b
Point4<T> Solve(const Point4<T> &b); // solve A <EFBFBD> x = b
///If you need to solve some equation you can use this function instead of Matrix44 one for speed.
T Determinant() const;
protected:

74
vcg/math/eigen.h
View File

@ -26,8 +26,40 @@
#define EIGEN_MATRIXBASE_PLUGIN <vcg/math/eigen_vcgaddons.h>
#include "../Eigen/LU"
#include "../Eigen/Geometry"
#include "../Eigen/Array"
#include "../Eigen/Core"
#include "base.h"
namespace Eigen {
template<> struct NumTraits<unsigned char>
{
typedef unsigned char Real;
typedef float FloatingPoint;
enum {
IsComplex = 0,
HasFloatingPoint = 0,
ReadCost = 1,
AddCost = 1,
MulCost = 1
};
};
template<> struct NumTraits<short int>
{
typedef short int Real;
typedef float FloatingPoint;
enum {
IsComplex = 0,
HasFloatingPoint = 0,
ReadCost = 1,
AddCost = 1,
MulCost = 1
};
};
}
#define VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATOR(Derived, Op) \
template<typename OtherDerived> \
@ -54,4 +86,46 @@
VCG_EIGEN_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, *=) \
VCG_EIGEN_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, /=)
namespace vcg {
template<typename Derived1, typename Derived2>
typename Eigen::ei_traits<Derived1>::Scalar
Angle(const Eigen::MatrixBase<Derived1>& p1, const Eigen::MatrixBase<Derived2> & p2)
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived1)
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived2)
EIGEN_STATIC_ASSERT_FIXED_SIZE(Derived1)
EIGEN_STATIC_ASSERT_FIXED_SIZE(Derived2)
EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(Derived1,Derived2)
typedef typename Eigen::ei_traits<Derived1>::Scalar Scalar;
Scalar w = p1.norm()*p2.norm();
if(w==0) return Scalar(-1);
Scalar t = (p1.dot(p2))/w;
if(t>1) t = 1;
else if(t<-1) t = -1;
return vcg::math::Acos(t);
}
template<typename Derived1>
inline typename Eigen::ei_traits<Derived1>::Scalar Norm( const Eigen::MatrixBase<Derived1>& p)
{ return p.norm(); }
template<typename Derived1>
inline typename Eigen::ei_traits<Derived1>::Scalar SquaredNorm( const Eigen::MatrixBase<Derived1>& p)
{ return p.norm2(); }
template<typename Derived1, typename Derived2>
inline typename Eigen::ei_traits<Derived1>::Scalar
Distance(const Eigen::MatrixBase<Derived1>& p1, const Eigen::MatrixBase<Derived2> & p2)
{ return (p1-p2).norm(); }
template<typename Derived1, typename Derived2>
inline typename Eigen::ei_traits<Derived1>::Scalar
SquaredDistance(const Eigen::MatrixBase<Derived1>& p1, const Eigen::MatrixBase<Derived2> & p2)
{ return (p1-p2).norm2(); }
}
#endif

22
vcg/math/eigen_vcgaddons.h
View File

@ -72,10 +72,15 @@ EIGEN_DEPRECATED inline unsigned int RowsNumber() const
* \param j the column index
* \return the element
*/
EIGEN_DEPRECATED inline Scalar ElementAt(unsigned int i, unsigned int j)
EIGEN_DEPRECATED inline Scalar ElementAt(unsigned int i, unsigned int j) const
{
return (*this)(i,j);
};
}
EIGEN_DEPRECATED inline Scalar& ElementAt(unsigned int i, unsigned int j)
{
return (*this)(i,j);
}
/*!
* \deprecated use *this.determinant() (or *this.lu().determinant() for large matrices)
@ -148,6 +153,8 @@ EIGEN_DEPRECATED void SwapRows(const unsigned int i, const unsigned int j)
row(i).swap(row(j));
};
Scalar* V() { return derived().data(); }
const Scalar* V() const { return derived().data(); }
/*!
* \deprecated use *this.cwise() += k
@ -288,6 +295,8 @@ EIGEN_DEPRECATED void SetDiagonal(Scalar *v)
diagonal() = Map<Matrix<Scalar,RowsAtCompileTime,1> >(v,cols(),1);
}
/** \deprecated use trace() */
EIGEN_DEPRECATED Scalar Trace() const { return trace(); }
/*!
* \deprecated use *this = *this.transpose()
@ -316,4 +325,13 @@ EIGEN_DEPRECATED void Dump()
printf("\n");
}
/** \deprecated use norm() */
EIGEN_DEPRECATED inline Scalar Norm() const { return norm(); }
/** \deprecated use squaredNorm() */
EIGEN_DEPRECATED inline Scalar SquaredNorm() const { return norm2(); }
/** \deprecated use normalize() or normalized() */
EIGEN_DEPRECATED inline Derived& Normalize() { normalize(); return derived(); }
/** \deprecated use .cross(p) */
inline EvalType operator ^ (const Derived& p ) const { return this->cross(p); }

4
vcg/math/gen_normal.h
View File

@ -73,7 +73,7 @@ static void UniformCone(int vn, std::vector<Point3<ScalarType > > &NN, ScalarTyp
ScalarType DotProd = cos(AngleRad);
for(vi=NNT.begin();vi!=NNT.end();++vi)
{
if(dir*(*vi) >= DotProd) NN.push_back(*vi);
if(dir.dot(*vi) >= DotProd) NN.push_back(*vi);
}
}
@ -124,7 +124,7 @@ static int BestMatchingNormal(const Point3x &n, std::vector<Point3x> &nv)
typename std::vector<Point3x>::iterator ni;
for(ni=nv.begin();ni!=nv.end();++ni)
{
cosang=(*ni)*n;
cosang=(*ni).dot(n);
if(cosang>bestang) {
bestang=cosang;
ret=ni-nv.begin();

18
vcg/math/matrix.h
View File

@ -23,8 +23,8 @@
#ifndef VCG_USE_EIGEN
#include "deprecated_matrix.h"
#endif
#else
#ifndef MATRIX_VCGLIB
#define MATRIX_VCGLIB
@ -69,10 +69,9 @@ public:
* @param Scalar (Templete Parameter) Specifies the ScalarType field.
*/
template<class _Scalar>
class Matrix : public Eigen::Matrix<_Scalar,Eigen::Dynamic,Eigen::Dynamic> // FIXME col or row major ?
class Matrix : public Eigen::Matrix<_Scalar,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> // FIXME col or row major ?
{
typedef Eigen::Matrix<_Scalar,Eigen::Dynamic,Eigen::Dynamic> _Base;
typedef Eigen::Matrix<_Scalar,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> _Base;
public:
@ -87,7 +86,7 @@ public:
* \param m the number of matrix rows
* \param n the number of matrix columns
*/
Matrix(unsigned int m, unsigned int n)
Matrix(int m, int n)
: Base(m,n)
{
memset(Base::data(), 0, m*n*sizeof(Scalar));
@ -100,10 +99,10 @@ public:
* \param n the number of matrix columns
* \param values the values of the matrix elements
*/
Matrix(unsigned int m, unsigned int n, Scalar *values)
: Base(m.n)
Matrix(int m, int n, Scalar *values)
: Base(m,n)
{
*this = Eigen::Map<Eigen::Matrix<Scalar,Dynamic,Dynamic,RowMajor> >(values, m , n);
*this = Eigen::Map<Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> >(values, m , n);
}
/*!
@ -197,3 +196,6 @@ void Invert(MatrixType & m)
} // end of namespace
#endif
#endif

75
vcg/math/matrix33.h
View File

@ -21,14 +21,15 @@
* *
****************************************************************************/
// #ifndef VCG_USE_EIGEN
#ifndef VCG_USE_EIGEN
#include "deprecated_matrix33.h"
// #endif
#else
#ifndef __VCGLIB_MATRIX33_H
#define __VCGLIB_MATRIX33_H
#include "eigen.h"
#include "matrix44.h"
namespace vcg{
template<class Scalar> class Matrix33;
@ -48,17 +49,20 @@ namespace vcg {
@param S (Templete Parameter) Specifies the ScalarType field.
*/
template<class _Scalar>
class Matrix33 : public Eigen::Matrix<_Scalar,3,3,RowMajor> // FIXME col or row major ?
class Matrix33 : public Eigen::Matrix<_Scalar,3,3,Eigen::RowMajor> // FIXME col or row major ?
{
typedef Eigen::Matrix<_Scalar,Eigen::Dynamic,Eigen::Dynamic> _Base;
typedef Eigen::Matrix<_Scalar,3,3,Eigen::RowMajor> _Base;
using _Base::coeff;
using _Base::coeffRef;
using _Base::setZero;
public:
_EIGEN_GENERIC_PUBLIC_INTERFACE(Matrix,_Base);
_EIGEN_GENERIC_PUBLIC_INTERFACE(Matrix33,_Base);
typedef _Scalar ScalarType;
VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Matrix)
VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Matrix33)
/// Default constructor
inline Matrix33() : Base() {}
@ -67,10 +71,13 @@ public:
Matrix33(const Matrix33& m ) : Base(m) {}
/// create from a \b row-major array
Matrix33(const Scalar * v ) : Base(Map<Matrix<Scalar,3,3,RowMajor>(v)) {}
Matrix33(const Scalar * v ) : Base(Eigen::Map<Eigen::Matrix<Scalar,3,3,Eigen::RowMajor> >(v)) {}
/// create from Matrix44 excluding row and column k
Matrix33(const Matrix44<S> & m, const int & k) : Base(m.minor(k,k)) {}
Matrix33(const Matrix44<Scalar> & m, const int & k) : Base(m.minor(k,k)) {}
template<typename OtherDerived>
Matrix33(const Eigen::MatrixBase<OtherDerived>& other) : Base(other) {}
/*!
* \deprecated use *this.row(i)
@ -86,13 +93,13 @@ public:
/** \deprecated */
Matrix33 & SetRotateRad(S angle, const Point3<S> & axis )
Matrix33 & SetRotateRad(Scalar angle, const Point3<Scalar> & axis )
{
*this = Eigen::AngleAxis<Scalar>(angle,axis).toRotationMatrix();
return (*this);
}
/** \deprecated */
Matrix33 & SetRotateDeg(S angle, const Point3<S> & axis ){
Matrix33 & SetRotateDeg(Scalar angle, const Point3<Scalar> & axis ){
return SetRotateRad(math::ToRad(angle),axis);
}
@ -100,7 +107,7 @@ public:
// Warning, this Inversion code can be HIGHLY NUMERICALLY UNSTABLE!
// In most case you are advised to use the Invert() method based on SVD decomposition.
/** \deprecated */
Matrix33 & FastInvert() { return *this = inverse(); }
Matrix33 & FastInvert() { return *this = Base::inverse(); }
void show(FILE * fp)
{
@ -111,7 +118,7 @@ public:
/*
compute the matrix generated by the product of a * b^T
*/
void ExternalProduct(const Point3<S> &a, const Point3<S> &b)
void ExternalProduct(const Point3<Scalar> &a, const Point3<Scalar> &b)
{
for(int i=0;i<3;++i)
for(int j=0;j<3;++j)
@ -120,15 +127,15 @@ public:
/* Compute the Frobenius Norm of the Matrix
*/
ScalarType Norm()
{
// FIXME looks like there is a bug: j is not used !!!
ScalarType SQsum=0;
for(int i=0;i<3;++i)
for(int j=0;j<3;++j)
SQsum += a[i]*a[i];
return (math::Sqrt(SQsum));
}
Scalar Norm() { return Base::cwise().abs2().sum(); }
// {
// // FIXME looks like there is a bug: j is not used !!!
// Scalar SQsum=0;
// for(int i=0;i<3;++i)
// for(int j=0;j<3;++j)
// SQsum += a[i]*a[i];
// return (math::Sqrt(SQsum));
// }
/**
@ -136,7 +143,7 @@ public:
*/
// FIXME should be outside Matrix
template <class STLPOINTCONTAINER >
void Covariance(const STLPOINTCONTAINER &points, Point3<S> &bp) {
void Covariance(const STLPOINTCONTAINER &points, Point3<Scalar> &bp) {
assert(!points.empty());
typedef typename STLPOINTCONTAINER::const_iterator PointIte;
// first cycle: compute the barycenter
@ -147,14 +154,12 @@ public:
this->setZero();
vcg::Matrix33<ScalarType> A;
for( PointIte pi = points.begin(); pi != points.end(); ++pi) {
Point3<S> p = (*pi)-bp;
Point3<Scalar> p = (*pi)-bp;
A.OuterProduct(p,p);
(*this)+= A;
}
}
/**
It computes the cross covariance matrix of two set of 3d points P and X;
it returns also the barycenters of P and X.
@ -168,14 +173,14 @@ public:
// FIXME should be outside Matrix
template <class STLPOINTCONTAINER >
void CrossCovariance(const STLPOINTCONTAINER &P, const STLPOINTCONTAINER &X,
Point3<S> &bp, Point3<S> &bx)
Point3<Scalar> &bp, Point3<Scalar> &bx)
{
setZero();
assert(P.size()==X.size());
bx.setZero();
bp.setZero();
Matrix33<S> tmp;
typename std::vector <Point3<S> >::const_iterator pi,xi;
Matrix33<Scalar> tmp;
typename std::vector <Point3<Scalar> >::const_iterator pi,xi;
for(pi=P.begin(),xi=X.begin();pi!=P.end();++pi,++xi){
bp+=*pi;
bx+=*xi;
@ -193,15 +198,15 @@ public:
void WeightedCrossCovariance(const STLREALCONTAINER & weights,
const STLPOINTCONTAINER &P,
const STLPOINTCONTAINER &X,
Point3<S> &bp,
Point3<S> &bx)
Point3<Scalar> &bp,
Point3<Scalar> &bx)
{
SetZero();
setZero();
assert(P.size()==X.size());
bx.SetZero();
bp.SetZero();
Matrix33<S> tmp;
typename std::vector <Point3<S> >::const_iterator pi,xi;
Matrix33<Scalar> tmp;
typename std::vector <Point3<Scalar> >::const_iterator pi,xi;
typename STLREALCONTAINER::const_iterator pw;
for(pi=P.begin(),xi=X.begin();pi!=P.end();++pi,++xi){
@ -317,10 +322,12 @@ template <class S>
///
typedef Matrix33<short> Matrix33s;
typedef Matrix33<int> Matrix33i;
typedef Matrix33<int> Matrix33i;
typedef Matrix33<float> Matrix33f;
typedef Matrix33<double> Matrix33d;
} // end of namespace
#endif
#endif

491
vcg/math/matrix44.h
View File

@ -8,7 +8,7 @@
* \ *
* All rights reserved. *
* *
* This program is free software; you can redistribute it and/or modify *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
@ -21,9 +21,9 @@
* *
****************************************************************************/
// #ifndef VCG_USE_EIGEN
#ifndef VCG_USE_EIGEN
#include "deprecated_matrix44.h"
// #endif
#else
#ifndef __VCGLIB_MATRIX44
#define __VCGLIB_MATRIX44
@ -83,17 +83,21 @@ for 'column' vectors.
/** This class represents a 4x4 matrix. T is the kind of element in the matrix.
*/
template<typename _Scalar>
class Matrix44 : public Eigen::Matrix<_Scalar,4,4,RowMajor> // FIXME col or row major !
class Matrix44 : public Eigen::Matrix<_Scalar,4,4,Eigen::RowMajor> // FIXME col or row major !
{
typedef Eigen::Matrix<_Scalar,Eigen::Dynamic,Eigen::Dynamic> _Base;
typedef Eigen::Matrix<_Scalar,4,4,Eigen::RowMajor> _Base;
using _Base::coeff;
using _Base::coeffRef;
using _Base::ElementAt;
using _Base::setZero;
public:
_EIGEN_GENERIC_PUBLIC_INTERFACE(Matrix,_Base);
_EIGEN_GENERIC_PUBLIC_INTERFACE(Matrix44,_Base);
typedef _Scalar ScalarType;
VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Matrix)
VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Matrix44)
/** \name Constructors
* No automatic casting and default constructor is empty
@ -101,50 +105,51 @@ public:
Matrix44() : Base() {}
~Matrix44() {}
Matrix44(const Matrix44 &m) : Base(m) {}
Matrix33(const Scalar * v ) : Base(Map<Matrix<Scalar,4,4,RowMajor>(v)) {}
Matrix44(const Scalar * v ) : Base(Eigen::Map<Eigen::Matrix<Scalar,4,4,Eigen::RowMajor> >(v)) {}
Scalar *V() { return Base::data(); }
const Scalar *V() const { return Base::data(); }
template<typename OtherDerived>
Matrix44(const Eigen::MatrixBase<OtherDerived>& other) : Base(other) {}
const typename Base::RowXpr operator[](int i) const { return Base::row(i); }
typename Base::RowXpr operator[](int i) { return Base::row(i); }
// return a copy of the i-th column
typename Base::ColXpr GetColumn4(const int& i) const { return col(i); }
Block<Matrix,3,1> GetColumn3(const int& i) const { return block<3,1>(0,i); }
typename Base::ColXpr GetColumn4(const int& i) const { return Base::col(i); }
const Eigen::Block<Base,3,1> GetColumn3(const int& i) const { return this->template block<3,1>(0,i); }
typename Base::RowXpr GetRow4(const int& i) const { return col(i); }
Block<Matrix,1,3> GetRow3(const int& i) const { return block<1,3>(i,0); }
typename Base::RowXpr GetRow4(const int& i) const { return Base::row(i); }
Eigen::Block<Base,1,3> GetRow3(const int& i) const { return this->template block<1,3>(i,0); }
template <class Matrix44Type>
void ToMatrix(Matrix44Type & m) const { m = (*this).cast<typename Matrix44Type::Scalar>(); }
void ToMatrix(Matrix44Type & m) const { m = (*this).template cast<typename Matrix44Type::Scalar>(); }
void ToEulerAngles(Scalar &alpha, Scalar &beta, Scalar &gamma);
template <class Matrix44Type>
void FromMatrix(const Matrix44Type & m) { for(int i = 0; i < 16; i++) data()[i] = m.data()[i]; }
void FromMatrix(const Matrix44Type & m) { for(int i = 0; i < 16; i++) Base::data()[i] = m.data()[i]; }
void FromEulerAngles(T alpha, T beta, T gamma);
void SetDiagonal(const T k);
Matrix44 &SetScale(const T sx, const T sy, const T sz);
Matrix44 &SetScale(const Point3<T> &t);
Matrix44 &SetTranslate(const Point3<T> &t);
Matrix44 &SetTranslate(const T sx, const T sy, const T sz);
Matrix44 &SetShearXY(const T sz);
Matrix44 &SetShearXZ(const T sy);
Matrix44 &SetShearYZ(const T sx);
void FromEulerAngles(Scalar alpha, Scalar beta, Scalar gamma);
void SetDiagonal(const Scalar k);
Matrix44 &SetScale(const Scalar sx, const Scalar sy, const Scalar sz);
Matrix44 &SetScale(const Point3<Scalar> &t);
Matrix44 &SetTranslate(const Point3<Scalar> &t);
Matrix44 &SetTranslate(const Scalar sx, const Scalar sy, const Scalar sz);
Matrix44 &SetShearXY(const Scalar sz);
Matrix44 &SetShearXZ(const Scalar sy);
Matrix44 &SetShearYZ(const Scalar sx);
///use radiants for angle.
Matrix44 &SetRotateDeg(T AngleDeg, const Point3<T> & axis);
Matrix44 &SetRotateRad(T AngleRad, const Point3<T> & axis);
T Determinant() const;
Matrix44 &SetRotateDeg(Scalar AngleDeg, const Point3<Scalar> & axis);
Matrix44 &SetRotateRad(Scalar AngleRad, const Point3<Scalar> & axis);
template <class Q> void Import(const Matrix44<Q> &m) {
for(int i = 0; i < 16; i++)
_a[i] = (T)(m.V()[i]);
Base::data()[i] = (Scalar)(m.data()[i]);
}
template <class Q>
static inline Matrix44 Construct( const Matrix44<Q> & b )
{
Matrix44<T> tmp; tmp.FromMatrix(b);
Matrix44 tmp; tmp.FromMatrix(b);
return tmp;
}
@ -173,7 +178,6 @@ protected:
///Premultiply
template <class T> Point3<T> operator*(const Matrix44<T> &m, const Point3<T> &p);
template <class T> Matrix44<T> &Transpose(Matrix44<T> &m);
//return NULL matrix if not invertible
template <class T> Matrix44<T> &Invert(Matrix44<T> &m);
template <class T> Matrix44<T> Inverse(const Matrix44<T> &m);
@ -184,27 +188,6 @@ typedef Matrix44<float> Matrix44f;
typedef Matrix44<double> Matrix44d;
template <class T> Matrix44<T>::Matrix44(const Matrix44<T> &m) {
memcpy((T *)_a, (T *)m._a, 16 * sizeof(T));
}
template <class T> Matrix44<T>::Matrix44(const T v[]) {
memcpy((T *)_a, v, 16 * sizeof(T));
}
template <class T> T &Matrix44<T>::ElementAt(const int row, const int col) {
assert(row >= 0 && row < 4);
assert(col >= 0 && col < 4);
return _a[(row<<2) + col];
}
template <class T> T Matrix44<T>::ElementAt(const int row, const int col) const {
assert(row >= 0 && row < 4);
assert(col >= 0 && col < 4);
return _a[(row<<2) + col];
}
//template <class T> T &Matrix44<T>::operator[](const int i) {
// assert(i >= 0 && i < 16);
// return ((T *)_a)[i];
@ -214,116 +197,7 @@ template <class T> T Matrix44<T>::ElementAt(const int row, const int col) const
// assert(i >= 0 && i < 16);
// return ((T *)_a)[i];
//}
template <class T> T *Matrix44<T>::operator[](const int i) {
assert(i >= 0 && i < 4);
return _a+i*4;
}
template <class T> const T *Matrix44<T>::operator[](const int i) const {
assert(i >= 0 && i < 4);
return _a+i*4;
}
template <class T> T *Matrix44<T>::V() { return _a;}
template <class T> const T *Matrix44<T>::V() const { return _a;}
template <class T> Matrix44<T> Matrix44<T>::operator+(const Matrix44 &m) const {
Matrix44<T> ret;
for(int i = 0; i < 16; i++)
ret.V()[i] = V()[i] + m.V()[i];
return ret;
}
template <class T> Matrix44<T> Matrix44<T>::operator-(const Matrix44 &m) const {
Matrix44<T> ret;
for(int i = 0; i < 16; i++)
ret.V()[i] = V()[i] - m.V()[i];
return ret;
}
template <class T> Matrix44<T> Matrix44<T>::operator*(const Matrix44 &m) const {
Matrix44 ret;
for(int i = 0; i < 4; i++)
for(int j = 0; j < 4; j++) {
T t = 0.0;
for(int k = 0; k < 4; k++)
t += ElementAt(i, k) * m.ElementAt(k, j);
ret.ElementAt(i, j) = t;
}
return ret;
}
template <class T> Matrix44<T> Matrix44<T>::operator*(const Matrix44Diag<T> &m) const {
Matrix44 ret = (*this);
for(int i = 0; i < 4; ++i)
for(int j = 0; j < 4; ++j)
ret[i][j]*=m[i];
return ret;
}
template <class T> Point4<T> Matrix44<T>::operator*(const Point4<T> &v) const {
Point4<T> ret;
for(int i = 0; i < 4; i++){
T t = 0.0;
for(int k = 0; k < 4; k++)
t += ElementAt(i,k) * v[k];
ret[i] = t;
}
return ret;
}
template <class T> bool Matrix44<T>::operator==(const Matrix44 &m) const {
for(int i = 0; i < 4; ++i)
for(int j = 0; j < 4; ++j)
if(ElementAt(i,j) != m.ElementAt(i,j))
return false;
return true;
}
template <class T> bool Matrix44<T>::operator!=(const Matrix44 &m) const {
for(int i = 0; i < 4; ++i)
for(int j = 0; j < 4; ++j)
if(ElementAt(i,j) != m.ElementAt(i,j))
return true;
return false;
}
template <class T> Matrix44<T> Matrix44<T>::operator-() const {
Matrix44<T> res;
for(int i = 0; i < 16; i++)
res.V()[i] = -V()[i];
return res;
}
template <class T> Matrix44<T> Matrix44<T>::operator*(const T k) const {
Matrix44<T> res;
for(int i = 0; i < 16; i++)
res.V()[i] =V()[i] * k;
return res;
}
template <class T> void Matrix44<T>::operator+=(const Matrix44 &m) {
for(int i = 0; i < 16; i++)
V()[i] += m.V()[i];
}
template <class T> void Matrix44<T>::operator-=(const Matrix44 &m) {
for(int i = 0; i < 16; i++)
V()[i] -= m.V()[i];
}
template <class T> void Matrix44<T>::operator*=( const Matrix44 & m ) {
*this = *this *m;
/*for(int i = 0; i < 4; i++) { //sbagliato
Point4<T> t(0, 0, 0, 0);
for(int k = 0; k < 4; k++) {
for(int j = 0; j < 4; j++) {
t[k] += ElementAt(i, k) * m.ElementAt(k, j);
}
}
for(int l = 0; l < 4; l++)
ElementAt(i, l) = t[l];
} */
}
template < class PointType , class T > void operator*=( std::vector<PointType> &vert, const Matrix44<T> & m ) {
typename std::vector<PointType>::iterator ii;
@ -331,21 +205,16 @@ template < class PointType , class T > void operator*=( std::vector<PointType> &
(*ii).P()=m * (*ii).P();
}
template <class T> void Matrix44<T>::operator*=( const T k ) {
for(int i = 0; i < 16; i++)
_a[i] *= k;
}
template <class T>
void Matrix44<T>::ToEulerAngles(T &alpha, T &beta, T &gamma)
void Matrix44<T>::ToEulerAngles(Scalar &alpha, Scalar &beta, Scalar &gamma)
{
alpha = atan2(ElementAt(1,2), ElementAt(2,2));
beta = asin(-ElementAt(0,2));
gamma = atan2(ElementAt(0,1), ElementAt(1,1));
alpha = atan2(coeff(1,2), coeff(2,2));
beta = asin(-coeff(0,2));
gamma = atan2(coeff(0,1), coeff(1,1));
}
template <class T>
void Matrix44<T>::FromEulerAngles(T alpha, T beta, T gamma)
void Matrix44<T>::FromEulerAngles(Scalar alpha, Scalar beta, Scalar gamma)
{
this->SetZero();
@ -371,28 +240,20 @@ void Matrix44<T>::FromEulerAngles(T alpha, T beta, T gamma)
ElementAt(3,3) = 1;
}
template <class T> void Matrix44<T>::SetZero() {
memset((T *)_a, 0, 16 * sizeof(T));
}
template <class T> void Matrix44<T>::SetIdentity() {
SetDiagonal(1);
}
template <class T> void Matrix44<T>::SetDiagonal(const T k) {
SetZero();
template <class T> void Matrix44<T>::SetDiagonal(const Scalar k) {
setZero();
ElementAt(0, 0) = k;
ElementAt(1, 1) = k;
ElementAt(2, 2) = k;
ElementAt(3, 3) = 1;
ElementAt(3, 3) = 1;
}
template <class T> Matrix44<T> &Matrix44<T>::SetScale(const Point3<T> &t) {
template <class T> Matrix44<T> &Matrix44<T>::SetScale(const Point3<Scalar> &t) {
SetScale(t[0], t[1], t[2]);
return *this;
}
template <class T> Matrix44<T> &Matrix44<T>::SetScale(const T sx, const T sy, const T sz) {
SetZero();
template <class T> Matrix44<T> &Matrix44<T>::SetScale(const Scalar sx, const Scalar sy, const Scalar sz) {
setZero();
ElementAt(0, 0) = sx;
ElementAt(1, 1) = sy;
ElementAt(2, 2) = sz;
@ -400,23 +261,23 @@ template <class T> Matrix44<T> &Matrix44<T>::SetScale(const T sx, const T sy, co
return *this;
}
template <class T> Matrix44<T> &Matrix44<T>::SetTranslate(const Point3<T> &t) {
template <class T> Matrix44<T> &Matrix44<T>::SetTranslate(const Point3<Scalar> &t) {
SetTranslate(t[0], t[1], t[2]);
return *this;
}
template <class T> Matrix44<T> &Matrix44<T>::SetTranslate(const T tx, const T ty, const T tz) {
SetIdentity();
template <class T> Matrix44<T> &Matrix44<T>::SetTranslate(const Scalar tx, const Scalar ty, const Scalar tz) {
Base::setIdentity();
ElementAt(0, 3) = tx;
ElementAt(1, 3) = ty;
ElementAt(2, 3) = tz;
return *this;
}
template <class T> Matrix44<T> &Matrix44<T>::SetRotateDeg(T AngleDeg, const Point3<T> & axis) {
template <class T> Matrix44<T> &Matrix44<T>::SetRotateDeg(Scalar AngleDeg, const Point3<Scalar> & axis) {
return SetRotateRad(math::ToRad(AngleDeg),axis);
}
template <class T> Matrix44<T> &Matrix44<T>::SetRotateRad(T AngleRad, const Point3<T> & axis) {
template <class T> Matrix44<T> &Matrix44<T>::SetRotateRad(Scalar AngleRad, const Point3<Scalar> & axis) {
//angle = angle*(T)3.14159265358979323846/180; e' in radianti!
T c = math::Cos(AngleRad);
T s = math::Sin(AngleRad);
@ -436,9 +297,9 @@ template <class T> Matrix44<T> &Matrix44<T>::SetRotateRad(T AngleRad, const Poin
ElementAt(2,2) = t[2]*t[2]*q +c;
ElementAt(2,3) = 0;
ElementAt(3,0) = 0;
ElementAt(3,1) = 0;
ElementAt(3,1) = 0;
ElementAt(3,2) = 0;
ElementAt(3,3) = 1;
ElementAt(3,3) = 1;
return *this;
}
@ -461,20 +322,20 @@ template <class T> Matrix44<T> &Matrix44<T>::SetRotateRad(T AngleRad, const Poin
*/
template <class T> Matrix44<T> & Matrix44<T>:: SetShearXY( const T sh) {// shear the X coordinate as the Y coordinate change
SetIdentity();
template <class T> Matrix44<T> & Matrix44<T>::SetShearXY( const Scalar sh) {// shear the X coordinate as the Y coordinate change
Base::setIdentity();
ElementAt(0,1) = sh;
return *this;
}
template <class T> Matrix44<T> & Matrix44<T>:: SetShearXZ( const T sh) {// shear the X coordinate as the Z coordinate change
SetIdentity();
template <class T> Matrix44<T> & Matrix44<T>::SetShearXZ( const Scalar sh) {// shear the X coordinate as the Z coordinate change
Base::setIdentity();
ElementAt(0,2) = sh;
return *this;
}
template <class T> Matrix44<T> &Matrix44<T>:: SetShearYZ( const T sh) {// shear the Y coordinate as the Z coordinate change
SetIdentity();
template <class T> Matrix44<T> &Matrix44<T>::SetShearYZ( const Scalar sh) {// shear the Y coordinate as the Z coordinate change
Base::setIdentity();
ElementAt(1,2) = sh;
return *this;
}
@ -546,28 +407,28 @@ bool Decompose(Matrix44<T> &M, Point3<T> &ScaleV, Point3<T> &ShearV, Point3<T> &
R[0]=M.GetColumn3(0);
R[0].Normalize();
ShearV[0]=R[0]*M.GetColumn3(1); // xy shearing
ShearV[0]=R[0].dot(M.GetColumn3(1)); // xy shearing
R[1]= M.GetColumn3(1)-R[0]*ShearV[0];
assert(math::Abs(R[1]*R[0])<1e-10);
assert(math::Abs(R[1].dot(R[0]))<1e-10);
ScaleV[1]=Norm(R[1]); // y scaling
R[1]=R[1]/ScaleV[1];
ShearV[0]=ShearV[0]/ScaleV[1];
ShearV[1]=R[0]*M.GetColumn3(2); // xz shearing
ShearV[1]=R[0].dot(M.GetColumn3(2)); // xz shearing
R[2]= M.GetColumn3(2)-R[0]*ShearV[1];
assert(math::Abs(R[2]*R[0])<1e-10);
assert(math::Abs(R[2].dot(R[0]))<1e-10);
R[2] = R[2]-R[1]*(R[2]*R[1]);
assert(math::Abs(R[2]*R[1])<1e-10);
assert(math::Abs(R[2]*R[0])<1e-10);
R[2] = R[2]-R[1]*(R[2].dot(R[1]));
assert(math::Abs(R[2].dot(R[1]))<1e-10);
assert(math::Abs(R[2].dot(R[0]))<1e-10);
ScaleV[2]=Norm(R[2]);
ShearV[1]=ShearV[1]/ScaleV[2];
R[2]=R[2]/ScaleV[2];
assert(math::Abs(R[2]*R[1])<1e-10);
assert(math::Abs(R[2]*R[0])<1e-10);
assert(math::Abs(R[2].dot(R[1]))<1e-10);
assert(math::Abs(R[2].dot(R[0]))<1e-10);
ShearV[2]=R[1]*M.GetColumn3(2); // yz shearing
ShearV[2]=R[1].dot(M.GetColumn3(2)); // yz shearing
ShearV[2]=ShearV[2]/ScaleV[2];
int i,j;
for(i=0;i<3;++i)
@ -608,15 +469,6 @@ bool Decompose(Matrix44<T> &M, Point3<T> &ScaleV, Point3<T> &ShearV, Point3<T> &
return true;
}
template <class T> T Matrix44<T>::Determinant() const {
LinearSolve<T> solve(*this);
return solve.Determinant();
}
template <class T> Point3<T> operator*(const Matrix44<T> &m, const Point3<T> &p) {
T w;
Point3<T> s;
@ -628,24 +480,19 @@ template <class T> Point3<T> operator*(const Matrix44<T> &m, const Point3<T> &p)
return s;
}
//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) {
// T w;
// Point3<T> s;
// s[0] = m.ElementAt(0, 0)*p[0] + m.ElementAt(1, 0)*p[1] + m.ElementAt(2, 0)*p[2] + m.ElementAt(3, 0);
// s[1] = m.ElementAt(0, 1)*p[0] + m.ElementAt(1, 1)*p[1] + m.ElementAt(2, 1)*p[2] + m.ElementAt(3, 1);
// s[2] = m.ElementAt(0, 2)*p[0] + m.ElementAt(1, 2)*p[1] + m.ElementAt(2, 2)*p[2] + m.ElementAt(3, 2);
// w = m.ElementAt(0, 3)*p[0] + m.ElementAt(1, 3)*p[1] + m.ElementAt(2, 3)*p[2] + m.ElementAt(3, 3);
// if(w != 0) s /= w;
// if(w != 0) s /= w;
// return s;
//}
// }
template <class T> Matrix44<T> &Transpose(Matrix44<T> &m) {
for(int i = 1; i < 4; i++)
for(int j = 0; j < i; j++) {
math::Swap(m.ElementAt(i, j), m.ElementAt(j, i));
}
return m;
}
/*
To invert a matrix you can
@ -658,194 +505,16 @@ template <class T> Matrix44<T> &Transpose(Matrix44<T> &m) {
invertedMatrix = vcg::Inverse(untouchedMatrix);
*/
template <class T> Matrix44<T> & Invert(Matrix44<T> &m) {
LinearSolve<T> solve(m);
for(int j = 0; j < 4; j++) { //Find inverse by columns.
Point4<T> col(0, 0, 0, 0);
col[j] = 1.0;
col = solve.Solve(col);
for(int i = 0; i < 4; i++)
m.ElementAt(i, j) = col[i];
}
return m;
template <class T> Matrix44<T> & Invert(Matrix44<T> &m) {
return m = m.lu().inverse();
}
template <class T> Matrix44<T> Inverse(const Matrix44<T> &m) {
LinearSolve<T> solve(m);
Matrix44<T> res;
for(int j = 0; j < 4; j++) { //Find inverse by columns.
Point4<T> col(0, 0, 0, 0);
col[j] = 1.0;
col = solve.Solve(col);
for(int i = 0; i < 4; i++)
res.ElementAt(i, j) = col[i];
}
return res;
template <class T> Matrix44<T> Inverse(const Matrix44<T> &m) {
return m.lu().inverse();
}
/* LINEAR SOLVE IMPLEMENTATION */
template <class T> LinearSolve<T>::LinearSolve(const Matrix44<T> &m): Matrix44<T>(m) {
if(!Decompose()) {
for(int i = 0; i < 4; i++)
index[i] = i;
Matrix44<T>::SetZero();
}
}
template <class T> T LinearSolve<T>::Determinant() const {
T det = d;
for(int j = 0; j < 4; j++)
det *= this-> ElementAt(j, j);
return det;
}
/*replaces a matrix by its LU decomposition of a rowwise permutation.
d is +or -1 depeneing of row permutation even or odd.*/
#define TINY 1e-100
template <class T> bool LinearSolve<T>::Decompose() {
/* Matrix44<T> A;
for(int i = 0; i < 16; i++)
A[i] = operator[](i);
SetIdentity();
Point4<T> scale;
// Set scale factor, scale(i) = max( |a(i,j)| ), for each row
for(int i = 0; i < 4; i++ ) {
index[i] = i; // Initialize row index list
T scalemax = (T)0.0;
for(int j = 0; j < 4; j++)
scalemax = (scalemax > math::Abs(A.ElementAt(i,j))) ? scalemax : math::Abs(A.ElementAt(i,j));
scale[i] = scalemax;
}
// Loop over rows k = 1, ..., (N-1)
int signDet = 1;
for(int k = 0; k < 3; k++ ) {
// Select pivot row from max( |a(j,k)/s(j)| )
T ratiomax = (T)0.0;
int jPivot = k;
for(int i = k; i < 4; i++ ) {
T ratio = math::Abs(A.ElementAt(index[i], k))/scale[index[i]];
if(ratio > ratiomax) {
jPivot = i;
ratiomax = ratio;
}
}
// Perform pivoting using row index list
int indexJ = index[k];
if( jPivot != k ) { // Pivot
indexJ = index[jPivot];
index[jPivot] = index[k]; // Swap index jPivot and k
index[k] = indexJ;
signDet *= -1; // Flip sign of determinant
}
// Perform forward elimination
for(int i=k+1; i < 4; i++ ) {
T coeff = A.ElementAt(index[i],k)/A.ElementAt(indexJ,k);
for(int j=k+1; j < 4; j++ )
A.ElementAt(index[i],j) -= coeff*A.ElementAt(indexJ,j);
A.ElementAt(index[i],k) = coeff;
//for( j=1; j< 4; j++ )
// ElementAt(index[i],j) -= A.ElementAt(index[i], k)*ElementAt(indexJ, j);
}
}
for(int i = 0; i < 16; i++)
operator[](i) = A[i];
d = signDet;
// this = A;
return true; */
d = 1; //no permutation still
T scaling[4];
int i,j,k;
//Saving the scvaling information per row
for(i = 0; i < 4; i++) {
T largest = 0.0;
for(j = 0; j < 4; j++) {
T t = math::Abs(this->ElementAt(i, j));
if (t > largest) largest = t;
}
if (largest == 0.0) { //oooppps there is a zero row!
return false;
}
scaling[i] = (T)1.0 / largest;
}
int imax = 0;
for(j = 0; j < 4; j++) {
for(i = 0; i < j; i++) {
T sum = this->ElementAt(i,j);
for(int k = 0; k < i; k++)
sum -= this->ElementAt(i,k)*this->ElementAt(k,j);
this->ElementAt(i,j) = sum;
}
T largest = 0.0;
for(i = j; i < 4; i++) {
T sum = this->ElementAt(i,j);
for(k = 0; k < j; k++)
sum -= this->ElementAt(i,k)*this->ElementAt(k,j);
this->ElementAt(i,j) = sum;
T t = scaling[i] * math::Abs(sum);
if(t >= largest) {
largest = t;
imax = i;
}
}
if (j != imax) {
for (int k = 0; k < 4; k++) {
T dum = this->ElementAt(imax,k);
this->ElementAt(imax,k) = this->ElementAt(j,k);
this->ElementAt(j,k) = dum;
}
d = -(d);
scaling[imax] = scaling[j];
}
index[j]=imax;
if (this->ElementAt(j,j) == 0.0) this->ElementAt(j,j) = (T)TINY;
if (j != 3) {
T dum = (T)1.0 / (this->ElementAt(j,j));
for (i = j+1; i < 4; i++)
this->ElementAt(i,j) *= dum;
}
}
return true;
}
template <class T> Point4<T> LinearSolve<T>::Solve(const Point4<T> &b) {
Point4<T> x(b);
int first = -1, ip;
for(int i = 0; i < 4; i++) {
ip = index[i];
T sum = x[ip];
x[ip] = x[i];
if(first!= -1)
for(int j = first; j <= i-1; j++)
sum -= this->ElementAt(i,j) * x[j];
else
if(sum) first = i;
x[i] = sum;
}
for (int i = 3; i >= 0; i--) {
T sum = x[i];
for (int j = i+1; j < 4; j++)
sum -= this->ElementAt(i, j) * x[j];
x[i] = sum / this->ElementAt(i, i);
}
return x;
}
} //namespace
#endif
#endif

22
vcg/math/point_matching.h
View File

@ -137,16 +137,12 @@ bool ComputeWeightedRigidMatchMatrix(Matrix44x &res,
)
{
Matrix33x tmp;
Matrix33x ccm;
Point3x bfix,bmov; // baricenter of src e trg
ccm.WeightedCrossCovariance(weights,Pmov,Pfix,bmov,bfix);
Matrix33x cyc; // the cyclic components of the cross covariance matrix.
cyc=ccm;
tmp=ccm;
tmp.Transpose();
cyc-=tmp;
cyc=ccm - ccm.transpose();
Matrix44x QQ;
QQ.SetZero();
@ -157,9 +153,7 @@ bool ComputeWeightedRigidMatchMatrix(Matrix44x &res,
RM[0][0]=-ccm.Trace();
RM[1][1]=-ccm.Trace();
RM[2][2]=-ccm.Trace();
RM+=ccm;
ccm.Transpose();
RM+=ccm;
RM += ccm + ccm.transpose();
QQ[0][0] = ccm.Trace();
QQ[0][1] = D[0]; QQ[0][2] = D[1]; QQ[0][3] = D[2];
@ -208,16 +202,12 @@ bool ComputeRigidMatchMatrix(Matrix44x &res,
Point3x &tr)
{
Matrix33x tmp;
Matrix33x ccm;
Point3x bfix,bmov; // baricenter of src e trg
ccm.CrossCovariance(Pmov,Pfix,bmov,bfix);
Matrix33x cyc; // the cyclic components of the cross covariance matrix.
cyc=ccm;
tmp=ccm;
tmp.Transpose();
cyc-=tmp;
cyc=ccm-ccm.transpose();
Matrix44x QQ;
QQ.SetZero();
@ -228,9 +218,7 @@ bool ComputeRigidMatchMatrix(Matrix44x &res,
RM[0][0]=-ccm.Trace();
RM[1][1]=-ccm.Trace();
RM[2][2]=-ccm.Trace();
RM+=ccm;
ccm.Transpose();
RM+=ccm;
RM += ccm + ccm.transpose();
QQ[0][0] = ccm.Trace();
QQ[0][1] = D[0]; QQ[0][2] = D[1]; QQ[0][3] = D[2];
@ -260,7 +248,7 @@ bool ComputeRigidMatchMatrix(Matrix44x &res,
maxv=d[i];
}
// The corresponding eigenvector define the searched rotation,
Matrix44x Rot;
Matrix44x Rot;
q.ToMatrix(Rot);
// the translation (last row) is simply the difference between the transformed src barycenter and the trg baricenter
tr= (bfix - Rot*bmov);

16
vcg/math/polar_decomposition.h
View File

@ -49,35 +49,31 @@ void RotationalPartByPolarDecomposition( const vcg::Matrix33<S> & m, vcg::Matrix
r.SetZero();
s.SetZero();
tmp= m;
tmp.Transpose();
tmp = m*tmp;
tmp = m*m.transpose();
Matrix33<S> res;
Point3<S> e;
bool ss = SingularValueDecomposition<vcg::Matrix33<S> >(tmp,&e[0],res);
res.Transpose();
e[0]=math::Sqrt(e[0]);
e[1]=math::Sqrt(e[1]);
e[2]=math::Sqrt(e[2]);
#ifdef VCG_USE_EIGEN
tmp = tmp*e.asDiagonal()*res;
tmp = tmp*e.asDiagonal()*res.transpose();
#else
tmp = tmp*Matrix33Diag<S>(e)*res;
tmp = tmp*Matrix33Diag<S>(e)*res.transpose();
#endif
bool s1 = SingularValueDecomposition<vcg::Matrix33<S> >(tmp,&e[0],res);
tmp.Transpose();
bool s1 = SingularValueDecomposition<vcg::Matrix33<S> >(tmp,&e[0],res.transpose());
e[0]=1/e[0];
e[1]=1/e[1];
e[2]=1/e[2];
#ifdef VCG_USE_EIGEN
tmp = res*e.asDiagonal()*tmp;
tmp = res*e.asDiagonal()*tmp.transpose();
#else
tmp = res*Matrix33Diag<S>(e)*tmp;
tmp = res*Matrix33Diag<S>(e)*tmp.transpose();
#endif
r = m*tmp;

2
vcg/math/quadric.h
View File

@ -88,7 +88,7 @@ template< class PlaneType >
c = (ScalarType)p.Offset()*p.Offset();
}
void Zero() // Azzera la quadrica
void SetZero() // Azzera la quadrica
{
a[0] = 0;
a[1] = 0;

4
vcg/math/quaternion.h
View File

@ -184,9 +184,9 @@ template <class S> Quaternion<S> Quaternion<S>::operator*(const S &s) const {
template <class S> Quaternion<S> Quaternion<S>::operator*(const Quaternion &q) const {
Point3<S> t1(V(1), V(2), V(3));
Point3<S> t2(q.V(1), q.V(2), q.V(3));
Point3<S> t2(q.V(1), q.V(2), q.V(3));
S d = t2 * t1;
S d = t2.dot(t1);
Point3<S> t3 = t1 ^ t2;
t1 *= q.V(0);

2
vcg/math/similarity.h
View File

@ -191,7 +191,7 @@ template <class S,class RotationType> Matrix44<S> Similarity<S,RotationType>::Ma
rot.ToMatrix(r);
Matrix44<S> s = Matrix44<S>().SetScale(sca, sca, sca);
Matrix44<S> t = Matrix44<S>().SetTranslate(tra[0], tra[1], tra[2]);
return s*r*t; // trans * scale * rot;
return Matrix44<S>(s*r*t); // trans * scale * rot;
}
template <class S,class RotationType> Matrix44<S> Similarity<S,RotationType>::InverseMatrix() const {

109
vcg/space/color4.h
View File

@ -89,6 +89,7 @@ namespace vcg {
template <class T>
class Color4 : public Point4<T>
{
typedef Point4<T> Base;
public:
/// Constant for storing standard colors.
/// Each color is stored in a simple in so that the bit pattern match with the one of Color4b.
@ -121,22 +122,30 @@ public:
inline Color4 ( const Point4<T> &c) :Point4<T>(c) {};
inline Color4 (){};
inline Color4 (ColorConstant cc);
#ifdef VCG_USE_EIGEN
template<typename OtherDerived>
inline Color4 (const Eigen::MatrixBase<OtherDerived>& other) : Base(other)
{
// TODO make sure the types are the same
}
#endif
template <class Q>
inline void Import(const Color4<Q> & b )
{
Point4<T>::_v[0] = T(b[0]);
Point4<T>::_v[1] = T(b[1]);
Point4<T>::_v[2] = T(b[2]);
Point4<T>::_v[3] = T(b[3]);
Point4<T>::V()[0] = T(b[0]);
Point4<T>::V()[1] = T(b[1]);
Point4<T>::V()[2] = T(b[2]);
Point4<T>::V()[3] = T(b[3]);
}
template <class Q>
inline void Import(const Point4<Q> & b )
{
Point4<T>::_v[0] = T(b[0]);
Point4<T>::_v[1] = T(b[1]);
Point4<T>::_v[2] = T(b[2]);
Point4<T>::_v[3] = T(b[3]);
Point4<T>::V()[0] = T(b[0]);
Point4<T>::V()[1] = T(b[1]);
Point4<T>::V()[2] = T(b[2]);
Point4<T>::V()[3] = T(b[3]);
}
template <class Q>
@ -150,7 +159,7 @@ public:
inline Color4 operator + ( const Color4 & p) const
{
return Color4( this->_v[0]+p._v[0], this->_v[1]+p._v[1], this->_v[2]+p._v[2], this->_v[3]+p._v[3] );
return Color4( (*this)[0]+p.V()[0], (*this)[1]+p.V()[1], (*this)[2]+p.V()[2], (*this)[3]+p.V()[3] );
}
@ -161,20 +170,20 @@ public:
inline void SetRGB( unsigned char r, unsigned char g, unsigned char b )
{
Point4<T>::_v[0] = r;
Point4<T>::_v[1] = g;
Point4<T>::_v[2] = b;
Point4<T>::_v[3] = 0;
Point4<T>::V()[0] = r;
Point4<T>::V()[1] = g;
Point4<T>::V()[2] = b;
Point4<T>::V()[3] = 0;
}
void SetHSVColor( float h, float s, float v){
float r,g,b;
if(s==0.0){ // gray color
r = g = b = v;
Point4<T>::_v[0]=(unsigned char)(255*r);
Point4<T>::_v[1]=(unsigned char)(255*g);
Point4<T>::_v[2]=(unsigned char)(255*b);
Point4<T>::_v[3]=255;
Point4<T>::V()[0]=(unsigned char)(255*r);
Point4<T>::V()[1]=(unsigned char)(255*g);
Point4<T>::V()[2]=(unsigned char)(255*b);
Point4<T>::V()[3]=255;
return;
}
if(h==1.0) h = 0.0;
@ -194,11 +203,11 @@ public:
case 4: r=t; g=p; b=v; break;
case 5: r=v; g=p; b=q; break;
}
Point4<T>::_v[0]=(unsigned char)(255*r);
Point4<T>::_v[1]=(unsigned char)(255*g);
Point4<T>::_v[2]=(unsigned char)(255*b);
Point4<T>::_v[3]=255;
// _v[0]=r*256;_v[1]=g*256;_v[2]=b*256;
Point4<T>::V()[0]=(unsigned char)(255*r);
Point4<T>::V()[1]=(unsigned char)(255*g);
Point4<T>::V()[2]=(unsigned char)(255*b);
Point4<T>::V()[3]=255;
// V()[0]=r*256;V()[1]=g*256;V()[2]=b*256;
}
inline static Color4 GrayShade(float f)
@ -246,10 +255,10 @@ inline void Color4<T>::lerp(const Color4<T> &c0, const Color4<T> &c1, const floa
assert(x>=0);
assert(x<=1);
Point4<T>::_v[0]=(T)(c1._v[0]*x + c0._v[0]*(1.0f-x));
Point4<T>::_v[1]=(T)(c1._v[1]*x + c0._v[1]*(1.0f-x));
Point4<T>::_v[2]=(T)(c1._v[2]*x + c0._v[2]*(1.0f-x));
Point4<T>::_v[3]=(T)(c1._v[3]*x + c0._v[3]*(1.0f-x));
Point4<T>::V()[0]=(T)(c1.V()[0]*x + c0.V()[0]*(1.0f-x));
Point4<T>::V()[1]=(T)(c1.V()[1]*x + c0.V()[1]*(1.0f-x));
Point4<T>::V()[2]=(T)(c1.V()[2]*x + c0.V()[2]*(1.0f-x));
Point4<T>::V()[3]=(T)(c1.V()[3]*x + c0.V()[3]*(1.0f-x));
}
template <class T>
@ -257,10 +266,10 @@ inline void Color4<T>::lerp(const Color4<T> &c0, const Color4<T> &c1, const Colo
{
assert(fabs(ip[0]+ip[1]+ip[2]-1)<0.00001);
Point4<T>::_v[0]=(T)(c0[0]*ip[0] + c1[0]*ip[1]+ c2[0]*ip[2]);
Point4<T>::_v[1]=(T)(c0[1]*ip[0] + c1[1]*ip[1]+ c2[1]*ip[2]);
Point4<T>::_v[2]=(T)(c0[2]*ip[0] + c1[2]*ip[1]+ c2[2]*ip[2]);
Point4<T>::_v[3]=(T)(c0[3]*ip[0] + c1[3]*ip[1]+ c2[3]*ip[2]);
Point4<T>::V()[0]=(T)(c0[0]*ip[0] + c1[0]*ip[1]+ c2[0]*ip[2]);
Point4<T>::V()[1]=(T)(c0[1]*ip[0] + c1[1]*ip[1]+ c2[1]*ip[2]);
Point4<T>::V()[2]=(T)(c0[2]*ip[0] + c1[2]*ip[1]+ c2[2]*ip[2]);
Point4<T>::V()[3]=(T)(c0[3]*ip[0] + c1[3]*ip[1]+ c2[3]*ip[2]);
}
@ -292,10 +301,10 @@ template <>
template <>
inline void Color4<float>::Import(const Color4<unsigned char> &b)
{
this->_v[0]=b[0]/255.0f;
this->_v[1]=b[1]/255.0f;
this->_v[2]=b[2]/255.0f;
this->_v[3]=b[3]/255.0f;
(*this)[0]=b[0]/255.0f;
(*this)[1]=b[1]/255.0f;
(*this)[2]=b[2]/255.0f;
(*this)[3]=b[3]/255.0f;
}
#if !defined(__GNUC__) || (__GNUC__ > 3)
@ -304,10 +313,10 @@ template <> // [Bug c++/14479] enum definition in template class with template m
template <>
inline void Color4<unsigned char>::Import(const Color4<float> &b)
{
this->_v[0]=(unsigned char)(b[0]*255.0f);
this->_v[1]=(unsigned char)(b[1]*255.0f);
this->_v[2]=(unsigned char)(b[2]*255.0f);
this->_v[3]=(unsigned char)(b[3]*255.0f);
(*this)[0]=(unsigned char)(b[0]*255.0f);
(*this)[1]=(unsigned char)(b[1]*255.0f);
(*this)[2]=(unsigned char)(b[2]*255.0f);
(*this)[3]=(unsigned char)(b[3]*255.0f);
}
#if !defined(__GNUC__) || (__GNUC__ > 3)
@ -316,10 +325,10 @@ template <> // [Bug c++/14479] enum definition in template class with template m
template <>
inline void Color4<unsigned char>::Import(const Point4<float> &b)
{
this->_v[0]=(unsigned char)(b[0]*255.0f);
this->_v[1]=(unsigned char)(b[1]*255.0f);
this->_v[2]=(unsigned char)(b[2]*255.0f);
this->_v[3]=(unsigned char)(b[3]*255.0f);
(*this)[0]=(unsigned char)(b[0]*255.0f);
(*this)[1]=(unsigned char)(b[1]*255.0f);
(*this)[2]=(unsigned char)(b[2]*255.0f);
(*this)[3]=(unsigned char)(b[3]*255.0f);
}
#if !defined(__GNUC__) || (__GNUC__ > 3)
@ -351,10 +360,10 @@ inline Color4<float> Color4<float>::Construct( const Color4<unsigned char> & b )
//template <class T,class S>
//inline void Color4<T>::Import(const Color4<S> &b)
//{
// _v[0] = T(b[0]);
// _v[1] = T(b[1]);
// _v[2] = T(b[2]);
// _v[3] = T(b[3]);
// V()[0] = T(b[0]);
// V()[1] = T(b[1]);
// V()[2] = T(b[2]);
// V()[3] = T(b[3]);
//}
//
template<>
@ -382,10 +391,10 @@ template<>
inline Color4<unsigned char> Color4<unsigned char>::operator + ( const Color4<unsigned char> & p) const
{
return Color4<unsigned char>(
math::Clamp(int(this->_v[0])+int(p._v[0]),0,255),
math::Clamp(int(this->_v[1])+int(p._v[1]),0,255),
math::Clamp(int(this->_v[2])+int(p._v[2]),0,255),
math::Clamp(int(this->_v[3])+int(p._v[3]),0,255)
math::Clamp(int((*this)[0])+int(p[0]),0,255),
math::Clamp(int((*this)[1])+int(p[1]),0,255),
math::Clamp(int((*this)[2])+int(p[2]),0,255),
math::Clamp(int((*this)[3])+int(p[3]),0,255)
);
}

572
vcg/space/deprecated_point3.h Normal file
View File

@ -0,0 +1,572 @@
/****************************************************************************
* VCGLib o o *
* Visual and Computer Graphics Library o o *
* _ O _ *
* Copyright(C) 2004 \/)\/ *
* Visual Computing Lab /\/| *
* ISTI - Italian National Research Council | *
* \ *
* All rights reserved. *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
* for more details. *
* *
****************************************************************************/
/****************************************************************************
History
$Log: not supported by cvs2svn $
Revision 1.26 2006/11/13 13:03:45 ponchio
Added GetBBox in Point3 (declaration) the body of the function is in box3.h
Revision 1.25 2006/10/13 12:59:24 cignoni
Added **explicit** constructor from three coords of a different scalartype
Revision 1.24 2006/09/28 13:37:35 m_di_benedetto
added non const * V()
Revision 1.23 2005/11/09 16:11:55 cignoni
Added Abs and LowClampToZero
Revision 1.22 2005/09/14 14:09:21 m_di_benedetto
Added specialized Convert() for the same scalar type.
Revision 1.21 2005/05/06 14:45:33 spinelli
cambiato parentesi nel costruttore di GetUV per rendere compatibile tale costruttore con MVC e borland
Revision 1.20 2005/04/27 16:05:19 callieri
line 466, added parentesis on default value creator getUV [borland]
Revision 1.19 2004/11/09 15:49:07 ganovelli
added GetUV
Revision 1.18 2004/10/13 12:45:51 cignoni
Better Doxygen documentation
Revision 1.17 2004/09/10 14:01:40 cignoni
Added polar to cartesian
Revision 1.16 2004/03/21 17:14:56 ponchio
Added a math::
Revision 1.15 2004/03/07 22:45:32 cignoni
Moved quality and normal functions to the triangle class.
Revision 1.14 2004/03/05 17:55:01 tarini
errorino: upper case in Zero()
Revision 1.13 2004/03/03 14:22:48 cignoni
Yet against cr lf mismatch
Revision 1.12 2004/02/23 23:42:26 cignoni
Translated comments, removed unusued stuff. corrected linefeed/cr
Revision 1.11 2004/02/19 16:12:28 cignoni
cr lf mismatch 2
Revision 1.10 2004/02/19 16:06:24 cignoni
cr lf mismatch
Revision 1.8 2004/02/19 15:13:40 cignoni
corrected sqrt and added doxygen groups
Revision 1.7 2004/02/17 02:08:47 cignoni
Di prova...
Revision 1.6 2004/02/15 23:35:47 cignoni
Cambiato nome type template in accordo alla styleguide
Revision 1.5 2004/02/10 01:07:15 cignoni
Edited Comments and GPL license
Revision 1.4 2004/02/09 13:48:02 cignoni
Edited doxygen comments
****************************************************************************/
#ifndef __VCGLIB_POINT3
#define __VCGLIB_POINT3
#include <assert.h>
#include <vcg/math/base.h>
namespace vcg {
/** \addtogroup space */
/*@{*/
/**
The templated class for representing a point in 3D space.
The class is templated over the ScalarType class that is used to represent coordinates. All the usual
operator overloading (* + - ...) is present.
*/
template <class T> class Box3;
template <class P3ScalarType> class Point3
{
protected:
/// The only data member. Hidden to user.
P3ScalarType _v[3];
public:
typedef P3ScalarType ScalarType;
enum {Dimension = 3};
//@{
/** @name Standard Constructors and Initializers
No casting operators have been introduced to avoid automatic unattended (and costly) conversion between different point types
**/
inline Point3 () { }
inline Point3 ( const P3ScalarType nx, const P3ScalarType ny, const P3ScalarType nz )
{
_v[0] = nx;
_v[1] = ny;
_v[2] = nz;
}
inline Point3 ( Point3 const & p )
{
_v[0]= p._v[0];
_v[1]= p._v[1];
_v[2]= p._v[2];
}
inline Point3 ( const P3ScalarType nv[3] )
{
_v[0] = nv[0];
_v[1] = nv[1];
_v[2] = nv[2];
}
inline Point3 & operator =( Point3 const & p )
{
_v[0]= p._v[0]; _v[1]= p._v[1]; _v[2]= p._v[2];
return *this;
}
inline void SetZero()
{
_v[0] = 0;
_v[1] = 0;
_v[2] = 0;
}
/// Padding function: give a default 0 value to all the elements that are not in the [0..2] range.
/// Useful for managing in a consistent way object that could have point2 / point3 / point4
inline P3ScalarType Ext( const int i ) const
{
if(i>=0 && i<=2) return _v[i];
else return 0;
}
template <class Q>
inline void Import( const Point3<Q> & b )
{
_v[0] = P3ScalarType(b[0]);
_v[1] = P3ScalarType(b[1]);
_v[2] = P3ScalarType(b[2]);
}
template <class Q>
static inline Point3 Construct( const Point3<Q> & b )
{
return Point3(P3ScalarType(b[0]),P3ScalarType(b[1]),P3ScalarType(b[2]));
}
template <class Q>
static inline Point3 Construct( const Q & P0, const Q & P1, const Q & P2)
{
return Point3(P3ScalarType(P0),P3ScalarType(P1),P3ScalarType(P2));
}
static inline Point3 Construct( const Point3<ScalarType> & b )
{
return b;
}
//@}
//@{
/** @name Data Access.
access to data is done by overloading of [] or explicit naming of coords (x,y,z)**/
inline P3ScalarType & operator [] ( const int i )
{
assert(i>=0 && i<3);
return _v[i];
}
inline const P3ScalarType & operator [] ( const int i ) const
{
assert(i>=0 && i<3);
return _v[i];
}
inline const P3ScalarType &X() const { return _v[0]; }
inline const P3ScalarType &Y() const { return _v[1]; }
inline const P3ScalarType &Z() const { return _v[2]; }
inline P3ScalarType &X() { return _v[0]; }
inline P3ScalarType &Y() { return _v[1]; }
inline P3ScalarType &Z() { return _v[2]; }
inline const P3ScalarType * V() const
{
return _v;
}
inline P3ScalarType * V()
{
return _v;
}
inline P3ScalarType & V( const int i )
{
assert(i>=0 && i<3);
return _v[i];
}
inline const P3ScalarType & V( const int i ) const
{
assert(i>=0 && i<3);
return _v[i];
}
//@}
//@{
/** @name Classical overloading of operators
Note
**/
inline Point3 operator + ( Point3 const & p) const
{
return Point3<P3ScalarType>( _v[0]+p._v[0], _v[1]+p._v[1], _v[2]+p._v[2] );
}
inline Point3 operator - ( Point3 const & p) const
{
return Point3<P3ScalarType>( _v[0]-p._v[0], _v[1]-p._v[1], _v[2]-p._v[2] );
}
inline Point3 operator * ( const P3ScalarType s ) const
{
return Point3<P3ScalarType>( _v[0]*s, _v[1]*s, _v[2]*s );
}
inline Point3 operator / ( const P3ScalarType s ) const
{
return Point3<P3ScalarType>( _v[0]/s, _v[1]/s, _v[2]/s );
}
/// Dot product
inline P3ScalarType operator * ( Point3 const & p ) const
{
return ( _v[0]*p._v[0] + _v[1]*p._v[1] + _v[2]*p._v[2] );
}
inline P3ScalarType dot( const Point3 & p ) const { return (*this) * p; }
/// Cross product
inline Point3 operator ^ ( Point3 const & p ) const
{
return Point3 <P3ScalarType>
(
_v[1]*p._v[2] - _v[2]*p._v[1],
_v[2]*p._v[0] - _v[0]*p._v[2],
_v[0]*p._v[1] - _v[1]*p._v[0]
);
}
inline Point3 & operator += ( Point3 const & p)
{
_v[0] += p._v[0];
_v[1] += p._v[1];
_v[2] += p._v[2];
return *this;
}
inline Point3 & operator -= ( Point3 const & p)
{
_v[0] -= p._v[0];
_v[1] -= p._v[1];
_v[2] -= p._v[2];
return *this;
}
inline Point3 & operator *= ( const P3ScalarType s )
{
_v[0] *= s;
_v[1] *= s;
_v[2] *= s;
return *this;
}
inline Point3 & operator /= ( const P3ScalarType s )
{
_v[0] /= s;
_v[1] /= s;
_v[2] /= s;
return *this;
}
// Norme
inline P3ScalarType Norm() const
{
return math::Sqrt( _v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] );
}
inline P3ScalarType SquaredNorm() const
{
return ( _v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] );
}
// Scalatura differenziata
inline Point3 & Scale( const P3ScalarType sx, const P3ScalarType sy, const P3ScalarType sz )
{
_v[0] *= sx;
_v[1] *= sy;
_v[2] *= sz;
return *this;
}
inline Point3 & Scale( const Point3 & p )
{
_v[0] *= p._v[0];
_v[1] *= p._v[1];
_v[2] *= p._v[2];
return *this;
}
// Normalizzazione
inline Point3 & Normalize()
{
P3ScalarType n = math::Sqrt(_v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2]);
if(n>0.0) { _v[0] /= n; _v[1] /= n; _v[2] /= n; }
return *this;
}
// for compatibility with eigen port
inline Point3 & normalized() { return Normalize(); }
/**
* Convert to polar coordinates from cartesian coordinates.
*
* Theta is the azimuth angle and ranges between [0, 360) degrees.
* Phi is the elevation angle (not the polar angle) and ranges between [-90, 90] degrees.
*
* /note Note that instead of the classical polar angle, which ranges between
* 0 and 180 degrees we opt for the elevation angle to obtain a more
* intuitive spherical coordinate system.
*/
void ToPolar(P3ScalarType &ro, P3ScalarType &theta, P3ScalarType &phi) const
{
ro = Norm();
theta = (P3ScalarType)atan2(_v[2], _v[0]);
phi = (P3ScalarType)asin(_v[1]/ro);
}
/**
* Convert from polar coordinates to cartesian coordinates.
*
* Theta is the azimuth angle and ranges between [0, 360) degrees.
* Phi is the elevation angle (not the polar angle) and ranges between [-90, 90] degrees.
*
* \note Note that instead of the classical polar angle, which ranges between
* 0 and 180 degrees, we opt for the elevation angle to obtain a more
* intuitive spherical coordinate system.
*/
void FromPolar(const P3ScalarType &ro, const P3ScalarType &theta, const P3ScalarType &phi)
{
_v[0]= ro*cos(theta)*cos(phi);
_v[1]= ro*sin(phi);
_v[2]= ro*sin(theta)*cos(phi);
}
Box3<P3ScalarType> GetBBox(Box3<P3ScalarType> &bb) const;
//@}
//@{
/** @name Comparison Operators.
Note that the reverse z prioritized ordering, useful in many situations.
**/
inline bool operator == ( Point3 const & p ) const
{
return _v[0]==p._v[0] && _v[1]==p._v[1] && _v[2]==p._v[2];
}
inline bool operator != ( Point3 const & p ) const
{
return _v[0]!=p._v[0] || _v[1]!=p._v[1] || _v[2]!=p._v[2];
}
inline bool operator < ( Point3 const & p ) const
{
return (_v[2]!=p._v[2])?(_v[2]<p._v[2]):
(_v[1]!=p._v[1])?(_v[1]<p._v[1]):
(_v[0]<p._v[0]);
}
inline bool operator > ( Point3 const & p ) const
{
return (_v[2]!=p._v[2])?(_v[2]>p._v[2]):
(_v[1]!=p._v[1])?(_v[1]>p._v[1]):
(_v[0]>p._v[0]);
}
inline bool operator <= ( Point3 const & p ) const
{
return (_v[2]!=p._v[2])?(_v[2]< p._v[2]):
(_v[1]!=p._v[1])?(_v[1]< p._v[1]):
(_v[0]<=p._v[0]);
}
inline bool operator >= ( Point3 const & p ) const
{
return (_v[2]!=p._v[2])?(_v[2]> p._v[2]):
(_v[1]!=p._v[1])?(_v[1]> p._v[1]):
(_v[0]>=p._v[0]);
}
inline Point3 operator - () const
{
return Point3<P3ScalarType> ( -_v[0], -_v[1], -_v[2] );
}
//@}
}; // end class definition
template <class P3ScalarType>
inline P3ScalarType Angle( Point3<P3ScalarType> const & p1, Point3<P3ScalarType> const & p2 )
{
P3ScalarType w = p1.Norm()*p2.Norm();
if(w==0) return -1;
P3ScalarType t = (p1*p2)/w;
if(t>1) t = 1;
else if(t<-1) t = -1;
return (P3ScalarType) acos(t);
}
// versione uguale alla precedente ma che assume che i due vettori sono unitari
template <class P3ScalarType>
inline P3ScalarType AngleN( Point3<P3ScalarType> const & p1, Point3<P3ScalarType> const & p2 )
{
P3ScalarType w = p1*p2;
if(w>1)
w = 1;
else if(w<-1)
w=-1;
return (P3ScalarType) acos(w);
}
template <class P3ScalarType>
inline P3ScalarType Norm( Point3<P3ScalarType> const & p )
{
return p.Norm();
}
template <class P3ScalarType>
inline P3ScalarType SquaredNorm( Point3<P3ScalarType> const & p )
{
return p.SquaredNorm();
}
template <class P3ScalarType>
inline Point3<P3ScalarType> & Normalize( Point3<P3ScalarType> & p )
{
p.Normalize();
return p;
}
template <class P3ScalarType>
inline P3ScalarType Distance( Point3<P3ScalarType> const & p1,Point3<P3ScalarType> const & p2 )
{
return (p1-p2).Norm();
}
template <class P3ScalarType>
inline P3ScalarType SquaredDistance( Point3<P3ScalarType> const & p1,Point3<P3ScalarType> const & p2 )
{
return (p1-p2).SquaredNorm();
}
// Dot product preciso numericamente (solo double!!)
// Implementazione: si sommano i prodotti per ordine di esponente
// (prima le piu' grandi)
template<class P3ScalarType>
double stable_dot ( Point3<P3ScalarType> const & p0, Point3<P3ScalarType> const & p1 )
{
P3ScalarType k0 = p0._v[0]*p1._v[0];
P3ScalarType k1 = p0._v[1]*p1._v[1];
P3ScalarType k2 = p0._v[2]*p1._v[2];
int exp0,exp1,exp2;
frexp( double(k0), &exp0 );
frexp( double(k1), &exp1 );
frexp( double(k2), &exp2 );
if( exp0<exp1 )
{
if(exp0<exp2)
return (k1+k2)+k0;
else
return (k0+k1)+k2;
}
else
{
if(exp1<exp2)
return(k0+k2)+k1;
else
return (k0+k1)+k2;
}
}
/// Point(p) Edge(v1-v2) dist, q is the point in v1-v2 with min dist
template<class P3ScalarType>
P3ScalarType PSDist( const Point3<P3ScalarType> & p,
const Point3<P3ScalarType> & v1,
const Point3<P3ScalarType> & v2,
Point3<P3ScalarType> & q )
{
Point3<P3ScalarType> e = v2-v1;
P3ScalarType t = ((p-v1)*e)/e.SquaredNorm();
if(t<0) t = 0;
else if(t>1) t = 1;
q = v1+e*t;
return Distance(p,q);
}
template <class P3ScalarType>
void GetUV( Point3<P3ScalarType> &n,Point3<P3ScalarType> &u, Point3<P3ScalarType> &v, Point3<P3ScalarType> up=(Point3<P3ScalarType>(0,1,0)) )
{
n.Normalize();
const double LocEps=double(1e-7);
u=n^up;
double len = u.Norm();
if(len < LocEps)
{
if(fabs(n[0])<fabs(n[1])){
if(fabs(n[0])<fabs(n[2])) up=Point3<P3ScalarType>(1,0,0); // x is the min
else up=Point3<P3ScalarType>(0,0,1); // z is the min
}else {
if(fabs(n[1])<fabs(n[2])) up=Point3<P3ScalarType>(0,1,0); // y is the min
else up=Point3<P3ScalarType>(0,0,1); // z is the min
}
u=n^up;
}
u.Normalize();
v=n^u;
v.Normalize();
Point3<P3ScalarType> 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<int> Point3i;
typedef Point3<float> Point3f;
typedef Point3<double> Point3d;
/*@}*/
} // end namespace
#endif

383
vcg/space/deprecated_point4.h Normal file
View File

@ -0,0 +1,383 @@
/****************************************************************************
* VCGLib o o *
* Visual and Computer Graphics Library o o *
* _ O _ *
* Copyright(C) 2004 \/)\/ *
* Visual Computing Lab /\/| *
* ISTI - Italian National Research Council | *
* \ *
* All rights reserved. *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
* for more details. *
* *
****************************************************************************/
/****************************************************************************
History
$Log: not supported by cvs2svn $
Revision 1.12 2006/06/21 11:06:16 ganovelli
changed return type of Zero() (to void)
Revision 1.11 2005/04/13 09:40:30 ponchio
Including math/bash.h
Revision 1.10 2005/03/18 16:34:42 fiorin
minor changes to comply gcc compiler
Revision 1.9 2005/01/21 18:02:11 ponchio
Removed dependence from matrix44 and changed VectProd
Revision 1.8 2005/01/12 11:25:02 ganovelli
added Dimension
Revision 1.7 2004/10/11 17:46:11 ganovelli
added definition of vector product (not implemented)
Revision 1.6 2004/05/10 11:16:19 ganovelli
include assert.h added
Revision 1.5 2004/03/31 10:09:58 cignoni
missing return value in zero()
Revision 1.4 2004/03/11 17:17:49 tarini
added commets (doxy), uniformed with new style, now using math::, ...
added HomoNormalize(), Zero()... remade StableDot() (hand made sort).
Revision 1.1 2004/02/10 01:11:28 cignoni
Edited Comments and GPL license
****************************************************************************/
#ifndef __VCGLIB_POINT4
#define __VCGLIB_POINT4
#include <assert.h>
#include <vcg/math/base.h>
namespace vcg {
/** \addtogroup space */
/*@{*/
/**
The templated class for representing a point in 4D space.
The class is templated over the ScalarType class that is used to represent coordinates.
All the usual operator (* + - ...) are defined.
*/
template <class T> class Point4
{
public:
/// The only data member. Hidden to user.
T _v[4];
public:
typedef T ScalarType;
enum {Dimension = 4};
//@{
/** @name Standard Constructors and Initializers
No casting operators have been introduced to avoid automatic unattended (and costly) conversion between different point types
**/
inline Point4 () { }
inline Point4 ( const T nx, const T ny, const T nz , const T nw )
{
_v[0] = nx; _v[1] = ny; _v[2] = nz; _v[3] = nw;
}
inline Point4 ( const T p[4] )
{
_v[0] = p[0]; _v[1]= p[1]; _v[2] = p[2]; _v[3]= p[3];
}
inline Point4 ( const Point4 & p )
{
_v[0]= p._v[0]; _v[1]= p._v[1]; _v[2]= p._v[2]; _v[3]= p._v[3];
}
inline void Zero()
{
_v[0] = _v[1] = _v[2] = _v[3]= 0;
}
template <class Q>
/// importer from different Point4 types
inline void Import( const Point4<Q> & b )
{
_v[0] = T(b[0]); _v[1] = T(b[1]);
_v[2] = T(b[2]);
_v[3] = T(b[3]);
}
/// constuctor that imports from different Point4 types
template <class Q>
static inline Point4 Construct( const Point4<Q> & b )
{
return Point4(T(b[0]),T(b[1]),T(b[2]),T(b[3]));
}
//@}
//@{
/** @name Data Access.
access to data is done by overloading of [] or explicit naming of coords (x,y,z,w)
**/
inline const T & operator [] ( const int i ) const
{
assert(i>=0 && i<4);
return _v[i];
}
inline T & operator [] ( const int i )
{
assert(i>=0 && i<4);
return _v[i];
}
inline T &X() {return _v[0];}
inline T &Y() {return _v[1];}
inline T &Z() {return _v[2];}
inline T &W() {return _v[3];}
inline T const * V() const
{
return _v;
}
inline const T & V ( const int i ) const
{
assert(i>=0 && i<4);
return _v[i];
}
inline T & V ( const int i )
{
assert(i>=0 && i<4);
return _v[i];
}
/// Padding function: give a default 0 value to all the elements that are not in the [0..2] range.
/// Useful for managing in a consistent way object that could have point2 / point3 / point4
inline T Ext( const int i ) const
{
if(i>=0 && i<=3) return _v[i];
else return 0;
}
//@}
//@{
/** @name Linear operators and the likes
**/
inline Point4 operator + ( const Point4 & p) const
{
return Point4( _v[0]+p._v[0], _v[1]+p._v[1], _v[2]+p._v[2], _v[3]+p._v[3] );
}
inline Point4 operator - ( const Point4 & p) const
{
return Point4( _v[0]-p._v[0], _v[1]-p._v[1], _v[2]-p._v[2], _v[3]-p._v[3] );
}
inline Point4 operator * ( const T s ) const
{
return Point4( _v[0]*s, _v[1]*s, _v[2]*s, _v[3]*s );
}
inline Point4 operator / ( const T s ) const
{
return Point4( _v[0]/s, _v[1]/s, _v[2]/s, _v[3]/s );
}
inline Point4 & operator += ( const Point4 & p)
{
_v[0] += p._v[0]; _v[1] += p._v[1]; _v[2] += p._v[2]; _v[3] += p._v[3];
return *this;
}
inline Point4 & operator -= ( const Point4 & p )
{
_v[0] -= p._v[0]; _v[1] -= p._v[1]; _v[2] -= p._v[2]; _v[3] -= p._v[3];
return *this;
}
inline Point4 & operator *= ( const T s )
{
_v[0] *= s; _v[1] *= s; _v[2] *= s; _v[3] *= s;
return *this;
}
inline Point4 & operator /= ( const T s )
{
_v[0] /= s; _v[1] /= s; _v[2] /= s; _v[3] /= s;
return *this;
}
inline Point4 operator - () const
{
return Point4( -_v[0], -_v[1], -_v[2], -_v[3] );
}
inline Point4 VectProd ( const Point4 &x, const Point4 &z ) const
{
Point4 res;
const Point4 &y = *this;
res[0] = y[1]*x[2]*z[3]-y[1]*x[3]*z[2]-x[1]*y[2]*z[3]+
x[1]*y[3]*z[2]+z[1]*y[2]*x[3]-z[1]*y[3]*x[2];
res[1] = y[0]*x[3]*z[2]-z[0]*y[2]*x[3]-y[0]*x[2]*
z[3]+z[0]*y[3]*x[2]+x[0]*y[2]*z[3]-x[0]*y[3]*z[2];
res[2] = -y[0]*z[1]*x[3]+x[0]*z[1]*y[3]+y[0]*x[1]*
z[3]-x[0]*y[1]*z[3]-z[0]*x[1]*y[3]+z[0]*y[1]*x[3];
res[3] = -z[0]*y[1]*x[2]-y[0]*x[1]*z[2]+x[0]*y[1]*
z[2]+y[0]*z[1]*x[2]-x[0]*z[1]*y[2]+z[0]*x[1]*y[2];
return res;
}
//@}
//@{
/** @name Norms and normalizations
**/
/// Euclidian normal
inline T Norm() const
{
return math::Sqrt( _v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] + _v[3]*_v[3] );
}
/// Squared euclidian normal
inline T SquaredNorm() const
{
return _v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] + _v[3]*_v[3];
}
/// Euclidian normalization
inline Point4 & Normalize()
{
T n = sqrt(_v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] + _v[3]*_v[3] );
if(n>0.0) { _v[0] /= n; _v[1] /= n; _v[2] /= n; _v[3] /= n; }
return *this;
}
/// Homogeneous normalization (division by W)
inline Point4 & HomoNormalize(){
if (_v[3]!=0.0) { _v[0] /= _v[3]; _v[1] /= _v[3]; _v[2] /= _v[3]; _v[3]=1.0; }
return *this;
};
//@}
//@{
/** @name Comparison operators (lexicographical order)
**/
inline bool operator == ( const Point4& p ) const
{
return _v[0]==p._v[0] && _v[1]==p._v[1] && _v[2]==p._v[2] && _v[3]==p._v[3];
}
inline bool operator != ( const Point4 & p ) const
{
return _v[0]!=p._v[0] || _v[1]!=p._v[1] || _v[2]!=p._v[2] || _v[3]!=p._v[3];
}
inline bool operator < ( Point4 const & p ) const
{
return (_v[3]!=p._v[3])?(_v[3]<p._v[3]):
(_v[2]!=p._v[2])?(_v[2]<p._v[2]):
(_v[1]!=p._v[1])?(_v[1]<p._v[1]):
(_v[0]<p._v[0]);
}
inline bool operator > ( const Point4 & p ) const
{
return (_v[3]!=p._v[3])?(_v[3]>p._v[3]):
(_v[2]!=p._v[2])?(_v[2]>p._v[2]):
(_v[1]!=p._v[1])?(_v[1]>p._v[1]):
(_v[0]>p._v[0]);
}
inline bool operator <= ( const Point4 & p ) const
{
return (_v[3]!=p._v[3])?(_v[3]< p._v[3]):
(_v[2]!=p._v[2])?(_v[2]< p._v[2]):
(_v[1]!=p._v[1])?(_v[1]< p._v[1]):
(_v[0]<=p._v[0]);
}
inline bool operator >= ( const Point4 & p ) const
{
return (_v[3]!=p._v[3])?(_v[3]> p._v[3]):
(_v[2]!=p._v[2])?(_v[2]> p._v[2]):
(_v[1]!=p._v[1])?(_v[1]> p._v[1]):
(_v[0]>=p._v[0]);
}
//@}
//@{
/** @name Dot products
**/
// dot product
inline T operator * ( const Point4 & p ) const
{
return _v[0]*p._v[0] + _v[1]*p._v[1] + _v[2]*p._v[2] + _v[3]*p._v[3];
}
inline T dot( const Point4 & p ) const { return (*this) * p; }
inline Point4 operator ^ ( const Point4& p ) const
{
assert(0);// not defined by two vectors (only put for metaprogramming)
return Point4();
}
/// slower version, more stable (double precision only)
T StableDot ( const Point4<T> & p ) const
{
T k0=_v[0]*p._v[0], k1=_v[1]*p._v[1], k2=_v[2]*p._v[2], k3=_v[3]*p._v[3];
int exp0,exp1,exp2,exp3;
frexp( double(k0), &exp0 );frexp( double(k1), &exp1 );
frexp( double(k2), &exp2 );frexp( double(k3), &exp3 );
if (exp0>exp1) { math::Swap(k0,k1); math::Swap(exp0,exp1); }
if (exp2>exp3) { math::Swap(k2,k3); math::Swap(exp2,exp3); }
if (exp0>exp2) { math::Swap(k0,k2); math::Swap(exp0,exp2); }
if (exp1>exp3) { math::Swap(k1,k3); math::Swap(exp1,exp3); }
if (exp2>exp3) { math::Swap(k2,k3); math::Swap(exp2,exp3); }
return ( (k0 + k1) + k2 ) +k3;
}
//@}
}; // end class definition
template <class T>
T Angle( const Point4<T>& p1, const Point4<T> & p2 )
{
T w = p1.Norm()*p2.Norm();
if(w==0) return -1;
T t = (p1*p2)/w;
if(t>1) t=1;
return T( math::Acos(t) );
}
template <class T>
inline T Norm( const Point4<T> & p )
{
return p.Norm();
}
template <class T>
inline T SquaredNorm( const Point4<T> & p )
{
return p.SquaredNorm();
}
template <class T>
inline T Distance( const Point4<T> & p1, const Point4<T> & p2 )
{
return Norm(p1-p2);
}
template <class T>
inline T SquaredDistance( const Point4<T> & p1, const Point4<T> & p2 )
{
return SquaredNorm(p1-p2);
}
/// slower version of dot product, more stable (double precision only)
template<class T>
double StableDot ( Point4<T> const & p0, Point4<T> const & p1 )
{
return p0.StableDot(p1);
}
typedef Point4<short> Point4s;
typedef Point4<int> Point4i;
typedef Point4<float> Point4f;
typedef Point4<double> Point4d;
/*@}*/
} // end namespace
#endif

4
vcg/space/intersection3.h
View File

@ -339,10 +339,10 @@ namespace vcg {
inline bool IntersectionLinePlane( const Plane3<T> & pl, const Line3<T> & li, Point3<T> & po){
const T epsilon = T(1e-8);
T k = pl.Direction() * li.Direction(); // Compute 'k' factor
T k = pl.Direction().dot(li.Direction()); // Compute 'k' factor
if( (k > -epsilon) && (k < epsilon))
return false;
T r = (pl.Offset() - pl.Direction()*li.Origin())/k; // Compute ray distance
T r = (pl.Offset() - pl.Direction().dot(li.Origin()))/k; // Compute ray distance
po = li.Origin() + li.Direction()*r;
return true;
}

4
vcg/space/line3.h
View File

@ -121,8 +121,8 @@ public:
{ return _ori!=p._ori || _dir!=p._dir; }
/// Projects a point on the line
inline ScalarType Projection( const PointType &p ) const
{ if (NORM) return ScalarType((p-_ori)*_dir);
else return ScalarType((p-_ori)*_dir/_dir.SquaredNorm());
{ if (NORM) return ScalarType((p-_ori).dot(_dir));
else return ScalarType((p-_ori).dot(_dir)/_dir.SquaredNorm());
}
/// returns wheter this type is normalized or not
static bool IsNormalized() {return NORM;};

10
vcg/space/plane3.h
View File

@ -141,7 +141,7 @@ public:
///Project a point on the plane
PointType Projection(const PointType &p) const {
ScalarType k = p * _dir - _offset;
ScalarType k = p.dot(_dir) - _offset;
return p - _dir * k;
}
@ -154,13 +154,13 @@ public:
void Init(const PointType &p0, const PointType &p1, const PointType &p2) {
_dir = (p2 - p0) ^ (p1 - p0);
if(NORM) Normalize();
_offset = p0 * _dir;
_offset = p0.dot(_dir);
}
/// Calculates the plane passing through a point and the normal (Rename this method
inline void Init(const PointType &p0, const PointType &norm) {
_dir = norm;
_offset = p0 * _dir;
_offset = p0.dot(_dir);
}
}; // end class Plane3
@ -169,12 +169,12 @@ typedef Plane3<double> Plane3d;
///Distance plane - point (Move this function to somewhere else)
template<class T> T Distance(const Plane3<T> &plane, const Point3<T> &point) {
return plane.Direction() * point - plane.Offset();
return plane.Direction().dot(point) - plane.Offset();
}
///Distance point-plane (Move this function to somewhere else)
template<class T> T Distance(const Point3<T> &point, const Plane3<T> &plane) {
return plane.Direction() * point - plane.Offset();
return plane.Direction().dot(point) - plane.Offset();
}
} // end namespace

431
vcg/space/point3.h
View File

@ -20,83 +20,26 @@
* for more details. *
* *
****************************************************************************/
/****************************************************************************
History
$Log: not supported by cvs2svn $
Revision 1.26 2006/11/13 13:03:45 ponchio
Added GetBBox in Point3 (declaration) the body of the function is in box3.h
Revision 1.25 2006/10/13 12:59:24 cignoni
Added **explicit** constructor from three coords of a different scalartype
Revision 1.24 2006/09/28 13:37:35 m_di_benedetto
added non const * V()
Revision 1.23 2005/11/09 16:11:55 cignoni
Added Abs and LowClampToZero
Revision 1.22 2005/09/14 14:09:21 m_di_benedetto
Added specialized Convert() for the same scalar type.
Revision 1.21 2005/05/06 14:45:33 spinelli
cambiato parentesi nel costruttore di GetUV per rendere compatibile tale costruttore con MVC e borland
Revision 1.20 2005/04/27 16:05:19 callieri
line 466, added parentesis on default value creator getUV [borland]
Revision 1.19 2004/11/09 15:49:07 ganovelli
added GetUV
Revision 1.18 2004/10/13 12:45:51 cignoni
Better Doxygen documentation
Revision 1.17 2004/09/10 14:01:40 cignoni
Added polar to cartesian
Revision 1.16 2004/03/21 17:14:56 ponchio
Added a math::
Revision 1.15 2004/03/07 22:45:32 cignoni
Moved quality and normal functions to the triangle class.
Revision 1.14 2004/03/05 17:55:01 tarini
errorino: upper case in Zero()
Revision 1.13 2004/03/03 14:22:48 cignoni
Yet against cr lf mismatch
Revision 1.12 2004/02/23 23:42:26 cignoni
Translated comments, removed unusued stuff. corrected linefeed/cr
Revision 1.11 2004/02/19 16:12:28 cignoni
cr lf mismatch 2
Revision 1.10 2004/02/19 16:06:24 cignoni
cr lf mismatch
Revision 1.8 2004/02/19 15:13:40 cignoni
corrected sqrt and added doxygen groups
Revision 1.7 2004/02/17 02:08:47 cignoni
Di prova...
Revision 1.6 2004/02/15 23:35:47 cignoni
Cambiato nome type template in accordo alla styleguide
Revision 1.5 2004/02/10 01:07:15 cignoni
Edited Comments and GPL license
Revision 1.4 2004/02/09 13:48:02 cignoni
Edited doxygen comments
****************************************************************************/
#ifndef VCG_USE_EIGEN
#include "deprecated_point3.h"
#else
#ifndef __VCGLIB_POINT3
#define __VCGLIB_POINT3
#include <assert.h>
#include "../math/eigen.h"
#include <vcg/math/base.h>
namespace vcg{
template<class Scalar> class Point3;
}
namespace Eigen{
template<typename Scalar>
struct ei_traits<vcg::Point3<Scalar> > : ei_traits<Eigen::Matrix<Scalar,3,1> > {};
}
namespace vcg {
/** \addtogroup space */
@ -108,14 +51,22 @@ namespace vcg {
*/
template <class T> class Box3;
template <class P3ScalarType> class Point3
template <class _Scalar> class Point3 : public Eigen::Matrix<_Scalar,3,1>
{
protected:
/// The only data member. Hidden to user.
P3ScalarType _v[3];
typedef Eigen::Matrix<_Scalar,3,1> _Base;
using _Base::coeff;
using _Base::coeffRef;
using _Base::setZero;
using _Base::data;
using _Base::V;
public:
typedef P3ScalarType ScalarType;
_EIGEN_GENERIC_PUBLIC_INTERFACE(Point3,_Base);
typedef Scalar ScalarType;
VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Point3)
enum {Dimension = 3};
@ -125,63 +76,40 @@ public:
No casting operators have been introduced to avoid automatic unattended (and costly) conversion between different point types
**/
inline Point3 () { }
inline Point3 ( const P3ScalarType nx, const P3ScalarType ny, const P3ScalarType nz )
{
_v[0] = nx;
_v[1] = ny;
_v[2] = nz;
}
inline Point3 ( Point3 const & p )
{
_v[0]= p._v[0];
_v[1]= p._v[1];
_v[2]= p._v[2];
}
inline Point3 ( const P3ScalarType nv[3] )
{
_v[0] = nv[0];
_v[1] = nv[1];
_v[2] = nv[2];
}
inline Point3 & operator =( Point3 const & p )
{
_v[0]= p._v[0]; _v[1]= p._v[1]; _v[2]= p._v[2];
return *this;
}
inline void SetZero()
{
_v[0] = 0;
_v[1] = 0;
_v[2] = 0;
}
inline Point3 () {}
inline Point3 ( const Scalar nx, const Scalar ny, const Scalar nz ) : Base(nx,ny,nz) {}
inline Point3 ( Point3 const & p ) : Base(p) {}
inline Point3 ( const Scalar nv[3] ) : Base(nv) {}
template<typename OtherDerived>
inline Point3(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 P3ScalarType Ext( const int i ) const
inline Scalar Ext( const int i ) const
{
if(i>=0 && i<=2) return _v[i];
if(i>=0 && i<=2) return data()[i];
else return 0;
}
template <class Q>
inline void Import( const Point3<Q> & b )
template<class OtherDerived>
inline void Import( const Eigen::MatrixBase<OtherDerived>& b )
{
_v[0] = P3ScalarType(b[0]);
_v[1] = P3ScalarType(b[1]);
_v[2] = P3ScalarType(b[2]);
EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(OtherDerived,3);
data()[0] = Scalar(b[0]);
data()[1] = Scalar(b[1]);
data()[2] = Scalar(b[2]);
}
template <class Q>
static inline Point3 Construct( const Point3<Q> & b )
{
return Point3(P3ScalarType(b[0]),P3ScalarType(b[1]),P3ScalarType(b[2]));
return Point3(Scalar(b[0]),Scalar(b[1]),Scalar(b[2]));
}
template <class Q>
static inline Point3 Construct( const Q & P0, const Q & P1, const Q & P2)
{
return Point3(P3ScalarType(P0),P3ScalarType(P1),P3ScalarType(P2));
return Point3(Scalar(P0),Scalar(P1),Scalar(P2));
}
static inline Point3 Construct( const Point3<ScalarType> & b )
@ -196,39 +124,21 @@ public:
/** @name Data Access.
access to data is done by overloading of [] or explicit naming of coords (x,y,z)**/
inline P3ScalarType & operator [] ( const int i )
inline const Scalar &X() const { return data()[0]; }
inline const Scalar &Y() const { return data()[1]; }
inline const Scalar &Z() const { return data()[2]; }
inline Scalar &X() { return data()[0]; }
inline Scalar &Y() { return data()[1]; }
inline Scalar &Z() { return data()[2]; }
inline Scalar & V( const int i )
{
assert(i>=0 && i<3);
return _v[i];
return data()[i];
}
inline const P3ScalarType & operator [] ( const int i ) const
inline const Scalar & V( const int i ) const
{
assert(i>=0 && i<3);
return _v[i];
}
inline const P3ScalarType &X() const { return _v[0]; }
inline const P3ScalarType &Y() const { return _v[1]; }
inline const P3ScalarType &Z() const { return _v[2]; }
inline P3ScalarType &X() { return _v[0]; }
inline P3ScalarType &Y() { return _v[1]; }
inline P3ScalarType &Z() { return _v[2]; }
inline const P3ScalarType * V() const
{
return _v;
}
inline P3ScalarType * V()
{
return _v;
}
inline P3ScalarType & V( const int i )
{
assert(i>=0 && i<3);
return _v[i];
}
inline const P3ScalarType & V( const int i ) const
{
assert(i>=0 && i<3);
return _v[i];
return data()[i];
}
//@}
//@{
@ -237,96 +147,19 @@ public:
Note
**/
inline Point3 operator + ( Point3 const & p) const
// Scalatura differenziata
inline Point3 & Scale( const Scalar sx, const Scalar sy, const Scalar sz )
{
return Point3<P3ScalarType>( _v[0]+p._v[0], _v[1]+p._v[1], _v[2]+p._v[2] );
}
inline Point3 operator - ( Point3 const & p) const
{
return Point3<P3ScalarType>( _v[0]-p._v[0], _v[1]-p._v[1], _v[2]-p._v[2] );
}
inline Point3 operator * ( const P3ScalarType s ) const
{
return Point3<P3ScalarType>( _v[0]*s, _v[1]*s, _v[2]*s );
}
inline Point3 operator / ( const P3ScalarType s ) const
{
return Point3<P3ScalarType>( _v[0]/s, _v[1]/s, _v[2]/s );
}
/// Dot product
inline P3ScalarType operator * ( Point3 const & p ) const
{
return ( _v[0]*p._v[0] + _v[1]*p._v[1] + _v[2]*p._v[2] );
}
/// Cross product
inline Point3 operator ^ ( Point3 const & p ) const
{
return Point3 <P3ScalarType>
(
_v[1]*p._v[2] - _v[2]*p._v[1],
_v[2]*p._v[0] - _v[0]*p._v[2],
_v[0]*p._v[1] - _v[1]*p._v[0]
);
}
inline Point3 & operator += ( Point3 const & p)
{
_v[0] += p._v[0];
_v[1] += p._v[1];
_v[2] += p._v[2];
return *this;
}
inline Point3 & operator -= ( Point3 const & p)
{
_v[0] -= p._v[0];
_v[1] -= p._v[1];
_v[2] -= p._v[2];
return *this;
}
inline Point3 & operator *= ( const P3ScalarType s )
{
_v[0] *= s;
_v[1] *= s;
_v[2] *= s;
return *this;
}
inline Point3 & operator /= ( const P3ScalarType s )
{
_v[0] /= s;
_v[1] /= s;
_v[2] /= s;
return *this;
}
// Norme
inline P3ScalarType Norm() const
{
return math::Sqrt( _v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] );
}
inline P3ScalarType SquaredNorm() const
{
return ( _v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] );
}
// Scalatura differenziata
inline Point3 & Scale( const P3ScalarType sx, const P3ScalarType sy, const P3ScalarType sz )
{
_v[0] *= sx;
_v[1] *= sy;
_v[2] *= sz;
data()[0] *= sx;
data()[1] *= sy;
data()[2] *= sz;
return *this;
}
inline Point3 & Scale( const Point3 & p )
{
_v[0] *= p._v[0];
_v[1] *= p._v[1];
_v[2] *= p._v[2];
return *this;
}
// Normalizzazione
inline Point3 & Normalize()
{
P3ScalarType n = math::Sqrt(_v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2]);
if(n>0.0) { _v[0] /= n; _v[1] /= n; _v[2] /= n; }
data()[0] *= p.data()[0];
data()[1] *= p.data()[1];
data()[2] *= p.data()[2];
return *this;
}
@ -340,11 +173,11 @@ public:
* 0 and 180 degrees we opt for the elevation angle to obtain a more
* intuitive spherical coordinate system.
*/
void ToPolar(P3ScalarType &ro, P3ScalarType &theta, P3ScalarType &phi) const
void ToPolar(Scalar &ro, Scalar &theta, Scalar &phi) const
{
ro = Norm();
theta = (P3ScalarType)atan2(_v[2], _v[0]);
phi = (P3ScalarType)asin(_v[1]/ro);
ro = this->norm();
theta = (Scalar)atan2(data()[2], data()[0]);
phi = (Scalar)asin(data()[1]/ro);
}
/**
@ -357,14 +190,14 @@ public:
* 0 and 180 degrees, we opt for the elevation angle to obtain a more
* intuitive spherical coordinate system.
*/
void FromPolar(const P3ScalarType &ro, const P3ScalarType &theta, const P3ScalarType &phi)
void FromPolar(const Scalar &ro, const Scalar &theta, const Scalar &phi)
{
_v[0]= ro*cos(theta)*cos(phi);
_v[1]= ro*sin(phi);
_v[2]= ro*sin(theta)*cos(phi);
data()[0]= ro*cos(theta)*cos(phi);
data()[1]= ro*sin(phi);
data()[2]= ro*sin(theta)*cos(phi);
}
Box3<P3ScalarType> GetBBox(Box3<P3ScalarType> &bb) const;
Box3<_Scalar> GetBBox(Box3<_Scalar> &bb) const;
//@}
//@{
@ -372,113 +205,64 @@ public:
Note that the reverse z prioritized ordering, useful in many situations.
**/
inline bool operator == ( Point3 const & p ) const
{
return _v[0]==p._v[0] && _v[1]==p._v[1] && _v[2]==p._v[2];
}
inline bool operator != ( Point3 const & p ) const
{
return _v[0]!=p._v[0] || _v[1]!=p._v[1] || _v[2]!=p._v[2];
}
inline bool operator < ( Point3 const & p ) const
{
return (_v[2]!=p._v[2])?(_v[2]<p._v[2]):
(_v[1]!=p._v[1])?(_v[1]<p._v[1]):
(_v[0]<p._v[0]);
return (data()[2]!=p.data()[2])?(data()[2]<p.data()[2]):
(data()[1]!=p.data()[1])?(data()[1]<p.data()[1]):
(data()[0]<p.data()[0]);
}
inline bool operator > ( Point3 const & p ) const
{
return (_v[2]!=p._v[2])?(_v[2]>p._v[2]):
(_v[1]!=p._v[1])?(_v[1]>p._v[1]):
(_v[0]>p._v[0]);
return (data()[2]!=p.data()[2])?(data()[2]>p.data()[2]):
(data()[1]!=p.data()[1])?(data()[1]>p.data()[1]):
(data()[0]>p.data()[0]);
}
inline bool operator <= ( Point3 const & p ) const
{
return (_v[2]!=p._v[2])?(_v[2]< p._v[2]):
(_v[1]!=p._v[1])?(_v[1]< p._v[1]):
(_v[0]<=p._v[0]);
return (data()[2]!=p.data()[2])?(data()[2]< p.data()[2]):
(data()[1]!=p.data()[1])?(data()[1]< p.data()[1]):
(data()[0]<=p.data()[0]);
}
inline bool operator >= ( Point3 const & p ) const
{
return (_v[2]!=p._v[2])?(_v[2]> p._v[2]):
(_v[1]!=p._v[1])?(_v[1]> p._v[1]):
(_v[0]>=p._v[0]);
}
inline Point3 operator - () const
{
return Point3<P3ScalarType> ( -_v[0], -_v[1], -_v[2] );
return (data()[2]!=p.data()[2])?(data()[2]> p.data()[2]):
(data()[1]!=p.data()[1])?(data()[1]> p.data()[1]):
(data()[0]>=p.data()[0]);
}
//@}
}; // end class definition
template <class P3ScalarType>
inline P3ScalarType Angle( Point3<P3ScalarType> const & p1, Point3<P3ScalarType> const & p2 )
{
P3ScalarType w = p1.Norm()*p2.Norm();
if(w==0) return -1;
P3ScalarType t = (p1*p2)/w;
if(t>1) t = 1;
else if(t<-1) t = -1;
return (P3ScalarType) acos(t);
}
// versione uguale alla precedente ma che assume che i due vettori sono unitari
template <class P3ScalarType>
inline P3ScalarType AngleN( Point3<P3ScalarType> const & p1, Point3<P3ScalarType> const & p2 )
template <class Scalar>
inline Scalar AngleN( Point3<Scalar> const & p1, Point3<Scalar> const & p2 )
{
P3ScalarType w = p1*p2;
Scalar w = p1*p2;
if(w>1)
w = 1;
else if(w<-1)
w=-1;
return (P3ScalarType) acos(w);
return (Scalar) acos(w);
}
template <class P3ScalarType>
inline P3ScalarType Norm( Point3<P3ScalarType> const & p )
{
return p.Norm();
}
template <class P3ScalarType>
inline P3ScalarType SquaredNorm( Point3<P3ScalarType> const & p )
{
return p.SquaredNorm();
}
template <class P3ScalarType>
inline Point3<P3ScalarType> & Normalize( Point3<P3ScalarType> & p )
template <class Scalar>
inline Point3<Scalar> & Normalize( Point3<Scalar> & p )
{
p.Normalize();
return p;
}
template <class P3ScalarType>
inline P3ScalarType Distance( Point3<P3ScalarType> const & p1,Point3<P3ScalarType> const & p2 )
{
return (p1-p2).Norm();
}
template <class P3ScalarType>
inline P3ScalarType SquaredDistance( Point3<P3ScalarType> const & p1,Point3<P3ScalarType> const & p2 )
{
return (p1-p2).SquaredNorm();
}
// Dot product preciso numericamente (solo double!!)
// Implementazione: si sommano i prodotti per ordine di esponente
// (prima le piu' grandi)
template<class P3ScalarType>
double stable_dot ( Point3<P3ScalarType> const & p0, Point3<P3ScalarType> const & p1 )
template<class Scalar>
double stable_dot ( Point3<Scalar> const & p0, Point3<Scalar> const & p1 )
{
P3ScalarType k0 = p0._v[0]*p1._v[0];
P3ScalarType k1 = p0._v[1]*p1._v[1];
P3ScalarType k2 = p0._v[2]*p1._v[2];
Scalar k0 = p0.data()[0]*p1.data()[0];
Scalar k1 = p0.data()[1]*p1.data()[1];
Scalar k2 = p0.data()[2]*p1.data()[2];
int exp0,exp1,exp2;
@ -505,14 +289,14 @@ double stable_dot ( Point3<P3ScalarType> const & p0, Point3<P3ScalarType> const
/// Point(p) Edge(v1-v2) dist, q is the point in v1-v2 with min dist
template<class P3ScalarType>
P3ScalarType PSDist( const Point3<P3ScalarType> & p,
const Point3<P3ScalarType> & v1,
const Point3<P3ScalarType> & v2,
Point3<P3ScalarType> & q )
template<class Scalar>
Scalar PSDist( const Point3<Scalar> & p,
const Point3<Scalar> & v1,
const Point3<Scalar> & v2,
Point3<Scalar> & q )
{
Point3<P3ScalarType> e = v2-v1;
P3ScalarType t = ((p-v1)*e)/e.SquaredNorm();
Point3<Scalar> e = v2-v1;
Scalar t = ((p-v1).dot(e))/e.SquaredNorm();
if(t<0) t = 0;
else if(t>1) t = 1;
q = v1+e*t;
@ -520,8 +304,8 @@ P3ScalarType PSDist( const Point3<P3ScalarType> & p,
}
template <class P3ScalarType>
void GetUV( Point3<P3ScalarType> &n,Point3<P3ScalarType> &u, Point3<P3ScalarType> &v, Point3<P3ScalarType> up=(Point3<P3ScalarType>(0,1,0)) )
template <class Scalar>
void GetUV( Point3<Scalar> &n,Point3<Scalar> &u, Point3<Scalar> &v, Point3<Scalar> up=(Point3<Scalar>(0,1,0)) )
{
n.Normalize();
const double LocEps=double(1e-7);
@ -530,18 +314,18 @@ void GetUV( Point3<P3ScalarType> &n,Point3<P3ScalarType> &u, Point3<P3ScalarType
if(len < LocEps)
{
if(fabs(n[0])<fabs(n[1])){
if(fabs(n[0])<fabs(n[2])) up=Point3<P3ScalarType>(1,0,0); // x is the min
else up=Point3<P3ScalarType>(0,0,1); // z is the min
if(fabs(n[0])<fabs(n[2])) up=Point3<Scalar>(1,0,0); // x is the min
else up=Point3<Scalar>(0,0,1); // z is the min
}else {
if(fabs(n[1])<fabs(n[2])) up=Point3<P3ScalarType>(0,1,0); // y is the min
else up=Point3<P3ScalarType>(0,0,1); // z is the min
if(fabs(n[1])<fabs(n[2])) up=Point3<Scalar>(0,1,0); // y is the min
else up=Point3<Scalar>(0,0,1); // z is the min
}
u=n^up;
}
u.Normalize();
v=n^u;
v.Normalize();
Point3<P3ScalarType> uv=u^v;
Point3<Scalar> uv=u^v;
}
@ -566,3 +350,4 @@ typedef Point3<double> Point3d;
#endif
#endif

320
vcg/space/point4.h
View File

@ -20,48 +20,24 @@
* for more details. *
* *
****************************************************************************/
/****************************************************************************
History
$Log: not supported by cvs2svn $
Revision 1.12 2006/06/21 11:06:16 ganovelli
changed return type of Zero() (to void)
Revision 1.11 2005/04/13 09:40:30 ponchio
Including math/bash.h
Revision 1.10 2005/03/18 16:34:42 fiorin
minor changes to comply gcc compiler
Revision 1.9 2005/01/21 18:02:11 ponchio
Removed dependence from matrix44 and changed VectProd
Revision 1.8 2005/01/12 11:25:02 ganovelli
added Dimension
Revision 1.7 2004/10/11 17:46:11 ganovelli
added definition of vector product (not implemented)
Revision 1.6 2004/05/10 11:16:19 ganovelli
include assert.h added
Revision 1.5 2004/03/31 10:09:58 cignoni
missing return value in zero()
Revision 1.4 2004/03/11 17:17:49 tarini
added commets (doxy), uniformed with new style, now using math::, ...
added HomoNormalize(), Zero()... remade StableDot() (hand made sort).
Revision 1.1 2004/02/10 01:11:28 cignoni
Edited Comments and GPL license
****************************************************************************/
#ifndef VCG_USE_EIGEN
#include "deprecated_point4.h"
#else
#ifndef __VCGLIB_POINT4
#define __VCGLIB_POINT4
#include <assert.h>
#include <vcg/math/base.h>
#include "../math/eigen.h"
namespace vcg{
template<class Scalar> class Point4;
}
namespace Eigen{
template<typename Scalar>
struct ei_traits<vcg::Point4<Scalar> > : ei_traits<Eigen::Matrix<Scalar,4,1> > {};
}
namespace vcg {
/** \addtogroup space */
@ -72,224 +48,119 @@ namespace vcg {
All the usual operator (* + - ...) are defined.
*/
template <class T> class Point4
template <class T> class Point4 : public Eigen::Matrix<T,4,1>
{
public:
/// The only data member. Hidden to user.
T _v[4];
typedef Eigen::Matrix<T,4,1> _Base;
using _Base::coeff;
using _Base::coeffRef;
using _Base::setZero;
using _Base::data;
using _Base::V;
public:
typedef T ScalarType;
_EIGEN_GENERIC_PUBLIC_INTERFACE(Point4,_Base);
typedef Scalar ScalarType;
VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Point4)
enum {Dimension = 4};
//@{
inline Point4() : Base() {}
inline Point4( const T nx, const T ny, const T nz , const T nw ) : Base(nx,ny,nz,nw) {}
inline Point4(const T p[4]) : Base(p) {}
inline Point4(const Point4& p) : Base(p) {}
template<typename OtherDerived>
inline Point4(const Eigen::MatrixBase<OtherDerived>& other) : Base(other) {}
/** @name Standard Constructors and Initializers
No casting operators have been introduced to avoid automatic unattended (and costly) conversion between different point types
**/
inline Point4 () { }
inline Point4 ( const T nx, const T ny, const T nz , const T nw )
{
_v[0] = nx; _v[1] = ny; _v[2] = nz; _v[3] = nw;
}
inline Point4 ( const T p[4] )
{
_v[0] = p[0]; _v[1]= p[1]; _v[2] = p[2]; _v[3]= p[3];
}
inline Point4 ( const Point4 & p )
{
_v[0]= p._v[0]; _v[1]= p._v[1]; _v[2]= p._v[2]; _v[3]= p._v[3];
}
inline void Zero()
{
_v[0] = _v[1] = _v[2] = _v[3]= 0;
}
template <class Q>
/// importer from different Point4 types
inline void Import( const Point4<Q> & b )
{
_v[0] = T(b[0]); _v[1] = T(b[1]);
_v[2] = T(b[2]);
_v[3] = T(b[3]);
}
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 Point4(T(b[0]),T(b[1]),T(b[2]),T(b[3]));
}
template <class Q>
static inline Point4 Construct( const Point4<Q> & b ) { return b.template cast<T>(); }
//@}
//@{
/** @name Data Access.
access to data is done by overloading of [] or explicit naming of coords (x,y,z,w)
**/
inline const T & operator [] ( const int i ) const
{
assert(i>=0 && i<4);
return _v[i];
}
inline T & operator [] ( const int i )
{
assert(i>=0 && i<4);
return _v[i];
}
inline T &X() {return _v[0];}
inline T &Y() {return _v[1];}
inline T &Z() {return _v[2];}
inline T &W() {return _v[3];}
inline T const * V() const
{
return _v;
}
inline T &X() {return Base::x();}
inline T &Y() {return Base::y();}
inline T &Z() {return Base::z();}
inline T &W() {return Base::w();}
inline const T & V ( const int i ) const
{
assert(i>=0 && i<4);
return _v[i];
return data()[i];
}
inline T & V ( const int i )
{
assert(i>=0 && i<4);
return _v[i];
return data()[i];
}
/// Padding function: give a default 0 value to all the elements that are not in the [0..2] range.
/// Useful for managing in a consistent way object that could have point2 / point3 / point4
/// Padding function: give a default 0 value to all the elements that are not in the [0..2] range.
/// Useful for managing in a consistent way object that could have point2 / point3 / point4
inline T Ext( const int i ) const
{
if(i>=0 && i<=3) return _v[i];
if(i>=0 && i<=3) return data()[i];
else return 0;
}
//@}
//@{
/** @name Linear operators and the likes
**/
inline Point4 operator + ( const Point4 & p) const
{
return Point4( _v[0]+p._v[0], _v[1]+p._v[1], _v[2]+p._v[2], _v[3]+p._v[3] );
}
inline Point4 operator - ( const Point4 & p) const
{
return Point4( _v[0]-p._v[0], _v[1]-p._v[1], _v[2]-p._v[2], _v[3]-p._v[3] );
}
inline Point4 operator * ( const T s ) const
{
return Point4( _v[0]*s, _v[1]*s, _v[2]*s, _v[3]*s );
}
inline Point4 operator / ( const T s ) const
{
return Point4( _v[0]/s, _v[1]/s, _v[2]/s, _v[3]/s );
}
inline Point4 & operator += ( const Point4 & p)
{
_v[0] += p._v[0]; _v[1] += p._v[1]; _v[2] += p._v[2]; _v[3] += p._v[3];
return *this;
}
inline Point4 & operator -= ( const Point4 & p )
{
_v[0] -= p._v[0]; _v[1] -= p._v[1]; _v[2] -= p._v[2]; _v[3] -= p._v[3];
return *this;
}
inline Point4 & operator *= ( const T s )
{
_v[0] *= s; _v[1] *= s; _v[2] *= s; _v[3] *= s;
return *this;
}
inline Point4 & operator /= ( const T s )
{
_v[0] /= s; _v[1] /= s; _v[2] /= s; _v[3] /= s;
return *this;
}
inline Point4 operator - () const
{
return Point4( -_v[0], -_v[1], -_v[2], -_v[3] );
}
inline Point4 VectProd ( const Point4 &x, const Point4 &z ) const
{
Point4 res;
const Point4 &y = *this;
res[0] = y[1]*x[2]*z[3]-y[1]*x[3]*z[2]-x[1]*y[2]*z[3]+
x[1]*y[3]*z[2]+z[1]*y[2]*x[3]-z[1]*y[3]*x[2];
res[1] = y[0]*x[3]*z[2]-z[0]*y[2]*x[3]-y[0]*x[2]*
z[3]+z[0]*y[3]*x[2]+x[0]*y[2]*z[3]-x[0]*y[3]*z[2];
res[0] = y[1]*x[2]*z[3]-y[1]*x[3]*z[2]-x[1]*y[2]*z[3]+
x[1]*y[3]*z[2]+z[1]*y[2]*x[3]-z[1]*y[3]*x[2];
res[1] = y[0]*x[3]*z[2]-z[0]*y[2]*x[3]-y[0]*x[2]*
z[3]+z[0]*y[3]*x[2]+x[0]*y[2]*z[3]-x[0]*y[3]*z[2];
res[2] = -y[0]*z[1]*x[3]+x[0]*z[1]*y[3]+y[0]*x[1]*
z[3]-x[0]*y[1]*z[3]-z[0]*x[1]*y[3]+z[0]*y[1]*x[3];
z[3]-x[0]*y[1]*z[3]-z[0]*x[1]*y[3]+z[0]*y[1]*x[3];
res[3] = -z[0]*y[1]*x[2]-y[0]*x[1]*z[2]+x[0]*y[1]*
z[2]+y[0]*z[1]*x[2]-x[0]*z[1]*y[2]+z[0]*x[1]*y[2];
z[2]+y[0]*z[1]*x[2]-x[0]*z[1]*y[2]+z[0]*x[1]*y[2];
return res;
}
//@}
//@{
/** @name Norms and normalizations
**/
/// Euclidian normal
inline T Norm() const
{
return math::Sqrt( _v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] + _v[3]*_v[3] );
}
/// Squared euclidian normal
inline T SquaredNorm() const
{
return _v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] + _v[3]*_v[3];
}
/// Euclidian normalization
inline Point4 & Normalize()
{
T n = sqrt(_v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] + _v[3]*_v[3] );
if(n>0.0) { _v[0] /= n; _v[1] /= n; _v[2] /= n; _v[3] /= n; }
return *this;
}
/// Homogeneous normalization (division by W)
inline Point4 & HomoNormalize(){
if (_v[3]!=0.0) { _v[0] /= _v[3]; _v[1] /= _v[3]; _v[2] /= _v[3]; _v[3]=1.0; }
inline Point4 & HomoNormalize() {
if (data()[3]!=0.0) { Base::template start<3>() /= coeff(3); coeffRef(3) = 1.0; }
return *this;
};
//@}
}
//@{
/** @name Comparison operators (lexicographical order)
**/
inline bool operator == ( const Point4& p ) const
{
return _v[0]==p._v[0] && _v[1]==p._v[1] && _v[2]==p._v[2] && _v[3]==p._v[3];
}
inline bool operator != ( const Point4 & p ) const
{
return _v[0]!=p._v[0] || _v[1]!=p._v[1] || _v[2]!=p._v[2] || _v[3]!=p._v[3];
}
inline bool operator < ( Point4 const & p ) const
{
return (_v[3]!=p._v[3])?(_v[3]<p._v[3]):
(_v[2]!=p._v[2])?(_v[2]<p._v[2]):
(_v[1]!=p._v[1])?(_v[1]<p._v[1]):
(_v[0]<p._v[0]);
return (data()[3]!=p.data()[3])?(data()[3]<p.data()[3]):
(data()[2]!=p.data()[2])?(data()[2]<p.data()[2]):
(data()[1]!=p.data()[1])?(data()[1]<p.data()[1]):
(data()[0]<p.data()[0]);
}
inline bool operator > ( const Point4 & p ) const
{
return (_v[3]!=p._v[3])?(_v[3]>p._v[3]):
(_v[2]!=p._v[2])?(_v[2]>p._v[2]):
(_v[1]!=p._v[1])?(_v[1]>p._v[1]):
(_v[0]>p._v[0]);
return (data()[3]!=p.data()[3])?(data()[3]>p.data()[3]):
(data()[2]!=p.data()[2])?(data()[2]>p.data()[2]):
(data()[1]!=p.data()[1])?(data()[1]>p.data()[1]):
(data()[0]>p.data()[0]);
}
inline bool operator <= ( const Point4 & p ) const
{
return (_v[3]!=p._v[3])?(_v[3]< p._v[3]):
(_v[2]!=p._v[2])?(_v[2]< p._v[2]):
(_v[1]!=p._v[1])?(_v[1]< p._v[1]):
(_v[0]<=p._v[0]);
return (data()[3]!=p.data()[3])?(data()[3]< p.data()[3]):
(data()[2]!=p.data()[2])?(data()[2]< p.data()[2]):
(data()[1]!=p.data()[1])?(data()[1]< p.data()[1]):
(data()[0]<=p.data()[0]);
}
inline bool operator >= ( const Point4 & p ) const
{
return (_v[3]!=p._v[3])?(_v[3]> p._v[3]):
(_v[2]!=p._v[2])?(_v[2]> p._v[2]):
(_v[1]!=p._v[1])?(_v[1]> p._v[1]):
(_v[0]>=p._v[0]);
return (data()[3]!=p.data()[3])?(data()[3]> p.data()[3]):
(data()[2]!=p.data()[2])?(data()[2]> p.data()[2]):
(data()[1]!=p.data()[1])?(data()[1]> p.data()[1]):
(data()[0]>=p.data()[0]);
}
//@}
@ -297,22 +168,16 @@ public:
/** @name Dot products
**/
// dot product
inline T operator * ( const Point4 & p ) const
{
return _v[0]*p._v[0] + _v[1]*p._v[1] + _v[2]*p._v[2] + _v[3]*p._v[3];
}
inline Point4 operator ^ ( const Point4& p ) const
{
assert(0);// not defined by two vectors (only put for metaprogramming)
assert(0 && "not defined by two vectors (only put for metaprogramming)");
return Point4();
}
/// slower version, more stable (double precision only)
T StableDot ( const Point4<T> & p ) const
{
T k0=_v[0]*p._v[0], k1=_v[1]*p._v[1], k2=_v[2]*p._v[2], k3=_v[3]*p._v[3];
T k0=data()[0]*p.data()[0], k1=data()[1]*p.data()[1], k2=data()[2]*p.data()[2], k3=data()[3]*p.data()[3];
int exp0,exp1,exp2,exp3;
frexp( double(k0), &exp0 );frexp( double(k1), &exp1 );
@ -331,39 +196,6 @@ public:
}; // end class definition
template <class T>
T Angle( const Point4<T>& p1, const Point4<T> & p2 )
{
T w = p1.Norm()*p2.Norm();
if(w==0) return -1;
T t = (p1*p2)/w;
if(t>1) t=1;
return T( math::Acos(t) );
}
template <class T>
inline T Norm( const Point4<T> & p )
{
return p.Norm();
}
template <class T>
inline T SquaredNorm( const Point4<T> & p )
{
return p.SquaredNorm();
}
template <class T>
inline T Distance( const Point4<T> & p1, const Point4<T> & p2 )
{
return Norm(p1-p2);
}
template <class T>
inline T SquaredDistance( const Point4<T> & p1, const Point4<T> & p2 )
{
return SquaredNorm(p1-p2);
}
/// slower version of dot product, more stable (double precision only)
template<class T>
@ -380,3 +212,5 @@ typedef Point4<double> Point4d;
/*@}*/
} // end namespace
#endif
#endif

4
vcg/space/ray3.h
View File

@ -112,8 +112,8 @@ public:
{ return _ori!=p._ori || _dir!=p._dir; }
/// Projects a point on the ray
inline ScalarType Projection( const PointType &p ) const
{ if (NORM) return ScalarType((p-_ori)*_dir);
else return ScalarType((p-_ori)*_dir/_dir.SquaredNorm());
{ if (NORM) return ScalarType((p-_ori).dot(_dir));
else return ScalarType((p-_ori).dot(_dir)/_dir.SquaredNorm());
}
/// returns wheter this type is normalized or not
static bool IsNormalized() {return NORM;};

38
wrap/gl/math.h
View File

@ -72,40 +72,32 @@ namespace vcg {
inline void glMultMatrixE(const Matrix44f &matrix) {
//glMultMatrixf((const GLfloat *)(matrix[0]));
if(glMultTransposeMatrixf) glMultTransposeMatrixf((const GLfloat *)(matrix[0]));
if(glMultTransposeMatrixf) glMultTransposeMatrixf((const GLfloat *)(matrix.V()));
else {
Matrix44f tmp(matrix);
Transpose(tmp);
glMultMatrixf((const GLfloat *)(tmp[0]));
glMultMatrixf((const GLfloat *)(matrix.transpose().eval().V()));
}
}
inline void glMultMatrix(const Matrix44f &matrix) {
Matrix44f tmp(matrix);
Transpose(tmp);
glMultMatrixf((const GLfloat *)(tmp[0]));
glMultMatrixf((const GLfloat *)(matrix.transpose().eval().V()));
}
inline void glMultMatrixE(const Matrix44d &matrix) {
if(glMultTransposeMatrixd) glMultTransposeMatrixd((const GLdouble *)(matrix[0]));
if(glMultTransposeMatrixd) glMultTransposeMatrixd((const GLdouble *)(matrix.V()));
else {
Matrix44d tmp(matrix);
Transpose(tmp);
glMultMatrixd((const GLdouble *)(tmp[0]));
glMultMatrixd((const GLdouble *)(matrix.transpose().eval().V()));
}
}
inline void glMultMatrix(const Matrix44d &matrix) {
Matrix44d tmp(matrix);
Transpose(tmp);
glMultMatrixd((const GLdouble *)(tmp[0]));
glMultMatrixd((const GLdouble *)(matrix.transpose().eval().V()));
}
inline void glMultMatrixDirect(const Matrix44f &matrix) {
glMultMatrixf((const GLfloat *)(matrix[0]));
glMultMatrixf((const GLfloat *)(matrix.V()));
}
inline void glMultMatrixDirect(const Matrix44d &matrix) {
glMultMatrixd((const GLdouble *)(matrix[0]));
glMultMatrixd((const GLdouble *)(matrix.V()));
}
inline void glMultMatrix(const Similarityf &s) {
@ -128,21 +120,23 @@ inline void glMultMatrix(const Similarityd &s) {
}
inline void glGetv(const GLenum pname, Matrix44f & m){
glGetFloatv(pname,&m[0][0]);
Transpose(m);
Matrix44f tmp;
glGetFloatv(pname,tmp.V());
m = tmp.transpose();
}
inline void glGetv(const GLenum pname, Matrix44d & m){
glGetDoublev(pname,&m[0][0]);
Transpose(m);
Matrix44d tmp;
glGetDoublev(pname,tmp.V());
m = tmp.transpose();
}
inline void glGetDirectv(const GLenum pname, Matrix44f & m){
glGetFloatv(pname,&m[0][0]);
glGetFloatv(pname,m.V());
}
inline void glGetDirecv(const GLenum pname, Matrix44d & m){
glGetDoublev(pname,&m[0][0]);
glGetDoublev(pname,m.V());
}

16
wrap/gl/space.h
View File

@ -243,9 +243,8 @@ inline void glBoxClip(const Box3<T> & b)
*(Point3<T>*)&m[2][0] = *(Point3<T>*)&v[0];m[2][3]=0;
*(Point3<T>*)&m[3][0] = *(Point3<T>*)&c1[0];m[3][3]=1;
vcg::Transpose(m);
glPushMatrix();
glMultMatrix(m);
glMultMatrix(m.transpose());
glBegin(GL_QUADS);
glNormal(Point3<T>(0,1,0));
@ -279,5 +278,18 @@ template <class TetraType>
glTriangle3(Triangle3<typename TetraType::ScalarType>(c.P(1),c.P(0),c.P(3)));
}
#ifdef VCG_USE_EIGEN
#define _WRAP_EIGEN_XPR(FUNC) template<typename Derived> \
inline void FUNC(const Eigen::MatrixBase<Derived>& p) { FUNC(p.eval()); }
_WRAP_EIGEN_XPR(glVertex)
_WRAP_EIGEN_XPR(glNormal)
_WRAP_EIGEN_XPR(glTexCoord)
_WRAP_EIGEN_XPR(glTranslate)
_WRAP_EIGEN_XPR(glScale)
#endif
}//namespace
#endif

8
wrap/gui/trackmode.cpp
View File

@ -486,8 +486,8 @@ int PathMode::Verse(Point3f reference_point,Point3f current_point,Point3f prev_p
prev_dir.Normalize();
next_dir.Normalize();
float prev_coeff,next_coeff;
prev_coeff = prev_dir * reference_dir;
next_coeff = next_dir * reference_dir;
prev_coeff = prev_dir.dot(reference_dir);
next_coeff = next_dir.dot(reference_dir);
if (prev_coeff < 0.0f)
prev_coeff = 0.0f;
if (next_coeff < 0.0f)
@ -573,8 +573,8 @@ void AreaMode::Init(const std::vector < Point3f > &pts)
bool pts_not_in_line=false;
Point3f a,b;
for(unsigned int i=0;i<onethird;i++){
a=(pts[(i+ onethird )%npts] - pts[i%npts]).Normalize();
b=(pts[(i+(2*onethird))%npts] - pts[i%npts]).Normalize();
a=(pts[(i+ onethird )%npts] - pts[i%npts]).normalized();
b=(pts[(i+(2*onethird))%npts] - pts[i%npts]).normalized();
pts_not_in_line = (a ^ b).Norm() > EPSILON;
if(pts_not_in_line){
plane.Init( pts[i%npts],

58
wrap/gui/trackutils.h
View File

@ -103,8 +103,8 @@ Plane3f GetViewPlane (const View < float >&camera, const Point3f & center)
Point3f vp = camera.ViewPoint ();
Plane3f pl;
Point3f plnorm = vp - center;
plnorm.Normalize ();
pl.Set (plnorm, plnorm * center);
plnorm.Normalize();
pl.Set(plnorm, plnorm.dot(center));
return pl;
}
@ -331,16 +331,16 @@ Point3f HitSphere (Trackball * tb, const Point3f & p)
std::pair< float, bool > LineLineDistance(const Line3f & P,const Line3f & Q,Point3f & P_s, Point3f & Q_t){
Point3f p0 = P.Origin (), Vp = P.Direction ();
Point3f q0 = Q.Origin (), Vq = Q.Direction ();
float VPVP = Vp * Vp;
float VQVQ = Vq * Vq;
float VPVQ = Vp * Vq;
float VPVP = Vp.dot(Vp);
float VQVQ = Vq.dot(Vq);
float VPVQ = Vp.dot(Vq);
const float det = ( VPVP * VQVQ ) - ( VPVQ * VPVQ );
const float EPSILON = 0.00001f;
if ( fabs(det) < EPSILON ) {
return std::make_pair(Distance(P,q0), true);
}
float b1= (q0 - p0) * Vp;
float b2= (p0 - q0) * Vq;
float b1= (q0 - p0).dot(Vp);
float b2= (p0 - q0).dot(Vq);
float s = ( (VQVQ * b1) + (VPVQ * b2) ) / det;
float t = ( (VPVQ * b1) + (VPVP * b2) ) / det;
P_s = p0 + (Vp * s);
@ -367,16 +367,16 @@ std::pair< float, bool > LineLineDistance(const Line3f & P,const Line3f & Q,Poin
std::pair< float, bool > RayLineDistance(const Ray3f & R,const Line3f & Q,Point3f & R_s, Point3f & Q_t){
Point3f r0 = R.Origin (), Vr = R.Direction ();
Point3f q0 = Q.Origin (), Vq = Q.Direction ();
float VRVR = Vr * Vr;
float VQVQ = Vq * Vq;
float VRVQ = Vr * Vq;
float VRVR = Vr.dot(Vr);
float VQVQ = Vq.dot(Vq);
float VRVQ = Vr.dot(Vq);
const float det = ( VRVR * VQVQ ) - ( VRVQ * VRVQ );
const float EPSILON = 0.00001f;
if ( ( det >= 0.0f ? det : -det) < EPSILON ) {
return std::make_pair(Distance(Q,r0), true);
}
float b1= (q0 - r0) * Vr;
float b2= (r0 - q0) * Vq;
float b1= (q0 - r0).dot(Vr);
float b2= (r0 - q0).dot(Vq);
float s = ( (VQVQ * b1) + (VRVQ * b2) ) / det;
float t = ( (VRVQ * b1) + (VRVR * b2) ) / det;
if(s<0){
@ -418,17 +418,17 @@ std::pair< float, bool > SegmentSegmentDistance(const Segment3f & R, const Segme
R_s=ClosestPoint(R,Q_t);
return std::make_pair(Distance(R_s,Q_t),true);
}
Point3f r0 = R.P0(), Vr = (R.P1()-R.P0()).Normalize();
Point3f q0 = Q.P0(), Vq = (Q.P1()-Q.P0()).Normalize();
float VRVR = Vr * Vr;
float VQVQ = Vq * Vq;
float VRVQ = Vr * Vq;
Point3f r0 = R.P0(), Vr = (R.P1()-R.P0()).normalized();
Point3f q0 = Q.P0(), Vq = (Q.P1()-Q.P0()).normalized();
float VRVR = Vr.dot(Vr);
float VQVQ = Vq.dot(Vq);
float VRVQ = Vr.dot(Vq);
const float det = ( VRVR * VQVQ ) - ( VRVQ * VRVQ );
const float EPSILON = 0.00001f;
if ( ( det >= 0.0f ? det : -det) < EPSILON ) {
Line3f lR(R.P0(),R.P1());
float qa=lR.Projection(Q.P0());
float qb=lR.Projection(Q.P1());
float qb=lR.Projection(Q.P1());
if( (qa<=0.0f) && qb<=(0.0f)){
R_s=R.P0();
Q_t=ClosestPoint(Q,R_s);
@ -453,8 +453,8 @@ std::pair< float, bool > SegmentSegmentDistance(const Segment3f & R, const Segme
}
return std::make_pair(Distance(R_s,Q_t),true);
}
float b1= (q0 - r0) * Vr;
float b2= (r0 - q0) * Vq;
float b1= (q0 - r0).dot(Vr);
float b2= (r0 - q0).dot(Vq);
float s = ( (VQVQ * b1) + (VRVQ * b2) ) / det;
float t = ( (VRVQ * b1) + (VRVR * b2) ) / det;
if( s < 0 ){
@ -529,7 +529,7 @@ Line3f ProjectLineOnPlane(const Line3f & ln, const Plane3f & pl)
*/
float signedDistance(Line3f line,Point3f pt,Point3f positive_dir)
{
return Distance(line,pt) * ((((pt-ClosestPoint(line,pt)) * positive_dir) >= 0.0f )? 1.0f: -1.0f);
return Distance(line,pt) * ((((pt-ClosestPoint(line,pt)).dot(positive_dir)) >= 0.0f )? 1.0f: -1.0f);
}
/*!
@ -557,10 +557,10 @@ template<class T>
inline bool IntersectionRayPlane( const Plane3<T> & pl, const Ray3<T> & ray, Point3<T> &po){
const T epsilon = T(1e-8);
T k = pl.Direction() * ray.Direction(); // Compute 'k' factor
T k = pl.Direction().dot(ray.Direction()); // Compute 'k' factor
if( (k > -epsilon) && (k < epsilon))
return false;
T r = (pl.Offset() - pl.Direction()*ray.Origin())/k; // Compute ray distance
T r = (pl.Offset() - pl.Direction().dot(ray.Origin()))/k; // Compute ray distance
if (r < 0)
return false;
po = ray.Origin() + ray.Direction()*r;
@ -887,8 +887,8 @@ void DrawUglyPlaneMode(Trackball * tb,Plane3f plane)
if(norm == d1 || norm == -d1)
d1 = Point3f(1,0,0);
d2=plane.Projection(d1);
d1=(d2 - p0).Normalize();
d2=(d1 ^ norm).Normalize();
d1=(d2 - p0).normalized();
d2=(d1 ^ norm).normalized();
glLineWidth(3.0);
glColor3f(0.2f, 0.2f, 0.9f);
glBegin(GL_LINES);
@ -945,8 +945,8 @@ void DrawUglyCylinderMode(Trackball * tb,Line3f axis)
if(norm == d1 || norm == -d1)
d1 = Point3f(1,0,0);
d2=plane.Projection(d1);
d1=(d2 - p0).Normalize();
d2=(d1 ^ norm).Normalize();
d1=(d2 - p0).normalized();
d2=(d1 ^ norm).normalized();
glLineWidth(1.0);
glColor3f(0.2f, 0.2f, 0.9f);
for(int i=-100;i<100;i++){
@ -1098,8 +1098,8 @@ void DrawUglyAreaMode(Trackball * tb,const std::vector < Point3f > &points,
if(norm == d1 || norm == -d1)
d1 = Point3f(1,0,0);
d2=plane.Projection(d1);
d1=(d2 - p0).Normalize();
d2=(d1 ^ norm).Normalize();
d1=(d2 - p0).normalized();
d2=(d1 ^ norm).normalized();
glLineWidth(3.0);
glColor3f(0.2f, 0.2f, 0.9f);
glBegin(GL_LINES);

8
wrap/gui/view.h
View File

@ -151,10 +151,10 @@ template <class T> void View<T>::GetView() {
glGetDoublev(GL_PROJECTION_MATRIX, m);
proj.Import(Matrix44d(m));
Transpose(proj);
proj = proj.transpose().eval();
glGetDoublev(GL_MODELVIEW_MATRIX, m);
model.Import(Matrix44d(m));
Transpose(model);
model = model.transpose().eval();
glGetIntegerv(GL_VIEWPORT, (GLint*)viewport);
@ -205,14 +205,14 @@ template <class T> Line3<T> View<T>::ViewLineFromWindow(const Point3<T> &p)
template <class T> Point3<T> View<T>::Project(const Point3<T> &p) const {
Point3<T> r;
r = matrix * p;
r = matrix * p;
return NormDevCoordToWindowCoord(r);
}
template <class T> Point3<T> View<T>::UnProject(const Point3<T> &p) const {
Point3<T> s = WindowCoordToNormDevCoord(p);
s = inverse * s ;
return s;
return s;
}
// Come spiegato nelle glspec

2
wrap/io_trimesh/import_ptx.h
View File

@ -230,7 +230,7 @@ namespace io {
else return false;
// PTX transformation matrix is transposed
Transpose(currtrasf);
currtrasf = currtrasf.transpose().eval();
// allocating vertex space
int vn = rownum*colnum;