Scipy计算Voronoi胞体体积结果异常问题求助

需要计算半径100、高度100的圆柱内封闭Voronoi胞体的体积，圆柱内是堆积率50%的非重叠球体中心坐标。使用scipy.spatial.Voronoi生成胞体，通过regions中是否包含-1判定开放胞体（体积设为0），封闭胞体用ConvexHull计算体积，但运行结果显示封闭胞体总体积远超圆柱体积，部分胞体体积异常庞大，需排查原因并修正。

• 未考虑圆柱边界约束：scipy.spatial.Voronoi生成的是无边界约束的全空间Voronoi图，对于靠近圆柱内壁/上下底面的点，其实际胞体的边界应被圆柱壁面截断，但当前代码直接用ConvexHull计算的是无限空间中该胞体的完整凸包体积，完全忽略了圆柱的边界限制，导致这些胞体的体积被过度计算。
• 开放胞体判定逻辑局限：仅通过regions包含-1判断开放胞体并不准确——部分靠近圆柱边界的点在全空间中属于封闭胞体，但实际在圆柱内的部分是有限的，代码却将其当作完整封闭胞体计算体积，这是体积异常的核心原因。

核心思路是将Voronoi胞体与圆柱空间做交集，只计算胞体在圆柱内部的体积，以下是两种可行实现方式：

方法1：用布尔运算裁剪胞体（推荐） #

借助trimesh库处理Voronoi胞体与圆柱的布尔交集，精准计算胞体在圆柱内的体积：

from scipy.spatial import Voronoi, ConvexHull
import numpy as np
import trimesh

# 几何参数
radius_spheres = 6
radius_cylinder = 100
height_cylinder = 100

# 读取点数据
points = np.genfromtxt('test_pos.txt', delimiter=',')
number_spheres = points.shape[0]

# 计算球体总体积和圆柱体积
volume_spheres_total = 4/3*np.pi*radius_spheres**3 * number_spheres
volume_cylinder = np.pi*(radius_cylinder**2)*height_cylinder

# 生成Voronoi图
v = Voronoi(points)
volume_Voronoi = np.zeros(v.npoints)

# 创建圆柱网格（轴线为z轴，z范围0~100）
cylinder = trimesh.creation.cylinder(
    radius=radius_cylinder, 
    height=height_cylinder, 
    transform=np.array([
        [1,0,0,0],
        [0,1,0,0],
        [0,0,1,height_cylinder/2],
        [0,0,0,1]
    ])
)

for idx, i_region in enumerate(v.point_region):
    region = v.regions[i_region]
    if -1 in region:
        volume_Voronoi[idx] = 0
        continue
    # 获取胞体顶点并生成凸包网格
    cell_vertices = v.vertices[region]
    try:
        hull = ConvexHull(cell_vertices)
        cell_mesh = trimesh.Trimesh(vertices=cell_vertices, faces=hull.simplices)
        # 计算胞体与圆柱的交集
        intersection = trimesh.boolean.intersection([cell_mesh, cylinder])
        if intersection is not None and len(intersection.vertices) >= 4:
            volume_Voronoi[idx] = intersection.volume
        else:
            volume_Voronoi[idx] = 0
    except:
        volume_Voronoi[idx] = 0

volume_Voronoi_total = sum(volume_Voronoi)

# 计算堆积率
packing_fraction_theoretical = volume_spheres_total/volume_cylinder
packing_fraction_Voronoi = volume_spheres_total/volume_Voronoi_total

print(f"理论堆积率: {packing_fraction_theoretical:.4f}")
print(f"Voronoi堆积率: {packing_fraction_Voronoi:.4f}")
print(f"封闭胞体总体积: {volume_Voronoi_total:.2f}")
print(f"圆柱体积: {volume_cylinder:.2f}")

方法2：手动约束胞体顶点到圆柱内 #

若不想引入额外依赖，可手动过滤/投影胞体顶点到圆柱范围内，再重新计算凸包体积：

1. 对每个胞体顶点，判断是否在圆柱内：满足x² + y² ≤ radius_cylinder²且0 ≤ z ≤ height_cylinder
2. 对超出圆柱的顶点，将其投影到圆柱边界（比如xy平面超出半径的点投影到圆柱侧壁，z超出范围的点投影到上下底面）
3. 用处理后的顶点重新构建凸包并计算体积

该方法实现繁琐，精度略低于布尔运算，但无需额外安装库。

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