泰森多边形
泰森多边形

泰森多边形

本节包含:基础知识、格林厄姆扫描法生成凸包(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 中通常需要用研究区边界或矩形框把它裁剪成有限区域。
1、Delaunay 三角网是泰森多边形的对偶图

如果两个点在 Delaunay 三角网中有边相连,那么它们的泰森多边形共享一条边。

2、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特性的三角网,广泛应用于地形建模、空间分析等场景。

2、求每个三角形的外心
3、把相邻三角形的外心连起来
4、得到泰森边
5、对凸包边向外延伸并裁剪
6、形成泰森多边形

发表回复

您的邮箱地址不会被公开。 必填项已用 * 标注