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

Question

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

如果我需要使用相同的点计算位于以下边界矩形内的Voronoi图会怎样？

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

有人曾解决过类似的问题吗？

有人可以帮忙吗？

- user3681833

Answer 1

1. 算法

我建议采用以下两步方法：

1. 首先，针对每个Voronoi区域创建一个凸多边形。在无限区域的情况下，通过将无穷远点分成两个足够远的点并以一条边相连来实现此操作。（“足够远”意味着额外的边完全经过边界多边形之外）。

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

这种方法比Ophir Cami的答案更好，因为它适用于非凸边界多边形，并且代码更简单。

2. 示例

让我们从Ophir Cami的答案中的点开始构建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节中的最后一个图形。

- Gareth Rees

---

我已经让这个例子运行起来了，但是使用 points = [[0, 6000042], [5, 6000060], [8, 6000076], [14, 6000016], [30, 6000069], [39, 6000000]]; boundary = [[-1, 5999999], [40, 5999999], [40, 6000077], [-1, 6000077], [-1, 5999999]] 失败了。它在 dir_j, dir_k = ridge_direction[i, j] 处引发了 ValueError: need more than 0 values to unpack。你有什么想法吗？然而，通过 shift = points.min(axis=0); points -= shift; boundary -= shift 将点和边界向原点移动，确实使它再次工作。这可能是为什么呢？ - anttikoo

我遇到了与@anttikoo相同的问题，涉及大量的地理坐标(WGS84)。在我的情况下，移动这些点是无效的。在本地系统(等距方位投影)中重新投影坐标也是无效的。但是，我仍然会沿着有效的最后一个山脊方向前进，该方向的 y 值接近于-1，我一直在想这是否会导致-1无限值的问题：ridge_direction[i, j]=[array([-0.08074507, -0.99673479])] - Arkeen

Answer 2

我使用了 voronoi_plot_2d 并进行了修改，如下所示。

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Voronoi
from shapely.geometry import Polygon, Point

# Voronoi - Compute exact boundaries of every region


def angle_between(v0, v1):
    return np.math.atan2(np.linalg.det([v0, v1]), np.dot(v0, v1))


def calc_angle(c0, c1, c2):
    return angle_between(np.array(c1) - np.array(c0), np.array(c2) - np.array(c1))


def is_convex(polygon):
    temp_coords = np.array(polygon.exterior.coords)
    temp_coords = np.vstack([temp_coords, temp_coords[1, :]])

    for i, (c0, c1, c2) in enumerate(zip(temp_coords, temp_coords[1:], temp_coords[2:])):
        if i == 0:
            first_angle_crit = calc_angle(c0, c1, c2) > 0
        elif (calc_angle(c0, c1, c2) > 0) != first_angle_crit:
            return False
    return True


def infinite_segments(vor_):
    line_segments = []
    center = vor_.points.mean(axis=0)
    for pointidx, simplex in zip(vor_.ridge_points, vor_.ridge_vertices):
        simplex = np.asarray(simplex)
        if np.any(simplex < 0):
            i = simplex[simplex >= 0][0]  # finite end Voronoi vertex

            t = vor_.points[pointidx[1]] - vor_.points[pointidx[0]]  # tangent
            t /= np.linalg.norm(t)
            n = np.array([-t[1], t[0]])  # normal

            midpoint = vor_.points[pointidx].mean(axis=0)
            direction = np.sign(np.dot(midpoint - center, n)) * n

            line_segments.append([(vor_.vertices[i, 0], vor_.vertices[i, 1]),
                                  (direction[0], direction[1])])
    return line_segments


class NotConvexException(Exception):
    def __str__(self):
        return 'The Polygon is not Convex!!!'


class NotAllPointsAreInException(Exception):
    def __str__(self):
        return 'Not all points are in the polygon!!!'


def intersect(p0, u, q0, q1):
    v = (q1 - q0)[np.newaxis].T
    A = np.hstack([u, -v])
    b = q0 - p0
    try:
        inv_A = np.linalg.inv(A)
    except np.linalg.LinAlgError:
        return np.nan, np.nan
    return np.dot(inv_A, b)


def _adjust_bounds(ax__, points_):
    ptp_bound = points_.ptp(axis=0)
    ax__.set_xlim(points_[:, 0].min() - 0.1*ptp_bound[0], points_[:, 0].max() + 0.1*ptp_bound[0])
    ax__.set_ylim(points_[:, 1].min() - 0.1*ptp_bound[1], points_[:, 1].max() + 0.1*ptp_bound[1])


def in_polygon(polygon, points_):
    return [polygon.contains(Point(x)) for x in points_]


def voronoi_plot_2d_inside_convex_polygon(vor_, polygon, ax__=None, **kw):
    from matplotlib.collections import LineCollection

    if not all(in_polygon(polygon, vor_.points_)):
        raise NotAllPointsAreInException()

    if not is_convex(polygon):
        raise NotConvexException()

    if vor_.points.shape[1] != 2:
        raise ValueError("Voronoi diagram is not 2-D")

    vor_inside_ind = np.array([i for i, x in enumerate(vor_.vertices) if polygon.contains(Point(x))])
    vor_outside_ind = np.array([i for i, x in enumerate(vor_.vertices) if not polygon.contains(Point(x))])
    ax__.plot(vor_.points[:, 0], vor_.points[:, 1], '.')
    if kw.get('show_vertices', True):
        ax__.plot(vor_.vertices[vor_inside_ind, 0], vor_.vertices[vor_inside_ind, 1], 'o')

    temp_coords = np.array(polygon.exterior.coords)
    line_segments = []
    for t0, t1 in zip(temp_coords, temp_coords[1:]):
        line_segments.append([t0, t1])
    ax__.add_collection(LineCollection(line_segments, colors='b', linestyle='solid'))
    line_segments = []
    for simplex in vor_.ridge_vertices:
        simplex = np.asarray(simplex)
        if np.all(simplex >= 0):
            if not all(in_polygon(polygon, vor_.vertices[simplex])):
                continue
            line_segments.append([(x, y) for x, y in vor_.vertices[simplex]])

    ax__.add_collection(LineCollection(line_segments, colors='k', linestyle='solid'))

    line_segments = infinite_segments(vor_)
    from_inside = np.array([x for x in line_segments if polygon.contains(Point(x[0]))])

    line_segments = []

    for f in from_inside:
        for coord0, coord1 in zip(temp_coords, temp_coords[1:]):
            s, t = intersect(f[0], f[1][np.newaxis].T, coord0, coord1)
            if 0 < t < 1 and s > 0:
                line_segments.append([f[0], f[0] + s * f[1]])
                break

    ax__.add_collection(LineCollection(np.array(line_segments), colors='k', linestyle='dashed'))

    line_segments = []

    for v_o_ind in vor_outside_ind:
        for simplex in vor_.ridge_vertices:
            simplex = np.asarray(simplex)
            if np.any(simplex < 0):
                continue
            if np.any(simplex == v_o_ind):
                i = simplex[simplex != v_o_ind][0]
                for coord0, coord1 in zip(temp_coords, temp_coords[1:]):
                    s, t = intersect(
                        vor_.vertices[i],
                        (vor_.vertices[v_o_ind] - vor_.vertices[i])[np.newaxis].T,
                        coord0,
                        coord1
                    )
                    if 0 < t < 1 and 0 < s < 1:
                        line_segments.append([
                            vor_.vertices[i],
                            vor_.vertices[i] + s * (vor_.vertices[v_o_ind] - vor_.vertices[i])
                        ])
                        break

    ax__.add_collection(LineCollection(np.array(line_segments), colors='r', linestyle='dashed'))

    _adjust_bounds(ax__, temp_coords)

    return ax__.figure

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

global_boundaries = Polygon([[-5, -2], [3.4, -2], [4.7, 4], [2.7, 5.7], [-1, 4]])
fig = plt.figure()
ax = fig.add_subplot(111)

vor = Voronoi(points)
voronoi_plot_2d_inside_convex_polygon(vor, global_boundaries, ax_=ax)
plt.show()

注意：

1. 多边形必须是凸的。
2. 所有的点都必须在多边形内。

颜色说明：

• 原始点为蓝色。
• Voronoi顶点为绿色。
• 有限内部Voronoi边缘为实线黑色。
• 有限外部Voronoi边缘为虚线红色。
• 无限Voronoi边缘为虚线黑色。

- Ophir Carmi