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
I think I've cracked it. My approach below is:
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.