搜索最优算法以在 2D 平面中找到一个点，该点的包含使新 voronoi 区域的面积最大化

Question

我正在使用仅限于平面的一部分的 Voronoi 图（以便我可以测量区域的面积，并且它们不会无穷大），并且在本练习中我只考虑整数坐标。

问题来了：
我有一个现有的 Voronoi 曲面细分，我想添加一个新点 (xi, yi)，以便由该点创建的新区域具有最大可能的面积。

现在，我想出了一个强力解决方案，它将针对每个可能的点运行（仅考虑整数），并计算新区域的面积。

我用 python 编写了代码，这是我测试过的一个示例：
我有 5 个顶点：

[[-5,  7], [ 2,  5], [-6, -2], [ 7,  2], [-4,  6]]

，它们的 Voronoi 细分为（红色数字表示该区域的面积）：

我尝试对 X 轴范围 [-5,5) 和 Y 轴范围 [-5,5) 范围内的顶点进行强力方法，并计算出使用顶点

(1, -5)

我们得到最大可能面积：

47.81

我尝试大量搜索如何进一步优化我的方法，但找不到任何方向。 有人可以指导我解决这个问题的数学和编程策略吗？任何见解将不胜感激。

Answer 1

我建议不要考虑每个可能的允许点作为新点的候选点并精确计算每个点的面积（这将是天真的蛮力方法，对于任何合理的“大”平面/网格来说绝对是不可行的），您可以简单地选择您认为运行时间和结果的准确性/最优性之间的权衡良好的粒度级别。

以下代码允许控制要考虑的新点的候选者的稀疏性（可通过参数

alpha

控制，这意味着仅考虑可被

alpha

整除的索引处的像素，例如

alpha=1

表示考虑每个像素，

 alpha=2

表示仅考虑四分之一的像素等）以及区域计算的粒度（可通过参数

beta

控制，意味着仅使用可被

beta

整除的索引处的像素来分配给簇以达到目的面积计算，例如

beta=1

会给出完整的精度）。 IE。值越高，结果越差，但速度越快。该算法只是选择最高区域中的候选点结果（近似为分配给它的像素份额）。

from random import randint, seed
import time
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import ConvexHull


def get_points(min_d, max_d, min_np, max_np):
    h = randint(min_d, max_d)
    w = randint(min_d, max_d)
    nd = randint(min_np, max_np)
    points = []
    while len(points) < nd:
        i = randint(0, h)
        j = randint(0, w)
        if (i, j) not in points:
            points.append((i, j))
    return h, w, points


def best_new_point(h, w, points, alpha, beta):
    cands = {(i, j) for i in range(0, h, alpha) for j in range(0, w, alpha)} - set(points)
    pixs = {(i, j) for i in range(0, h, beta) for j in range(0, w, beta)}
    
    best = None
    area = -1
    for cand in cands:
        areas = {k: set() for k in points}
        areas[cand] = set()
        for ij in pixs:
            f = lambda xy: (xy[0] - ij[0]) ** 2 + (xy[1] - ij[1]) ** 2
            k = min(areas.keys(), key=f)
            areas[k].add(ij)
        if len(areas[cand]) > area:
            best = cand
            area = len(areas[cand])
    return best

def draw_voronoi(h, w, points):
    img = np.zeros((h, w))
    areas = {k: set() for k in points}
    inds = {(i, j) for i in range(h) for j in range(w)}
    for ij in inds:
        f = lambda xy: (xy[0] - ij[0]) ** 2 + (xy[1] - ij[1]) ** 2
        k = min(points, key=f)
        areas[k].add(ij)
    for c, p in enumerate(points):
        area = areas[p]
        for i, j in area:
            img[i, j] = c
    return img

为了演示，忽略了“整数坐标”的附加约束（因为适应应该是直接的，这不是我首先想要关注的），以及少量的点以及小型飞机也被考虑。

seed(69)
dim = 128
nd = 4
h, w, points = get_points(dim, dim, nd, nd)
img = draw_voronoi(h, w, points)
for i, j in points:
    img[i, j] = nd+1
plt.imshow(img);
plt.scatter([y for x, y in points], [x for x, y in points], **{'c': 'black', 's': 12});

plt.subplots_adjust(left=0, right=1, top=1, bottom=0, wspace=0, hspace=0)
plt.axis('off')
plt.show()

原始 voronoi 图如下所示：

为了可视化此方法的两个“自由度”的效果，我选择可视化

alpha

（a，随列变化）和

beta

（b，随行变化）的 4 个不同值的每种组合（ 32、16、8 和 4)。红点显示考虑的点（alpha 的影响），红线显示用于面积计算的点的边界（beta 的影响）。值 p 是属于新点的像素的百分比，值 t 是计算新点所需的时间（以毫秒为单位）。

fig, ax = plt.subplots(nrows=4, ncols=4, figsize=(8, 8))
for i, a in enumerate([32, 16, 8, 4]):
    for j, b in enumerate([32, 16, 8, 4]):
        start = time.time()
        xy = best_new_point(h, w, points, a, b)
        end = time.time()
        t = end - start
        points2 = points + [xy]
        img2 = draw_voronoi(h, w, points2)
        area = np.count_nonzero(img2 == len(points)) / (h * w)
        ax[j, i].imshow(img2);
        ax[j, i].scatter([y for x, y in points2], [x for x, y in points2], **{'c': 'black', 's': 12});
        grid_a = [(i, j) for i in range(0, h, a) for j in range(0, w, a)]
        ax[j, i].scatter([y for x, y in grid_a], [x for x, y in grid_a], c='red', s=0.25);
        grid_b = [(i, j) for i in range(0, h, b) for j in range(0, w, b)]
        pts = [(i, j) for i, j in grid_b if img2[i, j] == nd]
        hull = ConvexHull(pts)
        pts_arr = np.array(pts)
        for simplex in hull.simplices:
            ax[j, i].plot(pts_arr[simplex, 1], pts_arr[simplex, 0], 'r-')
        
        ax[j, i].set_xticks([])
        ax[j, i].set_yticks([])
        ax[j, i].set_title(f'a: {a}, b: {b}\np: {area*100:.2f}% t: {t*1000:.0f}ms')
plt.tight_layout()
plt.show()