本节包含:基础知识、格林厄姆扫描法生成凸包(Python)、Bowyer-Watson法构造delaunay三角网(Python)。

声明:本文中所有机理图均由 F-PolyImage 图片生成工具(geocommunity.vip) 创作完成。
一、泰森多边形的定义

上图中,P1-P5为生成点,构成点集。根据泰森多边形的性质,两生成点的垂直平分线构成Voronoi边界。
(到这里,只需掌握泰森多边形是什么即可)
二、泰森多边形的基本性质
1、边界是两点连线的垂直平分线

2、每个泰森多边形是凸多边形
3、什么是凸包,如何构建?
所有离散点(样本点、站点、种子点)整体形成的最小凸边界。以下图直观解释(每个钉子看为每个生成点,橡皮筋即最小凸边界):

格林厄姆扫描法实现凸包构建:
- 选定原点:y最小,y相同取x最小。
- 按照极角逆时针排序,极角相同按距离排序。
- 栈过滤,叉积如果<=0,则删除顶点。(叉积<=0,表示右拐向。)
import math
# 格林厄姆扫描法生成凸壳
def cross_product(o, a, b):
'''
向量 OA = (a[0] - o[0], a[1] - o[1])
向量 OB = (b[0] - o[0], b[1] - o[1])
'''
return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])
def graham_scan(points):
if len(points) <= 2:
return
# 1、选定原点:y最小,y相同取x最小。
pivot = min(points, key = lambda p:(p[1], p[0]))
# 2、按照极角逆时针排序,极角相同按距离排序。
others = [p for p in points if p != pivot]
others.sort(key = lambda p:(math.atan2(p[1] - pivot[1], p[0] - pivot[0]), (p[0] - pivot) ** 2 + (p[1] - pivot[1]) ** 2))
# 3、栈过滤,叉积如果<=0,则删除顶点。
stack = []
for p in [pivot] + others:
while len(stack >= 2) and cross_product(stack[-2], stack[-1], p) <= 0:
stack.pop()
stack.append(p)
return stack
4、凸包上的点对应无界泰森多边形
如果某个点在点集的最外层,也就是在凸包上,那么它的势力范围会向外无限延伸。
所以:
- 内部点的泰森多边形通常是封闭有限多边形;
- 凸包点的泰森多边形是无界的;
- 在 GIS 中通常需要用研究区边界或矩形框把它裁剪成有限区域。
三、泰森多边形和 Delaunay 三角网的关系
1、Delaunay 三角网是泰森多边形的对偶图
如果两个点在 Delaunay 三角网中有边相连,那么它们的泰森多边形共享一条边。
2、Delaunay 三角形的外心就是泰森顶点
四、Delaunay 三角网的核心准则
1、Delaunay 三角网的关键是 空外接圆准则。
对于三角形 ABC,如果它的外接圆内部不包含其他点,那么这个三角形满足 Delaunay 条件。
2、最大最小角性质
它倾向于避免瘦长三角形,能最大化三角网中的最小角。
五、如何构造泰森多边形?
1、先生成 Delaunay 三角网(分割合并算法、三角网生成算法、Bowyer-Watson逐点插入法)
只介绍逐点插入法:
- (1)构建超级三角形
- (2)将超级三角形的顶点加入点集内,并加入临时三角形数组内
- (3)针对每个原始生成点进行讨论
- (4)拿到该点A,遍历所有临时三角形,如果该点在临时三角形内,则这个临时三角形是坏三角形
- (5)从临时三角形数组中移除坏三角形
- (6)如果坏三角形的边不重合,那么加入数组
- (7)构造新的临时三角形,并加入临时三角形数组
- (8)最后排除超级三角形的顶点构成的三角形
points中每个元素是一个点坐标的列表形式,triangles中每个元素存储三角形顶点在生成点列表中的索引。
def delaunay_(self,points):
# 1、构建超级三角形
n = len(points)
min_x = min([p[0] for p in points])
max_x = max([p[0] for p in points])
min_y = min([p[1] for p in points])
max_y = max([p[1] for p in points])
dx = max_x - min_x
dy = max_y - min_y
t = max(dx, dy)*10
mid_x = (min_x + max_x) / 2
mid_y = (min_y + max_y) / 2
v_p1 = [mid_x - t, mid_y - t] # 左下
v_p = [mid_x, mid_y + t] # 顶点
v_p2 = [mid_x + t, mid_y - t] # 右下
# 2、将超级三角形的顶点加入点集内,并加入临时三角形数组内
all_points = points + [v_p1, v_p, v_p2]
triangles = [[n, n+1, n+2]]
针对每个原始生成点进行讨论
# 3、针对每个原始生成点进行讨论
for i in range(n):
# 一、根据该点删除判断临时三角形数组中的坏三角形并从中删除
# 4、拿到该点,遍历所有临时三角形,如果该点在临时三角形内,则这个临时三角形是坏三角形
point = points[i]
bad_triangles = []
for triangle in triangles:
if self.point_in_circumcircle_extended(point, triangle, all_points):
bad_triangles.append(triangle)
# 5、从临时三角形数组中移除坏三角形
for triangle in bad_triangles:
triangles.remove(triangle)
# 二、准备基于该点构造新的三角形
# 6、如果坏三角形的边不重合,那么加入数组
polygon_edges = []
for triangle in bad_triangles:
for j in range(3):
edge = [triangle[j], triangle[(j+1) % 3]]
edge.sort()
flag = False
for o_triangle in bad_triangles:
if o_triangle == triangle:
continue
for k in range(3):
o_edge = [o_triangle[k], o_triangle[(k+1) % 3]]
o_edge.sort()
if o_edge == edge:
flag = True
break
if flag:
break
if not flag:
polygon_edges.append(edge)
# 7、构造新的临时三角形,并加入临时三角形数组
for edge in polygon_edges:
triangles.append([i, edge[0], edge[1]])
final_triangles = []
# 8、最后排除超级三角形的顶点构成的三角形
for triangle in triangles:
if all(t<n for t in triangle):
final_triangles.append(triangle)
return final_triangles
附录point_in_circumcircle_extended函数:
def point_in_circumcircle_extended(self, point, triangle, all_points):
"""判断点是否在三角形的外接圆内"""
px, py = point
# 获取三角形三个顶点
p1 = all_points[triangle[0]]
p2 = all_points[triangle[1]]
p3 = all_points[triangle[2]]
x1, y1 = p1
x2, y2 = p2
x3, y3 = p3
# 计算外接圆圆心和半径
d = 2 * ((x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1))
if abs(d) < 1e-10: # 三点共线
return False
# 克莱姆法则求解外接圆圆心
ux = ((x1**2 + y1**2) * (y2 - y3) + (x2**2 + y2**2) * (y3 - y1) + (x3**2 + y3**2) * (y1 - y2)) / d
uy = ((x1**2 + y1**2) * (x3 - x2) + (x2**2 + y2**2) * (x1 - x3) + (x3**2 + y3**2) * (x2 - x1)) / d
# 计算半径的平方
radius_sq = (x1 - ux)**2 + (y1 - uy)**2
# 计算点到圆心的距离的平方
dist_sq = (px - ux)**2 + (py - uy)**2
# 如果距离小于半径,则在圆内
return dist_sq < radius_sq
此外,分割合并算法的基本步骤为:
- 划分点集
- 逐个点集进行三角剖分(一般是划分到3个点直接出三角形)
- 合并三角剖分结果
合并三角剖分结果思路:

三角网生成算法:核心分为收缩生长与扩张生长两类,前者先构建点集的凸壳,再向内逐步生成三角网。后者则从初始三角形出发,利用空外接圆或张角最大准则向外层层扩展,通过局部优化保证三角网的“最圆润”特性,最终形成覆盖整个数据区域、满足Delaunay特性的三角网,广泛应用于地形建模、空间分析等场景。