彩色化 Voronoi 图表

Question

我正在尝试着对使用 scipy.spatial.Voronoi 创建的 Voronoi 图进行着色。以下是我的代码：

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Voronoi, voronoi_plot_2d

# make up data points
points = np.random.rand(15,2)

# compute Voronoi tesselation
vor = Voronoi(points)

# plot
voronoi_plot_2d(vor)

# colorize
for region in vor.regions:
    if not -1 in region:
        polygon = [vor.vertices[i] for i in region]
        plt.fill(*zip(*polygon))

plt.show()

生成的图像：

从图像边缘可见，一些Voronoi区域未被着色。这是因为这些区域的Voronoi顶点索引中有一些设置为-1，即对于那些在Voronoi图形之外的顶点。根据文档：

regions:（int列表的列表，形状为（nregions，*））形成每个Voronoi区域的Voronoi顶点的索引。-1表示Voronoi图外的顶点。

为了也可以为这些区域着色，我试图从多边形中删除这些“外部”顶点，但这并没有起作用。我想，我需要填充一些位于图像区域边缘的点，但我似乎无法合理地实现这一点。

谁能帮忙？

- moooeeeep

---

以下解决方案非常有效-但如果有人想要一个快速的不完美的实用解决方案，可以交换polygon = [vor.vertices[k] for k in（y for y in vor.regions [vor.point_region [i]] if y> -1）]并跳过if not -1 in region：检查。这将切掉任何延伸到无限远处的区域的角落，因此如果您不绘制单元格边界，则效果更好。 - zephyr

Answer 1

Voronoi数据结构包含构建“无限点”的位置所需的所有信息。Qhull还将它们报告为简单的-1索引，因此Scipy不会为您计算它们。

https://gist.github.com/pv/8036995

http://nbviewer.ipython.org/gist/pv/8037100

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Voronoi

def voronoi_finite_polygons_2d(vor, radius=None):
    """
    Reconstruct infinite voronoi regions in a 2D diagram to finite
    regions.

    Parameters
    ----------
    vor : Voronoi
        Input diagram
    radius : float, optional
        Distance to 'points at infinity'.

    Returns
    -------
    regions : list of tuples
        Indices of vertices in each revised Voronoi regions.
    vertices : list of tuples
        Coordinates for revised Voronoi vertices. Same as coordinates
        of input vertices, with 'points at infinity' appended to the
        end.

    """

    if vor.points.shape[1] != 2:
        raise ValueError("Requires 2D input")

    new_regions = []
    new_vertices = vor.vertices.tolist()

    center = vor.points.mean(axis=0)
    if radius is None:
        radius = vor.points.ptp().max()

    # Construct a map containing all ridges for a given point
    all_ridges = {}
    for (p1, p2), (v1, v2) in zip(vor.ridge_points, vor.ridge_vertices):
        all_ridges.setdefault(p1, []).append((p2, v1, v2))
        all_ridges.setdefault(p2, []).append((p1, v1, v2))

    # Reconstruct infinite regions
    for p1, region in enumerate(vor.point_region):
        vertices = vor.regions[region]

        if all(v >= 0 for v in vertices):
            # finite region
            new_regions.append(vertices)
            continue

        # reconstruct a non-finite region
        ridges = all_ridges[p1]
        new_region = [v for v in vertices if v >= 0]

        for p2, v1, v2 in ridges:
            if v2 < 0:
                v1, v2 = v2, v1
            if v1 >= 0:
                # finite ridge: already in the region
                continue

            # Compute the missing endpoint of an infinite ridge

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

            midpoint = vor.points[[p1, p2]].mean(axis=0)
            direction = np.sign(np.dot(midpoint - center, n)) * n
            far_point = vor.vertices[v2] + direction * radius

            new_region.append(len(new_vertices))
            new_vertices.append(far_point.tolist())

        # sort region counterclockwise
        vs = np.asarray([new_vertices[v] for v in new_region])
        c = vs.mean(axis=0)
        angles = np.arctan2(vs[:,1] - c[1], vs[:,0] - c[0])
        new_region = np.array(new_region)[np.argsort(angles)]

        # finish
        new_regions.append(new_region.tolist())

    return new_regions, np.asarray(new_vertices)

# make up data points
np.random.seed(1234)
points = np.random.rand(15, 2)

# compute Voronoi tesselation
vor = Voronoi(points)

# plot
regions, vertices = voronoi_finite_polygons_2d(vor)
print "--"
print regions
print "--"
print vertices

# colorize
for region in regions:
    polygon = vertices[region]
    plt.fill(*zip(*polygon), alpha=0.4)

plt.plot(points[:,0], points[:,1], 'ko')
plt.xlim(vor.min_bound[0] - 0.1, vor.max_bound[0] + 0.1)
plt.ylim(vor.min_bound[1] - 0.1, vor.max_bound[1] + 0.1)

plt.show()

- pv.

---

可能是一个小错误，不确定在新版本的numpy中是否已更改，但执行.ptp()会找到最大值和最小值之间的差异，然后.max()什么也不做。我认为你想要的是.ptp(axis=0).max()。 - Ehsan Kia

如果我提供x < 6，则某些区域保持白色。 我猜这种情况下端点未正确计算，或者我遗漏了什么？ - Karl Anka

这段代码有两个问题：1）半径可能需要任意大。2）你正在延伸/重构的半线脊的方向（direction = np.sign(np.dot(midpoint - center, n)) * n）并不总是正确的。我一直在努力开发一个能够始终正常工作的算法，但我还没有成功。 - Michele Piccolini

是的，这段代码肯定有问题。我正在使用投影数据点（墨卡托投影），代码中的一些远点变成了负数。 - n3rd

在这个解决方案中，有时会在以下位置出现键错误： # reconstruct a non-finite region ridges = all_ridges[p1] 不确定原因是什么。有什么见解吗？ - sandeepsign

Answer 2

我对这个问题有一个更简单的解决方案，那就是在调用 Voronoi 算法之前，将 4 个远离的虚拟点添加到您的点列表中。

根据您的代码，我添加了两行。

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Voronoi, voronoi_plot_2d

# make up data points
points = np.random.rand(15,2)

# add 4 distant dummy points
points = np.append(points, [[999,999], [-999,999], [999,-999], [-999,-999]], axis = 0)

# compute Voronoi tesselation
vor = Voronoi(points)

# plot
voronoi_plot_2d(vor)

# colorize
for region in vor.regions:
    if not -1 in region:
        polygon = [vor.vertices[i] for i in region]
        plt.fill(*zip(*polygon))

# fix the range of axes
plt.xlim([0,1]), plt.ylim([0,1])

plt.show()

然后得到的图形就像下面这样。

- Arrows

---

我有同样的图表，但是你有任何想法如何在没有橙色点的情况下绘制它吗？ - user2473664

只需使用voronoi_plot_2d(vor, show_vertices = False)。 - Arrows

这很引人入胜。但是，角度将会错误……您的坐标值与数据集中心点的比率越大，它就越精确，但仍然不准确。 - roberto tomás

@robertotomás 你确定吗？因为在我看来，似乎虚拟点只会影响兴趣区域外的Voronoi Ridge。 - sluki

@sluki非常感谢您在这里的评论和跟进。我原本假定这会是这样，但实际上当我去测试时并不是这样的。你可以找到边界大小和外部边界之间的静态比例，以确保角度正确。它有效地夹紧到正确的解决方案..所以例如在我的情况下，使用一个10个单位的正方形和一个1000个单位的边界（可能有些过度杀伤），这对于我来说似乎适用于所有情况。我在发现这个评论之后就把它丢了，否则我会早点跟进的。 - roberto tomás

Answer 3

我认为从vor结构中可用的数据中并没有足够的信息来解决这个问题，除非再次进行至少一些voronoi计算。既然如此，以下是原始的voronoi_plot_2d函数的相关部分，您应该能够使用它们提取与vor.max_bound或vor.min_bound相交的点，这些点是图表的左下角和右上角，以便找出多边形的其他坐标。

for simplex in vor.ridge_vertices:
    simplex = np.asarray(simplex)
    if np.all(simplex >= 0):
        ax.plot(vor.vertices[simplex,0], vor.vertices[simplex,1], 'k-')

ptp_bound = vor.points.ptp(axis=0)
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
        far_point = vor.vertices[i] + direction * ptp_bound.max()

        ax.plot([vor.vertices[i,0], far_point[0]],
                [vor.vertices[i,1], far_point[1]], 'k--')

- Troy Rockwood

---

我本来希望能够避免自己实现多边形点的计算。但是感谢指出vor.min_bound和vor.max_bound（之前没有注意到）。这些对于这个任务将会很有用，voronoi_plot_2d()的代码也是如此。 - moooeeeep