scipy.spatial.Voronoiの使い方

問題設定

点群(points cloud)が与えられた時に，各点のボロノイ多面体の体積を求めたい．

pythonの例

scipyにあるspatial.Voronoiを使うとボロノイ多面体を簡単に計算してくれる．spatial.VoronoiはQhullという凸解析ライブラリのラッパのようだ．
但し，spatial.Voronoiには各ボロノイ多面体の体積を計算するアルゴリズムは無いようだ．

多面体の体積の計算

ボロノイ多面体の体積を計算するためのアルゴリズムは無いが，各ボロノイ領域の頂点集合は計算されているため，それらの頂点集合を多面体の頂点集合として多面体の体積が求まればボロノイ多面体の体積を計算できる．
stackoverflowにstackoverflow.com同様の質問があったので参考にする．計算方法はドロネー三角分割で，四面体に分割して，各四面体の体積の和を取るというもの（たぶん）．

# 
import numpy as np
from scipy.spatial import Delaunay
from scipy.spatial import ConvexHull


def tetrahedron_volume(a, b, c, d):
    return np.abs(np.einsum('ij,ij->i', a-d, np.cross(b-d, c-d))) / 6

def convex_hull_volume(pts):
    ch = ConvexHull(pts)
    dt = Delaunay(pts[ch.vertices])
    tets = dt.points[dt.simplices]
    return np.sum(tetrahedron_volume(tets[:, 0], tets[:, 1],
                                     tets[:, 2], tets[:, 3]))


#pts = np.random.rand(10, 3)
pts = np.array([[0,0,0], [1,0,0], [1,1,0], [0,1,0],\
                [0,0,1], [1,0,1], [1,1,1], [0,1,1],\
                [0.5,0.5,2.0] ]           )

dt = Delaunay(pts)
tets = dt.points[dt.simplices]
print dt.simplices
print tets
vol = np.sum(tetrahedron_volume(tets[:, 0], tets[:, 1], 
                                tets[:, 2], tets[:, 3]))
print pts
print vol
print convex_hull_volume(pts)

同様の質問もあったけど，こっちはあっているのか分からない．テストしてみると違うような値がでてくる・・・もう少し調べよう． scipy - Volume of Voronoi cell (python) - Stack Overflow

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)
print points
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
 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)