# [ACCEPTED]-Voronoi - Compute exact boundaries of every region-voronoi

Accepted answer
Score: 14

### 1. Algorithm

I suggest the following two-step approach:

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

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

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

### 2. Example

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

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

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

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

### 3. Code

The 8 first step is to generate polygons corresponding 7 to the Voronoi regions. Like Ophir Cami, I 6 derived this from the implementation of 5 `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 4 polygons with the bounding polygon. We also 3 need to choose an appropriate diameter to 2 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 1 in §2 above.

Score: 4

I took the `voronoi_plot_2d` and modified it. see below.

``````import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Voronoi
from shapely.geometry import Polygon, Point

# Voronoi - Compute exact boundaries of every region

def angle_between(v0, v1):
return np.math.atan2(np.linalg.det([v0, v1]), np.dot(v0, v1))

def calc_angle(c0, c1, c2):
return angle_between(np.array(c1) - np.array(c0), np.array(c2) - np.array(c1))

def is_convex(polygon):
temp_coords = np.array(polygon.exterior.coords)
temp_coords = np.vstack([temp_coords, temp_coords[1, :]])

for i, (c0, c1, c2) in enumerate(zip(temp_coords, temp_coords[1:], temp_coords[2:])):
if i == 0:
first_angle_crit = calc_angle(c0, c1, c2) > 0
elif (calc_angle(c0, c1, c2) > 0) != first_angle_crit:
return False
return True

def infinite_segments(vor_):
line_segments = []
center = vor_.points.mean(axis=0)
for pointidx, simplex in zip(vor_.ridge_points, vor_.ridge_vertices):
simplex = np.asarray(simplex)
if np.any(simplex < 0):
i = simplex[simplex >= 0][0]  # finite end Voronoi vertex

t = vor_.points[pointidx[1]] - vor_.points[pointidx[0]]  # tangent
t /= np.linalg.norm(t)
n = np.array([-t[1], t[0]])  # normal

midpoint = vor_.points[pointidx].mean(axis=0)
direction = np.sign(np.dot(midpoint - center, n)) * n

line_segments.append([(vor_.vertices[i, 0], vor_.vertices[i, 1]),
(direction[0], direction[1])])
return line_segments

class NotConvexException(Exception):
def __str__(self):
return 'The Polygon is not Convex!!!'

class NotAllPointsAreInException(Exception):
def __str__(self):
return 'Not all points are in the polygon!!!'

def intersect(p0, u, q0, q1):
v = (q1 - q0)[np.newaxis].T
A = np.hstack([u, -v])
b = q0 - p0
try:
inv_A = np.linalg.inv(A)
except np.linalg.LinAlgError:
return np.nan, np.nan
return np.dot(inv_A, b)

def _adjust_bounds(ax__, points_):
ptp_bound = points_.ptp(axis=0)
ax__.set_xlim(points_[:, 0].min() - 0.1*ptp_bound[0], points_[:, 0].max() + 0.1*ptp_bound[0])
ax__.set_ylim(points_[:, 1].min() - 0.1*ptp_bound[1], points_[:, 1].max() + 0.1*ptp_bound[1])

def in_polygon(polygon, points_):
return [polygon.contains(Point(x)) for x in points_]

def voronoi_plot_2d_inside_convex_polygon(vor_, polygon, ax__=None, **kw):
from matplotlib.collections import LineCollection

if not all(in_polygon(polygon, vor_.points_)):
raise NotAllPointsAreInException()

if not is_convex(polygon):
raise NotConvexException()

if vor_.points.shape[1] != 2:
raise ValueError("Voronoi diagram is not 2-D")

vor_inside_ind = np.array([i for i, x in enumerate(vor_.vertices) if polygon.contains(Point(x))])
vor_outside_ind = np.array([i for i, x in enumerate(vor_.vertices) if not polygon.contains(Point(x))])
ax__.plot(vor_.points[:, 0], vor_.points[:, 1], '.')
if kw.get('show_vertices', True):
ax__.plot(vor_.vertices[vor_inside_ind, 0], vor_.vertices[vor_inside_ind, 1], 'o')

temp_coords = np.array(polygon.exterior.coords)
line_segments = []
for t0, t1 in zip(temp_coords, temp_coords[1:]):
line_segments.append([t0, t1])
ax__.add_collection(LineCollection(line_segments, colors='b', linestyle='solid'))
line_segments = []
for simplex in vor_.ridge_vertices:
simplex = np.asarray(simplex)
if np.all(simplex >= 0):
if not all(in_polygon(polygon, vor_.vertices[simplex])):
continue
line_segments.append([(x, y) for x, y in vor_.vertices[simplex]])

ax__.add_collection(LineCollection(line_segments, colors='k', linestyle='solid'))

line_segments = infinite_segments(vor_)
from_inside = np.array([x for x in line_segments if polygon.contains(Point(x[0]))])

line_segments = []

for f in from_inside:
for coord0, coord1 in zip(temp_coords, temp_coords[1:]):
s, t = intersect(f[0], f[1][np.newaxis].T, coord0, coord1)
if 0 < t < 1 and s > 0:
line_segments.append([f[0], f[0] + s * f[1]])
break

ax__.add_collection(LineCollection(np.array(line_segments), colors='k', linestyle='dashed'))

line_segments = []

for v_o_ind in vor_outside_ind:
for simplex in vor_.ridge_vertices:
simplex = np.asarray(simplex)
if np.any(simplex < 0):
continue
if np.any(simplex == v_o_ind):
i = simplex[simplex != v_o_ind][0]
for coord0, coord1 in zip(temp_coords, temp_coords[1:]):
s, t = intersect(
vor_.vertices[i],
(vor_.vertices[v_o_ind] - vor_.vertices[i])[np.newaxis].T,
coord0,
coord1
)
if 0 < t < 1 and 0 < s < 1:
line_segments.append([
vor_.vertices[i],
vor_.vertices[i] + s * (vor_.vertices[v_o_ind] - vor_.vertices[i])
])
break

ax__.add_collection(LineCollection(np.array(line_segments), colors='r', linestyle='dashed'))

_adjust_bounds(ax__, temp_coords)

return ax__.figure

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

global_boundaries = Polygon([[-5, -2], [3.4, -2], [4.7, 4], [2.7, 5.7], [-1, 4]])
fig = plt.figure()
ax = fig.add_subplot(111)

vor = Voronoi(points)
voronoi_plot_2d_inside_convex_polygon(vor, global_boundaries, ax_=ax)
plt.show()
``````

Note: There 1 are two simple constraints:

1. The polygon must be convex.
2. All the points must be inside the polygon.

colors:

• Original points are in blue.
• Voronoi vertices are in green.
• Finite inner Voronoi ridges are in solid black.
• Finite outer Voronoi ridges are in dashed red.
• Infinite Voronoi ridges are in dashed black.

More Related questions