import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Voronoi
def voronoi_finite_polygons_2d(vor, radius=None):
"""
Reconstruct infinite voronoi regions in a 2D diagram to finite
regions.
Parameters
----------
vor : Voronoi
Input diagram
radius : float, optional
Distance to 'points at infinity'.
Returns
-------
regions : list of tuples
Indices of vertices in each revised Voronoi regions.
vertices : list of tuples
Coordinates for revised Voronoi vertices. Same as coordinates
of input vertices, with 'points at infinity' appended to the
end.
"""
if vor.points.shape[1] != 2:
raise ValueError("Requires 2D input")
new_regions = []
new_vertices = vor.vertices.tolist()
center = vor.points.mean(axis=0)
if radius is None:
radius = vor.points.ptp().max()*2
# Construct a map containing all ridges for a given point
all_ridges = {}
for (p1, p2), (v1, v2) in zip(vor.ridge_points, vor.ridge_vertices):
all_ridges.setdefault(p1, []).append((p2, v1, v2))
all_ridges.setdefault(p2, []).append((p1, v1, v2))
# Reconstruct infinite regions
for p1, region in enumerate(vor.point_region):
vertices = vor.regions[region]
if all([v >= 0 for v in vertices]):
# finite region
new_regions.append(vertices)
continue
# reconstruct a non-finite region
ridges = all_ridges[p1]
new_region = [v for v in vertices if v >= 0]
for p2, v1, v2 in ridges:
if v2 < 0:
v1, v2 = v2, v1
if v1 >= 0:
# finite ridge: already in the region
continue
# Compute the missing endpoint of an infinite ridge
t = vor.points[p2] - vor.points[p1] # tangent
t /= np.linalg.norm(t)
n = np.array([-t[1], t[0]]) # normal
midpoint = vor.points[[p1, p2]].mean(axis=0)
direction = np.sign(np.dot(midpoint - center, n)) * n
far_point = vor.vertices[v2] + direction * radius
new_region.append(len(new_vertices))
new_vertices.append(far_point.tolist())
# sort region counterclockwise
vs = np.asarray([new_vertices[v] for v in new_region])
c = vs.mean(axis=0)
angles = np.arctan2(vs[:,1] - c[1], vs[:,0] - c[0])
new_region = np.array(new_region)[np.argsort(angles)]
# finish
new_regions.append(new_region.tolist())
return new_regions, np.asarray(new_vertices)
# make up data points
np.random.seed(1234)
points = np.random.rand(15, 2)
# compute Voronoi tesselation
vor = Voronoi(points)
# plot
regions, vertices = voronoi_finite_polygons_2d(vor)
print "--"
print regions
print "--"
print vertices
# colorize
for region in regions:
polygon = vertices[region]
plt.fill(*zip(*polygon), alpha=0.4)
plt.plot(points[:,0], points[:,1], 'ko')
plt.axis('equal')
plt.xlim(vor.min_bound[0] - 0.1, vor.max_bound[0] + 0.1)
plt.ylim(vor.min_bound[1] - 0.1, vor.max_bound[1] + 0.1)
-- [[10, 15, 14, 21, 20], [19, 7, 12, 17], [22, 23, 4, 3, 0], [24, 14, 16, 18, 25], [5, 26, 27], [9, 1, 0, 3, 2, 8], [12, 7, 6, 4, 3, 2, 11], [13, 8, 2, 11], [28, 29, 0, 1], [18, 19, 7, 6, 5, 30, 31], [17, 12, 11, 13, 15, 14, 16], [32, 1, 9, 10, 33], [34, 35, 5, 6, 4], [19, 17, 16, 18], [15, 10, 9, 8, 13]] -- [[ 0.51204146 0.28170809] [ 0.31568864 0.2231658 ] [ 0.60181497 0.48198612] [ 0.61195783 0.46638446] [ 0.76349576 0.49961565] [ 0.86163554 0.77300743] [ 0.82428491 0.74711566] [ 0.5955872 0.86737793] [ 0.3398531 0.5360898 ] [ 0.20717026 0.45505837] [ 0.19839983 0.46568522] [ 0.54203769 0.60556533] [ 0.52991347 0.64525971] [ 0.34632775 0.58619376] [ 0.28095359 0.68979933] [ 0.27956806 0.65401536] [ 0.32776593 0.71199644] [ 0.40994669 0.69667456] [ 0.41517255 1.51501742] [ 0.46233568 1.32637502] [-1.51846216 1.25291487] [-1.42673123 1.49674306] [ 1.80180548 -1.09809421] [ 2.5145059 -0.20841621] [-1.42673123 1.49674306] [ 0.11182282 3.37923966] [ 2.73878052 0.56402822] [ 0.92996562 2.66051282] [-0.89071967 -1.23008044] [ 1.80180548 -1.09809421] [ 0.92996562 2.66051282] [ 0.11182282 3.37923966] [-0.89071967 -1.23008044] [-1.51846216 1.25291487] [ 2.5145059 -0.20841621] [ 2.73878052 0.56402822]]
(-0.086231550409317764, 0.98264119063611655)