Skip to content

Instantly share code, notes, and snippets.

Last active June 9, 2022 22:16
Show Gist options
  • Save letmaik/8803860 to your computer and use it in GitHub Desktop.
Save letmaik/8803860 to your computer and use it in GitHub Desktop.
Creates a Voronoi diagram with cell polygons using scipy's Delaunay triangulation (scipy >= 0.9)
from __future__ import division
import collections
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from scipy.spatial import Delaunay, KDTree
# an adaptation of
# using ideas from
def voronoi(P):
Returns a list of all edges of the voronoi diagram for the given input points.
delauny = Delaunay(P)
triangles = delauny.points[delauny.vertices]
circum_centers = np.array([triangle_csc(tri) for tri in triangles])
long_lines_endpoints = []
lineIndices = []
for i, triangle in enumerate(triangles):
circum_center = circum_centers[i]
for j, neighbor in enumerate(delauny.neighbors[i]):
if neighbor != -1:
lineIndices.append((i, neighbor))
ps = triangle[(j+1)%3] - triangle[(j-1)%3]
ps = np.array((ps[1], -ps[0]))
middle = (triangle[(j+1)%3] + triangle[(j-1)%3]) * 0.5
di = middle - triangle[j]
ps /= np.linalg.norm(ps)
di /= np.linalg.norm(di)
if, ps) < 0.0:
ps *= -1000.0
ps *= 1000.0
long_lines_endpoints.append(circum_center + ps)
lineIndices.append((i, len(circum_centers) + len(long_lines_endpoints)-1))
vertices = np.vstack((circum_centers, long_lines_endpoints))
# filter out any duplicate lines
lineIndicesSorted = np.sort(lineIndices) # make (1,2) and (2,1) both (1,2)
lineIndicesTupled = [tuple(row) for row in lineIndicesSorted]
lineIndicesUnique = sorted(set(lineIndicesTupled))
return vertices, lineIndicesUnique
def triangle_csc(pts):
rows, cols = pts.shape
A = np.bmat([[2 *, pts.T), np.ones((rows, 1))],
[np.ones((1, rows)), np.zeros((1, 1))]])
b = np.hstack((np.sum(pts * pts, axis=1), np.ones((1))))
x = np.linalg.solve(A,b)
bary_coords = x[:-1]
return np.sum(pts * np.tile(bary_coords.reshape((pts.shape[0], 1)), (1, pts.shape[1])), axis=0)
def voronoi_cell_lines(points, vertices, lineIndices):
Returns a mapping from a voronoi cell to its edges.
:param points: shape (m,2)
:param vertices: shape (n,2)
:param lineIndices: shape (o,2)
:rtype: dict point index -> list of shape (n,2) with vertex indices
kd = KDTree(points)
cells = collections.defaultdict(list)
for i1,i2 in lineIndices:
v1,v2 = vertices[i1], vertices[i2]
mid = (v1+v2)/2
_, (p1Idx,p2Idx) = kd.query(mid, 2)
return cells
def voronoi_polygons(cells):
Transforms cell edges into polygons.
:param cells: as returned from voronoi_cell_lines
:rtype: dict point index -> list of vertex indices which form a polygon
# first, close the outer cells
for pIdx,lineIndices_ in cells.items():
dangling_lines = []
for i1,i2 in lineIndices_:
connections = list(filter(lambda i12_: (i1,i2) != (i12_[0],i12_[1]) and
(i1==i12_[0] or i1==i12_[1] or i2==i12_[0] or i2==i12_[1]),
assert 1 <= len(connections) <= 2
if len(connections) == 1:
assert len(dangling_lines) in [0,2]
if len(dangling_lines) == 2:
(i11,i12), (i21,i22) = dangling_lines
# determine which line ends are unconnected
connected = list(filter(lambda i12_: (i12_[0],i12_[1]) != (i11,i12) and (i12_[0] == i11 or i12_[1] == i11), lineIndices_))
i11Unconnected = len(connected) == 0
connected = list(filter(lambda i12_: (i12_[0],i12_[1]) != (i21,i22) and (i12_[0] == i21 or i12_[1] == i21), lineIndices_))
i21Unconnected = len(connected) == 0
startIdx = i11 if i11Unconnected else i12
endIdx = i21 if i21Unconnected else i22
cells[pIdx].append((startIdx, endIdx))
# then, form polygons by storing vertex indices in (counter-)clockwise order
polys = dict()
for pIdx,lineIndices_ in cells.items():
# get a directed graph which contains both directions and arbitrarily follow one of both
directedGraph = lineIndices_ + [(i2,i1) for (i1,i2) in lineIndices_]
directedGraphMap = collections.defaultdict(list)
for (i1,i2) in directedGraph:
orderedEdges = []
currentEdge = directedGraph[0]
while len(orderedEdges) < len(lineIndices_):
i1 = currentEdge[1]
i2 = directedGraphMap[i1][0] if directedGraphMap[i1][0] != currentEdge[0] else directedGraphMap[i1][1]
nextEdge = (i1, i2)
currentEdge = nextEdge
polys[pIdx] = [i1 for (i1,i2) in orderedEdges]
return polys
def polygons(points):
Returns the voronoi polygon for each input point.
:param points: shape (n,2)
:rtype: list of n polygons where each polygon is an array of vertices
vertices, lineIndices = voronoi(points)
cells = voronoi_cell_lines(points, vertices, lineIndices)
polys = voronoi_polygons(cells)
polylist = []
for i in range(len(points)):
poly = vertices[np.asarray(polys[i])]
return polylist
if __name__ == '__main__':
P = np.random.random((100,2))
fig = plt.figure(figsize=(4.5,4.5))
axes = plt.subplot(1,1,1)
polys = polygons(P)
for poly in polys:
p = matplotlib.patches.Polygon(poly, facecolor=np.random.rand(3,1))
X,Y = P[:,0],P[:,1]
plt.scatter(X, Y, marker='.', zorder=2)
Copy link

letmaik commented May 25, 2015

@BrianMacNamee This is because numpy's np.unique function changed slightly in behavior:

>>> np.__version__
>>> np.unique([(0, 187), (0, 20)])
array([[  0, 187],
       [  0,  20]])


>>> np.__version__
>>> np.unique([(0, 187), (0, 20)])
array([  0,  20, 187])

I guess you use numpy 1.9. I'll see if I can fix this.

EDIT: The code is fixed now and works with either version of numpy and also under Python 3.

Copy link

jwcarr commented Jul 22, 2015

Thanks! How easy would it be to get this working with Manhattan distance (as opposed to Euclidean)? Seems non-trivial.

Copy link

It should be noted that the points going into voronoi_polygons() have to be unique, otherwise the assertions will not hold.

Copy link

I have passed it a set of lat and long of shape (232, 2) and I get the following error:
File "", line 104, in voronoi_polygons
assert 1 <= len(connections) <= 2

Copy link

Hi, thank you for your great code ! I am getting the following error however when I try to run the code:

Traceback (most recent call last):
File "", line 169, in
p = matplotlib.patches.Polygon(poly, facecolor=np.random.rand(3,1))
File "G:\Anaconda3\lib\site-packages\matplotlib\", line 950, in init
Patch.init(self, **kwargs)
File "G:\Anaconda3\lib\site-packages\matplotlib\", line 126, in init
File "G:\Anaconda3\lib\site-packages\matplotlib\", line 334, in set_facecolor
File "G:\Anaconda3\lib\site-packages\matplotlib\", line 324, in _set_facecolor
self._facecolor = colors.to_rgba(color, alpha)
File "G:\Anaconda3\lib\site-packages\matplotlib\", line 143, in to_rgba
rgba = _to_rgba_no_colorcycle(c, alpha)
File "G:\Anaconda3\lib\site-packages\matplotlib\", line 194, in _to_rgba_no_colorcycle
raise ValueError("Invalid RGBA argument: {!r}".format(orig_c))
ValueError: Invalid RGBA argument: array([[ 0.16610242],
[ 0.15260443],
[ 0.91265765]])

I solved the error by replacing p = matplotlib.patches.Polygon(poly, facecolor=np.random.rand(3,1))
with p = matplotlib.patches.Polygon(poly, facecolor=(np.random.random((3,)))).And my matplotlib version is 2.0.2.

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