从Voronoi单元获取有界多边形坐标

Question

我有点（例如，电池塔位置的纬度，经度对），我需要获取它们形成的Voronoi细胞的多边形。

from scipy.spatial import Voronoi

tower = [[ 24.686 ,  46.7081],
       [ 24.686 ,  46.7081],
       [ 24.686 ,  46.7081]]

c = Voronoi(towers)

现在，我需要获取每个单元格的经度，经度坐标的多边形边界（以及该多边形所包围的质心是什么）。我也需要这个Voronoi。这意味着边界不会达到无穷大，而是在边界框内。

Answer 1

考虑到矩形边界框，我的第一个想法是定义此边界框与scipy.spatial.Voronoi生成的Voronoï图之间的相交运算。一个想法不一定很棒，因为这需要编写大量计算几何的基本函数。

但是，这是我想到的第二个主意（技巧？）：用于计算平面中一组scipy.spatial.Voronoi点的Voronoï图的算法的时间复杂度为n。如何添加点以约束初始点位于边界框中的Voronoï单元呢？

有界Voronoï图的解决方案

一张图片值得发表精彩演讲：

O(n ln(n))

我在这里做什么？那很简单！初始点（蓝色）位于中。然后，根据[0.0, 1.0] x [0.0, 1.0]（边界框的左边缘），通过反射对称性获得左侧（即[-1.0, 0.0] x [0.0, 1.0]）上的点（蓝色）。利用根据x = 0.0，x = 1.0和y = 0.0（边界框的其他边缘）的反射对称性，我得到了需要做的所有点（蓝色）。

然后我运行y = 1.0。上一张图片描绘了生成的Voronoï图（我使用scipy.spatial.Voronoi）。

下一步该怎么办？只需根据边界框过滤点，边或面。并根据众所周知的公式求出每张脸的质心，以计算出scipy.spatial.voronoi_plot_2d。这是结果的图像（质心为红色）：

scipy.spatial.voronoi_plot_2d

显示代码前的乐趣

太好了！它似乎有效。如果经过一次迭代后，我尝试对质心（红色）而不是初始点（蓝色）重新运行算法，该怎么办？如果我一次又一次尝试怎么办？

步骤2

centroid of polygon

步骤10

步骤25

很酷！ Voronoï细胞倾向于使其energy ...

最小化

这里是代码

Answer 2

我在使用scipy的voronoi函数和创建CVD时遇到了很多麻烦，因此这些精彩的帖子和评论对我们很有帮助。作为一名编程新手，我试图从Flabetvvibes的答案中理解代码，我将分享我对它如何与Energya以及我自己的修改一起工作的解释。我还在此答案的底部完整地发布了我的代码版本]

in_box函数使用numpy的logical_and方法返回一个布尔数组，该布尔数组表示塔中哪些坐标位于边界框中。

import matplotlib.pyplot as pl
import numpy as np
import scipy as sp
import scipy.spatial
import sys

eps = sys.float_info.epsilon

n_towers = 100
towers = np.random.rand(n_towers, 2)
bounding_box = np.array([0., 1., 0., 1.]) # [x_min, x_max, y_min, y_max]

def in_box(towers, bounding_box):
    return np.logical_and(np.logical_and(bounding_box[0] <= towers[:, 0],
                                         towers[:, 0] <= bounding_box[1]),
                          np.logical_and(bounding_box[2] <= towers[:, 1],
                                         towers[:, 1] <= bounding_box[3]))


def voronoi(towers, bounding_box):
    # Select towers inside the bounding box
    i = in_box(towers, bounding_box)
    # Mirror points
    points_center = towers[i, :]
    points_left = np.copy(points_center)
    points_left[:, 0] = bounding_box[0] - (points_left[:, 0] - bounding_box[0])
    points_right = np.copy(points_center)
    points_right[:, 0] = bounding_box[1] + (bounding_box[1] - points_right[:, 0])
    points_down = np.copy(points_center)
    points_down[:, 1] = bounding_box[2] - (points_down[:, 1] - bounding_box[2])
    points_up = np.copy(points_center)
    points_up[:, 1] = bounding_box[3] + (bounding_box[3] - points_up[:, 1])
    points = np.append(points_center,
                       np.append(np.append(points_left,
                                           points_right,
                                           axis=0),
                                 np.append(points_down,
                                           points_up,
                                           axis=0),
                                 axis=0),
                       axis=0)
    # Compute Voronoi
    vor = sp.spatial.Voronoi(points)
    # Filter regions
    regions = []
    for region in vor.regions:
        flag = True
        for index in region:
            if index == -1:
                flag = False
                break
            else:
                x = vor.vertices[index, 0]
                y = vor.vertices[index, 1]
                if not(bounding_box[0] - eps <= x and x <= bounding_box[1] + eps and
                       bounding_box[2] - eps <= y and y <= bounding_box[3] + eps):
                    flag = False
                    break
        if region != [] and flag:
            regions.append(region)
    vor.filtered_points = points_center
    vor.filtered_regions = regions
    return vor

def centroid_region(vertices):
    # Polygon's signed area
    A = 0
    # Centroid's x
    C_x = 0
    # Centroid's y
    C_y = 0
    for i in range(0, len(vertices) - 1):
        s = (vertices[i, 0] * vertices[i + 1, 1] - vertices[i + 1, 0] * vertices[i, 1])
        A = A + s
        C_x = C_x + (vertices[i, 0] + vertices[i + 1, 0]) * s
        C_y = C_y + (vertices[i, 1] + vertices[i + 1, 1]) * s
    A = 0.5 * A
    C_x = (1.0 / (6.0 * A)) * C_x
    C_y = (1.0 / (6.0 * A)) * C_y
    return np.array([[C_x, C_y]])

vor = voronoi(towers, bounding_box)

fig = pl.figure()
ax = fig.gca()
# Plot initial points
ax.plot(vor.filtered_points[:, 0], vor.filtered_points[:, 1], 'b.')
# Plot ridges points
for region in vor.filtered_regions:
    vertices = vor.vertices[region, :]
    ax.plot(vertices[:, 0], vertices[:, 1], 'go')
# Plot ridges
for region in vor.filtered_regions:
    vertices = vor.vertices[region + [region[0]], :]
    ax.plot(vertices[:, 0], vertices[:, 1], 'k-')
# Compute and plot centroids
centroids = []
for region in vor.filtered_regions:
    vertices = vor.vertices[region + [region[0]], :]
    centroid = centroid_region(vertices)
    centroids.append(list(centroid[0, :]))
    ax.plot(centroid[:, 0], centroid[:, 1], 'r.')

ax.set_xlim([-0.1, 1.1])
ax.set_ylim([-0.1, 1.1])
pl.savefig("bounded_voronoi.png")

sp.spatial.voronoi_plot_2d(vor)
pl.savefig("voronoi.png")

Flabetvvibes镜像这些点，以使边界框内边缘的区域有限。 Scipy的voronoi方法针对未定义的顶点返回-1，因此对这些点进行镜像可以使边界框内的所有区域都是有限的，并且所有无限区域都位于边界框外的镜像区域中，稍后将丢弃。

import matplotlib.pyplot as pl
import numpy as np
import scipy as sp
import scipy.spatial
import sys
import copy

eps = sys.float_info.epsilon

# Returns a new np.array of towers that within the bounding_box
def in_box(towers, bounding_box):
    return np.logical_and(np.logical_and(bounding_box[0] <= towers[:, 0],
                                         towers[:, 0] <= bounding_box[1]),
                          np.logical_and(bounding_box[2] <= towers[:, 1],
                                         towers[:, 1] <= bounding_box[3]))

bounded_voronoi方法的最后这一点调用scipy的voronoi函数，并为边界框中的已过滤点和区域添加新属性。 Energya建议删除Flabetvvibe的代码，该代码使用一个衬板手动找到边界框内的所有有限区域，该衬板获得区域的前五分之一，这是原始输入以及组成边界框的点。

# Generates a bounded vornoi diagram with finite regions in the bounding box
def bounded_voronoi(towers, bounding_box):
    # Select towers inside the bounding box
    i = in_box(towers, bounding_box)

    # Mirror points left, right, above, and under to provide finite regions for the
    # edge regions of the bounding box
    points_center = towers[i, :]

    points_left = np.copy(points_center)
    points_left[:, 0] = bounding_box[0] - (points_left[:, 0] - bounding_box[0])

    points_right = np.copy(points_center)
    points_right[:, 0] = bounding_box[1] + (bounding_box[1] - points_right[:, 0])

    points_down = np.copy(points_center)
    points_down[:, 1] = bounding_box[2] - (points_down[:, 1] - bounding_box[2])

    points_up = np.copy(points_center)
    points_up[:, 1] = bounding_box[3] + (bounding_box[3] - points_up[:, 1])

    points = np.append(points_center,
                       np.append(np.append(points_left,
                                           points_right,
                                           axis=0),
                                 np.append(points_down,
                                           points_up,
                                           axis=0),
                                 axis=0),
                       axis=0)

我采用了Flabetvvibe的代码，该代码执行了loyd算法的迭代，并将其形成一种易于使用的方法。对于每次迭代，将调用先前的bounded_voronoi函数，然后为每个单元找到质心，并且它们将成为下一次迭代的新点集。

# Compute Voronoi
vor = sp.spatial.Voronoi(points)

# creates a new attibute for points that form the diagram within the region
vor.filtered_points = points_center 
# grabs the first fifth of the regions, which are the original regions
vor.filtered_regions = np.array(vor.regions)[vor.point_region[:vor.npoints//5]]

return vor

只需获取当前区域的所有顶点，并将第一个顶点复制到末尾，以便第一个和最后一个顶点用于计算质心相同。

感谢Flabetvvibes和Energya。您的帖子/答案教会了我如何比它的文档更好地使用scipy的voronoi方法。我也将代码发布为一个单独的代码，放在其他任何要复制/粘贴的代码下面。

def generate_CVD(points, iterations, bounding_box):
    p = copy.copy(points)

    for i in range(iterations):
        vor = bounded_voronoi(p, bounding_box)
        centroids = []

        for region in vor.filtered_regions:
            # grabs vertices for the region and adds a duplicate
            # of the first one to the end
            vertices = vor.vertices[region + [region[0]], :] 
            centroid = centroid_region(vertices)
            centroids.append(list(centroid[0, :]))

        p = np.array(centroids)

    return bounded_voronoi(p, bounding_box)