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

创作时间:

作者:

@小白创作中心

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

引用

CSDN

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

Delaunay三角剖分是计算几何中的一个重要概念，在许多领域，如计算机图形学、有限元分析、地理信息系统等都有广泛应用。本文将详细介绍Delaunay三角剖分的基本概念、构建算法以及在GIS中的应用案例，帮助读者全面了解这一重要技术。

基本概念

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

构建算法

Bowyer-Watson算法

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

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

分治法

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

步骤1：排序
步骤2：划分
- 如果点集 P 中的点数较少（例如小于等于3个），则直接进行三角剖分（当点数为3时形成一个三角形，点数为2时为一条线段，点数为1时无三角形）。
- 否则，将点集 P 分成左右两个子集 PL 和 PR。
步骤3：递归求解
- 分别对 PL 和 PR 进行Delaunay三角剖分，得到两个子三角剖分 TL 和 TR。
步骤4：合并
- 找到连接 TL 和 TR 的“桥边”，并从下往上逐步合并两个子三角剖分，同时保证合并过程中所有三角形都满足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 函数将原始点绘制为圆形标记。最后设置图表的标题和坐标轴标签，并显示图表。