我在三维空间中有大约50,000个数据点，运行了新版scipy（版本号为0.10）的scipy.spatial.Delaunay函数获得了非常有用的三角剖分结果。
根据维基百科对Delaunay三角剖分的解释（链接如下： http://en.wikipedia.org/wiki/Delaunay_triangulation，"Relationship with the Voronoi diagram"一节），我想知道是否有一种简单的方法可以获取这个三角剖分的“对偶图”，也就是Voronoi图。
有任何线索吗？我的搜索似乎没有找到scipy中内置的相应功能，这让我感到相当奇怪！
谢谢， Edward
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()
由于我花费了相当多的时间，因此我想分享如何获取Voronoi 多边形而不仅仅是边缘的解决方案。
代码位于https://gist.github.com/letmaik/8803860，并扩展了tauran的解决方案。
首先，我更改了代码，使其分别给出顶点和（成对的）索引（=边缘），因为在处理索引而不是点坐标时可以简化许多计算。
然后，在voronoi_cell_lines方法中，我确定哪些边属于哪些单元格。为此，我使用了相关问题的Alink提出的解决方案。也就是说，对于每个边，找到两个最近的输入点（=单元格），并创建一个从那里映射过来的映射。
最后一步是创建实际的多边形（请参见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的答案和其他网上找到的代码片段，我制定了一个解决方案。该解决方案返回一个完整的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 重命名为更有意义的名称吗？我正在尝试理解无限部分。谢谢！ - letmaik我不知道有没有一个函数可以做到这一点，但这似乎不是一个过于复杂的任务。
Voronoi图是圆周圆的交汇处，正如维基百科文章所描述的那样。
因此，您可以从一个查找三角形圆周圆心的函数开始，这是基本数学知识（http://en.wikipedia.org/wiki/Circumscribed_circle）。
然后，只需连接相邻三角形的中心即可。
ncross2函数接受u和v作为参数，但计算的值仅取决于a和b。也许应该用u和v替换a和b？ - unutbu