Delaunay三角剖分：概念、算法与GIS应用

创作时间:

作者:

@小白创作中心

Delaunay三角剖分：概念、算法与GIS应用

引用

CSDN

https://blog.csdn.net/qq_42023999/article/details/145453247

Delaunay三角剖分是计算几何中的一个重要概念，在许多领域，如计算机图形学、有限元分析、地理信息系统等都有广泛应用。本文将详细介绍Delaunay三角剖分的基本概念、构建算法，并通过Python代码实现，最后结合GIS应用案例展示其实际应用场景。

基本概念

三角剖分：给定平面上的一组离散点集P = { p 1 , p 2 , ⋯ , p n } P = {p_1, p_2, \cdots, p_n}P={p1 ,p2 ,⋯,pn }，三角剖分是将这些点连接成互不重叠的三角形的过程，使得这些三角形的并集覆盖点集P PPP所构成的凸包。
Delaunay三角剖分：它是一种特殊的三角剖分，满足空圆特性。即在Delaunay三角剖分中，对于任意一个三角形，其外接圆内部不包含点集P PPP中的其他任何点。

构建算法

Bowyer-Watson算法

这是一种常用的增量式算法，其基本思想是逐步向已有的三角剖分中添加点，每次添加一个点后，对受影响的三角形进行调整以满足Delaunay条件。以下是具体步骤：

步骤1：初始化
构建一个超级三角形，该三角形要足够大，使得所有待剖分的点都在其内部。将这个超级三角形加入到初始的三角剖分集合中。
步骤2：逐点插入
从点集P PPP中依次取出一个点p pp。
找出所有外接圆包含点p pp的三角形，这些三角形构成一个“影响区域”。将这些三角形从当前的三角剖分中移除。
连接点p pp与“影响区域”的边界顶点，形成新的三角形。
检查新形成的三角形是否满足Delaunay条件（空圆特性），若不满足则进行局部调整（通过交换边的方式），直到所有三角形都满足Delaunay条件。
步骤3：移除超级三角形
当所有点都插入完毕后，移除与超级三角形顶点相连的所有三角形，剩下的就是点集P PPP的Delaunay三角剖分。

分治法

分治法的基本思想是将点集分成两个子集，分别对两个子集进行Delaunay三角剖分，然后将两个子三角剖分合并成一个完整的Delaunay三角剖分。具体步骤如下：

步骤1：排序
将点集P PPP按照x xx坐标进行排序。
步骤2：划分
如果点集P PPP中的点数较少（例如小于等于3个），则直接进行三角剖分（当点数为3时形成一个三角形，点数为2时为一条线段，点数为1时无三角形）。
否则，将点集P PPP分成左右两个子集P L P_LPL 和P R P_RPR 。
步骤3：递归求解
分别对P L P_LPL 和P R P_RPR 进行Delaunay三角剖分，得到两个子三角剖分T L T_LTL 和T R T_RTR 。
步骤4：合并
找到连接T L T_LTL 和T R T_RTR 的“桥边”，并从下往上逐步合并两个子三角剖分，同时保证合并过程中所有三角形都满足Delaunay条件。

Python代码实现（Bowyer-Watson算法）

import numpy as np

def circumcircle(p1, p2, p3):
    """
    计算三角形的外接圆
    """
    ax, ay = p1
    bx, by = p2
    cx, cy = p3
    d = 2 * (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by))
    if d == 0:
        return None
    ux = ((ax**2 + ay**2) * (by - cy) + (bx**2 + by**2) * (cy - ay) + (cx**2 + cy**2) * (ay - by)) / d
    uy = ((ax**2 + ay**2) * (cx - bx) + (bx**2 + by**2) * (ax - cx) + (cx**2 + cy**2) * (bx - ax)) / d
    r = np.sqrt((ax - ux)**2 + (ay - uy)**2)
    return (ux, uy), r

def is_point_inside_circle(point, center, radius):
    """
    判断点是否在圆内
    """
    px, py = point
    cx, cy = center
    return (px - cx)**2 + (py - cy)**2 < radius**2

def bowyer_watson(points):
    """
    Bowyer-Watson算法实现Delaunay三角剖分
    """
    # 构建超级三角形
    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
    delta_max = max(dx, dy)
    mid_x = (min_x + max_x) / 2
    mid_y = (min_y + max_y) / 2
    p1 = (mid_x - 20 * delta_max, mid_y - delta_max)
    p2 = (mid_x, mid_y + 20 * delta_max)
    p3 = (mid_x + 20 * delta_max, mid_y - delta_max)
    triangulation = [(p1, p2, p3)]
    for point in points:
        bad_triangles = []
        for triangle in triangulation:
            center, radius = circumcircle(*triangle)
            if center is not None and is_point_inside_circle(point, center, radius):
                bad_triangles.append(triangle)
        polygon = []
        for triangle in bad_triangles:
            for edge in [(triangle[0], triangle[1]), (triangle[1], triangle[2]), (triangle[2], triangle[0])]:
                is_shared = False
                for other_triangle in bad_triangles:
                    if triangle != other_triangle and edge in [(other_triangle[0], other_triangle[1]), (other_triangle[1], other_triangle[2]), (other_triangle[2], other_triangle[0])]:
                        is_shared = True
                        break
                if not is_shared:
                    polygon.append(edge)
        for triangle in bad_triangles:
            triangulation.remove(triangle)
        for edge in polygon:
            new_triangle = (edge[0], edge[1], point)
            triangulation.append(new_triangle)
    # 移除与超级三角形顶点相连的三角形
    final_triangulation = []
    for triangle in triangulation:
        if p1 not in triangle and p2 not in triangle and p3 not in triangle:
            final_triangulation.append(triangle)
    return final_triangulation

# 示例使用
points = [(0, 0), (1, 0), (0, 1), (1, 1)]
triangulation = bowyer_watson(points)
for triangle in triangulation:
    print(triangle)

应用场景

计算机图形学：用于生成三维模型的表面网格，进行曲面重建、纹理映射等操作。
有限元分析：将复杂的几何区域离散成三角形单元，便于进行力学、热学等物理问题的数值求解。
地理信息系统（GIS）：对地形数据进行三角剖分，用于地形建模、坡度分析、视线分析等。

GIS应用案例：地形数据处理

问题描述

假设我们有一组地形上的离散采样点，每个点包含其经纬度坐标以及对应的海拔高度信息。我们希望通过Delaunay三角剖分将这些离散点构建成三角网格，以便进行后续的地形建模和分析，比如计算坡度、进行视线分析等。

示例数据

为了简化问题，我们创建一个包含20个随机点的二维坐标数据，代表地形上的采样点。在实际的GIS应用中，这些点会带有经纬度和海拔信息。

import numpy as np

# 生成20个随机点作为示例地形采样点
np.random.seed(42)
points = np.random.rand(20, 2) * 100  # 坐标范围在0到100之间

代码实现

我们使用scipy库中的Delaunay函数来进行Delaunay三角剖分。该函数可以方便地对给定的点集进行三角剖分，并返回三角网格的相关信息。

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Delaunay

# 生成20个随机点作为示例地形采样点
np.random.seed(42)
points = np.random.rand(20, 2) * 100  # 坐标范围在0到100之间

# 进行Delaunay三角剖分
tri = Delaunay(points)

# 可视化三角网格
plt.triplot(points[:, 0], points[:, 1], tri.simplices)
plt.plot(points[:, 0], points[:, 1], 'o')
plt.title('Delaunay Triangulation of Terrain Points')
plt.xlabel('X Coordinate')
plt.ylabel('Y Coordinate')
plt.show()

代码解释

数据生成：使用numpy的random.rand函数生成20个二维随机点，每个点的坐标范围在0到100之间。
Delaunay三角剖分：调用scipy.spatial中的Delaunay函数对生成的点集进行三角剖分，返回一个Delaunay对象tri。该对象包含了三角网格的相关信息，如每个三角形的顶点索引（存储在tri.simplices中）。
可视化：使用matplotlib库的triplot函数绘制三角网格，同时使用plot函数将原始点绘制为圆形标记。最后设置图表的标题和坐标轴标签，并显示图表。