-
-
Save jeffmylife/ae31b7ff921c2a0b7285190b2df1a9b2 to your computer and use it in GitHub Desktop.
Colorized Voronoi diagram with Scipy, in 2D, including infinite regions
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 | |
try: | |
ridges = all_ridges[p1] | |
except KeyError: | |
raise KeyError("p1 not found in all_ridges. This usually means that there were duplicates in your input points. Try dropping duplicates and NaN's from the data you put into Voronoi().") | |
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) | |
plt.savefig('voro.png') | |
plt.show() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
I updated the print statements to be friendly to Python 3 and added a
try except
for letting the user know that having duplicates in the scipy.spatial.Voronoi might affect this code.