Voronoi 细胞的体积（python）答案

【问题标题】：Volume of Voronoi cell (python)Voronoi 细胞的体积（python）
【发布时间】：2013-11-07 05:41:12
【问题描述】：

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

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

【问题讨论】：

标签： python-2.7 scipy voronoi qhull

【解决方案1】：

我认为我已经破解了。我下面的方法是：

Voronoi 图的每个区域
对该区域的顶点执行 Delaunay 三角剖分
- 这将返回一组填充该区域的不规则四面体
可以计算四面体的体积easily (wikipedia)
- 将这些体积相加得到该区域的体积。

我敢肯定会有错误和糟糕的编码 - 我会寻找前者，欢迎 cmets 在后者 - 特别是因为我对 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.

【讨论】：

Voronoi 顶点位于数据凸包之外是预期行为。 Voronoi 顶点是同一点集的 Delaunay 三角剖分的三角形的外心，如果在边缘附近有细三角形，它们的外心可能会远远超出原始点集的范围。
谢谢@pv。我认为可能是这种情况，但由于我在 2d 测试用例中没有看到它可以绘制图表，所以我不确定。
Qhull 文档 recommends to compute the convex hull 而不是每个 voronoi 区域的 delaunay，但除此之外，那里的建议与您在上面所做的一致，因此它可能是最好的可用方式。原则上，Scipy Qhull 包装器也可以计算凸包体积，但 Qhull 似乎没有提供无需额外凸包计算即可直接获得 voronoi 区域体积的方法。
@pv。我终于找到了一些明确的说法，即 Qhull 不会使该卷可用，因此这种方式似乎接近最合适的方式。与我在此分析中必须担心的其他因素相比，与计算相关的因素现在很小，所以至少对我来说是排序的。

【解决方案2】：

正如 cmets 中提到的，您可以计算每个 Voronoi 单元的 ConvexHull。由于 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

此方法适用于任何维度

【讨论】：