Voronoi單元的體積（python）

我在Python 2.7中使用Scipy 0.13.0來計算3D中的一組Voronoi單元。我需要獲得每個單元的體積以用於專有仿真的（去）權重輸出。有沒有簡單的方法來做到這一點 - 這肯定是Voronoi單元的常見問題或常見用法，但我找不到任何東西。以下代碼將運行，並轉儲scipy.spatial.Voronoi manual知道的所有內容。Voronoi單元的體積（python）

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

來源

2013-10-28 Chris H

我認爲我已經破解了。下面我的做法是：

對於維諾圖的每個區域
執行的區域的頂點的Delaunay三角剖分
- 這將返回一組不規則的四面體，其填充區域的
可以計算四面體的體積easily (wikipedia)
- 將這些體積相加得到體積上。

我敢肯定，會有兩種錯誤和編碼差 - 我會找前者，評論對後者的歡迎 - 尤其是因爲我很新的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.

來源

2013-10-29 14:40:42

具有沃羅諾伊頂點位於外側的數據的凸包是預期行爲。沃羅諾伊頂點是相同點集的Delaunay三角剖分的三角形的外心，如果你有邊緣附近的薄三角形，他們的外部中心可能遠在博外原始點集的最後部分。 –

謝謝@pv。我認爲可能是這種情況，但由於我沒有在我可以繪製的2D測試案例中看到它，所以我不確定。 –

Qhull文檔[建議計算凸包]（http://www.qhull.org/html/qh-faq.htm#volume），而不是每個voronoi區域的delaunay，但否則其中的建議與您的做上面的，所以它可能是最好的可用方式。原則上，Scipy Qhull包裝器也可以計算凸包體積，但似乎Qhull似乎沒有提供直接獲得voronoi區域體積的方法，而無需額外的凸包計算。 –

Voronoi單元的體積（python）

回答

相關問題