基于C++(MFC)的二维Delaunay三角剖分与Voronoi图的算法及代码

一、 Delaunay三角网

Delaunay三角网的特性：

（1）空圆性，任一三角形外接圆内部不包含其他点。
（2）最接近：以最近临的三点形成三角形，且各线段(三角形的边)皆不相交。
（3）唯一性：不论从区域何处开始构建，最终都将得到一致的结果。
（4）最优性：任意两个相邻三角形形成的凸四边形的对角线如果可以互换的话，那么两个三角形六个内角中最小的角度不会变大。
（5）最规则：如果将三角网中的每个三角形的最小角进行升序排列，则Delaunay三角网的排列得到的数值最大。
（6）区域性：新增、删除、移动某一个顶点时只会影响临近的三角形。
（7）具有凸多边形的外壳：三角网最外层的边界形成一个凸多边形的外壳。

算法

Delaunay剖分是一种三角剖分的标准，实现它有多种算法。本次采用Bowyer-Watson算法，算法的基本步骤是：
（1）构造一个超级三角形，包含所有散点，放入三角形链表。
（2）将点集中的散点依次插入，在三角形链表中找出其外接圆包含
插入点的三角形（称为该点的影响三角形），删除影响三角形的公共边，将插入点同影响三角形的全部顶点连接起来，从而完成一个点在Delaunay三角形链表中的插入。
（3）根据优化准则对局部新形成的三角形进行优化。将形成的三角形放入Delaunay三角形链表。
（4）循环执行上述第2步，直到所有散点插入完毕。
（5）删除和超级三角形有关的三角形，形成凸包。

二、Voronoi Diagram

Voronoi图的定义
又叫泰森多边形或Dirichlet图，它是由一组由连接两邻点直线的垂直平分线组成的连续多边形组成。
Voronoi图的特点
（1）每个V多边形内有一个生成元；
（2）每个V多边形内点到该生成元距离短于到其它生成元距离；
（3）多边形边界上的点到生成此边界的生成元距离相等；
（4）邻接图形的Voronoi多边形界线以原邻接界线作为子集。
Voronoi的应用
在计算几何学科中的重要地位，由于其根据点集划分的区域到点的距离最近的特点，其在地理学、气象学、结晶学、航天、核物理学、机器人等领域具有广泛的应用。如在障碍物点集中，规避障碍寻找最佳路径。
建立Voronoi图的方法
V图有着按距离划分邻近区域的普遍特性，应用范围广。生成V图的方法很多，常见的有分治法、扫描线算法和Delaunay三角剖分算法。
本次实验采用的是Delaunay三角剖分算法。主要是指生成Voronoi图时先生成其对偶元Delaunay三角网，再找出三角网每一三角形的外接圆圆心，最后连接相邻三角形的外接圆圆心，形成以每一三角形顶点为生成元的多边形网。如下图所示（来源于网络，如有侵权，请联系删除）：

建立Voronoi图算法的关键是对离散数据点合理地连成三角网，即构建Delaunay三角网。
建立Voronoi图的步骤

（1）离散点自动构建三角网，即构建Delaunay三角网。对离散点和形成的三角形编号，记录每个三角形是由哪三个离散点构成的。
（2）计算每个三角形的外接圆圆心，并记录之。
（3）遍历三角形链表，寻找与当前三角形pTri三边共边的相邻三角形TriA，TriB和TriC。
（4）如果找到，则把寻找到的三角形的外心与pTri的外心连接，存入维诺边链表中。如果找不到，则求出最外边的中垂线射线存入维诺边链表中。
（5）遍历结束，所有维诺边被找到，根据边画出维诺图。

生成效果图：

相关代码下载链接：（等课程作业上交完，再分享吧，如有紧急需要，请联系我）
三、代码解析

新建项目
新建一个MFC程序项目，项目名Voronoi2。
资源视图的菜单设置
Voronoi图的生成主要分三步，鼠标点击屏幕获取点，根据点构造Delaunay三角形，根据三角形构造Voronoi图。为了重置视图，还需要增加一个清除屏幕按钮。
（步骤提示：资源视图—— 在菜单栏——编辑名称——属性——设置ID 。）

然后给每个菜单添加消息处理程序。在菜单按钮名右击添加。消息【鼠标点击屏幕获取点】是用来为数据开辟空间的，响应函数在Doc里面。因为要显示视图，所以消息【根据点构造Delaunay三角形】，【根据三角形构造Voronoi图】响应函数在View类里面。
封装类

新建类，本程序需要用到点，边和三角形的数据结构，放在自己新建的类头文件里，自己需要的函数，则是边编程边添加。
（步骤提示：解决方案资源管理器——右击——添加类——类名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 UniversityVec.h
Class for a constant-length vectorSupports the following operations:vec v1;         // Initialized to (0,0,0)vec v2(1,2,3);     // Initialized to (1,2,3)vec v3(v2);        // Copy constructorfloat farray[3];vec v4 = vec(farray);   // Explicit: "v4 = farray" won't workVec<3,double> vd;    // The "vec" used above is Vec<3,float>point p1, p2, p3;    // Same as vecv3 = 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/maxv1 = sin(v2);      // Componentwise - all the usual functions...swap(v1,v2);       // In-place swapv3 = v1 DOT v2;        // Actually operator^v3 = v1 CROSS v2; // Actually operator%float f = v1[0];  // Subscriptfloat *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 trianglecout << v1 << endl; // iostream output in the form (1,2,3)cin >> v2;      // iostream input using the same syntaxAlso 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 argumentsVec(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 checkingconst 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 operatorsVec<D,T> &operator += (const Vec<D,T> &x){for (int i = 0; i < D; i++)
#pragma omp atomicv[i] += x[i];return *this;}Vec<D,T> &operator -= (const Vec<D,T> &x){for (int i = 0; i < D; i++)
#pragma omp atomicv[i] -= x[i];return *this;}Vec<D,T> &operator *= (const Vec<D,T> &x){for (int i = 0; i < D; i++)
#pragma omp atomicv[i] *= x[i];return *this;}Vec<D,T> &operator *= (const T &x){for (int i = 0; i < D; i++)
#pragma omp atomicv[i] *= x;return *this;}Vec<D,T> &operator /= (const Vec<D,T> &x){for (int i = 0; i < D; i++)
#pragma omp atomicv[i] /= x[i];return *this;}Vec<D,T> &operator /= (const T &x){for (int i = 0; i < D; i++)
#pragma omp atomicv[i] /= x;return *this;}// Set each component to min/max of this and the other vectorVec<D,T> &min(const Vec<D,T> &x){
#pragma omp criticalfor (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 criticalfor (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 vectorstypedef 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;elseis.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_HANDLERSvirtual void InitializeSearchContent();virtual void OnDrawThumbnail(CDC& dc, LPRECT lprcBounds);
#endif // SHARED_HANDLERS// 实现
public:CVoronoidiagram * m_CVor;virtual ~CVoronoi2Doc();
#ifdef _DEBUGvirtual void AssertValid() const;virtual void Dump(CDumpContext& dc) const;
#endifprotected:// 生成的消息映射函数
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// CVoronoi2DocIMPLEMENT_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 onceclass 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 _DEBUGvirtual void AssertValid() const;virtual void Dump(CDumpContext& dc) const;
#endifprotected:// 生成的消息映射函数
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// CVoronoi2ViewIMPLEMENT_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;elsereturn 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));//新三角形添加进newTriListfor (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));//新三角形添加进newTriListfor (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));//新三角形添加进newTriListfor (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是自己定义的结构体，需要为该类型重载==操作符，再用findif (ret == ccommonEdges.end())return true;elsereturn false;}//判断某三角形是否存在于vector中bool CVoronoidiagram::triangleNotInVector(vector<DelaunayTriangle> Triangles, DelaunayTriangle oneTriangle){vector<DelaunayTriangle>::iterator ret;ret = find(Triangles.begin(), Triangles.end(), oneTriangle);//注：DelaunayTriangle是自己定义的结构体，需要为该类型重载==操作符，再用findif (ret == Triangles.end())return true;elsereturn 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);}

基于C++(MFC)的二维Delaunay三角剖分与Voronoi图的算法及代码相关推荐

matlab中的delaunay,基于MATLAB 实现二维delaunay 三角剖分
基于MATLAB 实现二维delaunay 三角剖分刘锋涛凡友华 (哈尔滨工业大学深圳研究生院深圳518055) [摘要]在已知凸多边形的顶点坐标的前提情况下,利用MATLAB 中的meshgrid ...
如何使用OpenCV进行Delaunay三角剖分和Voronoi图
图1.左:使用dlib检测到具有标志性建筑的奥巴马总统图像.中心:地标的Delaunay三角剖分.右:对应的Voronoi图. 俄国数学家鲍里斯·尼古拉耶维奇·德劳内(Boris Nikolaevic ...
数字图像处理 Delaunay三角剖分和Voronoi图
一.什么是Delaunay三角剖分? 给定平面中的一组点,三角测量是指将平面细分为三角形,以这些点为顶点.在下图1中,我们在左图像中看到一组地标,在中间图像中看到三角剖分.一组点可以有许多可能的三角剖 ...
计算机图形学 | 基于MFC和二维变换的画图软件
文章目录基于MFC和二维变换的画图软件摘要设计 1 程序总体结构 1.1 总体结构设计 1.1.1 绘图设计 1.1.2 变换设计 2 程序实现 2.1 鼠标绘图的消息映射 2.2 图形绘制 ...
二维Delaunay（德洛内）三角网剖分的matlab实现
二维Delaunay(德洛内)三角网的matlab实现二维Delaunay(德洛内)三角网的matlab实现 1.Delaunay三角网的概述 2.Delaunay三角网的算法 3.Delaunay ...
matlab二维谐振子,基于有限差分法求解的二维谐振子的MATLAB程序如下。哪位大神能帮我做个注明啊，完全看不懂啊，，急...
基于有限差分法求解的二维谐振子的MATLAB程序如下.哪位大神能帮我做个注明啊,完全看不懂啊,,急0 ____丿呆呆丶2017.04.15浏览20次分享举报 tic clc clear L=20; W ...
基于springboot的食品二维码溯源系统
1 简介今天向大家介绍一个帮助往届学生完成的毕业设计项目,基于springboot的食品二维码溯源系统. 计算机毕业生设计,课程设计需要帮助的可以找我 2 设计概要 21世纪是信息化时代,随着信息技 ...
【第 07 章基于主成分分析的人脸二维码识别MATLAB深度学习实战案例】
基于主成分分析的人脸二维码识别MATLAB深度学习实战案例人脸库全套文件资料目录下载链接–>传送门本文全文源码下载[链接–>传送门] 如下分析: 主文件 function varar ...
基于MATLAB的条码二维码识别系统
基于MATLAB的条码二维码识别系统课题介绍本设计研究的是基于数字图像处理的EAN-13条形码识别算法,通过工具平台MATLAB实现.其中图像处理部分是条码识别重要的前期工作,利用MATLAB强大 ...

基于C++(MFC)的二维Delaunay三角剖分与Voronoi图的算法及代码

基于C++(MFC)的二维Delaunay三角剖分与Voronoi图的算法及代码相关推荐

最新文章

热门文章