dpoint.hpp

最新delaunay 算法源代码

/*! \file dpoint.hpp
        \brief d-dimensional point class
	
	A d-dimensional point class which is written
	carefully using templates. It allows for basic
	operations on points in any dimension. Orientation
	tests for 2 and 3 dimensional points are supported
	using 
        <a href="http://www.cs.berkeley.edu/~jrs/">Jonathan's</a> 
        code. This class forms
	the building block of other classes like dplane, dsphere etc.

  \author <a href="www.compgeom.com/~piyush">Piyush Kumar</a>
  \bug    No known bugs.

 */

#ifndef REVIVER_POINT_HPP
#define REVIVER_POINT_HPP


#include "assert.hpp"
#include <iostream>
#include <valarray>
#include <stdio.h>
#include <limits>



//! The reviver namespace signifies the part of the code borrowed from reviver (dpoint.hpp). 
namespace reviver {


// Forward Declaration of the main Point Class
// Eucledian d-dimensional point. The distance
// is L_2

template<typename NumType, unsigned D>
class dpoint;


///////////////////////////////////////////////////////
// Internal number type traits for dpoint
///////////////////////////////////////////////////////

template<typename T> 
class InternalNumberType;

template<>
class InternalNumberType<float>{
public: 
  typedef double __INT;
};


template<>
class InternalNumberType<int>{
public: 
  typedef long long __INT;
};


template<>
class InternalNumberType<double>{
public: 
  typedef double __INT;
};

template<>
class InternalNumberType<long>{
public: 
  typedef long long __INT;
};




///////////////////////////////////////////////////////
// Origin of d-dimensional point
///////////////////////////////////////////////////////
template< typename NumType, unsigned D, unsigned I > struct origin
{
   static inline void eval( dpoint<NumType,D>& p )
   {
	  p[I] = 0.0;
      origin< NumType, D, I-1 >::eval( p );
   }
};


// Partial Template Specialization
template <typename NumType, unsigned D> struct origin<NumType, D, 0>
{
   static inline void eval( dpoint<NumType,D>& p )
   {
	  p[0] = 0.0;
   }
};


//! A structure to compute squared distances between points
/*!
    Uses unrolling of loops using templates.
*/
///////////////////////////////////////////////////////
// Squared Distance of d-dimensional point
///////////////////////////////////////////////////////
template< typename NumType, unsigned D, unsigned I > struct Distance
{
   static inline double eval( const dpoint<NumType,D>& p, const dpoint<NumType,D>& q )
   {
     double sum = ((double)p[I] -  (double)q[I] ) *( (double)p[I] - (double)q[I] );
	  return sum + Distance< NumType, D, I-1 >::eval( p,q );
   }
};


//!  Partial Template Specialization for distance calculations
template <typename NumType, unsigned D> struct Distance<NumType, D, 0>
{
   static inline double eval( const dpoint<NumType,D>& p, const dpoint<NumType,D>& q )
   {
     return ((double) p[0] - (double)q[0] )*( (double)p[0] - (double)q[0] );
   }
};


//! A structure to compute dot product between two points or associated vectors
/*!
    Uses unrolling of loops using templates.
*/
///////////////////////////////////////////////////////
// Dot Product of two d-dimensional points
///////////////////////////////////////////////////////
template< typename __INT, typename NumType, unsigned D, unsigned I > struct DotProd
{
   static inline __INT eval( const dpoint<NumType,D>& p, const dpoint<NumType,D>& q )
   {
	  __INT sum = ( ((__INT)p[I]) * ((__INT)q[I]) );
	  return sum + DotProd< __INT, NumType, D, I-1 >::eval( p,q );
   }
};


//! Partial Template Specialization for dot product calculations
template < typename __INT, typename NumType, unsigned D> struct DotProd<__INT,NumType, D, 0>
{
   static inline __INT eval( const dpoint<NumType,D>& p, const dpoint<NumType,D>& q )
   {
	  return ( ((__INT)p[0]) * ((__INT)q[0]) );
   }
};


///////////////////////////////////////////////////////
// Equality of two d-dimensional points
///////////////////////////////////////////////////////
template< typename NumType, unsigned D, unsigned I > struct IsEqual
{
   static inline bool eval( const dpoint<NumType,D>& p, const dpoint<NumType,D>& q )
   {
	  if( p[I]  != q[I] ) return false;
	  else return IsEqual< NumType, D, I-1 >::eval( p,q );
   }
};


// Partial Template Specialization
template <typename NumType, unsigned D> struct IsEqual<NumType, D, 0>
{
   static inline NumType eval( const dpoint<NumType,D>& p, const dpoint<NumType,D>& q )
   {
	   return (p[0] == q[0])?1:0;
   }
};



//!   Equate two d-dimensional points. 
/*!
    Uses unrolling of loops using templates. 
    A class used to implement operator= for points.
    This class also helps in automatic type conversions 
    of points (with explicit calls for conversion).
*/
template< 
  typename NumType1, 
  typename NumType2, 
  unsigned D, 
  unsigned I 
  > struct Equate
{
   static inline void eval( dpoint<NumType1,D>& p,const dpoint<NumType2,D>& q )
   {
	  p[I]  = q[I];
	  Equate< NumType1, NumType2, D, I-1 >::eval( p,q );
   }
};


//! Partial Template Specialization for Equate
template <
  typename NumType1, 
  typename NumType2,
  unsigned D
  > struct Equate<NumType1,NumType2, D, 0>
{
   static inline void  eval( dpoint<NumType1,D>& p,const dpoint<NumType2,D>& q )
   {
	   p[0] = q[0];
   }
};


//! A structure to add two points
/*!
    Uses unrolling of loops using templates.
*/
///////////////////////////////////////////////////////
// Add two d-dimensional points
///////////////////////////////////////////////////////
template< typename NumType, unsigned D, unsigned I > struct Add
{
   static inline void eval( dpoint<NumType,D>& result, const dpoint<NumType,D>& p, const dpoint<NumType,D>& q )
   {
	  result[I] = p[I]  + q[I];
	  Add< NumType, D, I-1 >::eval( result,p,q );
   }
};


//! Partial Template Specialization for Add structure
template <typename NumType, unsigned D> struct Add<NumType, D, 0>
{
   static inline void eval( dpoint<NumType,D>& result, const dpoint<NumType,D>& p, const dpoint<NumType,D>& q )
   {
	   result[0] = p[0] + q[0];
   }
};


///////////////////////////////////////////////////////
// Subtract two d-dimensional points
///////////////////////////////////////////////////////
// Could actually be done using scalar multiplication
// and addition


// What about unsigned types?
template< typename NumType > 
inline NumType Subtract_nums(const NumType& x, const NumType& y) {
  if(!std::numeric_limits<NumType>::is_signed) {std::cerr << "Exception: Can't subtract unsigned types."; exit(1);}
  return x - y;
}


//!   Subtract two d-dimensional vectors
/*!
      Caution: Do not use on unsigned types.
*/
template< typename NumType, unsigned D, unsigned I > struct Subtract
{
   static inline void eval( dpoint<NumType,D>& result, const dpoint<NumType,D>& p, const dpoint<NumType,D>& q )
   {
     
          result[I] = Subtract_nums(p[I] , q[I]);
	  Subtract< NumType, D, I-1 >::eval( result,p,q );
   }
};


//! Partial Template Specialization for subtraction of points (associated vectors)
template <typename NumType, unsigned D> struct Subtract<NumType, D, 0>
{
   static inline void eval( dpoint<NumType,D>& result, const dpoint<NumType,D>& p, const dpoint<NumType,D>& q )
   {
	   result[0] = Subtract_nums(p[0] , q[0]);
   }
};





//!   Mutiply scalar with d-dimensional point
/*!
      Scalar mulipltication of d-dimensional point
      with a number using template unrolling.
*/
template< typename NumType, unsigned D, unsigned I > struct Multiply
{
   static inline void eval( dpoint<NumType,D>& result, const dpoint<NumType,D>& p, NumType k)
   {
	  result[I] = p[I] * k;
	  Multiply< NumType, D, I-1 >::eval( result,p,k );
   }
};


//! Partial Template Specialization for scalar multiplication
template <typename NumType, unsigned D> struct Multiply<NumType, D, 0>
{
   static inline void eval( dpoint<NumType,D>& result, const dpoint<NumType,D>& p, NumType k )
   {
	   result[0] = p[0] * k;
   }
};



//!  Main d dimensional Point Class
/*!
	-  NumType = Floating Point Type
	-  D       = Dimension of Point
*/
template<typename NumType = double, unsigned D = 3>
class dpoint {

		// Makes Swap operation fast
		NumType  x[D];

public: