Bowyer-Watson算法

数学原理及算法过程

Delaunay 三角剖分是一种特殊的三角剖分方法，它满足以下两个重要性质：

• 最大化最小角性质：Delaunay 三角剖分通过避免细长的三角形来最大化所有三角形的最小角。
• 空外接圆性质：在 Delaunay 三角剖分中，每个三角形的外接圆不包含任何其他点。这意味着，对于三角剖分中的任意三角形，其外接圆内没有其他输入点。

基于这些性质，Delaunay 三角剖分算法的一种实现方式是 Bowyer-Watson 算法，这是一种增量算法。以下是具体的算法步骤：

算法过程

1. 初始化超级三角形：

• 创建一个足够大的超级三角形，包含所有输入点。这个三角形的三个顶点坐标远离实际输入点的范围，使其能够覆盖所有点。

2. 逐点插入：

• 对于每个输入点，找到所有包含该点的外接圆的三角形。这些三角形被称为“坏三角形”。

3. 构建多边形：

• 对于所有坏三角形，它们的每条边，如果只被一个坏三角形共享，则称其为边界边。这些边将形成一个多边形。

4. 删除坏三角形：

• 将所有坏三角形从三角剖分中删除。

5. 重新三角化多边形：

• 用新插入的点和多边形的边构成新的三角形，并将这些三角形加入三角剖分中。

6. 移除超级三角形的影响：

• 在所有点都插入后，移除包含超级三角形顶点的所有三角形，得到最终的 Delaunay 三角剖分。

数学原理

• 外接圆计算：

  • 对于每个三角形，计算其外接圆。外接圆的圆心（外心）和半径可以通过三角形顶点的坐标计算。
  • 设三角形顶点为 $(x_1,y_1),(x_2,y_2),(x_3,y_3)$。外接圆的圆心 $(u,v)$ 计算如下：
    $$d=2\left(x_1(y_2-y_3)+x_2(y_3-y_1)+x_3(y_1-y_2)\right)$$

  $$u=\frac{((x_1^2+y_1^2)(y_2-y_3)+(x_2^2+y_2^2)(y_3-y_1)+(x_3^2+y_3^2)(y_1-y_2))}{d}$$

  $$v=\frac{((x_1^2+y_1^2)(x_3-x_2)+(x_2^2+y_2^2)(x_1-x_3)+(x_3^2+y_3^2)(x_2-x_1))}{d}$$

  $$r=\sqrt{(x_1-u{)}^2+(y_1-v{)}^2}$$

import matplotlib.pyplot as plt
import numpy as npclass Point:def __init__(self, x, y):self.x = xself.y = yclass Triangle:def __init__(self, p1, p2, p3):self.p1 = p1self.p2 = p2self.p3 = p3self.circumcenter, self.circumradius = self.circumcircle()def circumcircle(self):"""Calculate the circumcenter and circumradius of the triangle."""ax, ay = self.p1.x, self.p1.ybx, by = self.p2.x, self.p2.ycx, cy = self.p3.x, self.p3.yd = 2 * (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by))ux = ((ax*ax + ay*ay) * (by - cy) + (bx*bx + by*by) * (cy - ay) + (cx*cx + cy*cy) * (ay - by)) / duy = ((ax*ax + ay*ay) * (cx - bx) + (bx*bx + by*by) * (ax - cx) + (cx*cx + cy*cy) * (bx - ax)) / dcircumcenter = Point(ux, uy)circumradius = np.sqrt((ax - ux)**2 + (ay - uy)**2)return circumcenter, circumradiusdef contains_point(self, p):"""Check if the point p is inside the circumcircle of the triangle."""return np.sqrt((p.x - self.circumcenter.x)**2 + (p.y - self.circumcenter.y)**2) < self.circumradiusdef delaunay_triangulation(points):"""Perform Delaunay triangulation on a set of points."""super_triangle = Triangle(Point(-1e5, -1e5), Point(1e5, -1e5), Point(0, 1e5))triangulation = [super_triangle]for p in points:bad_triangles = []for tri in triangulation:if tri.contains_point(p):bad_triangles.append(tri)polygon = []for tri in bad_triangles:for edge in [(tri.p1, tri.p2), (tri.p2, tri.p3), (tri.p3, tri.p1)]:is_shared = Falsefor other in bad_triangles:if other != tri and (edge in [(other.p1, other.p2), (other.p2, other.p3), (other.p3, other.p1)] or edge[::-1] in [(other.p1, other.p2), (other.p2, other.p3), (other.p3, other.p1)]):is_shared = Truebreakif not is_shared:polygon.append(edge)for tri in bad_triangles:triangulation.remove(tri)for edge in polygon:triangulation.append(Triangle(edge[0], edge[1], p))triangulation = [tri for tri in triangulation if not (super_triangle.p1 in [tri.p1, tri.p2, tri.p3] or super_triangle.p2 in [tri.p1, tri.p2, tri.p3] or super_triangle.p3 in [tri.p1, tri.p2, tri.p3])]return triangulationdef plot_triangulation(triangles, points):for tri in triangles:plt.plot([tri.p1.x, tri.p2.x], [tri.p1.y, tri.p2.y], 'b-')plt.plot([tri.p2.x, tri.p3.x], [tri.p2.y, tri.p3.y], 'b-')plt.plot([tri.p3.x, tri.p1.x], [tri.p3.y, tri.p1.y], 'b-')for p in points:plt.plot(p.x, p.y, 'ro')plt.show()
# Generate random points in the unit square
rectangle_corners = [Point(0, 0), Point(1, 0), Point(1, 1), Point(0, 1)]
random_points = [Point(np.random.rand(), np.random.rand()) for _ in range(20)]
points = rectangle_corners + random_pointstriangles = delaunay_triangulation(points)
plot_triangulation(triangles, points)