Python:从 3D 中 Scipy 的 Delaunay 三角剖分计算 Voronoi 细分 [英] Python: Calculate Voronoi Tesselation from Scipy's Delaunay Triangulation in 3D

问题描述

我有大约 50,000 个 3D 数据点，我在这些数据点上运行了来自新 scipy 的 scipy.spatial.Delaunay(我使用的是 0.10)，这给了我一个非常有用的三角剖分.

I have about 50,000 data points in 3D on which I have run scipy.spatial.Delaunay from the new scipy (I'm using 0.10) which gives me a very useful triangulation.

基于:http://en.wikipedia.org/wiki/Delaunay_triangulation(与 Voronoi 的关系"部分图")

Based on: http://en.wikipedia.org/wiki/Delaunay_triangulation (section "Relationship with the Voronoi diagram")

...我想知道是否有一种简单的方法可以得到这个三角剖分的双图"，即 Voronoi Tesselation.

...I was wondering if there is an easy way to get to the "dual graph" of this triangulation, which is the Voronoi Tesselation.

有什么线索吗?我在这方面的搜索似乎没有显示 scipy 中的任何预先构建的函数，我觉得这几乎很奇怪！

Any clues? My searching around on this seems to show no pre-built in scipy functions, which I find almost strange!

谢谢，爱德华

推荐答案

邻接信息可以在 Delaunay 对象的 neighbors 属性中找到.不幸的是，该代码目前并未向用户公开外心，因此您必须自己重新计算.

The adjacency information can be found in the neighbors attribute of the Delaunay object. Unfortunately, the code does not expose the circumcenters to the user at the moment, so you'll have to recompute those yourself.

另外，延伸到无穷远的Voronoi边也不是通过这种方式直接获得的.这仍然有可能，但需要更多思考.

Also, the Voronoi edges that extend to infinity are not directly obtained in this way. It's still probably possible, but needs some more thinking.

import numpy as np
from scipy.spatial import Delaunay

points = np.random.rand(30, 2)
tri = Delaunay(points)

p = tri.points[tri.vertices]

# Triangle vertices
A = p[:,0,:].T
B = p[:,1,:].T
C = p[:,2,:].T

# See http://en.wikipedia.org/wiki/Circumscribed_circle#Circumscribed_circles_of_triangles
# The following is just a direct transcription of the formula there
a = A - C
b = B - C

def dot2(u, v):
    return u[0]*v[0] + u[1]*v[1]

def cross2(u, v, w):
    """u x (v x w)"""
    return dot2(u, w)*v - dot2(u, v)*w

def ncross2(u, v):
    """|| u x v ||^2"""
    return sq2(u)*sq2(v) - dot2(u, v)**2

def sq2(u):
    return dot2(u, u)

cc = cross2(sq2(a) * b - sq2(b) * a, a, b) / (2*ncross2(a, b)) + C

# Grab the Voronoi edges
vc = cc[:,tri.neighbors]
vc[:,tri.neighbors == -1] = np.nan # edges at infinity, plotting those would need more work...

lines = []
lines.extend(zip(cc.T, vc[:,:,0].T))
lines.extend(zip(cc.T, vc[:,:,1].T))
lines.extend(zip(cc.T, vc[:,:,2].T))

# Plot it
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection

lines = LineCollection(lines, edgecolor='k')

plt.hold(1)
plt.plot(points[:,0], points[:,1], '.')
plt.plot(cc[0], cc[1], '*')
plt.gca().add_collection(lines)
plt.axis('equal')
plt.xlim(-0.1, 1.1)
plt.ylim(-0.1, 1.1)
plt.show()

这篇关于Python:从 3D 中 Scipy 的 Delaunay 三角剖分计算 Voronoi 细分的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！