我有大约 50,000 个 3D 数据点，我从新的 scipy（我使用的是 0.10）运行 scipy.spatial.Delaunay，这给了我一个非常有用的三角测量。
基于：http://en.wikipedia.org/wiki/Delaunay_triangulation（“与 Voronoi 图的关系”部分）
...我想知道是否有一种简单的方法可以得到这个三角剖分的“双图”，即 Voronoi Tesselation。
有什么线索吗？我在这方面的搜索似乎没有显示任何预建的 scipy 函数，我觉得这几乎很奇怪！
谢谢， 爱德华
【问题讨论】：
标签： python 3d scipy delaunay voronoi
邻接信息可以在 Delaunay 对象的neighbors 属性中找到。不幸的是，代码目前并未向用户公开外心，因此您必须自己重新计算。
此外，延伸到无穷大的 Voronoi 边也不是通过这种方式直接获得的。仍有可能，但需要更多思考。
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()
【讨论】：
ncross2 接受 u 和 v 是参数，但计算的值仅取决于 a 和 b。也许a 和b 应该替换为u 和v？
由于我在这方面花费了大量时间，我想分享我的解决方案，了解如何获得 Voronoi 多边形 而不仅仅是边缘。
代码在https://gist.github.com/letmaik/8803860，并在tauran的解决方案上扩展。
首先，我更改了代码，分别给出了顶点和（成对）索引（=边），因为在处理索引而不是点坐标时可以简化许多计算。
然后，在voronoi_cell_lines 方法中，我确定哪些边属于哪些单元格。为此，我使用相关问题中提出的Alink 解决方案。也就是说，为每条边找到两个最近的输入点（=cells）并从中创建一个映射。
最后一步是创建实际的多边形（参见voronoi_polygons 方法）。首先，需要关闭具有悬垂边缘的外部单元格。这就像查看所有边缘并检查哪些边缘只有一个相邻边缘一样简单。可以有零个或两个这样的边缘。如果有两个，然后我通过引入一个额外的边缘来连接它们。
最后，需要将每个单元格中的无序边按正确的顺序排列以从中导出多边形。
用法是：
P = np.random.random((100,2))
fig = plt.figure(figsize=(4.5,4.5))
axes = plt.subplot(1,1,1)
plt.axis([-0.05,1.05,-0.05,1.05])
vertices, lineIndices = voronoi(P)
cells = voronoi_cell_lines(P, vertices, lineIndices)
polys = voronoi_polygons(cells)
for pIdx, polyIndices in polys.items():
poly = vertices[np.asarray(polyIndices)]
p = matplotlib.patches.Polygon(poly, facecolor=np.random.rand(3,1))
axes.add_patch(p)
X,Y = P[:,0],P[:,1]
plt.scatter(X, Y, marker='.', zorder=2)
plt.axis([-0.05,1.05,-0.05,1.05])
plt.show()
哪个输出：
该代码可能不适合大量输入点，在某些方面可以改进。不过，它可能对遇到类似问题的其他人有所帮助。
【讨论】：
我遇到了同样的问题，并根据 pv. 的答案和我在网上找到的其他代码 sn-ps 构建了一个解决方案。该解决方案返回一个完整的 Voronoi 图，包括不存在三角形邻居的外线。
#!/usr/bin/env python
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from scipy.spatial import Delaunay
def voronoi(P):
delauny = Delaunay(P)
triangles = delauny.points[delauny.vertices]
lines = []
# Triangle vertices
A = triangles[:, 0]
B = triangles[:, 1]
C = triangles[:, 2]
lines.extend(zip(A, B))
lines.extend(zip(B, C))
lines.extend(zip(C, A))
lines = matplotlib.collections.LineCollection(lines, color='r')
plt.gca().add_collection(lines)
circum_centers = np.array([triangle_csc(tri) for tri in triangles])
segments = []
for i, triangle in enumerate(triangles):
circum_center = circum_centers[i]
for j, neighbor in enumerate(delauny.neighbors[i]):
if neighbor != -1:
segments.append((circum_center, circum_centers[neighbor]))
else:
ps = triangle[(j+1)%3] - triangle[(j-1)%3]
ps = np.array((ps[1], -ps[0]))
middle = (triangle[(j+1)%3] + triangle[(j-1)%3]) * 0.5
di = middle - triangle[j]
ps /= np.linalg.norm(ps)
di /= np.linalg.norm(di)
if np.dot(di, ps) < 0.0:
ps *= -1000.0
else:
ps *= 1000.0
segments.append((circum_center, circum_center + ps))
return segments
def triangle_csc(pts):
rows, cols = pts.shape
A = np.bmat([[2 * np.dot(pts, pts.T), np.ones((rows, 1))],
[np.ones((1, rows)), np.zeros((1, 1))]])
b = np.hstack((np.sum(pts * pts, axis=1), np.ones((1))))
x = np.linalg.solve(A,b)
bary_coords = x[:-1]
return np.sum(pts * np.tile(bary_coords.reshape((pts.shape[0], 1)), (1, pts.shape[1])), axis=0)
if __name__ == '__main__':
P = np.random.random((300,2))
X,Y = P[:,0],P[:,1]
fig = plt.figure(figsize=(4.5,4.5))
axes = plt.subplot(1,1,1)
plt.scatter(X, Y, marker='.')
plt.axis([-0.05,1.05,-0.05,1.05])
segments = voronoi(P)
lines = matplotlib.collections.LineCollection(segments, color='k')
axes.add_collection(lines)
plt.axis([-0.05,1.05,-0.05,1.05])
plt.show()
黑线 = Voronoi 图，红线 = Delauny 三角形
【讨论】：
di 和ps 重命名为更有意义的名称？我试图理解无穷大的部分。谢谢！
我不知道有什么功能可以做到这一点，但这似乎不是一项过于复杂的任务。
Voronoi 图是外接圆的交汇点，如维基百科文章中所述。
因此，您可以从找到三角形外接圆中心的函数开始，这是基本数学 (http://en.wikipedia.org/wiki/Circumscribed_circle)。
然后，只需连接相邻三角形的中心。
【讨论】：