基于C++(MFC)的二維Delaunay三角剖分與Voronoi圖的算法及代碼

一、 Delaunay三角網

1. Delaunay三角網的特性：

（1）空圓性，任一三角形外接圓内部不包含其他點。

（2）最接近：以最近臨的三點形成三角形，且各線段(三角形的邊)皆不相交。

（3）唯一性：不論從區域何處開始建構，最終都将得到一緻的結果。

（4）最優性：任意兩個相鄰三角形形成的凸四邊形的對角線如果可以互換的話，那麼兩個三角形六個内角中最小的角度不會變大。

（5）最規則：如果将三角網中的每個三角形的最小角進行升序排列，則Delaunay三角網的排列得到的數值最大。

（6）區域性：新增、删除、移動某一個頂點時隻會影響臨近的三角形。

（7）具有凸多邊形的外殼：三角網最外層的邊界形成一個凸多邊形的外殼。

1. 算法

Delaunay剖分是一種三角剖分的标準，實作它有多種算法。本次采用Bowyer-Watson算法，算法的基本步驟是：

（1） 構造一個超級三角形，包含所有散點，放入三角形連結清單。

（2）将點集中的散點依次插入，在三角形連結清單中找出其外接圓包含

插入點的三角形（稱為該點的影響三角形），删除影響三角形的公共邊，将插入點同影響三角形的全部頂點連接配接起來，進而完成一個點在Delaunay三角形連結清單中的插入。

（3）根據優化準則對局部新形成的三角形進行優化。将形成的三角形放入Delaunay三角形連結清單。

（4）循環執行上述第2步，直到所有散點插入完畢。

（5）删除和超級三角形有關的三角形，形成凸包。

二、Voronoi Diagram

1. Voronoi圖的定義

   又叫泰森多邊形或Dirichlet圖，它是由一組由連接配接兩鄰點直線的垂直平分線組成的連續多邊形組成。

2. Voronoi圖的特點

   （1）每個V多邊形内有一個生成元；

   （2）每個V多邊形内點到該生成元距離短于到其它生成元距離；

   （3）多邊形邊界上的點到生成此邊界的生成元距離相等；

   （4）鄰接圖形的Voronoi多邊形界線以原鄰接界線作為子集。

3. Voronoi的應用

   在計算幾何學科中的重要地位，由于其根據點集劃分的區域到點的距離最近的特點，其在地理學、氣象學、結晶學、航天、核實體學、機器人等領域具有廣泛的應用。如在障礙物點集中，規避障礙尋找最佳路徑。

4. 建立Voronoi圖的方法

   V圖有着按距離劃分鄰近區域的普遍特性，應用範圍廣。生成V圖的方法很多，常見的有分治法、掃描線算法和Delaunay三角剖分算法。

   本次實驗采用的是Delaunay三角剖分算法。主要是指生成Voronoi圖時先生成其對偶元Delaunay三角網，再找出三角網每一三角形的外接圓圓心，最後連接配接相鄰三角形的外接圓圓心，形成以每一三角形頂點為生成元的多邊形網。如下圖所示（來源于網絡，如有侵權，請聯系删除）：

   

   建立Voronoi圖算法的關鍵是對離散資料點合理地連成三角網，即建構Delaunay三角網。
5. 建立Voronoi圖的步驟

（1）離散點自動建構三角網，即建構Delaunay三角網。對離散點和形成的三角形編号，記錄每個三角形是由哪三個離散點構成的。

（2）計算每個三角形的外接圓圓心，并記錄之。

（3）周遊三角形連結清單，尋找與目前三角形pTri三邊共邊的相鄰三角形TriA，TriB和TriC。

（4）如果找到，則把尋找到的三角形的外心與pTri的外心連接配接，存入維諾邊連結清單中。如果找不到，則求出最外邊的中垂線射線存入維諾邊連結清單中。

（5）周遊結束，所有維諾邊被找到，根據邊畫出維諾圖。

生成效果圖：

相關代碼下載下傳連結：（等課程作業上交完，再分享吧，如有緊急需要，請聯系我）

三、代碼解析

1. 建立項目

   建立一個MFC程式項目，項目名Voronoi2。

2. 資源視圖的菜單設定

   Voronoi圖的生成主要分三步，滑鼠點選螢幕擷取點，根據點構造Delaunay三角形，根據三角形構造Voronoi圖。為了重置視圖，還需要增加一個清除螢幕按鈕。

   （步驟提示：資源視圖—— 在菜單欄——編輯名稱——屬性——設定ID 。）

   

   

   然後給每個菜單添加消息處理程式。在菜單按鈕名右擊添加。消息【滑鼠點選螢幕擷取點】是用來為資料開辟空間的，響應函數在Doc裡面。因為要顯示視圖，是以消息【根據點構造Delaunay三角形】，【根據三角形構造Voronoi圖】響應函數在View類裡面。

   

3. 封裝類

建立類，本程式需要用到點，邊和三角形的資料結構，放在自己建立的類頭檔案裡，自己需要的函數，則是邊程式設計邊添加。

（步驟提示：解決方案資料總管——右擊——添加類——類名CVoronoidiagram——編輯函數）

在Doc檔案裡，給自己建立的類執行個體化對象：

class CVoronoi2Doc : public CDocument
    {
        （省略.。。。。）
    // 實作
    public:
    	CVoronoidiagram * m_CVor;
    	（省略。。。。）
    	}
           

隻有觸發消息【滑鼠點選螢幕擷取點】，才可以使用類CVoronoidiagram，并為資料開辟空間。

在CPP檔案裡，對建立的指針對象，構造函數進行初始化，析構函數進行釋放

滑鼠點選螢幕，擷取點，并且在螢幕顯示出來。先添加一個滑鼠點選螢幕的觸發消息響應。我這篇文章有詳細操作，可以借鑒。

https://blog.csdn.net/weixin_44210987/article/details/90679214

//滑鼠點選螢幕擷取點

具體代碼如下

Vec.h

#ifndef VEC_H
#define VEC_H
/*
Szymon Rusinkiewicz
Princeton University

Vec.h
Class for a constant-length vector

Supports the following operations:
	vec v1;			// Initialized to (0,0,0)
	vec v2(1,2,3);		// Initialized to (1,2,3)
	vec v3(v2);		// Copy constructor
	float farray[3];
	vec v4 = vec(farray);	// Explicit: "v4 = farray" won't work
	Vec<3,double> vd;	// The "vec" used above is Vec<3,float>
	point p1, p2, p3;	// Same as vec

	v3 = v1 + v2;		// Also -, *, /  (all componentwise)
	v3 = 3.5f * v1;		// Also vec * scalar, vec / scalar
				// NOTE: scalar has to be the same type:
				// it won't work to do double * vec<float>
	v1 = min(v2,v3);	// Componentwise min/max
	v1 = sin(v2);		// Componentwise - all the usual functions...
	swap(v1,v2);		// In-place swap

	v3 = v1 DOT v2;		// Actually operator^
	v3 = v1 CROSS v2;	// Actually operator%

	float f = v1[0];	// Subscript
	float *fp = v1;		// Implicit conversion to float *

	f = len(v1);		// Length (also len2 == squared length)
	f = dist(p1, p2);	// Distance (also dist2 == squared distance)
	normalize(v1);		// Normalize (i.e., make it unit length)
				// normalize(vec(0,0,0)) => vec(1,0,0)
	v1 = trinorm(p1,p2,p3); // Normal of triangle

	cout << v1 << endl;	// iostream output in the form (1,2,3)
	cin >> v2;		// iostream input using the same syntax

Also defines the utility functions sqr, cube, sgn, fract, clamp, mix,
step, smoothstep, faceforward, reflect, and refract
*/


// Windows defines these as macros, which prevents us from using the
// type-safe versions from std::, as well as interfering with method defns
#undef min
#undef max


#include <cmath>
#include <iostream>
#include <algorithm>
using std::swap;
using std::sqrt;


// Let gcc optimize conditional branches a bit better...
#ifndef likely
#  if !defined(__GNUC__) || (__GNUC__ == 2 && __GNUC_MINOR__ < 96)
#    define likely(x) (x)
#    define unlikely(x) (x)
#  else
#    define likely(x)   (__builtin_expect((x), 1))
#    define unlikely(x) (__builtin_expect((x), 0))
#  endif
#endif


// Boost-like compile-time assertion checking
template <bool X> struct VEC_STATIC_ASSERTION_FAILURE;
template <> struct VEC_STATIC_ASSERTION_FAILURE<true>
	{ void operator () () {} };
#define VEC_STATIC_CHECK(expr) VEC_STATIC_ASSERTION_FAILURE<bool(expr)>()


template <int D, class T = float>
class Vec {
protected:
	T v[D];

public:
	// Constructor for no arguments.  Everything initialized to 0.
	Vec() { for (int i = 0; i < D; i++) v[i] = T(0); }

	// Uninitialized constructor - meant mostly for internal use
#define VEC_UNINITIALIZED ((void *) 0)
	Vec(void *) {}

	// Constructors for 2-4 arguments
	Vec(T x, T y)
		{ VEC_STATIC_CHECK(D == 2); v[0] = x; v[1] = y; }
	Vec(T x, T y, T z)
		{ VEC_STATIC_CHECK(D == 3); v[0] = x; v[1] = y; v[2] = z; }
	Vec(T x, T y, T z, T w)
		{ VEC_STATIC_CHECK(D == 4); v[0] = x; v[1] = y; v[2] = z; v[3] = w; }

	// Constructor from anything that can be accessed using []
	// Pretty aggressive, so marked as explicit.
	template <class S> explicit Vec(const S &x)
		{ for (int i = 0; i < D; i++) v[i] = T(x[i]); }

	// No destructor or assignment operator needed

	// Array reference and conversion to pointer - no bounds checking
	const T &operator [] (int i) const
		{ return v[i]; }
	T &operator [] (int i)
		{ return v[i]; }
	operator const T * () const
		{ return v; }
	operator const T * ()
		{ return v; }
	operator T * ()
		{ return v; }

	// Member operators
	Vec<D,T> &operator += (const Vec<D,T> &x)
	{
		for (int i = 0; i < D; i++)
#pragma omp atomic
			v[i] += x[i];
		return *this;
	}
	Vec<D,T> &operator -= (const Vec<D,T> &x)
	{
		for (int i = 0; i < D; i++)
#pragma omp atomic
			v[i] -= x[i];
		return *this;
	}
	Vec<D,T> &operator *= (const Vec<D,T> &x)
	{
		for (int i = 0; i < D; i++)
#pragma omp atomic
			v[i] *= x[i];
		return *this;
	}
	Vec<D,T> &operator *= (const T &x)
	{
		for (int i = 0; i < D; i++)
#pragma omp atomic
			v[i] *= x;
		return *this;
	}
	Vec<D,T> &operator /= (const Vec<D,T> &x)
	{
		for (int i = 0; i < D; i++)
#pragma omp atomic
			v[i] /= x[i];
		return *this;
	}
	Vec<D,T> &operator /= (const T &x)
	{
		for (int i = 0; i < D; i++)
#pragma omp atomic
			v[i] /= x;
		return *this;
	}

	// Set each component to min/max of this and the other vector
	Vec<D,T> &min(const Vec<D,T> &x)
	{
#pragma omp critical
		for (int i = 0; i < D; i++)
			if (x[i] < v[i]) v[i] = x[i];
		return *this;
	}
	Vec<D,T> &max(const Vec<D,T> &x)
	{
#pragma omp critical
		for (int i = 0; i < D; i++)
			if (x[i] > v[i]) v[i] = x[i];
		return *this;
	}

	// Outside of class: + - * / % ^ << >>

	// Some partial compatibility with valarrays and vectors
	typedef T value_type;
	size_t size() const
		{ return D; }
	T sum() const
		{ T total = v[0];
		  for (int i = 1; i < D; i++) total += v[i];
		  return total; }
	T avg() const
		{ return sum() / D; }
	T product() const
		{ T total = v[0];
		  for (int i = 1; i < D; i++) total *= v[i];
		  return total; }
	T min() const
		{ T m = v[0];
		  for (int i = 1; i < D; i++)
			if (v[i] < m)  m = v[i];
		  return m; }
	T max() const
		{ T m = v[0];
		  for (int i = 1; i < D; i++)
			if (v[i] > m)  m = v[i];
		  return m; }
	T *begin() { return &(v[0]); }
	const T *begin() const { return &(v[0]); }
	T *end() { return begin() + D; }
	const T *end() const { return begin() + D; }
	void clear() { for (int i = 0; i < D; i++) v[i] = T(0); }
	bool empty() const
		{ for (int i = 0; i < D; i++)
			if (v[i]) return false;
		  return true; }
	Vec<D,T> apply(T func(T)) const
		{ Vec<D,T> result(VEC_UNINITIALIZED);
		  for (int i = 0; i < D; i++) result[i] = func(v[i]);
		  return result; }
	Vec<D,T> apply(T func(const T&)) const
		{ Vec<D,T> result(VEC_UNINITIALIZED);
		  for (int i = 0; i < D; i++) result[i] = func(v[i]);
		  return result; }
};

typedef Vec<3,float> vec;
typedef Vec<3, float> point;
typedef Vec<2, float> point2;
typedef Vec<2, float> vec2;
typedef Vec<3,float> vec3;
typedef Vec<4,float> vec4;
typedef Vec<2,int> ivec2;
typedef Vec<3,int> ivec3;
typedef Vec<4,int> ivec4;


// Nonmember operators that take two Vecs
template <int D, class T>
static inline const Vec<D,T> operator + (const Vec<D,T> &v1, const Vec<D,T> &v2)
{
	Vec<D,T> result(VEC_UNINITIALIZED);
	for (int i = 0; i < D; i++)
		result[i] = v1[i] + v2[i];
	return result;
}

template <int D, class T>
static inline const Vec<D,T> operator - (const Vec<D,T> &v1, const Vec<D,T> &v2)
{
	Vec<D,T> result(VEC_UNINITIALIZED);
	for (int i = 0; i < D; i++)
		result[i] = v1[i] - v2[i];
	return result;
}

template <int D, class T>
static inline const Vec<D,T> operator * (const Vec<D,T> &v1, const Vec<D,T> &v2)
{
	Vec<D,T> result(VEC_UNINITIALIZED);
	for (int i = 0; i < D; i++)
		result[i] = v1[i] * v2[i];
	return result;
}

template <int D, class T>
static inline const Vec<D,T> operator / (const Vec<D,T> &v1, const Vec<D,T> &v2)
{
	Vec<D,T> result(VEC_UNINITIALIZED);
	for (int i = 0; i < D; i++)
		result[i] = v1[i] / v2[i];
	return result;
}

// Dot product in any dimension
template <int D, class T>
static inline const T operator ^ (const Vec<D,T> &v1, const Vec<D,T> &v2)
{
	T sum = v1[0] * v2[0];
	for (int i = 1; i < D; i++)
		sum += v1[i] * v2[i];
	return sum;
}
#define DOT ^

// Cross product - only in 3 dimensions
template <class T>
static inline const Vec<3,T> operator % (const Vec<3,T> &v1, const Vec<3,T> &v2)
{
	return Vec<3,T>(v1[1]*v2[2] - v1[2]*v2[1],
			v1[2]*v2[0] - v1[0]*v2[2],
			v1[0]*v2[1] - v1[1]*v2[0]);
}
// Cross product - only in 2 dimensions
template <class T>
static inline const Vec<2, T> operator % (const Vec<2, T> &v1, const Vec<2, T> &v2)
{
	return  Vec<3,T>(0,0,v1[0] * v2[1] - v1[1] * v2[0]);
}
#define CROSS %


// Component-wise equality and inequality (#include the usual caveats
// about comparing floats for equality...)
template <int D, class T>
static inline bool operator == (const Vec<D,T> &v1, const Vec<D,T> &v2)
{
	for (int i = 0; i < D; i++)
		if (v1[i] != v2[i])
			return false;
	return true;
}

template <int D, class T>
static inline bool operator != (const Vec<D,T> &v1, const Vec<D,T> &v2)
{
	for (int i = 0; i < D; i++)
		if (v1[i] != v2[i])
			return true;
	return false;
}


// Unary operators
template <int D, class T>
static inline const Vec<D,T> &operator + (const Vec<D,T> &v)
{
	return v;
}

template <int D, class T>
static inline const Vec<D,T> operator - (const Vec<D,T> &v)
{
	Vec<D,T> result(VEC_UNINITIALIZED);
	for (int i = 0; i < D; i++)
		result[i] = -v[i];
	return result;
}

template <int D, class T>
static inline bool operator ! (const Vec<D,T> &v)
{
	return v.empty();
}


// Vec/scalar operators
template <int D, class T>
static inline const Vec<D,T> operator * (const T &x, const Vec<D,T> &v)
{
	Vec<D,T> result(VEC_UNINITIALIZED);
	for (int i = 0; i < D; i++)
		result[i] = x * v[i];
	return result;
}

template <int D, class T>
static inline const Vec<D,T> operator * (const Vec<D,T> &v, const T &x)
{
	Vec<D,T> result(VEC_UNINITIALIZED);
	for (int i = 0; i < D; i++)
		result[i] = v[i] * x;
	return result;
}

template <int D, class T>
static inline const Vec<D,T> operator / (const T &x, const Vec<D,T> &v)
{
	Vec<D,T> result(VEC_UNINITIALIZED);
	for (int i = 0; i < D; i++)
		result[i] = x / v[i];
	return result;
}

template <int D, class T>
static inline const Vec<D,T> operator / (const Vec<D,T> &v, const T &x)
{
	Vec<D,T> result(VEC_UNINITIALIZED);
	for (int i = 0; i < D; i++)
		result[i] = v[i] / x;
	return result;
}


// iostream operators
template <int D, class T>
static inline std::ostream &operator << (std::ostream &os, const Vec<D,T> &v)

{
	os << "(";
	for (int i = 0; i < D-1; i++)
		os << v[i] << ", ";
	return os << v[D-1] << ")";
}

template <int D, class T>
static inline std::istream &operator >> (std::istream &is, Vec<D,T> &v)
{
	char c1 = 0, c2 = 0;

	is >> c1;
	if (c1 == '(' || c1 == '[') {
		is >> v[0] >> std::ws >> c2;
		for (int i = 1; i < D; i++) {
			if (c2 == ',')
				is >> v[i] >> std::ws >> c2;
			else
				is.setstate(std::ios::failbit);
		}
	}

	if (c1 == '(' && c2 != ')')
		is.setstate(std::ios::failbit);
	else if (c1 == '[' && c2 != ']')
		is.setstate(std::ios::failbit);

	return is;
}


// Functions on Vecs
template <int D, class T>
static inline void swap(const Vec<D,T> &v1, const Vec<D,T> &v2)
{
	for (int i = 0; i < D; i++)
		swap(v1[i], v2[i]);
}

template <int D, class T>
static inline const T len2(const Vec<D,T> &v)
{
	T l2 = v[0] * v[0];
	for (int i = 1; i < D; i++)
		l2 += v[i] * v[i];
	return l2;
}

template <int D, class T>
static inline const T len(const Vec<D,T> &v)
{
	return sqrt(len2(v));
}

template <int D, class T>
static inline const T dist2(const Vec<D,T> &v1, const Vec<D,T> &v2)
{
	T d2 = sqr(v2[0]-v1[0]);
	for (int i = 1; i < D; i++)
		d2 += sqr(v2[i]-v1[i]);
	return d2;
}

template <int D, class T>
static inline const T dist(const Vec<D,T> &v1, const Vec<D,T> &v2)
{
	return sqrt(dist2(v1,v2));
}

template <int D, class T>
static inline Vec<D,T> normalize(Vec<D,T> &v)
{
	T l = len(v);
	if (unlikely(l <= T(0))) {
		v[0] = T(1);
		for (int i = 1; i < D; i++)
			v[i] = T(0);
		return v;
	}

	l = T(1) / l;
	for (int i = 0; i < D; i++)
		v[i] *= l;

	return v;
}
/*static inline Vec<D,T> normalize(Vec<D,T> &v)
{
	T l = len2(v);
	T xhalf = 0.5*l;
	int f0 = *(int*)&l;
	f0 = 0x5f3759df - (f0 >> 1); // 計算第一個近似根
	l = *(T*)&f0;
	l = l*(1.5- xhalf*l*l); // 牛頓疊代法

	if (unlikely(l <= T(0))) {
		v[0] = T(1);
		for (int i = 1; i < D; i++)
			v[i] = T(0);
		return v;
	}

	for (int i = 0; i < D; i++)
		v[i] *= l;

	return v;
}*/


// Area-weighted triangle face normal
template <class T>
static inline T trinorm(const T &v0, const T &v1, const T &v2)
{
	return (typename T::value_type) 0.5 * ((v1 - v0) CROSS (v2 - v0));
}


// Utility functions for square and cube, to go along with sqrt and cbrt
template <class T>
static inline T sqr(const T &x)
{
	return x*x;
}

template <class T>
static inline T cube(const T &x)
{
	return x*x*x;
}


// Sign of a scalar
template <class T>
static inline T sgn(const T &x)
{
	return (x < T(0)) ? T(-1) : T(1);
}


// Utility functions based on GLSL
template <class T>
static inline T fract(const T &x)
{
	return x - floor(x);
}

template <class T>
static inline T clamp(const T &x, const T &a, const T &b)
{
	return x > a ? x < b ? x : b : a;  // returns a on NaN
}

template <class T, class S>
static inline T mix(const T &x, const T &y, const S &a)
{
	return (S(1)-a) * x + a * y;
}

template <class T>
static inline T step(const T &x, const T &a)
{
	return x < a ? T(0) : T(1);
}

template <class T>
static inline T smoothstep(const T &x, const T &a, const T &b)
{
	if (b <= a) return step(x,a);
	T t = (x - a) / (b - a);
	return t <= T(0) ? T(0) : t >= T(1) ? T(1) : t * t * (T(3) - T(2) * t);
}

template <int D, class T>
static inline T faceforward(const Vec<D,T> &N, const Vec<D,T> &I,
			    const Vec<D,T> &Nref)
{
	return ((Nref DOT I) < T(0)) ? N : -N;
}

template <int D, class T>
static inline T reflect(const Vec<D,T> &I, const Vec<D,T> &N)
{
	return I - (T(2) * (N DOT I)) * N;
}

template <int D, class T>
static inline T refract(const Vec<D,T> &I, const Vec<D,T> &N,
			const T &eta)
{
	T NdotI = N DOT I;
	T k = T(1) - sqr(eta) * (T(1) - sqr(NdotI));
	return (k < T(0)) ? T(0) : eta * I - (eta * NdotI * sqrt(k)) * N;
}


// Generic macros for declaring 1-, 2-, and 3- argument
// componentwise functions on vecs
#define VEC_DECLARE_ONEARG(name) \
 template <int D, class T> \
 static inline Vec<D,T> name(const Vec<D,T> &v) \
 { \
	Vec<D,T> result(VEC_UNINITIALIZED); \
	for (int i = 0; i < D; i++) \
		result[i] = name(v[i]); \
	return result; \
 }

#define VEC_DECLARE_TWOARG(name) \
 template <int D, class T> \
 static inline Vec<D,T> name(const Vec<D,T> &v, const T &w) \
 { \
	Vec<D,T> result(VEC_UNINITIALIZED); \
	for (int i = 0; i < D; i++) \
		result[i] = name(v[i], w); \
	return result; \
 } \
 template <int D, class T> \
 static inline Vec<D,T> name(const Vec<D,T> &v, const Vec<D,T> &w) \
 { \
	Vec<D,T> result(VEC_UNINITIALIZED); \
	for (int i = 0; i < D; i++) \
		result[i] = name(v[i], w[i]); \
	return result; \
 }
#define VEC_DECLARE_THREEARG(name) \
 template <int D, class T> \
 static inline Vec<D,T> name(const Vec<D,T> &v, const T &w, const T &x) \
 { \
	Vec<D,T> result(VEC_UNINITIALIZED); \
	for (int i = 0; i < D; i++) \
		result[i] = name(v[i], w, x); \
	return result; \
 } \
 template <int D, class T> \
 static inline Vec<D,T> name(const Vec<D,T> &v, const Vec<D,T> &w, const Vec<D,T> &x) \
 { \
	Vec<D,T> result(VEC_UNINITIALIZED); \
	for (int i = 0; i < D; i++) \
		result[i] = name(v[i], w[i], x[i]); \
	return result; \
 }

VEC_DECLARE_ONEARG(fabs)
VEC_DECLARE_ONEARG(floor)
VEC_DECLARE_ONEARG(ceil)
VEC_DECLARE_ONEARG(round)
VEC_DECLARE_ONEARG(trunc)
VEC_DECLARE_ONEARG(sin)
VEC_DECLARE_ONEARG(asin)
VEC_DECLARE_ONEARG(cos)
VEC_DECLARE_ONEARG(acos)
VEC_DECLARE_ONEARG(tan)
VEC_DECLARE_ONEARG(atan)
VEC_DECLARE_ONEARG(exp)
VEC_DECLARE_ONEARG(log)
VEC_DECLARE_ONEARG(sqrt)
VEC_DECLARE_ONEARG(sqr)
VEC_DECLARE_ONEARG(cbrt)
VEC_DECLARE_ONEARG(cube)
VEC_DECLARE_ONEARG(sgn)
VEC_DECLARE_TWOARG(min)
VEC_DECLARE_TWOARG(max)
VEC_DECLARE_TWOARG(atan2)
VEC_DECLARE_TWOARG(pow)
VEC_DECLARE_TWOARG(fmod)
VEC_DECLARE_TWOARG(step)
VEC_DECLARE_THREEARG(smoothstep)
VEC_DECLARE_THREEARG(clamp)

#undef VEC_DECLARE_ONEARG
#undef VEC_DECLARE_TWOARG
#undef VEC_DECLARE_THREEARG


// Both valarrays and GLSL use abs() on a vector to mean fabs().
// Let's be compatible...
template <int D, class T>
static inline Vec<D,T> abs(const Vec<D,T> &v)
{
	return fabs(v);
}

#endif
           

Voronoi2Doc.h

// Voronoi2Doc.h : CVoronoi2Doc 類的接口
//


#pragma once
#include "Voronoidiagram.h"

class CVoronoi2Doc : public CDocument
{
protected: // 僅從序列化建立
	CVoronoi2Doc();
	DECLARE_DYNCREATE(CVoronoi2Doc)

// 特性
public:
	
// 操作
public:

// 重寫
public:
	virtual BOOL OnNewDocument();
	virtual void Serialize(CArchive& ar);
#ifdef SHARED_HANDLERS
	virtual void InitializeSearchContent();
	virtual void OnDrawThumbnail(CDC& dc, LPRECT lprcBounds);
#endif // SHARED_HANDLERS

// 實作
public:
	CVoronoidiagram * m_CVor;

	virtual ~CVoronoi2Doc();
#ifdef _DEBUG
	virtual void AssertValid() const;
	virtual void Dump(CDumpContext& dc) const;
#endif

protected:

// 生成的消息映射函數
protected:
	DECLARE_MESSAGE_MAP()

#ifdef SHARED_HANDLERS
	// 用于為搜尋處理程式設定搜尋内容的 Helper 函數
	void SetSearchContent(const CString& value);
#endif // SHARED_HANDLERS
public:
	afx_msg void OnGenerateRandom();
};
           

Doc.cpp

// Voronoi2Doc.cpp : CVoronoi2Doc 類的實作
//

#include "stdafx.h"
// SHARED_HANDLERS 可以在實作預覽、縮略圖和搜尋篩選器句柄的
// ATL 項目中進行定義，并允許與該項目共享文檔代碼。
#ifndef SHARED_HANDLERS
#include "Voronoi2.h"
#endif

#include "Voronoi2Doc.h"

#include <propkey.h>

#ifdef _DEBUG
#define new DEBUG_NEW
#endif

// CVoronoi2Doc

IMPLEMENT_DYNCREATE(CVoronoi2Doc, CDocument)

BEGIN_MESSAGE_MAP(CVoronoi2Doc, CDocument)
	ON_COMMAND(ID_GENERATE_RANDOM, &CVoronoi2Doc::OnGenerateRandom)
END_MESSAGE_MAP()


// CVoronoi2Doc 構造/析構

CVoronoi2Doc::CVoronoi2Doc()
{
	// TODO:  在此添加一次性構造代碼
	m_CVor = NULL;
}

CVoronoi2Doc::~CVoronoi2Doc()
{
	if (m_CVor)
	{
		delete m_CVor;
		m_CVor = NULL;

	}
}

BOOL CVoronoi2Doc::OnNewDocument()
{
	if (!CDocument::OnNewDocument())
		return FALSE;

	// TODO:  在此添加重新初始化代碼
	// (SDI 文檔将重用該文檔)

	return TRUE;
}




// CVoronoi2Doc 序列化

void CVoronoi2Doc::Serialize(CArchive& ar)
{
	if (ar.IsStoring())
	{
		// TODO:  在此添加存儲代碼
	}
	else
	{
		// TODO:  在此添加加載代碼
	}
}

#ifdef SHARED_HANDLERS

// 縮略圖的支援
void CVoronoi2Doc::OnDrawThumbnail(CDC& dc, LPRECT lprcBounds)
{
	// 修改此代碼以繪制文檔資料
	dc.FillSolidRect(lprcBounds, RGB(255, 255, 255));

	CString strText = _T("TODO: implement thumbnail drawing here");
	LOGFONT lf;

	CFont* pDefaultGUIFont = CFont::FromHandle((HFONT) GetStockObject(DEFAULT_GUI_FONT));
	pDefaultGUIFont->GetLogFont(&lf);
	lf.lfHeight = 36;

	CFont fontDraw;
	fontDraw.CreateFontIndirect(&lf);

	CFont* pOldFont = dc.SelectObject(&fontDraw);
	dc.DrawText(strText, lprcBounds, DT_CENTER | DT_WORDBREAK);
	dc.SelectObject(pOldFont);
}

// 搜尋處理程式的支援
void CVoronoi2Doc::InitializeSearchContent()
{
	CString strSearchContent;
	// 從文檔資料設定搜尋内容。
	// 内容部分應由“;”分隔

	// 例如:     strSearchContent = _T("point;rectangle;circle;ole object;")；
	SetSearchContent(strSearchContent);
}

void CVoronoi2Doc::SetSearchContent(const CString& value)
{
	if (value.IsEmpty())
	{
		RemoveChunk(PKEY_Search_Contents.fmtid, PKEY_Search_Contents.pid);
	}
	else
	{
		CMFCFilterChunkValueImpl *pChunk = NULL;
		ATLTRY(pChunk = new CMFCFilterChunkValueImpl);
		if (pChunk != NULL)
		{
			pChunk->SetTextValue(PKEY_Search_Contents, value, CHUNK_TEXT);
			SetChunkValue(pChunk);
		}
	}
}

#endif // SHARED_HANDLERS

// CVoronoi2Doc 診斷

#ifdef _DEBUG
void CVoronoi2Doc::AssertValid() const
{
	CDocument::AssertValid();
}

void CVoronoi2Doc::Dump(CDumpContext& dc) const
{
	CDocument::Dump(dc);
}
#endif //_DEBUG


// CVoronoi2Doc 指令


void CVoronoi2Doc::OnGenerateRandom()
{
	// TODO:  在此添加指令處理程式代碼
	m_CVor = new CVoronoidiagram();

}
           

Vew.h

// Voronoi2View.h : CVoronoi2View 類的接口
//

#pragma once


class CVoronoi2View : public CView
{
protected: // 僅從序列化建立
	CVoronoi2View();
	DECLARE_DYNCREATE(CVoronoi2View)

// 特性
public:
	CVoronoi2Doc* GetDocument() const;

// 操作
public:

// 重寫
public:
	virtual void OnDraw(CDC* pDC);  // 重寫以繪制該視圖
	virtual BOOL PreCreateWindow(CREATESTRUCT& cs);
protected:
	virtual BOOL OnPreparePrinting(CPrintInfo* pInfo);
	virtual void OnBeginPrinting(CDC* pDC, CPrintInfo* pInfo);
	virtual void OnEndPrinting(CDC* pDC, CPrintInfo* pInfo);

// 實作
public:
	virtual ~CVoronoi2View();
#ifdef _DEBUG
	virtual void AssertValid() const;
	virtual void Dump(CDumpContext& dc) const;
#endif

protected:

// 生成的消息映射函數
protected:
	DECLARE_MESSAGE_MAP()
public:
	afx_msg void OnLButtonDown(UINT nFlags, CPoint point);
	afx_msg void OnCleanWindow();
	afx_msg void OnDelaunayTriangle();
	afx_msg void OnGenerateVoronoi();
};

#ifndef _DEBUG  // Voronoi2View.cpp 中的調試版本
inline CVoronoi2Doc* CVoronoi2View::GetDocument() const
   { return reinterpret_cast<CVoronoi2Doc*>(m_pDocument); }
#endif
           

Vew,cpp

// Voronoi2View.cpp : CVoronoi2View 類的實作
//

#include "stdafx.h"
// SHARED_HANDLERS 可以在實作預覽、縮略圖和搜尋篩選器句柄的
// ATL 項目中進行定義，并允許與該項目共享文檔代碼。
#ifndef SHARED_HANDLERS
#include "Voronoi2.h"
#endif

#include "Voronoi2Doc.h"
#include "Voronoi2View.h"
#include "Voronoidiagram.h"

#ifdef _DEBUG
#define new DEBUG_NEW
#endif


// CVoronoi2View

IMPLEMENT_DYNCREATE(CVoronoi2View, CView)

BEGIN_MESSAGE_MAP(CVoronoi2View, CView)
	// 标準列印指令
	ON_COMMAND(ID_FILE_PRINT, &CView::OnFilePrint)
	ON_COMMAND(ID_FILE_PRINT_DIRECT, &CView::OnFilePrint)
	ON_COMMAND(ID_FILE_PRINT_PREVIEW, &CView::OnFilePrintPreview)
	ON_WM_LBUTTONDOWN()
	ON_COMMAND(ID_CLEAN_WINDOW, &CVoronoi2View::OnCleanWindow)
	ON_COMMAND(ID_DELAUNAY_TRIANGLE, &CVoronoi2View::OnDelaunayTriangle)
	ON_COMMAND(ID_GENERATE_VORONOI, &CVoronoi2View::OnGenerateVoronoi)
END_MESSAGE_MAP()

// CVoronoi2View 構造/析構

CVoronoi2View::CVoronoi2View()
{
	// TODO:  在此處添加構造代碼

}

CVoronoi2View::~CVoronoi2View()
{
	
	//CWnd::ReleaseDC();
}

BOOL CVoronoi2View::PreCreateWindow(CREATESTRUCT& cs)
{
	// TODO:  在此處通過修改
	//  CREATESTRUCT cs 來修改視窗類或樣式

	return CView::PreCreateWindow(cs);
}

// CVoronoi2View 繪制

void CVoronoi2View::OnDraw(CDC* /*pDC*/)
{
	CVoronoi2Doc* pDoc = GetDocument();
	ASSERT_VALID(pDoc);
	if (!pDoc)
		return;

	// TODO:  在此處為本機資料添加繪制代碼
	//if (pDoc->m_CVor != nullptr)
	//{
	//	CDC* pDC = GetDC();
	//	//pDoc->m_CVor->drawPoint(pDC);//繪制螢幕上點選的點
	//	pDoc->m_CVor->drawTriangle(pDC);//繪制delaunay三角形
	//}
}


// CVoronoi2View 列印

BOOL CVoronoi2View::OnPreparePrinting(CPrintInfo* pInfo)
{
	// 預設準備
	return DoPreparePrinting(pInfo);
}

void CVoronoi2View::OnBeginPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/)
{
	// TODO:  添加額外的列印前進行的初始化過程
}

void CVoronoi2View::OnEndPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/)
{
	// TODO:  添加列印後進行的清理過程
}


// CVoronoi2View 診斷

#ifdef _DEBUG
void CVoronoi2View::AssertValid() const
{
	CView::AssertValid();
}

void CVoronoi2View::Dump(CDumpContext& dc) const
{
	CView::Dump(dc);
}

CVoronoi2Doc* CVoronoi2View::GetDocument() const // 非調試版本是内聯的
{
	ASSERT(m_pDocument->IsKindOf(RUNTIME_CLASS(CVoronoi2Doc)));
	return (CVoronoi2Doc*)m_pDocument;
}
#endif //_DEBUG


// CVoronoi2View 消息處理程式



//滑鼠點選螢幕擷取點
void CVoronoi2View::OnLButtonDown(UINT nFlags, CPoint point)
{
	// TODO:  在此添加消息處理程式代碼和/或調用預設值
	//擷取點
	CVoronoi2Doc* pDoc = GetDocument();
	ASSERT_VALID(pDoc);
	if (pDoc->m_CVor != nullptr)
	{
		CDC* pDC = GetDC();
		pDoc->m_CVor->getpoint(point); //擷取螢幕點選的點
		pDoc->m_CVor->drawPoint(pDC);//繪制螢幕上點選的點
	}
	CView::OnLButtonDown(nFlags, point);
}

//清理螢幕
void CVoronoi2View::OnCleanWindow()
{
	// TODO:  在此添加指令處理程式代碼
	CVoronoi2Doc* pDoc = GetDocument();
	ASSERT_VALID(pDoc);
	CDC* pDC = GetDC();
	if (pDoc->m_CVor != NULL)
	{
		delete pDoc->m_CVor;
		pDoc->m_CVor = NULL;
		pDoc->UpdateAllViews(NULL);
	}
}

//生成delaunay三角形
void CVoronoi2View::OnDelaunayTriangle()
{
	// TODO:  在此添加指令處理程式代碼
	//Invalidate(false);
	CVoronoi2Doc* pDoc = GetDocument();
	ASSERT_VALID(pDoc);
	if (pDoc->m_CVor != NULL)
	{
		CDC* pDC = GetDC();
		pDoc->m_CVor->setDelaunayTriangle(pDoc->m_CVor->m_pt);//生成delaunay三角形,存在“allTriangle”全局變量中
		pDoc->m_CVor->drawTriangle(pDC);//繪制delaunay三角形
		//pDoc->UpdateAllViews(NULL);
	}

}

// 繪制Voronoi圖
void CVoronoi2View::OnGenerateVoronoi()
{
	// TODO:  在此添加指令處理程式代碼
		CVoronoi2Doc* pDoc = GetDocument();
	ASSERT_VALID(pDoc);
	if (pDoc->m_CVor != NULL)
	{
		CDC* pDC = GetDC();
		pDoc->m_CVor->drawVoronoi(pDC);//繪制Voronoi圖
		
	}
}
           

Voronoidiagram.h

#pragma once
#include "Vec.h"
#include <vector>
#include <algorithm>
using namespace std;

class CVoronoidiagram;
struct Edge;
struct DelaunayTriangle;


//類外定義結構體
struct Edge
{
	point2 pt1, pt2;
	Edge(){ pt1 = NULL; pt2 = NULL; }//預設構造函數
	Edge(point2 pt_1, point2 pt_2)
	{
		pt1 = pt_1;
		pt2 = pt_2;
	}
	
};

//delaunay三角形的資料結構
struct DelaunayTriangle
{
	point2 pts1, pts2, pts3;//三角形三點
	Edge  edge1, edge2, edge3;//三邊
	point2 centerPoint;//外接圓圓心
	double radius;    //外接圓半徑
	int obliqueAngle1, obliqueAngle2, obliqueAngle3;  //鈍角或者銳角
	//三角形外心和外接圓半徑

	void circle_center(point2& center, const point2 &pt1, const point2 &pt2, const point2 &pt3, double& radius)
	{
		double x1, x2, x3, y1, y2, y3;
		double x = 0;
		double y = 0;

		x1 = pt1[0];
		x2 = pt2[0];
		x3 = pt3[0];
		y1 = pt1[1];
		y2 = pt2[1];
		y3 = pt3[1];
		x = ((y2 - y1)*(y3*y3 - y1*y1 + x3*x3 - x1*x1) - (y3 - y1)*(y2*y2 - y1*y1 + x2*x2 - x1*x1)) / (2 * (x3 - x1)*(y2 - y1) - 2 * ((x2 - x1)*(y3 - y1)));
		y = ((x2 - x1)*(x3*x3 - x1*x1 + y3*y3 - y1*y1) - (x3 - x1)*(x2*x2 - x1*x1 + y2*y2 - y1*y1)) / (2 * (y3 - y1)*(x2 - x1) - 2 * ((y2 - y1)*(x3 - x1)));

		center[0] = x;
		center[1] = y;
		radius = sqrt(abs(pt1[0] - x)*abs(pt1[0] - x) + abs(pt1[1] - y) * abs(pt1[1] - y));
	}
	int ObliqueAngle(point2 site1, point2 site2, point2 site3)
	{
		vec2 a, b;
		int obliqueAngle;
		a = site2 - site1;
		b = site3 - site1;
		return (obliqueAngle = a DOT b);
	}
	//構造函數
	DelaunayTriangle(point2 &site1, point2 &site2, point2 &site3)
	{
		pts1 = site1;
		pts2 = site2;
		pts3 = site3;
		circle_center(centerPoint, pts1, pts2, pts3, radius);//構造外接圓圓心以及半徑
		edge1 = Edge(pts2, pts3);
		edge2 = Edge(pts1, pts3);
		edge3 = Edge(pts1, pts2);
		obliqueAngle1 = ObliqueAngle(pts1, pts2, pts3);
		obliqueAngle2 = ObliqueAngle(pts2, pts3, pts1);
		obliqueAngle3 = ObliqueAngle(pts3, pts1, pts2);
	}


};

  //重載==運算符
static inline bool operator == (const Edge & edge1, const Edge & edge2)
{
	if (((edge1.pt1 == edge2.pt1) && (edge1.pt2 == edge2.pt2)) || ((edge1.pt1 == edge2.pt2) && (edge1.pt2 == edge2.pt1)))
		return true;
	return false;
}

static inline bool operator == (const DelaunayTriangle & Triangle1, const DelaunayTriangle & Triangle2)
{
	if ((Triangle1.edge1 == Triangle2.edge1) || (Triangle1.edge2 == Triangle2.edge2) || (Triangle1.edge3 == Triangle2.edge3))
		return true;
	else if ((Triangle1.edge1 == Triangle2.edge1) || (Triangle1.edge2 == Triangle2.edge3) || (Triangle1.edge3 == Triangle2.edge2))
		return true;
	else if ((Triangle1.edge1 == Triangle2.edge2) || (Triangle1.edge2 == Triangle2.edge1) || (Triangle1.edge3 == Triangle2.edge3))
		return true;
	else if ((Triangle1.edge1 == Triangle2.edge2) || (Triangle1.edge2 == Triangle2.edge3) || (Triangle1.edge3 == Triangle2.edge1))
		return true;
	else if ((Triangle1.edge1 == Triangle2.edge3) || (Triangle1.edge2 == Triangle2.edge1) || (Triangle1.edge3 == Triangle2.edge2))
		return true;
	else if ((Triangle1.edge1 == Triangle2.edge3) || (Triangle1.edge2 == Triangle2.edge2) || (Triangle1.edge3 == Triangle2.edge1))
		return true;
	else
	 return false;
}

 
     // 類的定義
class CVoronoidiagram
{
public:
	CVoronoidiagram();
	~CVoronoidiagram();

	//函數聲明
public:
	void getpoint(CPoint _pt);
	static bool sortFun(const point2 &p1, const point2 &p2);//兩點排序

	double distance2Point(const point2 &p1, const point2 &p2);//兩點之間距離
	static void circle_center(point2& center, const point2 &pt1, const point2 &pt2, const point2 &pt3, double& radius);//求三角形的外接圓心
	void SetBigTriangle(vector<DelaunayTriangle> &allTriangles);//設定超級三角形

	void setDelaunayTriangle(vector<point2> &m_pt);//根據點集構造Delaunay三角網
	void drawPoint(CDC* pDC);//繪制螢幕點選的點
	void drawTriangle(CDC* pDC);// 繪圖函數
	void drawVoronoi(CDC* pDC);// 繪制Voronoi圖
	void addNewDelaunayTriangle(vector<DelaunayTriangle> &allTriangles, DelaunayTriangle influenedTri, point2 pointS); //從受影響的三角形連結清單中，形成新的三角形連結清單
	void findCommonEdges(vector<DelaunayTriangle> &influenedTriangles, vector<Edge>& coomonEdges);//找出受影響的三角形的公共邊
	void findCommonEdge(DelaunayTriangle chgTri1, DelaunayTriangle chgTri2, vector<Edge> &edge);//找出兩個三角形的公共邊
	bool siteIsEqual(point2 a, point2 b);//判斷兩點是否相同
	
	void remmoveTrianglesByEdges(vector<DelaunayTriangle> &allTriangless, vector<Edge> &edges);//移除所有公共邊所在的三角形
	
	bool isEdgeOnTriangle(DelaunayTriangle triangel, Edge edge);//判斷邊是否屬于三角形
	void LOP(vector<DelaunayTriangle> &newTriList);//對新形成的三角形進行局部優化
	bool isInCircle(DelaunayTriangle triangle,point2 site);//判斷點是否在三角形外接圓的内部
	void returnVoronoiEdgesFromDelaunayTriangles(vector<Edge> triedges, vector<Edge> &voredges);//根據Delaunay三角形網構造Voronoi圖的邊
	void returnEdgesofTriangleList(const vector<DelaunayTriangle> allTriangle, vector<Edge> &allEdges);//根據三角形容器傳回三角形所有的邊
	void DeleteBigTriEdge(vector<Edge> &allEdge); //删除和超級三角形相連的邊
	void DeleteConnectedBigTriangleByPoint(vector<DelaunayTriangle> &allTriangle);// 删除和超級三角形相連的三角形
	bool notInvector(vector<Edge> ccommonEdges, Edge  edgee); //判斷某邊是否存在于vector中
	bool triangleNotInVector(vector<DelaunayTriangle> Triangles, DelaunayTriangle oneTriangle); //判斷某三角形是否存在于vector中
	bool isPointOnEdge(Edge edge, point2 site);//判斷點是否在邊上
	bool isPointOnTriangle(DelaunayTriangle allTriangle, point2 site);//判斷點是否在三角形上
	point2 findMidPoint(point2 a, point2 b);//找出兩點的中點
	Edge produceRayEdge(point2 start, point2 direction);//根據兩點求以第一個點為起點的射線邊
	//容器對象聲明
public:
	vector<point2> m_pt;     //存儲點
	vector<DelaunayTriangle>::iterator intit;//疊代器
	vector<Edge>::iterator intitS;//疊代器
	vector<DelaunayTriangle> allTriangle;  //存貯三角形

	
	
};
           

.cpp

#include "stdafx.h"
#include "Voronoidiagram.h"
using namespace std;

CVoronoidiagram::CVoronoidiagram(){}
CVoronoidiagram::~CVoronoidiagram(){}



   //根據點集構造Delaunay三角網
   void CVoronoidiagram::setDelaunayTriangle( vector<point2> &m_pt)
   {
	   sort(m_pt.begin(),m_pt.end(),sortFun);
	   allTriangle.clear();
	   SetBigTriangle(allTriangle);//設定超級三角形
	   
	   vector<DelaunayTriangle> influenedTriangles;//受影響的三角形
	   
	   vector<DelaunayTriangle> newTriangles;//新形成的三角形
	   
	   vector<Edge> commonEdgess;//受影響三角形的公共邊
	   commonEdgess.clear();
	   //點排序

	   for (int i = 0; i < m_pt.size(); i++)
	   {
		   influenedTriangles.clear();
		   for (int j = 0; j < allTriangle.size(); j++)
		   {
			   double lengthToCenter;//該點到圓心距離
			   lengthToCenter = distance2Point(allTriangle[j].centerPoint, m_pt[i]);
			   if (lengthToCenter < allTriangle[j].radius)
			   {
				   influenedTriangles.push_back(allTriangle[j]);//添加到受影響的三角形連結清單
				   
				   intit = allTriangle.begin() + j;  //用疊代器查找目前三角形位址
				   allTriangle.erase(intit);//移除目前三角形
				   j--;
			   }
		   }
		   newTriangles.clear();
		   //從受影響的三角形容器中，形成新的三角形連結清單
		   for (int k = 0; k<influenedTriangles.size(); k++)
		   {
			   addNewDelaunayTriangle(newTriangles, influenedTriangles[k], m_pt[i]);
		   }
		   if (influenedTriangles.size() ==1)
		   {
			   intit = influenedTriangles.begin();
			   influenedTriangles.erase(intit);//移除目前三角形
		   }
		   //查找受影響三角形的公共邊
		   if (influenedTriangles.size() > 1)
		   {
			  findCommonEdges(influenedTriangles, commonEdgess);
		   }
		   //将受影響三角形容器中的公共邊所在的新形成三角形全部排除
		   if (commonEdgess.size()> 0)
		   {
			   remmoveTrianglesByEdges(newTriangles, commonEdgess);
		   }
		   commonEdgess.clear();
		   // 對新形成的三角形進行局部優化
		      //LOP(newTriangles);          

		   //将新形成的三角形添加到三角形連結清單中
		   for (int k = 0; k < newTriangles.size(); k++)
		   {
			   allTriangle.push_back(newTriangles[k]);
		   }
	   }
   }


   //将點與受影響的三角形三點連接配接，形成新的三個三角形添加到三角形鍊中
   void CVoronoidiagram::addNewDelaunayTriangle(vector<DelaunayTriangle> &allTriangles, DelaunayTriangle influenedTri, point2 pointS)
   {
	   allTriangles.push_back(DelaunayTriangle(influenedTri.pts1, influenedTri.pts2, pointS));
	   allTriangles.push_back(DelaunayTriangle(influenedTri.pts1, influenedTri.pts3, pointS));
	   allTriangles.push_back(DelaunayTriangle(influenedTri.pts2, influenedTri.pts3, pointS));
   }

   //找出受影響的三角形的公共邊
   void  CVoronoidiagram::findCommonEdges(vector<DelaunayTriangle> &influenedTriangles, vector<Edge>& coomonEdges)
   {
	   vector<Edge> cooomonEdges;
	   for (int i = 0; i < influenedTriangles.size(); i++)
	   {
		   for (int j = i + 1; j < influenedTriangles.size(); j++)
		   {
			   findCommonEdge(influenedTriangles[i], influenedTriangles[j], cooomonEdges);
		   }
	   }
	   for (int k = 0; k < cooomonEdges.size(); k++)
	   {
		   coomonEdges.push_back(cooomonEdges[k]);
	   }

   }

   //找出兩個三角形的公共邊
   void CVoronoidiagram::findCommonEdge(DelaunayTriangle chgTri1, DelaunayTriangle chgTri2, vector<Edge> &edge)
   {
	 
	   vector<point2> commonSites;
	   if (siteIsEqual(chgTri1.pts1, chgTri2.pts1) || siteIsEqual(chgTri1.pts1, chgTri2.pts2) || siteIsEqual(chgTri1.pts1, chgTri2.pts3))
	   {
		   commonSites.push_back(chgTri1.pts1);
	   }
	   if (siteIsEqual(chgTri1.pts2, chgTri2.pts1) || siteIsEqual(chgTri1.pts2, chgTri2.pts2) || siteIsEqual(chgTri1.pts2, chgTri2.pts3))
	   {
		   commonSites.push_back(chgTri1.pts2);
	   }
	   if (siteIsEqual(chgTri1.pts3, chgTri2.pts1) || siteIsEqual(chgTri1.pts3, chgTri2.pts2) || siteIsEqual(chgTri1.pts3, chgTri2.pts3))
	   {
		   commonSites.push_back(chgTri1.pts3);
	   }
	   if (commonSites.size()== 2)
	   {
		   Edge m_edge(commonSites[0], commonSites[1]);
		   edge.push_back(m_edge);
	   }
   }

   //判斷兩點是否相同
   bool CVoronoidiagram::siteIsEqual(point2 a, point2 b)
   {
	   if (a[0] == b[0] && a[1] == b[1])
		   return true;
	   return false;
   }

   //将受影響三角形中的公共邊所在的新三角形從新形成的三角形容器中排除
   void CVoronoidiagram::remmoveTrianglesByEdges(vector<DelaunayTriangle> &allTriangless, vector<Edge> &edges)
   {
	 
	   for (int i = 0; i < allTriangless.size(); i++)
	   {
		   for (int j = 0; j < edges.size(); j++)
		   {
			   if (isEdgeOnTriangle(allTriangless[i], edges[j]))
			   {
				     
					   intit = allTriangless.begin() + i;  //用疊代器查找目前三角形位址
					   allTriangless.erase(intit);//移除目前三角形
					   i--;
					 
			   }
		   }
	   }
	  

   }

   //判斷邊是否屬于三角形
   bool CVoronoidiagram::isEdgeOnTriangle(DelaunayTriangle triangel, Edge edge)
   {
	   int samePointNum = 0;
	   if (siteIsEqual(edge.pt1, triangel.pts1) || siteIsEqual(edge.pt1, triangel.pts2) || siteIsEqual(edge.pt1, triangel.pts3))
		   samePointNum++;
	   if (siteIsEqual(edge.pt2, triangel.pts1) || siteIsEqual(edge.pt2, triangel.pts2) || siteIsEqual(edge.pt2, triangel.pts3))
		   samePointNum++;
	   if (samePointNum == 2)
		   return true;
	   return false;
   }

   //對新形成的三角形進行局部優化
   void CVoronoidiagram::LOP(vector<DelaunayTriangle> &newTriList)
   {
	   vector<Edge> CopyCommonEdges;//拷貝所有找到的公共邊
	   vector<DelaunayTriangle> resultTriList;
	   for (int i = 0; i < newTriList.size(); i++)
	   {
		   for (int j = i + 1; j < newTriList.size(); j++)
		   {
			   resultTriList.clear();
			   vector<Edge>  commonEdge;//需要調整對角線的兩三角形的公共邊
			   commonEdge.clear();
			   point2 anotherPoint;//新對角線的另一點
			   if (isInCircle(newTriList[j], newTriList[i].pts1))//三角形點在外接圓内
			   {
				   //找出兩個三角形的公共邊
				   findCommonEdge(newTriList[i], newTriList[j], commonEdge);
				   //拷貝所有找到的公共邊
				   for (int k = 0; k < commonEdge.size(); k++)
				   {
					   CopyCommonEdges.push_back(commonEdge[k]);
				   }
				   //用換邊法則替換,形成新的三角形
				   if (commonEdge.size()>0)
				   {
					  
					   //找出對角線的另一點
					   if ((siteIsEqual(newTriList[j].pts1, commonEdge[0].pt1) == false) && (siteIsEqual(newTriList[j].pts1, commonEdge[0].pt2) == false))
						   anotherPoint = newTriList[j].pts1;
					   if ((siteIsEqual(newTriList[j].pts2, commonEdge[0].pt1) == false) && (siteIsEqual(newTriList[j].pts2, commonEdge[0].pt2) == false))
						   anotherPoint = newTriList[j].pts2;
					   if ((siteIsEqual(newTriList[j].pts3, commonEdge[0].pt1) == false) && (siteIsEqual(newTriList[j].pts3, commonEdge[0].pt2) == false))
						   anotherPoint = newTriList[j].pts3;
					   //形成兩個新的三角形
					   resultTriList.push_back(DelaunayTriangle(newTriList[i].pts1, anotherPoint, commonEdge[0].pt1));
					   resultTriList.push_back(DelaunayTriangle(newTriList[i].pts1, anotherPoint, commonEdge[0].pt2));
					   //新三角形添加進newTriList
					   for (int m = 0; m < resultTriList.size(); m++)
					   {
						   newTriList.push_back(resultTriList[m]);
					   }
					   resultTriList.clear();
					   //舊三角形從newTriList清除
					   remmoveTrianglesByEdges(newTriList, CopyCommonEdges);
				   }
			   }

			   if (isInCircle(newTriList[j], newTriList[i].pts2))//三角形點在外接圓内
			   {
				   //找出兩個三角形的公共邊
				   findCommonEdge(newTriList[i], newTriList[j], commonEdge);
				   //拷貝所有找到的公共邊
				   for (int k = 0; k < commonEdge.size(); k++)
				   {
					   CopyCommonEdges.push_back(commonEdge[k]);
				   }
				   if (commonEdge.size() > 0)
				   {
					  
					   //找出對角線的另一點
					   if ((siteIsEqual(newTriList[j].pts1, commonEdge[0].pt1) == false) && (siteIsEqual(newTriList[j].pts1, commonEdge[0].pt2) == false))
						   anotherPoint = newTriList[j].pts1;
					   if ((siteIsEqual(newTriList[j].pts2, commonEdge[0].pt1) == false) && (siteIsEqual(newTriList[j].pts2, commonEdge[0].pt2) == false))
						   anotherPoint = newTriList[j].pts2;
					   if ((siteIsEqual(newTriList[j].pts3, commonEdge[0].pt1) == false) && (siteIsEqual(newTriList[j].pts3, commonEdge[0].pt2) == false))
						   anotherPoint = newTriList[j].pts3;
					   //形成兩個新的三角形
					   resultTriList.push_back(DelaunayTriangle(newTriList[i].pts2, anotherPoint, commonEdge[0].pt1));
					   resultTriList.push_back(DelaunayTriangle(newTriList[i].pts2, anotherPoint, commonEdge[0].pt2));
					   //新三角形添加進newTriList
					   for (int m = 0; m < resultTriList.size(); m++)
					   {
						   newTriList.push_back(resultTriList[m]);
					   }
					   resultTriList.clear();
					   //舊三角形從newTriList清除
					   remmoveTrianglesByEdges(newTriList, CopyCommonEdges);
				   }
			   }

			   if (isInCircle(newTriList[j], newTriList[i].pts3))//三角形點在外接圓内
			   {
				   //找出兩個三角形的公共邊
				   findCommonEdge(newTriList[i], newTriList[j], commonEdge);
				   //拷貝所有找到的公共邊
				   for (int k = 0; k < commonEdge.size(); k++)
				   {
					   CopyCommonEdges.push_back(commonEdge[k]);
				   }
				   if (commonEdge.size()>0)
				   {
					   //找出對角線的另一點
					   if ((siteIsEqual(newTriList[j].pts1, commonEdge[0].pt1) == false) && (siteIsEqual(newTriList[j].pts1, commonEdge[0].pt2) == false))
						   anotherPoint = newTriList[j].pts1;
					   if ((siteIsEqual(newTriList[j].pts2, commonEdge[0].pt1) == false) && (siteIsEqual(newTriList[j].pts2, commonEdge[0].pt2) == false))
						   anotherPoint = newTriList[j].pts2;
					   if ((siteIsEqual(newTriList[j].pts3, commonEdge[0].pt1) == false) && (siteIsEqual(newTriList[j].pts3, commonEdge[0].pt2) == false))
						   anotherPoint = newTriList[j].pts3;
					   //形成兩個新的三角形
					   resultTriList.push_back(DelaunayTriangle(newTriList[i].pts3, anotherPoint, commonEdge[0].pt1));
					   resultTriList.push_back(DelaunayTriangle(newTriList[i].pts3, anotherPoint, commonEdge[0].pt2));
					   //新三角形添加進newTriList
					   for (int m = 0; m < resultTriList.size(); m++)
					   {
						   newTriList.push_back(resultTriList[m]);
					   }
					   resultTriList.clear();
					   //舊三角形從newTriList清除
					   remmoveTrianglesByEdges(newTriList, CopyCommonEdges);
				   }
			   }
		   }
	   }
	   
	   
   }

   //判斷點是否在三角形外接圓的内部
   bool CVoronoidiagram::isInCircle(DelaunayTriangle triangle, point2 site)
   {
	   double lengthToCenter;//該點到圓心距離
	   lengthToCenter = distance2Point(triangle.centerPoint, site);
	   if (lengthToCenter < triangle.radius)
	   {
		   return true;
	   }
	   return false;
   }

   // 螢幕上擷取點
   void CVoronoidiagram::getpoint(CPoint _pt)
   {
	   point2 pt(_pt.x, _pt.y);
	   m_pt.push_back(pt);
   }

   //排序函數
   bool CVoronoidiagram::sortFun(const point2 &p1, const point2 &p2)
   {
	   if (p1[1] < p2[1]) return true;
	   else if (p1[1] == p2[1])
	   if (p1[0] < p2[0])return true;
	   else return false;
	   else return false;
   }

   //求兩點之間距離
   double CVoronoidiagram::distance2Point(const point2 &p1, const point2 &p2)
   {
	   double value = sqrt(abs(p1[0] - p2[0])*abs(p1[0] - p2[0]) + abs(p1[1] - p2[1])*abs(p1[1] - p2[1]));
	   return value;
   }


   //三角形的外接圓心
   void CVoronoidiagram::circle_center(point2& center, const point2 &pt1, const point2 &pt2, const point2 &pt3, double& radius)
   {
	   double x1, x2, x3, y1, y2, y3;
	   double x = 0;
	   double y = 0;

	   x1 = pt1[0];
	   x2 = pt2[0];
	   x3 = pt3[0];
	   y1 = pt1[1];
	   y2 = pt2[1];
	   y3 = pt3[1];
	   x = ((y2 - y1)*(y3*y3 - y1*y1 + x3*x3 - x1*x1) - (y3 - y1)*(y2*y2 - y1*y1 + x2*x2 - x1*x1)) / (2 * (x3 - x1)*(y2 - y1) - 2 * ((x2 - x1)*(y3 - y1)));
	   y = ((x2 - x1)*(x3*x3 - x1*x1 + y3*y3 - y1*y1) - (x3 - x1)*(x2*x2 - x1*x1 + y2*y2 - y1*y1)) / (2 * (y3 - y1)*(x2 - x1) - 2 * ((y2 - y1)*(x3 - x1)));

	   center[0] = x;
	   center[1] = y;
	   radius = sqrt(abs(pt1[0] - x)*abs(pt1[0] - x) + abs(pt1[1] - y) * abs(pt1[1] - y));
   }

   //設定超級三角形
   void CVoronoidiagram::SetBigTriangle(vector<DelaunayTriangle> &allTriangles)
   {
	   //将超級三角形的三點添加到三角形網中
	   point2 A(250, -5000);
	   point2 B(-5000, 4000);
	   point2 C(5000, 4000);
	   DelaunayTriangle dt(A, B, C);
	   allTriangles.push_back(dt); //将超級三角形放入容器中
   }
   //繪制螢幕點選的點
   void CVoronoidiagram::drawPoint(CDC* pDC)
   {
	   for (unsigned int i = 0; i < m_pt.size(); i++)
	   {
		   pDC->SetPixel(CPoint(m_pt[i][0], m_pt[i][1]), RGB(100, 100, 100));
		   CBrush brush;
		   brush.CreateSysColorBrush(HS_BDIAGONAL);
		   CBrush* oldBr = pDC->SelectObject(&brush);
		   pDC->Ellipse(m_pt[i][0] - 5, m_pt[i][1] + 5, m_pt[i][0] + 5, m_pt[i][1] - 5);
		   pDC->SelectObject(oldBr);
	   }
   }

   //繪制delaunay三角形
   void CVoronoidiagram::drawTriangle(CDC* pDC)
   {
	   
	   vector<Edge> edgess;
	   edgess.clear();
	   DeleteConnectedBigTriangleByPoint(allTriangle);//删除超級三角形的有關的三角形
	   returnEdgesofTriangleList(allTriangle, edgess);//獲得所有的邊
	   //DeleteBigTriEdge(edgess);//删除超級三角形的有關的邊
	   CPen NewPen, *pOldPen;
	   NewPen.CreatePen(PS_DOT, 2, RGB(0, 255, 0));
	   pOldPen = pDC->SelectObject(&NewPen);
	   for (int i = 0; i < edgess.size(); i++)
	   {
		   
		   pDC->MoveTo(edgess[i].pt1[0], edgess[i].pt1[1]);
		   pDC->LineTo(edgess[i].pt2[0], edgess[i].pt2[1]);
	   }
	   pDC->SelectObject(pOldPen);
	   NewPen.DeleteObject();
   }

   // 繪制Voronoi圖
   void CVoronoidiagram::drawVoronoi(CDC* pDC)
   {
	   vector<Edge> triedges;//三角形邊
	   triedges.clear();
	   vector<Edge> voredges;//voronoi邊
	   voredges.clear();
	   //DeleteConnectedBigTriangleByPoint(allTriangle);//删除超級三角形的有關的三角形
	   returnEdgesofTriangleList(allTriangle, triedges);//獲得所有的邊
	   //DeleteBigTriEdge(triedges);//删除超級三角形的有關的邊
	   returnVoronoiEdgesFromDelaunayTriangles(triedges, voredges);//生成voronoi圖的邊

	   CPen NewPen, *pOldPen;
	   NewPen.CreatePen(PS_SOLID, 4, RGB(255, 0, 255));
	   pOldPen = pDC->SelectObject(&NewPen);
	   for (int i = 0; i < voredges.size(); i++)
	   {
		   pDC->MoveTo(voredges[i].pt1[0], voredges[i].pt1[1]);
		   pDC->LineTo(voredges[i].pt2[0], voredges[i].pt2[1]);
	   }
	   pDC->SelectObject(pOldPen);
	   NewPen.DeleteObject();
   }

   //根據Delaunay三角形網構造Voronoi圖的邊
   void CVoronoidiagram::returnVoronoiEdgesFromDelaunayTriangles(vector<Edge> triedges, vector<Edge> &voronoiEdgeList)
   {
	   vector<Edge> voronoiRayEdgeList; //voronoi圖外圍射線邊
	   voronoiRayEdgeList.clear();
	   voronoiEdgeList.clear();//voronoiEdgeList voronoi圖的邊
	   vector<Edge> neighborEdgeList;   //三角形所有鄰接邊集合
	   vector<Edge> neighboroneEdgeList;    //三角形一個鄰接邊
	   neighborEdgeList.clear();
	   point2  midpoint;
	   Edge rayEdge;
	   for (int i = 0; i < allTriangle.size(); i++)
	   {
		   neighborEdgeList.clear();
		   for (int j = 0; j < allTriangle.size(); j++)//為了找出鄰接三角形數為2的三角形，即最外邊的三角形，循環隻能從0開始
		   {
			   if (j != i)//不與自身比較
			   {
				   neighboroneEdgeList.clear();
				   findCommonEdge(allTriangle[i], allTriangle[j], neighboroneEdgeList);
				   if (neighboroneEdgeList.size()>0)
				   {
					   //構造Voronoi邊
					   Edge voronoiEdge = Edge(allTriangle[i].centerPoint, allTriangle[j].centerPoint);
					   if (notInvector(voronoiEdgeList, voronoiEdge))
					   {
						   voronoiEdgeList.push_back(voronoiEdge);
					   }
					   neighborEdgeList.push_back(neighboroneEdgeList[0]);
				   }
			   }
		    }
		   if (neighborEdgeList.size() <3)//表示此三角形是獨立三角形，Voronoi邊需要3條射線
		     {
			    
			   if (notInvector(neighborEdgeList, allTriangle[i].edge1))
			   {
				   //這條邊是邊界邊 
				   if (allTriangle[i].obliqueAngle1>0)
				   {
					   midpoint = findMidPoint(allTriangle[i].pts2, allTriangle[i].pts3);
					   rayEdge = produceRayEdge(allTriangle[i].centerPoint, midpoint);
					   voronoiEdgeList.push_back(rayEdge);
				   }
				   else if (allTriangle[i].obliqueAngle1 < 0)
				   {
					   midpoint = findMidPoint(allTriangle[i].pts2, allTriangle[i].pts3);
					   rayEdge = produceRayEdge(allTriangle[i].centerPoint, allTriangle[i].centerPoint - midpoint + allTriangle[i].centerPoint);
					   voronoiEdgeList.push_back(rayEdge);
				   }
			   }
			   if (notInvector(neighborEdgeList, allTriangle[i].edge2))
			   {
				   //這條邊是邊界邊 
				   if (allTriangle[i].obliqueAngle2>0)
				   {
					   midpoint = findMidPoint(allTriangle[i].pts1, allTriangle[i].pts3);
					   rayEdge = produceRayEdge(allTriangle[i].centerPoint, midpoint);
					   voronoiEdgeList.push_back(rayEdge);
				   }
				   else if (allTriangle[i].obliqueAngle2 < 0)
				   {
					   midpoint = findMidPoint(allTriangle[i].pts1, allTriangle[i].pts3);
					   rayEdge = produceRayEdge(allTriangle[i].centerPoint, allTriangle[i].centerPoint - midpoint + allTriangle[i].centerPoint);
					   voronoiEdgeList.push_back(rayEdge);
				   }
			   }
			   if (notInvector(neighborEdgeList, allTriangle[i].edge3))
			   {
				   //這條邊是邊界邊 
				   if (allTriangle[i].obliqueAngle3>0)
				   {
					   midpoint = findMidPoint(allTriangle[i].pts1, allTriangle[i].pts2);
					   rayEdge = produceRayEdge(allTriangle[i].centerPoint, midpoint);
					   voronoiEdgeList.push_back(rayEdge);
				   }
				   else if (allTriangle[i].obliqueAngle3 < 0)
				   {
					   midpoint = findMidPoint(allTriangle[i].pts1, allTriangle[i].pts2);
					   rayEdge = produceRayEdge(allTriangle[i].centerPoint, allTriangle[i].centerPoint - midpoint + allTriangle[i].centerPoint);
					   voronoiEdgeList.push_back(rayEdge);
				   }
			   }
			 
			 }
	   }
	   
   }
   //根據三角形容器傳回三角形所有的邊
   void CVoronoidiagram::returnEdgesofTriangleList( const vector<DelaunayTriangle> allTriangle, vector<Edge> &allEdges)
   {
	  
	   for (int i = 0; i < allTriangle.size(); i++)
	   {
		   Edge edge1 = Edge(allTriangle[i].pts1, allTriangle[i].pts2);
		   Edge edge2 = Edge(allTriangle[i].pts1, allTriangle[i].pts3);
		   Edge edge3 = Edge(allTriangle[i].pts2, allTriangle[i].pts3);
		   if (notInvector(allEdges, edge1))
			   allEdges.push_back(edge1);
		   if (notInvector(allEdges, edge2))
			   allEdges.push_back(edge2);
		   if (notInvector(allEdges, edge3))
			   allEdges.push_back(edge3);
	   }
   }

   //删除和超級三角形相連的邊
   void CVoronoidiagram::DeleteBigTriEdge(vector<Edge> &allEdge)
   {
	   point2 A(250, -5000);
	   point2 B(-5000, 4000);
	   point2 C(5000, 4000);
	   vector <point2> bigpoint{A,B,C};
	   for (int j = 0; j < bigpoint.size();j++)
	   {
		   for (int i = 0; i < allEdge.size(); i++)
		   {
			   if (isPointOnEdge(allEdge[i], bigpoint[j]))//找到含有超級三角形的點的邊
			   {
				   intitS = allEdge.begin() + i;  //用疊代器查找目前三角形位址
				   allEdge.erase(intitS);//移除目前三角形
				   i--;
			   }

		   }
	   }
   }

   // 删除和超級三角形相連的三角形
   void CVoronoidiagram::DeleteConnectedBigTriangleByPoint(vector<DelaunayTriangle>& allTriangle)
   {
		   point2 A(250, -5000);
		   point2 B(-5000, 4000);
		   point2 C(5000, 4000);
		   vector <point2> bigpoint{ A, B, C };
		   for (int j = 0; j < bigpoint.size(); j++)
		   {
			   for (int i = 0; i < allTriangle.size(); i++)
			   {
				   if (isPointOnTriangle(allTriangle[i], bigpoint[j]))//找到含有超級三角形的點的邊
				   {
					   intit = allTriangle.begin() + i;  //用疊代器查找目前三角形位址
					   allTriangle.erase(intit);//移除目前三角形
					   i--;
				   }

			   }
		   }
	   }

   //判斷某邊是否存在于vector中
   bool CVoronoidiagram::notInvector(vector<Edge> ccommonEdges,  Edge  edgee)
   {
	   vector<Edge>::iterator ret;
	   ret = find(ccommonEdges.begin(), ccommonEdges.end(), edgee);//注：Edge是自己定義的結構體，需要為該類型重載==操作符，再用find
	   if (ret == ccommonEdges.end())
		   return true;
	   else
		   return false;
	    
   }

   //判斷某三角形是否存在于vector中
   bool CVoronoidiagram::triangleNotInVector(vector<DelaunayTriangle> Triangles, DelaunayTriangle oneTriangle)
   {
	   vector<DelaunayTriangle>::iterator ret;
	   ret = find(Triangles.begin(), Triangles.end(), oneTriangle);//注：DelaunayTriangle是自己定義的結構體，需要為該類型重載==操作符，再用find
	   if (ret == Triangles.end())
		   return true;
	   else
		   return false;

   }

  

   //判斷點是否在邊上
   bool CVoronoidiagram::isPointOnEdge(Edge edge, point2 site)
   {
	   if (siteIsEqual(site, edge.pt1) || siteIsEqual(site, edge.pt2))
		   return true;
	   return false;
   }
   //判斷點是否在三角形上
   bool CVoronoidiagram::isPointOnTriangle(DelaunayTriangle allTriangle, point2 site)
   {
	   if (siteIsEqual(site, allTriangle.pts1) || siteIsEqual(site, allTriangle.pts2) || siteIsEqual(site, allTriangle.pts3))
		   return true;
	   return false;
   }
  
   //找出兩點的中點
   point2 CVoronoidiagram::findMidPoint(point2 a, point2 b)
   {
	   return  point2((a[0] + b[0]) / 2.0, (a[1] + b[1]) / 2.0);
	   
   }

   //根據兩點求以第一個點為起點的射線邊
   Edge CVoronoidiagram::produceRayEdge(point2 start, point2 direction)//produceRayEdge(allTriangle[i].centerPoint, midpoint);
   {
	   point2 end;
	   end[0]= 100 * (direction[0] - start[0]) + start[0];//找出射線方向的較大的x終點
	   end[1] = (direction[1] - start[1]) * (end[0] - start[0]) / (direction[0] - start[0]) + start[1];
	   return  Edge(start, end);
   }