Voronoi - Compute exact boundaries of every region

computational-geometry, numpy, python, voronoi

Solution

1. Algorithm

I suggest the following two-step approach:

First, make a convex polygon for each of the Voronoi regions. In the case of the infinite regions, do this by splitting the point at infinity into two points that are far enough away, joined by an edge. ("Far enough away" means that the extra edge passes entirely outside the bounding polygon.)

Intersect each of the polygons from step (1) with the bounding polygon, using shapely's `intersection` method.

The benefit of this approach over Ophir Cami's answer is that it works with non-convex bounding polygons, and the code is a little simpler.

2. Example

Let's start with the Voronoi diagram for the points from Ophir Cami's answer. The infinite ridges are shown as dashed lines by `scipy.spatial.voronoi_plot_2d`:

Then for each Voronoi region we construct a convex polygon. This is easy for the finite regions, but we have to zoom out a long way to see what happens to the infinite Voronoi regions. The polygons corresponding to these regions have an extra edge that is far enough away that it lies entirely outside the bounding polygon:

Now we can intersect the polygons for each Voronoi region with the bounding polygon:

In this case, all the Voronoi polygons have non-empty intersection with the bounding polygon, but in the general case some of them might vanish.

3. Code

The first step is to generate polygons corresponding to the Voronoi regions. Like Ophir Cami, I derived this from the implementation of `scipy.spatial.voronoi_plot_2d`.

from collections import defaultdict

from shapely.geometry import Polygon

def voronoi_polygons(voronoi, diameter):
    """Generate shapely.geometry.Polygon objects corresponding to the
    regions of a scipy.spatial.Voronoi object, in the order of the
    input points. The polygons for the infinite regions are large
    enough that all points within a distance 'diameter' of a Voronoi
    vertex are contained in one of the infinite polygons.

    """
    centroid = voronoi.points.mean(axis=0)

    # Mapping from (input point index, Voronoi point index) to list of
    # unit vectors in the directions of the infinite ridges starting
    # at the Voronoi point and neighbouring the input point.
    ridge_direction = defaultdict(list)
    for (p, q), rv in zip(voronoi.ridge_points, voronoi.ridge_vertices):
        u, v = sorted(rv)
        if u == -1:
            # Infinite ridge starting at ridge point with index v,
            # equidistant from input points with indexes p and q.
            t = voronoi.points[q] - voronoi.points[p] # tangent
            n = np.array([-t[1], t[0]]) / np.linalg.norm(t) # normal
            midpoint = voronoi.points[[p, q]].mean(axis=0)
            direction = np.sign(np.dot(midpoint - centroid, n)) * n
            ridge_direction[p, v].append(direction)
            ridge_direction[q, v].append(direction)

    for i, r in enumerate(voronoi.point_region):
        region = voronoi.regions[r]
        if -1 not in region:
            # Finite region.
            yield Polygon(voronoi.vertices[region])
            continue
        # Infinite region.
        inf = region.index(-1)              # Index of vertex at infinity.
        j = region[(inf - 1) % len(region)] # Index of previous vertex.
        k = region[(inf + 1) % len(region)] # Index of next vertex.
        if j == k:
            # Region has one Voronoi vertex with two ridges.
            dir_j, dir_k = ridge_direction[i, j]
        else:
            # Region has two Voronoi vertices, each with one ridge.
            dir_j, = ridge_direction[i, j]
            dir_k, = ridge_direction[i, k]

        # Length of ridges needed for the extra edge to lie at least
        # 'diameter' away from all Voronoi vertices.
        length = 2 * diameter / np.linalg.norm(dir_j + dir_k)

        # Polygon consists of finite part plus an extra edge.
        finite_part = voronoi.vertices[region[inf + 1:] + region[:inf]]
        extra_edge = [voronoi.vertices[j] + dir_j * length,
                      voronoi.vertices[k] + dir_k * length]
        yield Polygon(np.concatenate((finite_part, extra_edge)))

The second step is to intersect the Voronoi polygons with the bounding polygon. We also need to choose an appropriate diameter to pass to `voronoi_polygons`.

import matplotlib.pyplot as plt
from scipy.spatial import Voronoi

points = np.array([[0.1, -0.4], [0, 1.5], [0, 2.25], [1, 0], [1, 1], [1, 2],
                   [2, 0], [2.5, 1], [2, 2], [2.3, 2.3], [-0.5, -1.3], [-1.5, 3]])
boundary = np.array([[-5, -2], [3.4, -2], [4.7, 4], [2.7, 5.7], [-1, 4]])

x, y = boundary.T
plt.xlim(round(x.min() - 1), round(x.max() + 1))
plt.ylim(round(y.min() - 1), round(y.max() + 1))
plt.plot(*points.T, 'b.')

diameter = np.linalg.norm(boundary.ptp(axis=0))
boundary_polygon = Polygon(boundary)
for p in voronoi_polygons(Voronoi(points), diameter):
    x, y = zip(*p.intersection(boundary_polygon).exterior.coords)
    plt.plot(x, y, 'r-')

plt.show()

This plots the last of the figures in §2 above.

Problem

Original source

Related problems

• Colorize Voronoi Diagram