Search of optimal algorithm to find a point in 2D plane, whose inclusion maximizes the area of new voronoi region

Question

I am working with Voronoi diagrams limited to a section of plane (so that I can measure areas of regions, and they don't go off to infinity), and I'm considering only integer co-ordintes in this exercise.

Here's the problem:
I have an existing Voronoi tessellation and I want to add a new point (xi, yi) such that the new region created by this point has the maximum possible area.

Now, what I came up with was a brute force solution, which will run for each possible point (considering integers only), and calculate area of new region.

I wrote the code in python and here is one example I tested:
I have 5 vertices: [[-5, 7], [ 2, 5], [-6, -2], [ 7, 2], [-4, 6]] and their Voronoi tessellation is (numbers in red denote area of that region):

I tried brute force approach for vertices in range [-5,5) along X-axis and [-5,5) along Y-axis, and calculated that with vertex (1, -5) we get maximum possible area: 47.81

I tried to search a lot on how to optimize my approach further from brute force, but couldn't find any direction. Could someone kindly guide me on the mathematical and programming strategies to address this problem? Any insights would be greatly appreciated.

Answer 1

I suggest to - instead of considering each possible allowed point as a candidate for the new point and calculating the area for each of those precisely (which would be the naive brute-force approach, most definitely infeasible for any reasonably "large" planes/grids) - simply choose granularity levels for which you deem the trade-off between runtime and accuracy/optimality of results to be good.

The following code allows to control both the sparsity of the candidates to consider for the new point (controllable via the parameter alpha, meaning considering only pixels at indices divisible by alpha, e.g. alpha=1 means every pixel is considered, alpha=2 means only a quarter of the pixels are considered, etc.) as well as the granularity of the area calculation (controllable via the parameter beta, meaning only using pixels at indices divisible by beta to assign to clusters for the purpose of area calculation, e.g. beta=1 would give full precision). I.e. higher values give worse but faster results. The algorithm simply picks whichever candidate point results in the highest area (approximated as the share of pixels assigned to it).

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

For the sake of demonstration, the additional constraint of "integer coordinates" is ignored (as the adaptation should be straight-forward & that wasn't what I wanted to focus on here), and only a small number of points as well as only a small plane is considered.

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

The original voronoi diagram looks as follows:

To visualize the effects of the two "degrees of freedom" for this approach, I chose to visualize each combination of 4 different values for alpha (a, varying across columns) and beta (b, varying across rows) (32, 16, 8 and 4). The red dots show the points considered (effect of alpha), and the red line shows the boundary of the points used for the area calculation (effect of beta). The value p is percentage of the pixels that belong to the new point, and the value t is the time it took to compute the new point, in miliseconds. As you can see, the percentage is overall increasing as alpha and beta are decreased, as expected.

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