如何在正方形边界内绘制排除边界顶点的Voronoi胞元

核心修改点 #

• 保留原有添加四个正方形顶点作为辅助采样点的逻辑，仅需要在输出阶段跳过这四个点对应的Voronoi胞元即可
• 裁剪蒙版直接使用目标正方形边界，替代原有点集凸包，裁剪结果更稳定
• 直接调用shapely对象的area属性即可获取每个胞元的精确面积

可运行完整代码 #

import matplotlib.pyplot as plt
import numpy as np
from shapely.geometry import Polygon
from scipy.spatial import Voronoi

def voronoi_finite_polygons_2d(vor, radius=None):
    """
    重构2D Voronoi图中的无限区域为有限区域
    参数
    ----------
    vor : Voronoi
        输入的Voronoi图对象
    radius : float, 可选
        到"无穷远点"的距离
    返回值
    -------
    regions : 元组列表
        修正后的每个Voronoi区域的顶点索引
    vertices : 元组列表
        修正后的Voronoi顶点坐标，在输入顶点基础上末尾追加了"无穷远点"对应的坐标
    """

    if vor.points.shape[1] != 2:
        raise ValueError("仅支持2D输入")

    new_regions = []
    new_vertices = vor.vertices.tolist()

    center = vor.points.mean(axis=0)
    if radius is None:
        radius = vor.points.ptp().max()*2

    # 构建每个点对应的所有脊的映射
    all_ridges = {}
    for (p1, p2), (v1, v2) in zip(vor.ridge_points, vor.ridge_vertices):
        all_ridges.setdefault(p1, []).append((p2, v1, v2))
        all_ridges.setdefault(p2, []).append((p1, v1, v2))

    # 重构无限区域
    for p1, region in enumerate(vor.point_region):
        vertices = vor.regions[region]

        if all(v >= 0 for v in vertices):
            # 有限区域
            new_regions.append(vertices)
            continue

        # 重构非有限区域
        ridges = all_ridges[p1]
        new_region = [v for v in vertices if v >= 0]

        for p2, v1, v2 in ridges:
            if v2 < 0:
                v1, v2 = v2, v1
            if v1 >= 0:
                # 有限脊：已经在区域中
                continue

            # 计算无限脊的缺失端点
            t = vor.points[p2] - vor.points[p1] # 切向量
            t /= np.linalg.norm(t)
            n = np.array([-t[1], t[0]])  # 法向量

            midpoint = vor.points[[p1, p2]].mean(axis=0)
            direction = np.sign(np.dot(midpoint - center, n)) * n
            far_point = vor.vertices[v2] + direction * radius

            new_region.append(len(new_vertices))
            new_vertices.append(far_point.tolist())

        # 按逆时针排序区域顶点
        vs = np.asarray([new_vertices[v] for v in new_region])
        c = vs.mean(axis=0)
        angles = np.arctan2(vs[:,1] - c[1], vs[:,0] - c[0])
        new_region = np.array(new_region)[np.argsort(angles)]

        # 完成重构
        new_regions.append(new_region.tolist())

    return new_regions, np.asarray(new_vertices)

# 生成随机采样点
point_count = 10
a = np.random.random(point_count)
b = np.random.random(point_count)
data = list(zip(a,b))
# 正方形边界顶点，作为辅助采样点
bound = [[0,0], [1,0], [1,1], [0,1]]
data.extend(bound)
points = np.array(data)

# 生成Voronoi图并重构有限区域
vor = Voronoi(points)
regions, vertices = voronoi_finite_polygons_2d(vor)

# 定义正方形裁剪蒙版
mask = Polygon(bound)
cell_areas = []
valid_points = points[:point_count]

# 仅遍历前point_count个采样点对应的胞元，跳过最后4个边界顶点的胞元
for idx in range(point_count):
    region_idx = vor.point_region[idx]
    region = regions[region_idx]
    polygon_vertices = vertices[region]
    # 构建多边形并和正方形边界求交
    voronoi_poly = Polygon(polygon_vertices)
    clipped_poly = voronoi_poly.intersection(mask)
    # 保存面积和绘制
    cell_areas.append(clipped_poly.area)
    poly_coords = np.array(clipped_poly.exterior.coords)
    plt.fill(poly_coords[:,0], poly_coords[:,1], "white", edgecolor = 'black', zorder = 1)

# 输出每个胞元的面积
for i, area in enumerate(cell_areas):
    print(f"第{i+1}个胞元面积：{area:.4f}")

# 绘制采样点
plt.scatter(valid_points[:,0], valid_points[:,1], color = "blue", marker = "o", zorder = 2)
plt.xlim(0,1)
plt.ylim(0,1)
plt.show()

效果说明 #

运行代码后只会输出生成的10个随机点对应的Voronoi胞元，四个角落的辅助顶点胞元会被完全过滤，每个胞元的面积会直接打印输出，所有胞元都会被限制在[0,1]的正方形范围内。

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