获取 Voronoi 图中特定点对应区域的面积

编程问答 2022-06-01

问题描述

使用 this answer，我可以创建有界 Voronoi 图（此代码归功于 @Flabetvibes）：

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

eps = sys.float_info.epsilon

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

points = np.array([[0.17488374,0.36498964],[0.94904866,0.80085891],[0.89265224,0.4160692 ],[0.17035869,0.82769497],[0.30274931,0.04572908],[0.40515272,0.1445514 ],[0.23191921,0.08250689],[0.48713553,0.94806717],[0.77714412,0.46517511],[0.25945989,0.76444964]])

vor = voronoi(points,(0,1,1))

fig = plt.figure()
ax = fig.gca()
# Plot initial points
ax.plot(vor.filtered_points[:,vor.filtered_points[:,'b.')
# Plot ridges points
for region in vor.filtered_regions:
    vertices = vor.vertices[region,:]
    ax.plot(vertices[:,vertices[:,'go')
# Plot ridges
for region in vor.filtered_regions:
    vertices = vor.vertices[region + [region[0]],'k-')

enter image description here

现在，我想获取包含蓝色原始点之一的区域的面积，例如点 [0]。在这个例子中，points[0] 是点 (0.17488374,0.36498964)。我以为我可以使用以下代码找到这一点的区域：

area = ConvexHull(vor.vertices[vor.filtered_regions[0],:]).volume

因为我认为 points[0] 中的 0 索引将对应 vor.filtered_regions[0] 中的 0 索引。但它没有 - vor.filtered_regions[9] 实际上是我正在寻找的（我手动计算出来的，但我希望它是自动化的）。在另一个示例中，索引为 2 的区域正是我要查找的区域，因此它看起来也不一致。

有没有办法找到 vor.filtered_regions 的索引，它会给我我想要的区域？或者还有其他方法可以解决这个问题吗？即使我用所有 10 个点创建整个 Voronoi 图，点 [0] 区域的面积就是我真正要寻找的（虽然仍然有界），所以我假设可能有一个更快的这样做的方法，但我不知道那可能是什么。

解决方法

scipy Voronoi 图的 point_region 属性告诉您哪个区域与哪个点相关联。因此，您可以使用该数据来查找关联的区域。

这是您的 voronoi 函数的简化版本，它使用该属性来确保 filted_points 和 filters_regions 的构造一致，即第一个区域是与第一个点相关联的区域。

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)
    # Compute Voronoi
    vor = sp.spatial.Voronoi(points)
    # Filter regions
    regions = []
    [vor.point_region[i] for i in range(10)]

    vor.filtered_points = points_center
    vor.filtered_regions = [vor.regions[vor.point_region[i]] for i in range(len(points_center))]
    return vor
pythonscipyvoronoi