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()
相邻信息可以在Delaunay对象的
neighbors
属性中找到。不幸的是,代码目前没有向用户公开环心,因此您必须自己重新计算这些环心。同样,扩展到无穷远的Voronoi边也不是用这种方法直接得到的。可能还是有可能的,但还需要更多的思考。
由于我花了大量的时间在这个问题上,我想分享我的解决方案,如何得到Voronoi多边形而不仅仅是边。
代码位于https://gist.github.com/letmaik/8803860,并扩展到tauran的解。
首先,我修改了代码,分别给出了顶点和(对)索引(=边),因为在处理索引而不是点坐标时,许多计算都可以简化。
然后,在
voronoi_cell_lines
方法中,我确定哪些边属于哪个单元格。为此,我使用一个相关问题的Alink的建议解。也就是说,对于每条边,找到两个最近的输入点(=单元格),并从中创建映射。最后一步是创建实际的多边形(请参见
voronoi_polygons
方法)。首先,需要关闭具有悬挂边缘的外部单元。这就像查看所有边并检查哪些边只有一个相邻边一样简单。这样的边可以是零,也可以是两条。在两个例子中,我通过引入一个额外的边来连接它们。最后,每个单元中的无序边需要按正确的顺序排列,以便从中派生多边形。
用法是:
哪些输出:
该代码可能不适合大量的输入点,在某些方面可以改进。不过,这可能对其他有类似问题的人有所帮助。
我遇到了同样的问题,并根据pv的答案和我在网上找到的其他代码片段构建了一个解决方案。该解返回一个完整的Voronoi图,包括不存在三角形邻域的外线。
黑线=Voronoi图,红线=Delauny三角形
相关问题 更多 >
编程相关推荐