Voronoi细胞(蟒蛇)的体积

Question

我在Python 2.7中使用Scipy 0.13.0来计算一组3d的Voronoi单元格。 我需要获取每个单元格的体积，以便对专有模拟进行（去）加权输出。 有没有简单的方法可以做到这一点 - 肯定这是一个常见问题或 Voronoi 细胞的常见用途，但我找不到任何东西。 下面的代码运行，并转储scipy.spatial.Voronoi手册所知道的所有内容。

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

Answer 1

如注释中所述，您可以计算每个 Voronoi 单元的凸壳。由于Voronoi细胞是凸的，您将获得适当的体积。

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

此方法适用于任何维度

Answer 2

我想我已经破解了它。 我下面的方法是：

• 对于沃罗诺伊图的每个区域
• 对区域的顶点执行 Delaunay 三角测量
  • 这将返回一组不规则的四面体，填充该区域
• 可以轻松计算四面体的体积（维基百科）
  • 将这些卷相加以得到该区域的卷。

我相信会有错误和糟糕的编码 - 我会寻找前者，欢迎评论后者 - 特别是因为我对 Python 很陌生。 我仍在检查一些事情 - 有时会给出 -1 的顶点索引，根据 scipy 手册"指示 Voronoi 图之外的顶点"，但此外，顶点是用远在原始数据之外的坐标生成的（插入numpy.random.seed(42)并查看点 7 的区域坐标，它们转到 ~（7，-14,6），第 49 点类似。 所以我需要弄清楚为什么有时会发生这种情况，有时我会得到索引 -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.