Removed all references to the useless old vcg::math::Swap

This commit is contained in:
Paolo Cignoni 2013-07-26 07:01:21 +00:00
parent f9db54196b
commit 04268b170f
7 changed files with 280 additions and 282 deletions

View File

@ -90,6 +90,7 @@ Edited Comments and GPL license
#include <math.h>
#include <assert.h>
#include <limits>
#include <algorithm>
#ifdef __BORLANDC__
@ -153,15 +154,12 @@ namespace math {
}
}
template<class T> inline void Swap(T &a, T &b){
T tmp=a; a=b; b=tmp;
}
template<class T> inline void Sort(T &a, T &b){
if (a>b) Swap(a,b);
if (a>b) std::swap(a,b);
}
template<class T> inline void Sort(T &a, T &b, T &c){
if (a>b) Swap(a,b);
if (b>c) {Swap(b,c); if (a>b) Swap(a,b);}
if (a>b) std::swap(a,b);
if (b>c) {std::swap(b,c); if (a>b) std::swap(a,b);}
}
/* Some <math.h> files do not define M_PI... */

View File

@ -30,18 +30,18 @@ namespace vcg {
namespace math {
/*
* Contrary to their definition, the Associated Legendre Polynomials presented here are
* Contrary to their definition, the Associated Legendre Polynomials presented here are
* only computed for positive m values.
*
*
*/
template <typename ScalarType>
class Legendre {
protected :
/**
* Legendre Polynomial three term Recurrence Relation
*/
*/
static inline ScalarType legendre_next(unsigned l, ScalarType P_lm1, ScalarType P_lm2, ScalarType x)
{
return ((2 * l + 1) * x * P_lm1 - l * P_lm2) / (l + 1);
@ -50,113 +50,113 @@ protected :
/**
* Associated Legendre Polynomial three term Recurrence Relation.
* Raises the band index.
*/
*/
static inline double legendre_next(unsigned l, unsigned m, ScalarType P_lm1, ScalarType P_lm2, ScalarType x)
{
return ((2 * l + 1) * x * P_lm1 - (l + m) * P_lm2) / (l + 1 - m);
}
/**
* Recurrence relation to compute P_m_(m+1) given P_m_m at the same x
*/
*/
static inline double legendre_P_m_mplusone(unsigned m, ScalarType p_m_m, ScalarType x)
{
return x * (2.0 * m + 1.0) * p_m_m;
}
/**
* Starting relation to compute P_m_m according to the formula:
*
*
* pow(-1, m) * double_factorial(2 * m - 1) * pow(1 - x*x, abs(m)/2)
*
*
* which becomes, if x = cos(theta) :
*
*
* pow(-1, m) * double_factorial(2 * m - 1) * pow(sin(theta), abs(m)/2)
*/
static inline double legendre_P_m_m(unsigned m, ScalarType sin_theta)
{
ScalarType p_m_m = 1.0;
if (m > 0)
{
ScalarType fact = 1.0; //Double factorial here
ScalarType fact = 1.0; //Double factorial here
for (unsigned i = 1; i <= m; ++i)
{
p_m_m *= fact * sin_theta; //raising sin_theta to the power of m/2
fact += 2.0;
}
if (m&1) //odd m
{
// Condon-Shortley Phase term
p_m_m *= -1;
}
}
return p_m_m;
}
static inline double legendre_P_l(unsigned l, ScalarType x)
{
ScalarType p0 = 1;
ScalarType p1 = x;
if (l == 0) return p0;
for (unsigned n = 1; n < l; ++n)
{
Swap(p0, p1);
std::swap(p0, p1);
p1 = legendre_next(n, p0, p1, x);
}
return p1;
}
/**
* Computes the Associated Legendre Polynomial for any given
* positive m and l, with m <= l and -1 <= x <= 1.
*/
static inline double legendre_P_l_m(unsigned l, unsigned m, ScalarType cos_theta, ScalarType sin_theta)
{
{
if(m > l) return 0;
if(m == 0) return legendre_P_l(l, cos_theta); //OK
else
{
ScalarType p_m_m = legendre_P_m_m(m, sin_theta); //OK
if (l == m) return p_m_m;
ScalarType p_m_mplusone = legendre_P_m_mplusone(m, p_m_m, cos_theta); //OK
if (l == m + 1) return p_m_mplusone;
if (l == m + 1) return p_m_mplusone;
unsigned n = m + 1;
while(n < l)
{
Swap(p_m_m, p_m_mplusone);
p_m_mplusone = legendre_next(n, m, p_m_m, p_m_mplusone, cos_theta);
++n;
std::swap(p_m_m, p_m_mplusone);
p_m_mplusone = legendre_next(n, m, p_m_m, p_m_mplusone, cos_theta);
++n;
}
return p_m_mplusone;
}
}
}
public :
static double Polynomial(unsigned l, ScalarType x)
{
assert (x <= 1 && x >= -1);
return legendre_P_l(l, x);
}
static double AssociatedPolynomial(unsigned l, unsigned m, ScalarType x)
{
assert (m <= l && x <= 1 && x >= -1);
return legendre_P_l_m(l, m, x, Sqrt(1.0 - x * x) );
}
static double AssociatedPolynomial(unsigned l, unsigned m, ScalarType cos_theta, ScalarType sin_theta)
{
assert (m <= l && cos_theta <= 1 && cos_theta >= -1 && sin_theta <= 1 && sin_theta >= -1);
@ -175,7 +175,7 @@ private:
ScalarType _sin_theta;
void generate(ScalarType cos_theta, ScalarType sin_theta)
{
{
//generate all 'l's with m = 0
matrix[0][0] = 1;

View File

@ -261,9 +261,9 @@ public:
/// Funzione per eseguire la trasposta della matrice
Matrix33 & Transpose()
{
math::Swap(a[1],a[3]);
math::Swap(a[2],a[6]);
math::Swap(a[5],a[7]);
std::swap(a[1],a[3]);
std::swap(a[2],a[6]);
std::swap(a[5],a[7]);
return *this;
}

View File

@ -516,7 +516,7 @@ This procedure decompose it in a sequence of
- ScaleV and Tranv are obiviously scaling and translation.
- ShearV contains three scalars with, respectively,
ShearXY, ShearXZ and ShearYZ
ShearXY, ShearXZ and ShearYZ
- RotateV contains the rotations (in degree!) around the x,y,z axis
The input matrix is modified leaving inside it a simple roto translation.
@ -661,7 +661,7 @@ template <class T> Point3<T> operator*(const Matrix44<T> &m, const Point3<T> &p)
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));
std::swap(m.ElementAt(i, j), m.ElementAt(j, i));
}
return m;
}

View File

@ -61,7 +61,7 @@ Revision 1.1 2004/03/16 03:07:38 tarini
namespace vcg {
namespace ndim{
namespace ndim{
//template <int N, class S>
@ -72,8 +72,8 @@ namespace vcg {
/**
The templated class for representing a point in R^N space.
The class is templated over the ScalarType class that is used to represent coordinates.
PointBase provides the interface and the common operators for points
of any dimensionality.
PointBase provides the interface and the common operators for points
of any dimensionality.
*/
template <int N, class S>
@ -88,7 +88,7 @@ public:
protected:
/// The only data member. Hidden to user.
S _v[N];
S _v[N];
public:
@ -103,11 +103,11 @@ public:
/// Padding function: give a default 0 value to all the elements that are not in the [0..2] range.
/// Useful for managing in a consistent way object that could have point2 / point3 / point4
inline S Ext( const int i ) const
{
if(i>=0 && i<=N) return _v[i];
else return 0;
}
inline S Ext( const int i ) const
{
if(i>=0 && i<=N) return _v[i];
else return 0;
}
/// importer for points with different scalar type and-or dimensionality
template <int N2, class S2>
@ -123,7 +123,7 @@ public:
template <int N2, class S2>
static inline PointType Construct( const Point<N2,S2> & b )
{
PointType p; p.Import(b);
PointType p; p.Import(b);
return p;
}
@ -141,7 +141,7 @@ public:
template <int N2, class S2>
static inline PointType Construct( const Point<N-1,S2> & b )
{
PointType p; p.ImportHomo(b);
PointType p; p.ImportHomo(b);
return p;
}
@ -163,30 +163,30 @@ public:
return _v[i];
}
inline const S &X() const { return _v[0]; }
inline const S &Y() const { return _v[1]; }
inline const S &Z() const { static_assert(N>2); return _v[2]; }
/// W is in any case the last coordinate.
/// (in a 2D point, W() == Y(). In a 3D point, W()==Z()
/// in a 4D point, W() is a separate component)
inline const S &W() const { return _v[N-1]; }
inline S &X() { return _v[0]; }
inline S &Y() { return _v[1]; }
inline S &Z() { static_assert(N>2); return _v[2]; }
inline S &W() { return _v[N-1]; }
inline const S * V() const
{
return _v;
}
inline S & V( const int i )
{
assert(i>=0 && i<N);
return _v[i];
}
inline const S & V( const int i ) const
{
assert(i>=0 && i<N);
return _v[i];
}
inline const S &Y() const { return _v[1]; }
inline const S &Z() const { static_assert(N>2); return _v[2]; }
/// W is in any case the last coordinate.
/// (in a 2D point, W() == Y(). In a 3D point, W()==Z()
/// in a 4D point, W() is a separate component)
inline const S &W() const { return _v[N-1]; }
inline S &X() { return _v[0]; }
inline S &Y() { return _v[1]; }
inline S &Z() { static_assert(N>2); return _v[2]; }
inline S &W() { return _v[N-1]; }
inline const S * V() const
{
return _v;
}
inline S & V( const int i )
{
assert(i>=0 && i<N);
return _v[i];
}
inline const S & V( const int i ) const
{
assert(i>=0 && i<N);
return _v[i];
}
//@}
//@{
@ -272,8 +272,8 @@ public:
/** @name Dot products (cross product "%" is defined olny for 3D points)
**/
/// Dot product
inline S operator * ( PointType const & p ) const;
/// Dot product
inline S operator * ( PointType const & p ) const;
/// slower version, more stable (double precision only)
inline S StableDot ( const PointType & p ) const;
@ -351,8 +351,8 @@ public:
//@{
/** @name
Glocal to Local and viceversa
(provided for uniformity with other spatial classes. trivial for points)
Glocal to Local and viceversa
(provided for uniformity with other spatial classes. trivial for points)
**/
inline PointType LocalToGlobal(ParamType p) const{
@ -381,7 +381,7 @@ public:
using Point<2,S>::V;
using Point<2,S>::W;
//@{
//@{
/** @name Special members for 2D points. **/
/// default
@ -426,7 +426,7 @@ public:
//@}
//@{
//@{
/** @name Implementation of standard functions for 3D points **/
inline void Zero(){
@ -475,7 +475,7 @@ public:
template <class PT> static S SquaredNorm(const PT &p ) {
return ( p.V(0)*p.V(0) + p.V(1)*p.V(1) );}
inline S operator * ( Point2 const & p ) const {
inline S operator * ( Point2 const & p ) const {
return ( _v[0]*p._v[0] + _v[1]*p._v[1]) ; }
inline bool operator == ( Point2 const & p ) const {
@ -545,7 +545,7 @@ public:
using Point<3,S>::W;
//@{
//@{
/** @name Special members for 3D points. **/
/// default
@ -563,7 +563,7 @@ public:
}
//@}
//@{
//@{
/** @name Implementation of standard functions for 3D points **/
inline void Zero(){
@ -612,7 +612,7 @@ public:
template <class PT> static S SquaredNorm(const PT &p ) {
return ( p.V(0)*p.V(0) + p.V(1)*p.V(1) + p.V(2)*p.V(2) );}
inline S operator * ( PointType const & p ) const {
inline S operator * ( PointType const & p ) const {
return ( _v[0]*p._v[0] + _v[1]*p._v[1] + _v[2]*p._v[2]) ; }
inline bool operator == ( PointType const & p ) const {
@ -623,19 +623,19 @@ public:
inline bool operator < ( PointType 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]); }
(_v[1]!=p._v[1])?(_v[1]< p._v[1]) : (_v[0]<p._v[0]); }
inline bool operator > ( PointType 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]); }
(_v[1]!=p._v[1])?(_v[1]> p._v[1]) : (_v[0]>p._v[0]); }
inline bool operator <= ( PointType const & p ) {
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]); }
(_v[1]!=p._v[1])?(_v[1]< p._v[1]) : (_v[0]<=p._v[0]); }
inline bool operator >= ( PointType 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]); }
(_v[1]!=p._v[1])?(_v[1]> p._v[1]) : (_v[0]>=p._v[0]); }
inline PointType & Normalize() {
S n = Norm();
@ -650,7 +650,7 @@ public:
inline S NormInfinity() const {
return math::Max( math::Max( math::Abs(_v[0]), math::Abs(_v[1]) ),
math::Abs(_v[3]) ); }
math::Abs(_v[3]) ); }
inline S NormOne() const {
return math::Abs(_v[0])+ math::Abs(_v[1])+math::Max(math::Abs(_v[2]));}
@ -701,10 +701,10 @@ public:
using Point<3,S>::V;
using Point<3,S>::W;
//@{
//@{
/** @name Special members for 4D points. **/
/// default
inline Point4 (){}
/// default
inline Point4 (){}
/// xyzw constructor
//@}
@ -764,7 +764,7 @@ public:
template <class PT> static S SquaredNorm(const PT &p ) {
return ( p.V(0)*p.V(0) + p.V(1)*p.V(1) + p.V(2)*p.V(2) + p.V(3)*p.V(3) );}
inline S operator * ( PointType const & p ) const {
inline S operator * ( PointType const & 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 == ( PointType const & p ) const {
@ -775,19 +775,19 @@ public:
inline bool operator < ( PointType 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]); }
(_v[1]!=p._v[1])?(_v[1]< p._v[1]) : (_v[0]<p._v[0]); }
inline bool operator > ( PointType 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]); }
(_v[1]!=p._v[1])?(_v[1]> p._v[1]) : (_v[0]>p._v[0]); }
inline bool operator <= ( PointType const & p ) {
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]); }
(_v[1]!=p._v[1])?(_v[1]< p._v[1]) : (_v[0]<=p._v[0]); }
inline bool operator >= ( PointType 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]); }
(_v[1]!=p._v[1])?(_v[1]> p._v[1]) : (_v[0]>=p._v[0]); }
inline PointType & Normalize() {
PointType n = Norm(); if(n!=0.0) { n=1.0/n; _v[0]*=n; _v[1]*=n; _v[2]*=n; _v[3]*=n; }
@ -803,7 +803,7 @@ public:
inline S NormInfinity() const {
return math::Max( math::Max( math::Abs(_v[0]), math::Abs(_v[1]) ),
math::Max( math::Abs(_v[2]), math::Abs(_v[3]) ) ); }
math::Max( math::Abs(_v[2]), math::Abs(_v[3]) ) ); }
inline S NormOne() const {
return math::Abs(_v[0])+ math::Abs(_v[1])+math::Max(math::Abs(_v[2]),math::Abs(_v[3]));}
@ -837,11 +837,11 @@ public:
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); }
if (exp0>exp1) { std::swap(k0,k1); std::swap(exp0,exp1); }
if (exp2>exp3) { std::swap(k2,k3); std::swap(exp2,exp3); }
if (exp0>exp2) { std::swap(k0,k2); std::swap(exp0,exp2); }
if (exp1>exp3) { std::swap(k1,k3); std::swap(exp1,exp3); }
if (exp2>exp3) { std::swap(k2,k3); std::swap(exp2,exp3); }
return ( (k0 + k1) + k2 ) +k3; }
//@}
@ -856,7 +856,7 @@ inline S Angle( Point3<S> const & p1, Point3<S> const & p2 )
S t = (p1*p2)/w;
if(t>1) t = 1;
else if(t<-1) t = -1;
return (S) acos(t);
return (S) acos(t);
}
// versione uguale alla precedente ma che assume che i due vettori siano unitari
@ -875,7 +875,7 @@ inline S AngleN( Point3<S> const & p1, Point3<S> const & p2 )
template <int N,class S>
inline S Norm( Point<N,S> const & p )
{
return p.Norm();
return p.Norm();
}
template <int N,class S>

View File

@ -69,14 +69,14 @@ namespace vcg {
/**
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.
All the usual operator (* + - ...) are defined.
*/
template <class T> class Point4
{
public:
/// The only data member. Hidden to user.
T _v[4];
T _v[4];
public:
typedef T ScalarType;
@ -134,104 +134,104 @@ public:
/** @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 * V()
{
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;
}
**/
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 * V()
{
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;
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];
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];
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;
}
//@}
@ -239,70 +239,70 @@ public:
//@{
/** @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
/// 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;
};
{
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]);
}
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]);
}
//@}
//@{
@ -316,10 +316,10 @@ public:
}
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();
}
{
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
@ -331,11 +331,11 @@ public:
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); }
if (exp0>exp1) { std::swap(k0,k1); std::swap(exp0,exp1); }
if (exp2>exp3) { std::swap(k2,k3); std::swap(exp2,exp3); }
if (exp0>exp2) { std::swap(k0,k2); std::swap(exp0,exp2); }
if (exp1>exp3) { std::swap(k1,k3); std::swap(exp1,exp3); }
if (exp2>exp3) { std::swap(k2,k3); std::swap(exp2,exp3); }
return ( (k0 + k1) + k2 ) +k3;
}
@ -363,19 +363,19 @@ inline T Norm( const Point4<T> & p )
template <class T>
inline T SquaredNorm( const Point4<T> & p )
{
return p.SquaredNorm();
return p.SquaredNorm();
}
template <class T>
inline T Distance( const Point4<T> & p1, const Point4<T> & p2 )
{
return Norm(p1-p2);
return Norm(p1-p2);
}
template <class T>
inline T SquaredDistance( const Point4<T> & p1, const Point4<T> & p2 )
{
return SquaredNorm(p1-p2);
return SquaredNorm(p1-p2);
}
/// slower version of dot product, more stable (double precision only)

View File

@ -88,13 +88,13 @@ public:
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];
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];
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;
}
@ -117,11 +117,11 @@ public:
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); }
if (exp0>exp1) { std::swap(k0,k1); std::swap(exp0,exp1); }
if (exp2>exp3) { std::swap(k2,k3); std::swap(exp2,exp3); }
if (exp0>exp2) { std::swap(k0,k2); std::swap(exp0,exp2); }
if (exp1>exp3) { std::swap(k1,k3); std::swap(exp1,exp3); }
if (exp2>exp3) { std::swap(k2,k3); std::swap(exp2,exp3); }
return ( (k0 + k1) + k2 ) +k3;
}