Voronoi cells volumes are off with Scipy

Question

I would like to calculate the volume of (closed) Voronoi cells. There are points which correspond to non-overlapping sphere centers, and which have a packing fraction within a cylinder of radius 100, height 100, of 50%. These points can be found here.

Voronoi cells can easily be generated with spicy.spatial. When, in a region contains -1, it means that it is open, and therefore I cannot really calculate the cell volumes (can I?), so it is set to 0. Otherwise, with ConvexHull, it is possible to calculate the volume of a polyhedron. This is done with the following script:

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

With this implementation, I would expect the total volume of the Voronoi cells to be roughly equal to the volume of the cylinder in which the spheres are. It should even be a bit underestimated since the open cells where not considered. However, when printing the results with:

# 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]))

I obtained the following:

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

We can see that the total volume is way overestimated by the Voronoi technique. Looking at the distribution, it seems that some cells have huge volumes, and I do not really understand why. Would you have any hints of what might be happening or what to change?

Thank you.

Answer 1

The resulting Voronoi cells are not limited to the cylinder volume and can protrude well beyond its surface.

Think of three points just on the flat surface of the cylinder, arranged in an equilateral triangle. Then, place a fourth point in the center of that triangle, just below the surface. This will make the cell associated with this fourth point to protrude well beyond the surface.

A trivial way to address this would be to check if all vertices of a voronoi cell are inside the cylinder, and ignore them otherwise. Another approach would be to compute intersection of voronoi cell with the cyclinder.