Skip to content

Instantly share code, notes, and snippets.

Forked from pv/
Created March 26, 2017 00:36
Show Gist options
  • Star 7 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
  • Save Sklavit/e05f0b61cb12ac781c93442fbea4fb55 to your computer and use it in GitHub Desktop.
Save Sklavit/e05f0b61cb12ac781c93442fbea4fb55 to your computer and use it in GitHub Desktop.
Colorized Voronoi diagram with Scipy, in 2D, including infinite regions which are clipped to given box.
# coding=utf-8
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Voronoi
from shapely.geometry import Polygon
def voronoi_finite_polygons_2d(vor, radius=None):
Reconstruct infinite voronoi regions in a 2D diagram to finite
vor : Voronoi
Input diagram
radius : float, optional
Distance to 'points at infinity'.
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
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
# 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
# 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( - center, n)) * n
far_point = vor.vertices[v2] + direction * radius
# 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
return new_regions, np.asarray(new_vertices)
# make up data points
points = np.random.rand(15, 2)
# compute Voronoi tesselation
vor = Voronoi(points)
# plot
regions, vertices = voronoi_finite_polygons_2d(vor)
min_x = vor.min_bound[0] - 0.1
max_x = vor.max_bound[0] + 0.1
min_y = vor.min_bound[1] - 0.1
max_y = vor.max_bound[1] + 0.1
mins = np.tile((min_x, min_y), (vertices.shape[0], 1))
bounded_vertices = np.max((vertices, mins), axis=0)
maxs = np.tile((max_x, max_y), (vertices.shape[0], 1))
bounded_vertices = np.min((bounded_vertices, maxs), axis=0)
box = Polygon([[min_x, min_y], [min_x, max_y], [max_x, max_y], [max_x, min_y]])
# colorize
for region in regions:
polygon = vertices[region]
# Clipping polygon
poly = Polygon(polygon)
poly = poly.intersection(box)
polygon = [p for p in poly.exterior.coords]
plt.fill(*zip(*polygon), alpha=0.4)
plt.plot(points[:, 0], points[:, 1], 'ko')
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)
Copy link

I might be wrong, but the line 53 can throw an error if p1 not in all_ridges. Adding a line if p1 in all_ridges: simply solves the problem I think

Copy link

Hi, In which list are the generated vertices stored? I managed to graph such points by modifying some attributes, but I would like to know the position (x, y) of each generated vertices.

Copy link

Sklavit commented Nov 29, 2020

@Javier1314 It is hard to recall and test now for me, but I guess that generated vertices (which are for infinity) is in the end of new_vertices.
So you should take new_vertices[len(vor.vertices):]. In this list tuples (or lists) are stored, so I suppose point[0] will be for x and point[1] will be for y.

Copy link

Vosatorp commented Mar 2, 2021

What is poly.exterior. coords? My code don't work because <<AttributeError: 'GeometryCollection' object has no attribute 'exterior'>>

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment