From 223902e61cdb15ea2725b560fb106bd273414b6a Mon Sep 17 00:00:00 2001 From: mtarini Date: Tue, 16 Mar 2004 03:07:38 +0000 Subject: [PATCH] "dimensionally unified" version: first commit --- vcg/space/box.h | 391 ++++++++++++++++++++++++++++++ vcg/space/point.h | 588 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 979 insertions(+) create mode 100644 vcg/space/box.h create mode 100644 vcg/space/point.h diff --git a/vcg/space/box.h b/vcg/space/box.h new file mode 100644 index 00000000..9afa8955 --- /dev/null +++ b/vcg/space/box.h @@ -0,0 +1,391 @@ +/**************************************************************************** +* 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.5 2004/03/05 17:51:28 tarini +Errorino "ScalarType" -> "S" + +Revision 1.4 2004/03/03 14:32:13 cignoni +Yet another cr lf mismatch + +Revision 1.3 2004/02/23 23:44:21 cignoni +cr lf mismatch + +Revision 1.2 2004/02/19 15:40:56 cignoni +Added doxygen groups + +Revision 1.1 2004/02/13 02:16:22 cignoni +First working release. + + +****************************************************************************/ + + +#ifndef __VCGLIB_BOX +#define __VCGLIB_BOX + +#include +#include +#include + +namespace vcg { + +/** \addtogroup space */ +/*@{*/ +/** +Templated class for 3D boxes. + This is the class for definition of a axis aligned box in 2D or 3D space. + Typically used as bounding boxes. + It is stored just as two Points (at the opposite vertices). + @param S (template parameter) Specifies the type of scalar used to represent coords. +*/ +template +class Box : public Space , Linear +{ +public: + typedef S ScalarType; + typedef Point ParamType; + typedef Point PointType; + enum {Dimension=N}; + + /// The scalar type +protected: + /// _min coordinate point + Point<3,S> _min; + /// _max coordinate point + Point<3,S> _max; + +public: + + inline const PointType &Max() const { return _max; } + inline PointType &Max() { return _max; } + inline const PointType &Min() const { return _min; } + inline PointType &Min() { return _min; } + + /// The box constructor + inline Box() { + _min.X()= 1;_max.X()= -1; + _min.Y()= 1;_max.Y()= -1; + if (N>2) {_min.Z()= 1;_max.Z()= -1;} + } + /// Min Max constructor + inline Box( const PointType & mi, const PointType & ma ) { _min = mi; _max = ma; } + /// The box distructor + inline ~Box() { } + /// Operator to compare two boxes + inline bool operator == ( Box const & p ) const + { + return _min==p._min && _max==p._max; + } + /// Operator to dispare two boxes + inline bool operator != ( Box const & p ) const + { + return _min!=p._min || _max!=p._max; + } + /** Infaltes the box of a percentage.. + @param s Scalar value. E.g if s=0.1 the box enlarges of 10% in every direction + if S==0.5 box doubles (+50% in every direction) + if S < 0 box shrinks + if S==0.5 box reduces to a point + */ + void Inflate( const S s ) + { + Inflate( (_max-_min)*s ); + } + /** Enlarges the box dimensions by k in every direction, with k = bbox.diag*s + */ + void InflateFix( const S s ) + { + S k = Diag()*s; + if (N==2) Inflate( PointType (k,k)); + if (N==3) Inflate( PointType (k,k,k)); + } + /** Enlarges the box dimensions by a fixed delta. + @param delta Point in D space. If delta > 0 box enlarges. If delta < 0 box reduces. + */ + void Inflate( const PointType & delta ) + { + _min -= delta; + _max += delta; + } + /// Initializing the box + void Set( const PointType & p ) + { + _min = _max = p; + } + /// Set the box to a null value + void SetNull() + { + _min.X()= 1; _max.X()= -1; + _min.Y()= 1; _max.Y()= -1; + _min.Z()= 1; _max.Z()= -1; + } + /** Add two boxex: + Returns minimal box that contains both operands. + @param b The box to add + */ + void Add( Box const & b ) + { + if(IsNull()) *this=b; + else + Add(_min); Add(_max); + } + /** Add a point to a box. + The box is modified is the added point is aoutside it. + @param p The point to add + */ + void Add( const PointType & p ) + { + if(IsNull()) Set(p); + else + { + if(_min.X() > p.X()) _min.X() = p.X(); + else if(_max.X() < p.X()) _max.X() = p.X(); + if(_min.Y() > p.Y()) _min.Y() = p.Y(); + else if(_max.Y() < p.Y()) _max.Y() = p.Y(); + if (N>2) { + if(_min.Z() > p.Z()) _min.Z() = p.Z(); + else if(_max.Z() < p.Z()) _max.Z() = p.Z(); + }; + } + } + /** Coputes intersection of Boxes: the minimal box containing both operands. + @param b The other operand + */ + void Intersect( const Box & b ) + { + if(_min.X() < b._min.X()) _min.X() = b._min.X(); + if(_min.Y() < b._min.Y()) _min.Y() = b._min.Y(); + if (N>2) if(_min.Z() < b._min.Z()) _min.Z() = b._min.Z(); + + if(_max.X() > b._max.X()) _max.X() = b._max.X(); + if(_max.Y() > b._max.Y()) _max.Y() = b._max.Y(); + if (N>2) if(_max.Z() > b._max.Z()) _max.Z() = b._max.Z(); + + if(_min.X()>_max.X() || _min.Y()>_max.Y() ) SetNull(); + else if (N>2) if (_min.Z()>_max.Z()) SetNull(); + } + /** Traslalate the box. + @param p: the translation vector + */ + void Translate( const PointType & p ) + { + _min += p; + _max += p; + } + /** Check wheter a point is inside box. + @param p The point + @returns True if inside, false otherwise + */ + bool IsIn( PointType const & p ) const + { + if (N==2) return ( + _min.X() <= p.X() && p.X() <= _max.X() && + _min.Y() <= p.Y() && p.Y() <= _max.Y() + ); + if (N==3) return ( + _min.X() <= p.X() && p.X() <= _max.X() && + _min.Y() <= p.Y() && p.Y() <= _max.Y() && + _min.Z() <= p.Z() && p.Z() <= _max.Z() + ); + } + /** Check wheter a point is inside box, open at left and closed at right [min..max) + @param p The point 3D + @returns True if inside, false otherwise + */ + bool IsInEx( PointType const & p ) const + { + if (N==2) return ( + _min.X() <= p.X() && p.X() < _max.X() && + _min.Y() <= p.Y() && p.Y() < _max.Y() + ); + if (N==3) return ( + _min.X() <= p.X() && p.X() < _max.X() && + _min.Y() <= p.Y() && p.Y() < _max.Y() && + _min.Z() <= p.Z() && p.Z() < _max.Z() + ); + } + /** + TODO: Move TO COLLIDE!!! + Verifica se due box collidono cioe' se hanno una intersezione non vuota. Per esempio + due box adiacenti non collidono. + @param b A box + @return True se collidoo, false altrimenti + */ + + bool Collide(Box const &b) + { + return b._min.X()<_max.X() && b._max.X()>_min.X() && + b._min.Y()<_max.Y() && b._max.Y()>_min.Y() && + b._min.Z()<_max.Z() && b._max.Z()>_min.Z() ; + } + /** Controlla se il box e' nullo. + @return True se il box e' nullo, false altrimenti + */ + bool IsNull() const { return _min.X()>_max.X() || _min.Y()>_max.Y() || _min.Z()>_max.Z(); } + /** Controlla se il box e' vuoto. + @return True se il box e' vuoto, false altrimenti + */ + bool IsEmpty() const { return _min==_max; } + /// Restituisce la lunghezza della diagonale del box. + S Diag() const + { + return Distance(_min,_max); + } + /// Calcola il quadrato della diagonale del box. + S SquaredDiag() const + { + return SquaredDistance(_min,_max); + } + /// Calcola il centro del box. + PointType Center() const + { + return (_min+_max)/2; + } + + /// Returns global coords of a local point expressed in [0..1]^3 + PointType LocalToGlobal(PointType const & p) const{ + return PointType( + _min[0] + p[0]*(_max[0]-_min[0]), + _min[1] + p[1]*(_max[1]-_min[1]), + _min[2] + p[2]*(_max[2]-_min[2])); + } + /// Returns local coords expressed in [0..1]^3 of a point in 3D + PointType GlobalToLocal(PointType const & p) const{ + return PointType( + (p[0]-_min[0])/(_max[0]-_min[0]), + (p[1]-_min[1])/(_max[1]-_min[1]), + (p[2]-_min[2])/(_max[2]-_min[2]) + ); + } + /// Computes the Volume for the box. + inline S Volume() const + { + if (N==2) return (_max.X()-_min.X())*(_max.Y()-_min.Y()); + if (N==3) return (_max.X()-_min.X())*(_max.Y()-_min.Y())*(_max.Z()-_min.Z()); + } + /// Area() and Volume() are sinonims (a + inline S Area() const { + return Volume(); + }; + /// Compute box size. + PointType Size() const + { + return (_max-_min); + } + /// Compute box size X. + inline S SizeX() const { return _max.X()-_min.X();} + /// Compute box size Y. + inline S SizeY() const { return _max.Y()-_min.Y();} + /// Compute box size Z. + inline S SizeZ() const { static_assert(N>2); return _max.Z()-_min.Z();} + + /** @name Linearity for boxes + **/ + + /// sets a point to Zero + inline void Zero() + { + _min.Zero(); + _max.Zero(); + } + inline Box operator + ( Box const & p) const + { + return Box(_min+p._min,_max+p._max); + } + inline Box operator - ( Box const & p) const + { + return Box(_min-p._min,_max-p._max); + } + inline Box operator * ( const S s ) const + { + return Box(_min*s,_max*s); + } + inline Box operator / ( const S s ) const + { + S inv=S(1.0)/s; + return Box(_min*inv,_max*inv); + } + inline Box & operator += ( Box const & p) + { + _min+=p._min; _max+=p._max; return *this; + } + inline Box & operator -= ( Box const & p) + { + _min-=p._min; _max-=p._max; return *this; + } + inline Box & operator *= ( const S s ) + { + _min*=s; _max*=s; return *this; + } + inline Box & operator /= ( const S s ) + { + S inv=S(1.0)/s; + _min*=s; _max*=s; return *this; + return *this; + } + inline Box operator - () const + { + return Box(-_min,-_max); + } +//@} + +//@{ + /** @name Iporters (for boxes in different spaces and with different scalar types) + **/ + + /// imports the box + template + inline void Import( const Box & b ) + { _max.Import( b._max );_min.Import( b._min ); + } + template + /// constructs a new ray importing it from an existing one + static Box Construct( const Box & b ) + { + return Box(PointType::Construct(b._min),PointType::Construct(b._max)); + } + + +}; // end class definition + + + + +typedef Box<3,short> Box3s; +typedef Box<3,int> Box3i; +typedef Box<3,float> Box3f; +typedef Box<3,double> Box3d; + +typedef Box<2,short> Box2s; +typedef Box<2,int> Box2i; +typedef Box<2,float> Box2f; +typedef Box<2,double> Box2d; + + +/*@}*/ + +} // end namespace +#endif diff --git a/vcg/space/point.h b/vcg/space/point.h new file mode 100644 index 00000000..e872385c --- /dev/null +++ b/vcg/space/point.h @@ -0,0 +1,588 @@ +/**************************************************************************** +* 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.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_POINT +#define __VCGLIB_POINT + +#include +#include +#include + +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 Point : public Space, Linear +{ +public: + typedef S ScalarType; + typedef VoidType ParamType; + typedef Point PointType; + enum {Dimension=N}; + +protected: + /// The only data member. Hidden to user. + S _v[N]; + +public: + +//@{ + + /** @name Standard Constructors and Initializers + No casting operators have been introduced to avoid automatic unattended (and costly) conversion between different point types + **/ + + inline Point () { } + inline Point ( const S nx, const S ny, const S nz, const S nw ) + { + static_assert(N==4); + _v[0] = nx; + _v[1] = ny; + _v[2] = nz; + _v[3] = nw; + } + inline Point ( const S nx, const S ny, const S nz) + { + static_assert(N==3); + _v[0] = nx; + _v[1] = ny; + _v[2] = nz; + } + inline Point ( const S nx, const S ny) + { + static_assert(N==2); + _v[0] = nx; + _v[1] = ny; + } + inline Point ( const S nv[N] ) + { + _v[0] = nv[0]; + _v[1] = nv[1]; + if (N>2) _v[2] = nv[2]; + if (N>3) _v[3] = nv[3]; + } + + /// 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; + } + + template + inline void Import( const Point & b ) + { + _v[0] = ScalarType(b[0]); + _v[1] = ScalarType(b[1]); + if (N>2) { if (N2>2) _v[2] = ScalarType(b[2]); else _v[2] = 0}; + if (N>3) { if (N2>3) _v[3] = ScalarType(b[3]); else _v[3] = 0}; + } + + static inline Point Construct( const PointType & b ) + { + PointType p; p.Import(b); + return p; + } + +//@} + +//@{ + + /** @name Data Access. + access to data is done by overloading of [] or explicit naming of coords (x,y,z)**/ + + inline S & operator [] ( const int i ) + { + assert(i>=0 && i=0 && i<3); + 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]; } + inline const S &W() const { static_assert(N>3); return _v[3]; } + 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() { static_assert(N>3); return _v[3]; } + inline const S * V() const + { + return _v; + } + inline S & V( const int i ) + { + assert(i>=0 && i=0 && i2) _v[2] = 0; + if (N>3) _v[3] = 0; + } + inline Point operator + ( Point const & p) const + { + if (N==2) return Point( _v[0]+p._v[0], _v[1]+p._v[1] ); + if (N==3) return Point( _v[0]+p._v[0], _v[1]+p._v[1], _v[2]+p._v[2] ); + if (N==4) return Point( _v[0]+p._v[0], _v[1]+p._v[1], _v[2]+p._v[2], _v[3]+p._v[3] ); + } + inline Point operator - ( Point const & p) const + { + if (N==2) return Point( _v[0]-p._v[0], _v[1]-p._v[1] ); + if (N==3) return Point( _v[0]-p._v[0], _v[1]-p._v[1], _v[2]-p._v[2] ); + if (N==4) return Point( _v[0]-p._v[0], _v[1]-p._v[1], _v[2]-p._v[2], _v[3]-p._v[3] ); + } + inline Point operator * ( const S s ) const + { + if (N==2) return Point( _v[0]*s, _v[1]*s ); + if (N==3) return Point( _v[0]*s, _v[1]*s, _v[2]*s ); + if (N==4) return Point( _v[0]*s, _v[1]*s, _v[2]*s, _v[3]*s ); + } + inline Point operator / ( const S s ) const + { + if (N==2) return Point( _v[0]/s, _v[1]/s ); + if (N==3) return Point( _v[0]/s, _v[1]/s, _v[2]/s ); + if (N==4) return Point( _v[0]/s, _v[1]/s, _v[2]/s, _v[3]/s ); + } + inline Point & operator += ( Point const & p) + { + _v[0] += p._v[0]; + _v[1] += p._v[1]; + if (N>2) _v[2] += p._v[2]; + if (N>3) _v[3] += p._v[3]; + return *this; + } + inline Point & operator -= ( Point const & p) + { + _v[0] -= p._v[0]; + _v[1] -= p._v[1]; + if (N>2) _v[2] -= p._v[2]; + if (N>3) _v[3] -= p._v[3]; + return *this; + } + inline Point & operator *= ( const S s ) + { + _v[0] *= s; + _v[1] *= s; + if (N>2) _v[2] *= s; + if (N>3) _v[3] *= s; + return *this; + } + inline Point & operator /= ( const S s ) + { + _v[0] /= s; + _v[1] /= s; + if (N>2) _v[2] /= s; + if (N>3) _v[3] /= s; + return *this; + } + inline Point operator - () const + { + if (N==2) return Point ( -_v[0], -_v[1] ); + if (N==3) return Point ( -_v[0], -_v[1], -_v[2] ); + if (N==4) return Point ( -_v[0], -_v[1], -_v[2] , -_v[3] ); + } +//@} +//@{ + + /** @name Dot products + **/ + /// Dot product + inline S operator * ( Point const & p ) const + { + if (N==2) return ( _v[0]*p._v[0] + _v[1]*p._v[1] ); + if (N==3) return ( _v[0]*p._v[0] + _v[1]*p._v[1] + _v[2]*p._v[2] ); + if (N==4) return ( _v[0]*p._v[0] + _v[1]*p._v[1] + _v[2]*p._v[2] + _v[2]*p._v[2] ); + }; + /// slower version, more stable (double precision only) + inline S StableDot ( const Point & p ) const + { + if (N==2) return _v[0]*p._v[0] + _v[1]*p._v[1]; + if (N==4) { + + S 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; + }; + if (N==3) { + T k0=_v[0]*p._v[0], k1=_v[1]*p._v[1], k2=_v[2]*p._v[2]; + int exp0,exp1,exp2; + + frexp( double(k0), &exp0 ); + frexp( double(k1), &exp1 ); + frexp( double(k2), &exp2 ); + + if( exp0 + ( + _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] + ); + } + /// Cross product for 2D Point + /// if called from a 3D or 4D points, returns the z component of the cross prod. + inline S operator % ( Point const & p ) const + { + return _v[0]*p._v[1] - _v[1]*p._v[0]; + } +//@} +//@{ + + /** @name Norms + **/ + + + // Euclidean norm + inline S Norm() const + { + if (N==2) return math::Sqrt( _v[0]*_v[0] + _v[1]*_v[1] ); + if (N==3) return math::Sqrt( _v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] ); + if (N==4) return math::Sqrt( _v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] + _v[4]*_v[4] ); + } + // Squared Euclidean norm + inline S SquaredNorm() const + { + if (N==2) return ( _v[0]*_v[0] + _v[1]*_v[1] ); + if (N==3) return ( _v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] ); + if (N==4) return ( _v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] + _v[3]*_v[3] ); + } + // Normalization (division by norm) + inline Point & Normalize() + { + S n = Norm(); + if(n>0.0) (*this)/=n; + return *this; + } + /// Homogeneous normalization (division by W) + inline Point & HomoNormalize(){ + if (_v[3]!=0.0) { _v[0] /= _v[3]; _v[1] /= _v[3]; _v[2] /= _v[3]; _v[3]=1.0; } + return *this; + }; + +//@} + // Per component scaling + inline Point & Scale( const Point & p ) + { + _v[0] *= p._v[0]; + _v[1] *= p._v[1]; + if (N>2) _v[2] *= p._v[2]; + if (N>3) _v[3] *= p._v[3]; + return *this; + } + + + // Convert to polar coordinates + void ToPolar( S & ro, S & tetha, S & fi ) const + { + ro = Norm(); + tetha = (S)atan2( _v[1], _v[0] ); + fi = (S)acos( _v[2]/ro ); + } + +//@{ + + /** @name Comparison Operators. + Lexicographical order. + **/ + + inline bool operator == ( Point const & p ) const + { + if (N==2) return _v[0]==p._v[0] && _v[1]==p._v[1]; + if (N==3) return _v[0]==p._v[0] && _v[1]==p._v[1] && _v[2]==p._v[2]; + if (N==4) return _v[0]==p._v[0] && _v[1]==p._v[1] && _v[2]==p._v[2] && _v[3]==p._v[3]; + } + inline bool operator != ( Point const & p ) const + { + if (N==2) return _v[0]!=p._v[0] || _v[1]!=p._v[1] ; + if (N==3) return _v[0]!=p._v[0] || _v[1]!=p._v[1] || _v[2]!=p._v[2]; + if (N==4) return _v[0]!=p._v[0] || _v[1]!=p._v[1] || _v[2]!=p._v[2] || _v[3]!=p._v[3]; + } + inline bool operator < ( Point const & p ) const + { + return (_v[2]!=p._v[2])?(_v[2] ( Point 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 <= ( Point 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 >= ( Point 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 PointType LocalToGlobal(ParamType p) const{ + return *this; + }; + + //@} +}; // end class definition + + +template +inline S Angle( Point<3,S> const & p1, Point<3,S> const & p2 ) +{ + S w = p1.Norm()*p2.Norm(); + if(w==0) return -1; + S t = (p1*p2)/w; + if(t>1) t = 1; + else if(t<-1) t = -1; + return (S) acos(t); +} + +// versione uguale alla precedente ma che assume che i due vettori sono unitari +template +inline S AngleN( Point<3,S> const & p1, Point<3,S> const & p2 ) +{ + S w = p1*p2; + if(w>1) + w = 1; + else if(w<-1) + w=-1; + return (S) acos(w); +} + + +template +inline S Norm( Point const & p ) +{ + return p.Norm(); +} + +template +inline S SquaredNorm( Point const & p ) +{ + return p.SquaredNorm(); +} + +template +inline Point & Normalize( Point & p ) +{ + p.Normalize(); + return p; +} + +template +inline S Distance( Point const & p1,Point const & p2 ) +{ + return (p1-p2).Norm(); +} + +template +inline S SquaredDistance( Point const & p1,Point 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 +double StableDot ( Point<3,S> const & p0, Point<3,S> const & p1 ) +{ + +} + +/// Computes a shape quality measure of the triangle composed by points p0,p1,p2 +/// It Returns 2*AreaTri/(MaxEdge^2), +/// the range is range [0.0, 0.866] +/// e.g. Equilateral triangle sqrt(3)/2, halfsquare: 1/2, ... up to a line that has zero quality. +template +S Quality( Point<3,S> const &p0, Point<3,S> const & p1, Point<3,S> const & p2) +{ + PointType d10=p1-p0; + PointType d20=p2-p0; + PointType d12=p1-p2; + PointType x = d10^d20; + + S a = Norm( x ); + if(a==0) return 0; // Area zero triangles have surely quality==0; + S b = SquaredNorm( d10 ); + S t = b; + t = SquaredNorm( d20 ); if ( b +Point<3,S> Normal(const Point<3,S> & p0, const Point<3,S> & p1, const Point<3,S> & p2) +{ + return ((p1 - p0) ^ (p2 - p0)); +} + +/// Like the above, it returns the normal to the plane passing through p0,p1,p2, but normalized. +template +Point<3,S> NormalizedNormal(const Point<3,S> & p0, const Point<3,S> & p1, const Point<3,S> & p2) +{ + return ((p1 - p0) ^ (p2 - p0)).Normalize(); +} + + +/// Point(p) Edge(v1-v2) dist, q is the point in v1-v2 with min dist +template +S PSDist( const Point<3,S> & p, + const Point<3,S> & v1, + const Point<3,S> & v2, + Point<3,S> & q ) +{ + Point<3,S> e = v2-v1; + S 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 +inline Point<2,S>::Point ( const S nx, const S ny ) +{_v[0]=nx;_v[1]=ny;};*/ + + +/*template +inline Point<4,S>::Point ( const S nx, const S ny , const S nz , const S nw ) +{_v[0]=nx;_v[1]=ny;_v[2]=nz;_v[3]=nw;};*/ + +/*template < class S> + Point<3,S> Point<3,S>::operator * ( const S s ) const + { + return Point<3,S>( _v[0]*s, _v[1]*s , _v[2]*s ); + }*/ + +//template < class S> + /*Point<3,double> Point<2,double>::operator * ( const double s ) const + { + return Point<2,double>( _v[0]*s, _v[1]*s ); + }*/ + +typedef Point<3,short> Point3s; +typedef Point<3,int> Point3i; +typedef Point<3,float> Point3f; +typedef Point<3,double> Point3d; +/*@}*/ + + + +} // end namespace +#endif +