我想计算(封闭的)Voronoi单元的体积。有一些点对应于不重叠的球心,并且在半径为100、高度为100的圆柱体内具有50%的填充分数。这些点可以被找到here。
使用spicy.spatial
可以很容易地生成Voronoi单元格。当区域中包含-1时,意味着它是开放的,因此我无法真正计算单元格的体积(可以吗?),因此将其设置为0。否则,使用ConvexHull
,可以计算多面体的体积。这是通过以下脚本完成的:
from scipy.spatial import Voronoi, ConvexHull
import numpy as np
# These geometrical values were verified
radius_spheres = 6
radius_cylinder = 100
height_cylinder = 100
# Read positions file
points = np.genfromtxt('test_pos.txt', delimiter=',')
number_spheres = points.shape[0]
# Calculate total volume of the spheres
volume_spheres_total = 4/3*np.pi*radius_spheres**3 * number_spheres
# Calculate real cylinder volume
volume_cylinder = np.pi*(radius_cylinder**2)*height_cylinder
# Calculate Voronoi cells volumes
v = Voronoi(points)
volume_Voronoi = np.zeros(v.npoints)
for idx, i_region in enumerate(v.point_region):
# Open cells are removed from the calculation of the total volume
if -1 not in v.regions[i_region]:
volume_Voronoi[idx] = ConvexHull(v.vertices[v.regions[i_region]]).volume
volume_Voronoi_total = sum(volume_Voronoi)
# Analyze distribution of Voronoi volume values
counts, values = np.histogram(volume_Voronoi, bins=np.logspace(np.log10(np.min(volume_Voronoi[volume_Voronoi!=0])), np.log10(np.max(volume_Voronoi)), 10))
# Calculate theorical and Voronoi packing fraction
packing_fraction_theoretical = volume_spheres_total/volume_cylinder
packing_fraction_Voronoi = volume_spheres_total/volume_Voronoi_total
通过这种实现,我希望Voronoi单元的总体积大致等于球体所在的圆柱体的体积。由于未考虑开放单元,因此该值甚至可能被低估。但是,当打印结果时,请使用:
# Summarize
print('Total volume of the spheres: {:.0f}\n'.format(volume_spheres_total))
print('Volume of the cylinder: {:.0f}'.format(volume_cylinder))
print('-> Theoretical packing fraction: {:.0f}%\n'.format(packing_fraction_theoretical*100))
print('Total volume of the (closed) Voronoi cells: {:.0f}'.format(volume_Voronoi_total))
print('Difference between Voronoi volume and cylinder volume: {:.0f}%'.format((volume_Voronoi_total-volume_cylinder)/volume_cylinder))
print('-> Voronoi packing fraction: {:.0f}%\n'.format(packing_fraction_Voronoi*100))
print('Volumes distribution:')
for i in range(len(values)-1):
if counts[i]!=0:
print('\t{} Voronoi cell volumes are between {:.2E} and {:.2E}'.format(counts[i], values[i], values[i+1]))
我获得了以下信息:
Total volume of the spheres: 1570696
Volume of the cylinder: 3141593
-> Theoretical packing fraction: 50%
Total volume of the (closed) Voronoi cells: 5547909930
Difference between Voronoi volume and cylinder volume: 1765%
-> Voronoi packing fraction: 0%
Volumes distribution:
1181 Voronoi cell volumes are between 1.46E+03 and 7.47E+03
136 Voronoi cell volumes are between 7.47E+03 and 3.82E+04
96 Voronoi cell volumes are between 3.82E+04 and 1.95E+05
63 Voronoi cell volumes are between 1.95E+05 and 9.99E+05
36 Voronoi cell volumes are between 9.99E+05 and 5.11E+06
15 Voronoi cell volumes are between 5.11E+06 and 2.61E+07
12 Voronoi cell volumes are between 2.61E+07 and 1.33E+08
2 Voronoi cell volumes are between 1.33E+08 and 6.82E+08
我们可以看到总体积被Voronoi技术高估了很多。从分布来看,似乎有些细胞的体积很大,我真的不明白为什么。你能告诉我可能发生了什么或者需要改变什么吗?
- 谢谢-谢谢
1条答案
按热度按时间rxztt3cl1#
所得到的Voronoi单元不限于圆柱体体积,并且可以突出到其表面之外。
假设圆柱体的平面上有三个点,排列成一个等边三角形,然后在三角形的中心放置第四个点,使与第四个点相关的细胞突出到平面之外。
解决这个问题的一个简单方法是检查Voronoi单元的所有顶点是否都在圆柱体内,否则忽略它们。另一种方法是计算Voronoi单元与圆柱体的交集。