text: /****************************************************************************
  * 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.9  2006/12/20 15:23:52  ganovelli
  using of locally defined variable removed
  
  Revision 1.8  2006/04/11 08:10:05  zifnab1974
  changes necessary for gcc 3.4.5 on linux 64bit.
  
  Revision 1.7  2005/12/12 11:22:32  ganovelli
  compiled with gcc
  
  Revision 1.6  2005/01/12 11:25:52  ganovelli
  corrected Point<3
  
  Revision 1.5  2004/10/20 16:45:21  ganovelli
  first compiling version (MC,INtel,gcc)
  
  Revision 1.4  2004/04/29 10:47:06  ganovelli
  some siyntax error corrected
  
  Revision 1.3  2004/04/05 12:36:43  tarini
  unified version: PointBase version, with no guards "(N==3)"
  
  
  
  Revision 1.1  2004/03/16 03:07:38  tarini
  "dimensionally unified" version: first commit
  
  ****************************************************************************/
  
  #ifndef __VCGLIB_POINT
  #define __VCGLIB_POINT
  
  #include <assert.h>
  #include <vcg/math/base.h>
  #include <vcg/space/space.h>
  
  namespace vcg {
  
      namespace ndim{
  
  
  //template <int N, class S>
  //class Point;
  
  /** \addtogroup space */
  /*@{*/
      /**
          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.
       */
  
  template <int N, class S>
  class Point
  {
  public:
  	typedef S          ScalarType;
  	typedef VoidType   ParamType;
  	typedef Point<N,S> 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 PointType types
     **/
  
    inline Point () { };
  //	inline Point ( const S nv[N] );
  
    /// 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;
      }
  
  	/// importer for points with different scalar type and-or dimensionality
  	template <int N2, class S2>
  	inline void Import( const Point<N2,S2> & 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;};
  	}
  
  	/// constructor for points with different scalar type and-or dimensionality
  	template <int N2, class S2>
    static inline PointType Construct( const Point<N2,S2> & b )
    {
          PointType p; p.Import(b);
      return p;
    }
  
  	  /// importer for homogeneous points
  	template <class S2>
  	inline void ImportHomo( const Point<N-1,S2> & b )
  	{
  		_v[0] = ScalarType(b[0]);
  		_v[1] = ScalarType(b[1]);
  		if (N>2) { _v[2] = ScalarType(_v[2]); };
  		_v[N-1] = 1.0;
  	}
  
  		/// constructor for homogeneus point.
  	template <int N2, class S2>
    static inline PointType Construct( const Point<N-1,S2> & b )
    {
          PointType p; p.ImportHomo(b);
      return p;
    }
  
  //@}
  
  //@{
  
    /** @name Data Access.
     access to data is done by overloading of [] or explicit naming of coords (x,y,z)**/
  
  	inline S & operator [] ( const int i )
  	{
  		assert(i>=0 && i<N);
  		return _v[i];
  	}
  	inline const S & operator [] ( const int i ) const
  	{
  		assert(i>=0 && i<N);
  		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()
  	{
  		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];
      }
  //@}
  //@{
  
    /** @name Linearity for points
    **/
  
  	/// sets a PointType to Zero
  	inline void SetZero()
  	{
  		for(unsigned int ii = 0; ii < Dimension;++ii)
  			V(ii) = S();
  	}
  
  	inline PointType operator + ( PointType const & p) const
  	{
  		PointType res;
  		for(unsigned int ii = 0; ii < Dimension;++ii)
  			res[ii] = V(ii) + p[ii];
  		return res;
  	}
  
  	inline PointType operator - ( PointType const & p) const
  	{
  		PointType res;
  		for(unsigned int ii = 0; ii < Dimension;++ii)
  			res[ii] = V(ii) - p[ii];
  		return res;
  	}
  
  	inline PointType operator * ( const S s ) const
  	{
  		PointType res;
  		for(unsigned int ii = 0; ii < Dimension;++ii)
  			res[ii] = V(ii) * s;
  		return res;
  	}
  
  	inline PointType operator / ( const S s ) const
  	{
  		PointType res;
  		for(unsigned int ii = 0; ii < Dimension;++ii)
  			res[ii] = V(ii) / s;
  		return res;
  	}
  
  	inline PointType & operator += ( PointType const & p)
  	{
  		for(unsigned int ii = 0; ii < Dimension;++ii)
  			V(ii) += p[ii];
  		return *this;
  	}
  
  	inline PointType & operator -= ( PointType const & p)
  	{
  		for(unsigned int ii = 0; ii < Dimension;++ii)
  			V(ii) -= p[ii];
  		return *this;
  	}
  
  	inline PointType & operator *= ( const S s )
  	{
  		for(unsigned int ii = 0; ii < Dimension;++ii)
  			V(ii) *= s;
  		return *this;
  	}
  
  	inline PointType & operator /= ( const S s )
  	{
  		for(unsigned int ii = 0; ii < Dimension;++ii)
  			V(ii) *= s;
  		return *this;
  	}
  
  	inline PointType operator - () const
  	{
  		PointType res;
  		for(unsigned int ii = 0; ii < Dimension;++ii)
  			res[ii] = - V(ii);
  		return res;
  	}
  //@}
  //@{
  
    /** @name Dot products (cross product "%" is defined olny for 3D points)
    **/
          /// Dot product
      inline S operator * ( PointType const & p ) const;
  
  	/// slower version, more stable (double precision only)
  	inline S StableDot ( const PointType & p ) const;
  
  //@}
  
  //@{
  
    /** @name Norms
    **/
  
  		/// Euclidean norm
  	inline S Norm() const;
  		/// Euclidean norm, static version
  	template <class PT> static S Norm(const PT &p );
  		/// Squared Euclidean norm
  	inline S SquaredNorm() const;
  		/// Squared Euclidean norm, static version
  	template <class PT> static S SquaredNorm(const PT &p );
  		/// Normalization (division by norm)
  	inline PointType & Normalize();
  		/// Normalization (division by norm), static version
  	template <class PT> static PointType & Normalize(const PT &p);
  		/// Homogeneous normalization (division by W)
  	inline PointType & HomoNormalize();
  
  	  /// norm infinity: largest absolute value of compoenet
  	inline S NormInfinity() const;
  		/// norm 1: sum of absolute values of components
  	inline S NormOne() const;
  
  //@}
  		/// Signed area operator
  	  /// a % b returns the signed area of the parallelogram inside a and b
  	inline S operator % ( PointType const & p ) const;
  
  	  /// the sum of the components
  	inline S Sum() const;
  	  /// returns the biggest component
  	inline S Max() const;
  	  /// returns the smallest component
  	inline S Min() const;
  	  /// returns the index of the biggest component
  	inline int MaxI() const;
  	  /// returns the index of the smallest component
  	inline int MinI() const;
  
  
  	 /// Per component scaling
  	inline PointType & Scale( const PointType & p );
  
  
  	 /// 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.
     Lexicographic order.
     **/
  
  	inline bool operator == ( PointType const & p ) const;
  	inline bool operator != ( PointType const & p ) const;
  	inline bool operator <  ( PointType const & p ) const;
  	inline bool operator >  ( PointType const & p ) const;
  	inline bool operator <= ( PointType const & p ) const;
  	inline bool operator >= ( PointType const & p ) const;
   //@}
  
  //@{
  
    /** @name
      Glocal to Local and viceversa
      (provided for uniformity with other spatial classes. trivial for points)
     **/
  
  	inline PointType LocalToGlobal(ParamType p) const{
  		return *this; }
  
    inline ParamType GlobalToLocal(PointType /*p*/) const{
        return ParamType(); }
  //@}
  
  }; // end class definition
  
  
  
  
  
  
  
  
  
  template <class S>
  class Point2 : public Point<2,S> {
  public:
  	typedef S ScalarType;
  	typedef Point2 PointType;
  	using Point<2,S>::_v;
  	using Point<2,S>::V;
  	using Point<2,S>::W;
  
      //@{
    /** @name Special members for 2D points. **/
  
  	/// default
  	inline Point2 (){}
  
  	  /// yx constructor
  	inline Point2 ( const S a, const S b){
  		_v[0]=a; _v[1]=b;  };
  
  	  /// unary orthogonal operator (2D equivalent of cross product)
  	  /// returns orthogonal vector (90 deg left)
  	inline Point2 operator ~ () const {
  		return Point2 ( -_v[2], _v[1] );
  	}
  
  	  /// returns the angle with X axis (radiants, in [-PI, +PI] )
  	inline ScalarType &Angle(){
  		return math::Atan2(_v[1],_v[0]);}
  
  		/// transform the point in cartesian coords into polar coords
  	inline Point2 & ToPolar(){
  		ScalarType t = Angle();
  		_v[0] = Norm();
  		_v[1] = t;
  		return *this;}
  
  		/// transform the point in polar coords into cartesian coords
  	inline Point2 & ToCartesian() {
  		ScalarType l = _v[0];
  		_v[0] = (ScalarType)(l*math::Cos(_v[1]));
  		_v[1] = (ScalarType)(l*math::Sin(_v[1]));
  		return *this;}
  
  			/// rotates the point of an angle (radiants, counterclockwise)
  	inline Point2 & Rotate( const ScalarType rad ){
  		ScalarType t = _v[0];
  		ScalarType s = math::Sin(rad);
  		ScalarType c = math::Cos(rad);
  		_v[0] = _v[0]*c - _v[1]*s;
  		_v[1] =   t *s + _v[1]*c;
  		return *this;}
  
  	 //@}
  
      //@{
    /** @name Implementation of standard functions for 3D points **/
  
  	inline void Zero(){
  		_v[0]=0; _v[1]=0;  };
  
  	inline Point2 ( const S nv[2] ){
  		_v[0]=nv[0]; _v[1]=nv[1];   };
  
  	inline Point2 operator + ( Point2 const & p) const {
  		return Point2( _v[0]+p._v[0], _v[1]+p._v[1]); }
  
  	inline Point2 operator - ( Point2 const & p) const {
  		return Point2( _v[0]-p._v[0], _v[1]-p._v[1]); }
  
  	inline Point2 operator * ( const S s ) const {
  		return Point2( _v[0]*s, _v[1]*s  ); }
  
  	inline Point2 operator / ( const S s ) const {
  		S t=1.0/s;
  		return Point2( _v[0]*t, _v[1]*t   ); }
  
  	inline Point2 operator - () const {
  		return Point2 ( -_v[0], -_v[1]  ); }
  
  	inline Point2 & operator += ( Point2 const & p ) {
  		_v[0] += p._v[0]; _v[1] += p._v[1];   return *this; }
  
  	inline Point2 & operator -= ( Point2 const & p ) {
  		_v[0] -= p._v[0]; _v[1] -= p._v[1];   return *this; }
  
  	inline Point2 & operator *= ( const S s ) {
  		_v[0] *= s; _v[1] *= s;  return *this; }
  
  	inline Point2 & operator /= ( const S s ) {
  		S t=1.0/s; _v[0] *= t; _v[1] *= t;   return *this; }
  
  	inline S Norm() const {
  		return math::Sqrt( _v[0]*_v[0] + _v[1]*_v[1]   );}
  
  	template <class PT> static S Norm(const PT &p ) {
  		return math::Sqrt( p.V(0)*p.V(0) + p.V(1)*p.V(1)  );}
  
  	inline S SquaredNorm() const {
  		return ( _v[0]*_v[0] + _v[1]*_v[1] );}
  
  	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 {
      return ( _v[0]*p._v[0] + _v[1]*p._v[1]) ; }
  
  	inline bool operator == ( Point2 const & p ) const {
  		return _v[0]==p._v[0] && _v[1]==p._v[1] ;}
  
  	inline bool operator != ( Point2 const & p ) const {
  		return _v[0]!=p._v[0] || _v[1]!=p._v[1] ;}
  
  	inline bool operator <  ( Point2 const & p ) const{
  		return (_v[1]!=p._v[1])?(_v[1]< p._v[1]) : (_v[0]<p._v[0]); }
  
  	inline bool operator >  ( Point2 const & p ) const	{
  		return (_v[1]!=p._v[1])?(_v[1]> p._v[1]) : (_v[0]>p._v[0]); }
  
  	inline bool operator <= ( Point2 const & p ) {
  		return (_v[1]!=p._v[1])?(_v[1]< p._v[1]) : (_v[0]<=p._v[0]); }
  
  	inline bool operator >= ( Point2 const & p ) const {
  		return (_v[1]!=p._v[1])?(_v[1]> p._v[1]) : (_v[0]>=p._v[0]); }
  
  	inline Point2 & Normalize() {
  		PointType n = Norm(); if(n!=0.0) { n=1.0/n;	_v[0]*=n;	_v[1]*=n;} return *this;}
  
  	inline Point2 & HomoNormalize(){
  		if (_v[2]!=0.0) {	_v[0] /= W(); W()=1.0; } return *this;}
  
  	inline S NormInfinity() const {
  		return math::Max( math::Abs(_v[0]), math::Abs(_v[1]) ); }
  
  	inline S NormOne() const {
  		return math::Abs(_v[0])+ math::Abs(_v[1]);}
  
  	inline S operator % ( Point2 const & p ) const {
  		return _v[0] * p._v[1] - _v[1] * p._v[0]; }
  
  	inline S Sum() const {
  		return _v[0]+_v[1];}
  
  	inline S Max() const {
  		return math::Max( _v[0], _v[1] ); }
  
  	inline S Min() const {
  		return math::Min( _v[0], _v[1] ); }
  
  	inline int MaxI() const {
  		return (_v[0] < _v[1]) ? 1:0; };
  
  	inline int MinI() const {
  	  return (_v[0] > _v[1]) ? 1:0; };
  
  	inline PointType & Scale( const PointType & p ) {
  		_v[0] *= p._v[0]; _v[1] *= p._v[1];  return *this; }
  
  	inline S StableDot ( const PointType & p ) const {
  		return _v[0]*p._v[0] +_v[1]*p._v[1]; }
  	 //@}
  
  };
  
  template <typename S>
  class Point3 : public Point<3,S> {
  public:
  	typedef S ScalarType;
  	typedef Point3<S>  PointType;
  	using Point<3,S>::_v;
  	using Point<3,S>::V;
  	using Point<3,S>::W;
  
  
      //@{
    /** @name Special members for 3D points. **/
  
  	/// default
  	inline Point3 ():Point<3,S>(){}
  	  /// yxz constructor
  	inline Point3 ( const S a, const S b,  const S c){
  		_v[0]=a; _v[1]=b; _v[2]=c; };
  
  		/// Cross product for 3D points
  	inline PointType operator ^ ( PointType const & p ) const {
  		return Point3 (
  			_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] );
  	}
  	 //@}
  
      //@{
    /** @name Implementation of standard functions for 3D points **/
  
  	inline void Zero(){
  		_v[0]=0; _v[1]=0; _v[2]=0; };
  
  	inline Point3 ( const S nv[3] ){
  		_v[0]=nv[0]; _v[1]=nv[1]; _v[2]=nv[2];  };
  
  	inline Point3 operator + ( Point3 const & p)  const{
  		return Point3( _v[0]+p._v[0], _v[1]+p._v[1], _v[2]+p._v[2]); }
  
  	inline Point3 operator - ( Point3 const & p) const {
  		return Point3( _v[0]-p._v[0], _v[1]-p._v[1], _v[2]-p._v[2]); }
  
  	inline Point3 operator * ( const S s ) const {
  		return Point3( _v[0]*s, _v[1]*s , _v[2]*s  ); }
  
  	inline Point3 operator / ( const S s ) const {
  		S t=1.0/s;
  		return Point3( _v[0]*t, _v[1]*t , _v[2]*t  ); }
  
  	inline Point3 operator - () const {
  		return Point3 ( -_v[0], -_v[1] , -_v[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 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 S s ) {
  		_v[0] *= s; _v[1] *= s; _v[2] *= s; return *this; }
  
  	inline Point3 & operator /= ( const S s ) {
  		S t=1.0/s; _v[0] *= t; _v[1] *= t;  _v[2] *= t;   return *this; }
  
  	inline S Norm() const {
  		return math::Sqrt( _v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2]  );}
  
  	template <class PT> static S Norm(const PT &p ) {
  		return math::Sqrt( p.V(0)*p.V(0) + p.V(1)*p.V(1) + p.V(2)*p.V(2)  );}
  
  	inline S SquaredNorm() const {
  		return ( _v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2]  );}
  
  	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 {
      return ( _v[0]*p._v[0] + _v[1]*p._v[1] + _v[2]*p._v[2]) ; }
  
  	inline bool 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 {
  		return _v[0]!=p._v[0] || _v[1]!=p._v[1] || _v[2]!=p._v[2] ;}
  
  	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]); }
  
  	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]); }
  
  	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]); }
  
  	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]); }
  
  	inline PointType & Normalize() {
  		S n = Norm();
  		if(n!=0.0) {
  			n=S(1.0)/n;
  		_v[0]*=n;	_v[1]*=n;	_v[2]*=n; }
  		return *this;}
  
  	inline PointType & HomoNormalize(){
  		if (_v[2]!=0.0) {	_v[0] /= W();	_v[1] /= W(); W()=1.0; }
  		return *this;}
  
  	inline S NormInfinity() const {
  		return math::Max( math::Max( math::Abs(_v[0]), math::Abs(_v[1]) ),
  							math::Abs(_v[3])  ); }
  
  	inline S NormOne() const {
  		return math::Abs(_v[0])+ math::Abs(_v[1])+math::Max(math::Abs(_v[2]));}
  
  	inline S operator % ( PointType const & p ) const {
  		S t = (*this)*p; /* Area, general formula */
  		return math::Sqrt( SquaredNorm() * p.SquaredNorm() - (t*t) );};
  
  	inline S Sum() const {
  		return _v[0]+_v[1]+_v[2];}
  
  	inline S Max() const {
  		return math::Max( math::Max( _v[0], _v[1] ),  _v[2] ); }
  
  	inline S Min() const {
  		return math::Min( math::Min( _v[0], _v[1] ),  _v[2] ); }
  
  	inline int MaxI() const {
  		int i= (_v[0] < _v[1]) ? 1:0; if (_v[i] < _v[2]) i=2; return i;};
  
  	inline int MinI() const {
  		int i= (_v[0] > _v[1]) ? 1:0; if (_v[i] > _v[2]) i=2; return i;};
  
  	inline PointType & Scale( const PointType & p ) {
  		_v[0] *= p._v[0]; _v[1] *= p._v[1]; _v[2] *= p._v[2]; return *this; }
  
  	inline S StableDot ( const PointType & p ) const {
  		S 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<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;
  	}
  	 //@}
  
  };
  
  template <typename S>
  class Point4 : public Point<4,S> {
  public:
  	typedef S ScalarType;
  	typedef Point4<S> PointType;
  	using Point<3,S>::_v;
  	using Point<3,S>::V;
  	using Point<3,S>::W;
  
      //@{
    /** @name Special members for 4D points. **/
          /// default
          inline Point4 (){}
  
  	  /// xyzw constructor
  	//@}
  	inline Point4 ( const S a, const S b,  const S c, const S d){
  		_v[0]=a; _v[1]=b; _v[2]=c; _v[3]=d; };
  	//@{
    /** @name Implementation of standard functions for 3D points **/
  
  	inline void Zero(){
  		_v[0]=0; _v[1]=0; _v[2]=0; _v[3]=0; };
  
  	inline Point4 ( const S nv[4] ){
  		_v[0]=nv[0]; _v[1]=nv[1]; _v[2]=nv[2]; _v[3]=nv[3]; };
  
  	inline Point4 operator + ( Point4 const & 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 - ( Point4 const & 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 S s ) const {
  		return Point4( _v[0]*s, _v[1]*s , _v[2]*s , _v[3]*s ); }
  
  	inline PointType operator ^ ( PointType const & p ) const {
  		assert(0);
  		return *this;
  	}
  
  	inline Point4 operator / ( const S s ) const {
  		S t=1.0/s;
  		return Point4( _v[0]*t, _v[1]*t , _v[2]*t , _v[3]*t ); }
  
  	inline Point4 operator - () const {
  		return Point4 ( -_v[0], -_v[1] , -_v[2] , -_v[3] ); }
  
  	inline Point4 & operator += ( Point4 const & 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 -= ( Point4 const & 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 S s ) {
  		_v[0] *= s; _v[1] *= s; _v[2] *= s; _v[3] *= s; return *this; }
  
  	inline Point4 & operator /= ( const S s ) {
  		S t=1.0/s; _v[0] *= t; _v[1] *= t;  _v[2] *= t;  _v[3] *= t; return *this; }
  
  	inline S Norm() const {
  		return math::Sqrt( _v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] + _v[3]*_v[3] );}
  
  	template <class PT> static S Norm(const PT &p ) {
  		return math::Sqrt( 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 SquaredNorm() const {
  		return ( _v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] + _v[3]*_v[3] );}
  
  	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 {
      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 {
  		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 {
  		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{
  		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 >  ( 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]); }
  
  	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]); }
  
  	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]); }
  
  	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; }
  		return *this;};
  
  	template <class PT>  PointType & Normalize(const PT &p){
  		PointType n = Norm(); if(n!=0.0) { n=1.0/n;	V(0)*=n;	V(1)*=n;	V(2)*=n; V(3)*=n; }
  		return *this;};
  
  	inline PointType & HomoNormalize(){
  		if (_v[3]!=0.0) {	_v[0] /= W();	_v[1] /= W();	_v[2] /= W(); W()=1.0; }
  		return *this;};
  
  	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]) ) ); }
  
  	inline S NormOne() const {
  		return math::Abs(_v[0])+ math::Abs(_v[1])+math::Max(math::Abs(_v[2]),math::Abs(_v[3]));}
  
  	inline S operator % ( PointType const & p ) const {
  		S t = (*this)*p; /* Area, general formula */
  		return math::Sqrt( SquaredNorm() * p.SquaredNorm() - (t*t) );};
  
  	inline S Sum() const {
  		return _v[0]+_v[1]+_v[2]+_v[3];}
  
  	inline S Max() const {
  		return math::Max( math::Max( _v[0], _v[1] ), math::Max( _v[2], _v[3] )); }
  
  	inline S Min() const {
  		return math::Min( math::Min( _v[0], _v[1] ), math::Min( _v[2], _v[3] )); }
  
  	inline int MaxI() const {
  		int i= (_v[0] < _v[1]) ? 1:0; if (_v[i] < _v[2]) i=2; if (_v[i] < _v[3]) i=3;
  		return i;};
  
  	inline int MinI() const {
  		int i= (_v[0] > _v[1]) ? 1:0; if (_v[i] > _v[2]) i=2; if (_v[i] > _v[3]) i=3;
  		return i;};
  
  	inline PointType & Scale( const PointType & p ) {
  		_v[0] *= p._v[0]; _v[1] *= p._v[1]; _v[2] *= p._v[2]; _v[3] *= p._v[3]; return *this; }
  
  	inline S StableDot ( const PointType & p ) const {
  		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) { 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; }
  
  	 //@}
  };
  
  
  template <class S>
  inline S Angle( Point3<S>  const & p1, Point3<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 siano unitari
  template <class S>
  inline S AngleN( Point3<S>  const & p1, Point3<S>  const & p2 )
  {
  	S w = p1*p2;
  	if(w>1)
  		w = 1;
  	else if(w<-1)
  		w=-1;
    return (S) acos(w);
  }
  
  
  template <int N,class S>
  inline S Norm( Point<N,S> const & p )
  {
      return p.Norm();
  }
  
  template <int N,class S>
  inline S SquaredNorm( Point<N,S> const & p )
  {
    return p.SquaredNorm();
  }
  
  template <int N,class S>
  inline Point<N,S> & Normalize( Point<N,S> & p )
  {
    p.Normalize();
    return p;
  }
  
  template <int N, class S>
  inline S Distance( Point<N,S> const & p1,Point<N,S> const & p2 )
  {
      return (p1-p2).Norm();
  }
  
  template <int N, class S>
  inline S SquaredDistance( Point<N,S> const & p1,Point<N,S> const & p2 )
  {
      return (p1-p2).SquaredNorm();
  }
  
  
  //template <typename S>
  //struct Point2:public Point<2,S>{
  //	inline Point2(){};
  //	inline Point2(Point<2,S> const & p):Point<2,S>(p){} ;
  //	inline Point2( const S a, const S b):Point<2,S>(a,b){};
  //};
  //
  //template <typename S>
  //struct Point3:public Point3<S> {
  //	inline Point3(){};
  //	inline Point3(Point3<S>  const & p):Point3<S> (p){}
  //	inline Point3( const S a, const S b,  const S c):Point3<S> (a,b,c){};
  //};
  //
  //
  //template <typename S>
  //struct Point4:public Point4<S>{
  //	inline Point4(){};
  //	inline Point4(Point4<S> const & p):Point4<S>(p){}
  //	inline Point4( const S a, const S b,  const S c, const S d):Point4<S>(a,b,c,d){};
  //};
  
  typedef Point2<short>  Point2s;
  typedef Point2<int>	  Point2i;
  typedef Point2<float>  Point2f;
  typedef Point2<double> Point2d;
  typedef Point2<short>  Vector2s;
  typedef Point2<int>	  Vector2i;
  typedef Point2<float>  Vector2f;
  typedef Point2<double> Vector2d;
  
  typedef Point3<short>  Point3s;
  typedef Point3<int>	  Point3i;
  typedef Point3<float>  Point3f;
  typedef Point3<double> Point3d;
  typedef Point3<short>  Vector3s;
  typedef Point3<int>	  Vector3i;
  typedef Point3<float>  Vector3f;
  typedef Point3<double> Vector3d;
  
  
  typedef Point4<short>  Point4s;
  typedef Point4<int>	  Point4i;
  typedef Point4<float>  Point4f;
  typedef Point4<double> Point4d;
  typedef Point4<short>  Vector4s;
  typedef Point4<int>	  Vector4i;
  typedef Point4<float>  Vector4f;
  typedef Point4<double> Vector4d;
  
  /*@}*/
  
  //added only for backward compatibility
  template<unsigned int N,typename S>
  struct PointBase : Point<N,S>
  {
  	PointBase()
  		:Point<N,S>()
  	{
  	}
  };
  
  } // end namespace ndim
  } // end namespace vcg
  #endif
  
  