Volume of Voronoi cell (python)

Question

I'm using Scipy 0.13.0 in Python 2.7 to calculate a set of Voronoi cells in 3d. I need to get the volume of each cell for (de)weighting output of a proprietary simulation. Is there any simple way of doing this - surely it's a common problem or a common use of Voronoi cells but I can't find anything. The following code runs, and dumps everything that the scipy.spatial.Voronoi manual knows about.

from scipy.spatial import Voronoi
x=[0,1,0,1,0,1,0,1,0,1]
y=[0,0,1,1,2,2,3,3.5,4,4.5]
z=[0,0,0,0,0,1,1,1,1,1]
points=zip(x,y,z)
print points
vor=Voronoi(points)
print vor.regions
print vor.vertices
print vor.ridge_points
print vor.ridge_vertices
print vor.points
print vor.point_region

Hypogeous

Answer 1

I think I've cracked it. My approach below is:

• For each region of the Voronoi diagram
• perform a Delaunay triangulation of the vertices of the region
  • this will return a set of irregular tetrahedra which fill the region
• The volume of a tetrahedron can be calculated easily (wikipedia)
  • sum these volumes to get the volume of the region.

I'm sure there will be both bugs and poor coding - I'll be looking for the former, comments welcome on the latter - especially as I'm quite new to Python. I'm still checking a couple of things - sometimes a vertex index of -1 is given, which according to the scipy manual "indicates vertex outside the Voronoi diagram", but in addition, vertices are generated with coordinates well outside the original data (insert numpy.random.seed(42) and check out the coordinates of the region for point 7, they go to ~(7,-14,6), point 49 is similar. So I need to figure out why sometimes this happens, and sometimes I get index -1.

from scipy.spatial import Voronoi,Delaunay
import numpy as np
import matplotlib.pyplot as plt

def tetravol(a,b,c,d):
 '''Calculates the volume of a tetrahedron, given vertices a,b,c and d (triplets)'''
 tetravol=abs(np.dot((a-d),np.cross((b-d),(c-d))))/6
 return tetravol

def vol(vor,p):
 '''Calculate volume of 3d Voronoi cell based on point p. Voronoi diagram is passed in v.'''
 dpoints=[]
 vol=0
 for v in vor.regions[vor.point_region[p]]:
  dpoints.append(list(vor.vertices[v]))
 tri=Delaunay(np.array(dpoints))
 for simplex in tri.simplices:
  vol+=tetravol(np.array(dpoints[simplex[0]]),np.array(dpoints[simplex[1]]),np.array(dpoints[simplex[2]]),np.array(dpoints[simplex[3]]))
 return vol

x= [np.random.random() for i in xrange(50)]
y= [np.random.random() for i in xrange(50)]
z= [np.random.random() for i in xrange(50)]
dpoints=[]
points=zip(x,y,z)
vor=Voronoi(points)
vtot=0


for i,p in enumerate(vor.points):
 out=False
 for v in vor.regions[vor.point_region[i]]:
  if v<=-1: #a point index of -1 is returned if the vertex is outside the Vornoi diagram, in this application these should be ignorable edge-cases
   out=True
  else:
 if not out:
  pvol=vol(vor,i)
  vtot+=pvol
  print "point "+str(i)+" with coordinates "+str(p)+" has volume "+str(pvol)

print "total volume= "+str(vtot)

#oddly, some vertices outside the boundary of the original data are returned, meaning that the total volume can be greater than the volume of the original.

Hypogeous

Having Voronoi vertices lie outside the convex hull of the data is expected behavior. Voronoi vertices are circumcenters of the triangles of the Delaunay triangulation of the same point set, and if you have thin triangles near the edge, their circumcenters can be far outside the bounds of the original point set. – Tasset

Thanks @Tasset I thought that might be the case but as I didn't see it in my 2d test case that I could graph, I wasn't sure. – Hypogeous

The Qhull documentation recommends to compute the convex hull instead of delaunay for each voronoi region, but otherwise the recommendation there agrees with what you do above, so it's probably the best available way. In principle, the Scipy Qhull wrappers could also compute convex hull volumes, but it doesn't seem that Qhull offers ways to directly obtain the voronoi region volumes without additional convex hull computations. – Tasset

@Tasset I finally found something explicitly saying Qhull doesn't make the volume available, so this way appears to be close to the most appropriate. Factors related to the computation are now small compared to the other factors I have to worry about in this analysis, so for me at least it's sorted. – Hypogeous

Answer 2

As was mentioned in comments, you can compute ConvexHull of each Voronoi cell. Since Voronoi cells are convex, you will get the proper volumes.

def voronoi_volumes(points):
    v = Voronoi(points)
    vol = np.zeros(v.npoints)
    for i, reg_num in enumerate(v.point_region):
        indices = v.regions[reg_num]
        if -1 in indices: # some regions can be opened
            vol[i] = np.inf
        else:
            vol[i] = ConvexHull(v.vertices[indices]).volume
    return vol

This method works in any dimensions

Sible

Is it possible to compute the volume of the unbounded Voronoi cells considering the bounded domain (considering bounded axes of the domain)? – Tenno