从Voronoi图到Delaunay三角剖分用Python手把手实现Bowyer-Watson算法计算几何领域中Voronoi图和Delaunay三角剖分是两种紧密相关的空间分割方法。它们不仅在理论上具有优美的数学性质更在计算机图形学、地理信息系统、有限元分析等众多领域有着广泛应用。本文将带您深入理解这两种结构的本质关联并重点讲解如何用Python实现经典的Bowyer-Watson算法来构建Delaunay三角剖分。1. Voronoi图自然界的空间分割想象一下在荒野中分布的几处水源每处水源自然会影响其周边区域——这就是Voronoi图在自然界中的直观体现。给定平面上一组离散点称为基点Voronoi图将平面划分为若干区域每个区域包含距离对应基点最近的所有点。Voronoi图的关键特性每个Voronoi单元都是凸多边形Voronoi顶点是三个基点外接圆的圆心与Delaunay三角剖分存在对偶关系import numpy as np import matplotlib.pyplot as plt def voronoi_plot(points): from scipy.spatial import Voronoi vor Voronoi(points) fig plt.figure() ax fig.add_subplot(111) voronoi_plot_2d(vor, axax) ax.plot(points[:,0], points[:,1], ko) plt.show() # 示例使用 points np.random.rand(15, 2) voronoi_plot(points)2. Delaunay三角剖分的数学之美Delaunay三角剖分是Voronoi图的对偶图具有以下重要性质空圆特性任意三角形的外接圆内不包含其他基点最大化最小角所有可能的三角剖分中Delaunay剖分的最小内角最大这些特性使得Delaunay三角剖分成为许多应用的首选有限元分析中的网格生成三维重建中的表面建模路径规划中的空间划分3. Bowyer-Watson算法详解Bowyer-Watson算法是一种增量式Delaunay三角剖分算法其核心思想是逐点插入并维护Delaunay性质。让我们深入剖析其实现步骤3.1 算法流程初始化超级三角形创建一个包含所有输入点的巨大三角形逐点插入定位影响三角形外接圆包含新点的三角形构建空穴多边形连接新点与空穴边界形成新三角形清理阶段移除与超级三角形相关的所有三角形3.2 关键实现细节超级三角形构造需要确保包含所有输入点。一个可靠的方法是计算点集的包围盒然后向外扩展一定比例def get_super_triangle(points): xs [p.x for p in points] ys [p.y for p in points] xmin, xmax min(xs), max(xs) ymin, ymax min(ys), max(ys) dx, dy xmax - xmin, ymax - ymin dmax max(dx, dy) * 2 # 扩展系数确保包含所有点 return Triangle( Point(xmin - dmax, ymin - dmax), Point(xmin dx/2, ymax dmax), Point(xmax dmax, ymin - dmax) )空圆检测是算法的核心操作。我们可以通过计算行列式来高效判断点是否在三角形外接圆内def point_in_circumcircle(triangle, point): 判断点是否在三角形外接圆内 a, b, c triangle.vertices d point # 转换为齐次坐标行列式 def det(a, b, c): return (a.x - d.x) * (b.y - d.y) * (c.x**2 c.y**2 - d.x**2 - d.y**2) \ (a.y - d.y) * (b.x**2 b.y**2 - d.x**2 - d.y**2) * (c.x - d.x) \ (a.x**2 a.y**2 - d.x**2 - d.y**2) * (b.x - d.x) * (c.y - d.y) - \ (a.x**2 a.y**2 - d.x**2 - d.y**2) * (b.y - d.y) * (c.x - d.x) - \ (a.y - d.y) * (b.x - d.x) * (c.x**2 c.y**2 - d.x**2 - d.y**2) - \ (a.x - d.x) * (b.x**2 b.y**2 - d.x**2 - d.y**2) * (c.y - d.y) return det(a, b, c) 04. 完整Python实现下面给出完整的Bowyer-Watson算法实现包含详细的类型注解和注释from typing import List, Set, Tuple import math class Point: def __init__(self, x: float, y: float): self.x x self.y y def __eq__(self, other): return self.x other.x and self.y other.y def __hash__(self): return hash((self.x, self.y)) class Edge: def __init__(self, p1: Point, p2: Point): self.p1 p1 self.p2 p2 def __eq__(self, other): return (self.p1 other.p1 and self.p2 other.p2) or \ (self.p1 other.p2 and self.p2 other.p1) def __hash__(self): return hash((self.p1, self.p2)) hash((self.p2, self.p1)) class Triangle: def __init__(self, a: Point, b: Point, c: Point): self.vertices (a, b, c) self.edges (Edge(a, b), Edge(b, c), Edge(c, a)) self.circumcircle self.calculate_circumcircle() def calculate_circumcircle(self) - Tuple[Point, float]: 计算三角形的外接圆(圆心和半径) a, b, c self.vertices d 2 * (a.x*(b.y - c.y) b.x*(c.y - a.y) c.x*(a.y - b.y)) if d 0: # 共线情况 return (Point(0,0), float(inf)) ux ((a.x**2 a.y**2)*(b.y - c.y) (b.x**2 b.y**2)*(c.y - a.y) (c.x**2 c.y**2)*(a.y - b.y)) / d uy ((a.x**2 a.y**2)*(c.x - b.x) (b.x**2 b.y**2)*(a.x - c.x) (c.x**2 c.y**2)*(b.x - a.x)) / d center Point(ux, uy) radius math.sqrt((a.x - ux)**2 (a.y - uy)**2) return (center, radius) def contains_point(self, point: Point) - bool: 判断点是否在三角形外接圆内 center, radius self.circumcircle if radius float(inf): return False return (point.x - center.x)**2 (point.y - center.y)**2 radius**2 1e-8 class DelaunayTriangulation: def __init__(self, points: List[Point]): self.points points self.triangles [] self.triangulate() def triangulate(self): if len(self.points) 3: return # 步骤1创建超级三角形 super_triangle self.create_super_triangle() self.triangles.append(super_triangle) # 步骤2逐点插入 for point in self.points: self.add_point(point) # 步骤3移除超级三角形相关的三角形 self.remove_super_triangle(super_triangle) def create_super_triangle(self) - Triangle: 创建包含所有点的超级三角形 min_x min(p.x for p in self.points) max_x max(p.x for p in self.points) min_y min(p.y for p in self.points) max_y max(p.y for p in self.points) dx max_x - min_x dy max_y - min_y delta max(dx, dy) * 2 p1 Point(min_x - delta, min_y - delta) p2 Point(min_x delta, min_y delta) p3 Point(max_x delta, min_y - delta) return Triangle(p1, p2, p3) def add_point(self, point: Point): 向三角剖分中添加一个新点 bad_triangles [] # 查找所有外接圆包含该点的三角形 for triangle in self.triangles: if triangle.contains_point(point): bad_triangles.append(triangle) polygon [] # 查找空穴边界 for triangle in bad_triangles: for edge in triangle.edges: # 统计共享该边的坏三角形数量 shared sum(1 for t in bad_triangles if edge in t.edges) if shared 1: # 只被一个坏三角形包含的边就是边界边 polygon.append(edge) # 移除坏三角形 for triangle in bad_triangles: self.triangles.remove(triangle) # 用新点连接空穴边界形成新三角形 for edge in polygon: new_triangle Triangle(edge.p1, edge.p2, point) self.triangles.append(new_triangle) def remove_super_triangle(self, super_triangle: Triangle): 移除与超级三角形顶点相关的所有三角形 super_points super_triangle.vertices self.triangles [ t for t in self.triangles if not any(p in super_points for p in t.vertices) ] def get_edges(self) - Set[Edge]: 获取所有唯一的边 edges set() for triangle in self.triangles: for edge in triangle.edges: edges.add(edge) return edges5. 可视化与验证为了验证我们的实现我们可以将结果与scipy的Delaunay实现进行对比def plot_triangulation(points, triangulation): plt.figure(figsize(12, 6)) # 绘制我们的实现结果 plt.subplot(121) for edge in triangulation.get_edges(): plt.plot([edge.p1.x, edge.p2.x], [edge.p1.y, edge.p2.y], b-) plt.plot([p.x for p in points], [p.y for p in points], ro) plt.title(Our Implementation) # 绘制scipy的结果 plt.subplot(122) from scipy.spatial import Delaunay tri Delaunay(np.array([[p.x, p.y] for p in points])) plt.triplot(np.array([p.x for p in points]), np.array([p.y for p in points]), tri.simplices, b-) plt.plot(np.array([p.x for p in points]), np.array([p.y for p in points]), ro) plt.title(SciPy Implementation) plt.show() # 测试用例 test_points [Point(x, y) for x, y in np.random.rand(20, 2)] dt DelaunayTriangulation(test_points) plot_triangulation(test_points, dt)6. 性能分析与优化Bowyer-Watson算法的时间复杂度在最坏情况下为O(n²)但通过一些优化可以显著提高实际性能优化策略点定位加速使用空间索引结构如KD树来快速定位影响三角形增量排序按x坐标排序输入点减少搜索范围内存管理重用已删除的三角形数据结构from scipy.spatial import KDTree class OptimizedDelaunay(DelaunayTriangulation): def __init__(self, points: List[Point]): self.kdtree KDTree([[p.x, p.y] for p in points]) super().__init__(points) def find_bad_triangles(self, point: Point) - List[Triangle]: 使用KD树加速坏三角形查找 # 先找到最近的几个三角形进行初步筛选 _, indices self.kdtree.query([point.x, point.y], k5) nearby_points [self.points[i] for i in indices] # 查找包含这些点的三角形 candidate_triangles set() for triangle in self.triangles: if any(p in triangle.vertices for p in nearby_points): candidate_triangles.add(triangle) # 精确检查 return [t for t in candidate_triangles if t.contains_point(point)]在实际项目中当处理大规模点集时可以考虑分治算法或并行计算来进一步提高性能。对于实时应用还可以研究增量式更新算法在已有剖分的基础上高效处理新增或删除的点。理解Bowyer-Watson算法不仅有助于掌握计算几何的基础知识更能为各种空间分析问题提供强大的工具。本文实现的完整代码库包含了更多实用功能和测试用例读者可以在此基础上继续扩展和优化。