python - Voronoi - 计算每个区域的确切边界

Question

我正在尝试使用 scipy.spatial.Voronoi 计算 Voronoi 图每个区域的确切边界，前提是所有点都在预定义的多边形内。例如，使用此 documentation 中的示例.

如果我需要计算具有相同点但位于具有以下边界的矩形内的 Voroni 怎么办

global_boundaries = np.array([[-2, -2], [4, -2], [4, 4], [-2, 4], [-2, -2]])

我需要像那样计算每个 Voronoi 区域的精确边界吗？

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 个区域，而不是

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]]

如何计算无限脊的缺失端点？

我已尝试修改此代码 http://nbviewer.ipython.org/gist/pv/8037100

与此问题相关Colorize Voronoi Diagram

但它仅适用于圆形边界。我修改了它，考虑了一个半径，使我的区域完全在圆内，然后计算连接点和圆周的线与边界之间的交点。它有效，但仅适用于第一点，之后我得到“GEOMETRYCOLLECTION EMPTY”结果。

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())

有没有人解决过类似的问题？

有人能帮忙吗？

Answer 1

1。算法

我建议采用以下两步法:

1. 首先，为每个 Voronoi 区域制作一个凸多边形。在无限区域的情况下，通过将无限远的点分成两个足够远的点，由一条边连接起来。 (“足够远”意味着额外的边缘完全通过边界多边形之外。)

2. 使用 shapely 的 intersection 将步骤 (1) 中的每个多边形与边界多边形相交方法。

这种方法相对于 Ophir Cami's answer 的好处是它适用于非凸边界多边形，代码更简单一些。

2。例子

让我们从 Ophir Cami's answer 中的点的 Voronoi 图开始.无限脊由 scipy.spatial.voronoi_plot_2d 显示为虚线:

然后，我们为每个 Voronoi 区域构造一个凸多边形。这对于有限区域来说很容易，但我们必须拉远很远才能看到无限 Voronoi 区域会发生什么。对应于这些区域的多边形有一个额外的边缘，它距离足够远，完全位于边界多边形之外:

现在我们可以将每个 Voronoi 区域的多边形与边界多边形相交:

在这种情况下，所有 Voronoi 多边形都与边界多边形有非空交集，但在一般情况下，其中一些可能会消失。

3。代码

第一步是生成对应于 Voronoi 区域的多边形。与 Ophir Cami 一样，我从 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。

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 中的最后一个数字。