问题描述
在所有点都在预定义多边形内的情况下，我正在尝试使用scipy.spatial.Voronoi计算Voronoi图每个区域的确切边界.
I'm trying to compute the exact boundaries of every region of a Voronoi Diagram using scipy.spatial.Voronoi, in the case when all the points are inside a pre-defined polygon.
例如，使用文档中的示例，
For example, using the example in the documentation,
global_boundaries = np.array([[-2, -2], [4, -2], [4, 4], [-2, 4], [-2, -2]])
我需要计算每个voronoi地区的精确边界吗?
and I need to compute the precise boundaries of every voronoi region, like that?
voronoi_region_1_boundaries = [[-2, -2], [0.5, -2], [0.5, 0.5], [-2, 0-5], [-2, -2]]
voronoi_region_2_boundaries = [[-2, 1.5], [0.5, 1.5], [0.5, 4], [-2, 4], [-2, 1.5]]
voronoi_region_3_boundaries = [[-2, 0.5], [0.5, 0.5], [0.5, 1.5], [-2, 1.5], [-2, 0.5]]
对于所有9个区域，依次为
等等，
and so on for all the 9 regions, instead of
vor.regions
[[], [-1, 0], [-1, 1], [1, -1, 0], [3, -1, 2], [-1, 3], [-1, 2], [3, 2, 0, 1], [2, -1, 0], [3, -1, 1]]
如何计算无限脊的终点?
How do I compute the missing endpoint of an infinite ridge?
我尝试修改此代码 http://nbviewer.ipython.org/gist /pv/8037100
与此问题相关的为Voronoi图着色
，但仅适用于圆形边界.我已经考虑到半径来修改了它，以使我的区域完全在圆内，然后计算了连接点与圆周和边界的线之间的交点.它有效，但仅适用于第一点，其后，我得到了"GEOMETRYCOLLECTION EMPTY".
but it's working only for rounded boundaries.I've modified it considering a radius such that my area is completely inside the circle, and then computing the intersection between the line connecting the points and the circumference and the boundaries. It works, but only for the first point and after that I have "GEOMETRYCOLLECTION EMPTY" as a result.
direction = np.sign(np.dot(midpoint - center, n)) * n
super_far_point = vor.vertices[v2] + direction * radius
line_0 = LineString([midpoint, super_far_point])
for i in range(0, len(map_boundaries)-1):
i += 1
line_i = LineString([(map_boundaries[i-1]), (map_boundaries[i])])
if line_0.intersection(line_i) != 0:
far_point = line_0.intersection(line_i)
new_region.append(len(new_vertices))
new_vertices.append(far_point.tolist())
有人解决过类似的问题吗?
Has anyone ever solved a similar problem?
任何人都可以帮忙吗?
推荐答案
1.算法
我建议采用以下两步方法:
1. Algorithm
I suggest the following two-step approach:
-
首先，为每个Voronoi区域制作一个凸多边形.对于无限区域，可通过将无限远处的点分成两个足够远的点(由边连接)来实现. (足够远"表示多余的边缘完全通过边界多边形的外部.)
First, make a convex polygon for each of the Voronoi regions. In the case of the infinite regions, do this by splitting the point at infinity into two points that are far enough away, joined by an edge. ("Far enough away" means that the extra edge passes entirely outside the bounding polygon.)
使用shapely的 intersection 方法.
Intersect each of the polygons from step (1) with the bounding polygon, using shapely's intersection method.
此方法相对于 Ophir Cami的答案的好处在于，它可以与非凸边界多边形配合使用，并且代码比较简单.
The benefit of this approach over Ophir Cami's answer is that it works with non-convex bounding polygons, and the code is a little simpler.
让我们从Voronoi图开始，了解 Ophir Cami的答案. scipy.spatial.voronoi_plot_2d :
Let's start with the Voronoi diagram for the points from Ophir Cami's answer. The infinite ridges are shown as dashed lines by scipy.spatial.voronoi_plot_2d:
然后为每个Voronoi区域构造一个凸多边形.这对于有限区域很容易，但是我们必须缩小很长的距离才能看到无限Voronoi地区发生了什么.与这些区域相对应的多边形具有一个额外的边缘，该边缘足够远，以至于它完全位于边界多边形的外部:
Then for each Voronoi region we construct a convex polygon. This is easy for the finite regions, but we have to zoom out a long way to see what happens to the infinite Voronoi regions. The polygons corresponding to these regions have an extra edge that is far enough away that it lies entirely outside the bounding polygon:
现在，我们可以将每个Voronoi区域的多边形与边界多边形相交:
Now we can intersect the polygons for each Voronoi region with the bounding polygon:
在这种情况下，所有Voronoi多边形都与边界多边形具有非空相交，但是在一般情况下，其中一些可能会消失.
In this case, all the Voronoi polygons have non-empty intersection with the bounding polygon, but in the general case some of them might vanish.
第一步是生成与Voronoi区域相对应的多边形.像Ophir Cami一样，我从 .
The first step is to generate polygons corresponding to the Voronoi regions. Like Ophir Cami, I derived this from the implementation of scipy.spatial.voronoi_plot_2d.
from collections import defaultdict
from shapely.geometry import Polygon
def voronoi_polygons(voronoi, diameter):
"""Generate shapely.geometry.Polygon objects corresponding to the
regions of a scipy.spatial.Voronoi object, in the order of the
input points. The polygons for the infinite regions are large
enough that all points within a distance 'diameter' of a Voronoi
vertex are contained in one of the infinite polygons.
"""
centroid = voronoi.points.mean(axis=0)
# Mapping from (input point index, Voronoi point index) to list of
# unit vectors in the directions of the infinite ridges starting
# at the Voronoi point and neighbouring the input point.
ridge_direction = defaultdict(list)
for (p, q), rv in zip(voronoi.ridge_points, voronoi.ridge_vertices):
u, v = sorted(rv)
if u == -1:
# Infinite ridge starting at ridge point with index v,
# equidistant from input points with indexes p and q.
t = voronoi.points[q] - voronoi.points[p] # tangent
n = np.array([-t[1], t[0]]) / np.linalg.norm(t) # normal
midpoint = voronoi.points[[p, q]].mean(axis=0)
direction = np.sign(np.dot(midpoint - centroid, n)) * n
ridge_direction[p, v].append(direction)
ridge_direction[q, v].append(direction)
for i, r in enumerate(voronoi.point_region):
region = voronoi.regions[r]
if -1 not in region:
# Finite region.
yield Polygon(voronoi.vertices[region])
continue
# Infinite region.
inf = region.index(-1) # Index of vertex at infinity.
j = region[(inf - 1) % len(region)] # Index of previous vertex.
k = region[(inf + 1) % len(region)] # Index of next vertex.
if j == k:
# Region has one Voronoi vertex with two ridges.
dir_j, dir_k = ridge_direction[i, j]
else:
# Region has two Voronoi vertices, each with one ridge.
dir_j, = ridge_direction[i, j]
dir_k, = ridge_direction[i, k]
# Length of ridges needed for the extra edge to lie at least
# 'diameter' away from all Voronoi vertices.
length = 2 * diameter / np.linalg.norm(dir_j + dir_k)
# Polygon consists of finite part plus an extra edge.
finite_part = voronoi.vertices[region[inf + 1:] + region[:inf]]
extra_edge = [voronoi.vertices[j] + dir_j * length,
voronoi.vertices[k] + dir_k * length]
yield Polygon(np.concatenate((finite_part, extra_edge)))
第二步是将Voronoi多边形与边界多边形相交.我们还需要选择一个合适的直径以传递到voronoi_polygons.
The second step is to intersect the Voronoi polygons with the bounding polygon. We also need to choose an appropriate diameter to pass to voronoi_polygons.
import matplotlib.pyplot as plt
from scipy.spatial import Voronoi
points = np.array([[0.1, -0.4], [0, 1.5], [0, 2.25], [1, 0], [1, 1], [1, 2],
[2, 0], [2.5, 1], [2, 2], [2.3, 2.3], [-0.5, -1.3], [-1.5, 3]])
boundary = np.array([[-5, -2], [3.4, -2], [4.7, 4], [2.7, 5.7], [-1, 4]])
x, y = boundary.T
plt.xlim(round(x.min() - 1), round(x.max() + 1))
plt.ylim(round(y.min() - 1), round(y.max() + 1))
plt.plot(*points.T, 'b.')
diameter = np.linalg.norm(boundary.ptp(axis=0))
boundary_polygon = Polygon(boundary)
for p in voronoi_polygons(Voronoi(points), diameter):
x, y = zip(*p.intersection(boundary_polygon).exterior.coords)
plt.plot(x, y, 'r-')
plt.show()
这绘制了上面§2中的最后一个数字.
This plots the last of the figures in §2 above.
这篇关于Voronoi-计算每个区域的确切边界的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！