Объем клетки Вороного (питон)

Я использую Scipy 0.13.0 в Python 2.7 для расчета набора ячеек Вороного в 3D. Мне нужно получить объем каждой ячейки для (де)взвешивания вывода проприетарной симуляции. Есть ли простой способ сделать это - конечно, это общая проблема или обычное использование ячеек Вороного, но я ничего не могу найти. Следующий код запускается и выводит все, что scipy .spatial.Вороной руководство об этом знает.

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
person Chris H    schedule 28.10.2013    source источник

Ответы (2)


arrow_upward
11
arrow_downward

Как упоминалось в комментариях, вы можете вычислить ConvexHull каждой ячейки Вороного. Поскольку клетки Вороного выпуклые, вы получите правильные объемы.

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

Этот метод работает в любых измерениях

person ybeltukov    schedule 19.01.2019

arrow_upward
4
arrow_downward

Я думаю, что взломал его. Мой подход ниже:

  • Для каждой области диаграммы Вороного
  • perform a Delaunay triangulation of the vertices of the region
    • this will return a set of irregular tetrahedra which fill the region
  • The volume of a tetrahedron can be calculated easily (wikipedia)
    • sum these volumes to get the volume of the region.

Я уверен, что будут как ошибки, так и плохое кодирование - я буду искать первое, комментарии приветствуются по последнему - тем более, что я новичок в Python. Я все еще проверяю пару вещей - иногда дается индекс вершины -1, который согласно руководству по scipy "указывает на вершину вне диаграммы Вороного", но кроме того, вершины генерируются с координаты далеко за пределами исходных данных (вставьте 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.
person Chris H    schedule 29.10.2013
comment
Наличие вершин Вороного вне выпуклой оболочки данных является ожидаемым поведением. Вершины Вороного являются центрами описанных окружностей треугольников триангуляции Делоне одного и того же набора точек, и если у вас есть тонкие треугольники рядом с краем, их центры описанных окружностей могут находиться далеко за пределами исходного набора точек. - person pv.; 29.10.2013
comment
Спасибо @pv. Я думал, что это может быть так, но, поскольку я не видел в своем 2-мерном тестовом примере, что я могу построить график, я не был уверен. - person Chris H; 29.10.2013
comment
Документация Qhull рекомендует вычислять выпуклую оболочку вместо делене для каждого Вороной области, но в остальном рекомендация там совпадает с тем, что вы делаете выше, так что это, вероятно, лучший доступный способ. В принципе, оболочки Scipy Qhull также могут вычислять объемы выпуклой оболочки, но не похоже, что Qhull предлагает способы прямого получения объемов области вороного без дополнительных вычислений выпуклой оболочки. - person pv.; 04.11.2013
comment
@pv. Наконец-то я нашел что-то, прямо говорящее, что Qhull не делает том доступным, поэтому этот способ кажется наиболее подходящим. Факторы, связанные с вычислением, теперь малы по сравнению с другими факторами, о которых мне приходится беспокоиться в этом анализе, так что, по крайней мере, для меня это отсортировано. - person Chris H; 04.11.2013