Voronoi 细胞的体积(python) [英] Volume of Voronoi cell (python)

问题描述

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

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:

• 对于 Voronoi 图的每个区域
• 对该区域的顶点执行 Delaunay 三角剖分
  • 这将返回一组填充该区域的不规则四面体
  • 将这些体积相加以获得该区域的体积.

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

  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.
  

  这篇关于Voronoi 细胞的体积(python)的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！