基于scipy.spatial计算Voronoi图与站点凸包所有交点的高效方案

scipy和numpy目前没有现成的一步式API来直接计算Voronoi图与站点集合凸包的所有交点，但也不需要逐条暴力校验所有边，借助scipy.spatial的内置属性可以大幅降低计算量，高效完成需求。

实现思路 #

• 第一步：分别用scipy.spatial.Voronoi生成站点的Voronoi图，用scipy.spatial.ConvexHull生成站点的凸包
• 第二步：过滤不需要校验的Voronoi边：两个顶点都在凸包内部的有限边、半无限边（ridge_vertices中包含-1的边）的有限端点在凸包外的，都不会和凸包边界相交，可以直接跳过
• 第三步：对剩余的Voronoi边，仅校验和凸包边界的相交情况，收集交点后去重即可（多条边可能交于凸包同一个顶点，会产生重复结果）

代码示例 #

import numpy as np
from scipy.spatial import Voronoi, ConvexHull
from shapely.geometry import LineString, Polygon

# 生成测试站点
points = np.random.rand(20, 2)
# 生成Voronoi图
vor = Voronoi(points)
# 生成站点凸包，转成Polygon方便相交计算
hull = ConvexHull(points)
hull_poly = Polygon(points[hull.vertices])

intersections = []
for ridge_idx, ridge in enumerate(vor.ridge_vertices):
    # 处理有限Voronoi边
    if -1 not in ridge:
        line = LineString(vor.vertices[ridge])
        inter = line.intersection(hull_poly.boundary)
        if not inter.is_empty:
            intersections.extend(np.array(inter.coords))
    # 处理半无限Voronoi边
    else:
        # 提取半无限边的有限端点
        finite_vert_idx = ridge[0] if ridge[1] == -1 else ridge[1]
        finite_vert = vor.vertices[finite_vert_idx]
        # 计算半无限边的方向向量
        p1, p2 = vor.points[vor.ridge_points[ridge_idx]]
        direction = np.array([-(p2[1] - p1[1]), p2[0] - p1[0]])
        # 延伸足够长的线段模拟无限延伸的边
        extended_line = LineString([finite_vert, finite_vert + direction * 1e6])
        inter = extended_line.intersection(hull_poly.boundary)
        if not inter.is_empty:
            intersections.extend(np.array(inter.coords))

# 去重得到最终交点列表
intersections = np.unique(np.array(intersections), axis=0)

注意事项 #

• 示例中用shapely做相交计算是为了代码简洁，你也可以替换为自己实现的线段相交算法，不需要额外引入依赖
• 凸包顶点如果刚好落在Voronoi边上的情况会被自动去重，不需要额外处理
• 对于站点数在1e4以内的场景，这个方法的运行效率远高于全量边暴力校验的方案

内容的提问来源于stack exchange，提问作者jeongbyulji