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

Question

我有点(例如,lat,lon对的细胞塔位置),我需要得到它们形成的Voronoi细胞的多边形.

from scipy.spatial import Voronoi

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

c = Voronoi(towers)

现在,我需要为每个单元格的lat,lon坐标获取多边形边界(以及该多边形所包围的质心).我需要这个Voronoi也是有限的.意味着边界不会变为无穷大,而是在边界框内.

Answer 1

给定一个矩形边界框,我的第一个想法是在这个边界框和Voronoï图之间定义一种交叉操作scipy.spatial.Voronoi.一个想法不一定很好,因为这需要编码计算几何的大量基本功能.

然而,这是我想到的第二个想法(黑客？):计算平面中一组n点的Voronoï图的算法具有时间复杂度O(n ln(n)).如何添加点来限制初始点的Voronoï单元格位于边界框中？

有界Voronoï图的解决方案

一张图片值得一个伟大的演讲:

我在这做了什么？这很简单!初始点(蓝色)在于[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(边框的其它边缘),我得到的所有的点(蓝色),我需要做的工作.

然后我跑了scipy.spatial.Voronoi.上图描绘了生成的Voronoï图(我使用scipy.spatial.voronoi_plot_2d).

接下来做什么？只需根据边界框过滤点,边或面.并根据众所周知的公式计算每个面的质心,以计算多边形的质心.这是结果的图像(质心是红色的):

显示代码之前有些乐趣

大!它似乎工作.如果在一次迭代后我尝试在质心(红色)而不是初始点(蓝色)上重新运行算法怎么办？如果我一次又一次地尝试怎么办？

第2步

第10步

第25步

凉!Voronoï细胞倾向于最小化它们的能量 ......

这是代码

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