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

Question

我正在尝试使用 scipy.spatial.Voronoi 计算 Voronoi 图的每个区域的确切边界，以防所有点都在预定义的多边形内。

例如，使用文档中的示例，

http://docs.scipy.org/doc/scipy-dev/reference/generated/scipy.spatial.Voronoi.html 如果我需要使用相同的点但在具有以下边界的矩形内计算 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 节中的最后一个数字。

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 脊以黑色虚线显示。