Instantly share code, notes, and snippets.

# letmaik/voronoi_polygons.py

Last active June 9, 2022 22:16
Show Gist options
• 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)
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
 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 https://stackoverflow.com/a/15783581/60982 # using ideas from https://stackoverflow.com/a/9471601/60982 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)) else: 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 np.dot(di, ps) < 0.0: ps *= -1000.0 else: 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 * np.dot(pts, 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) cells[p1Idx].append((i1,i2)) cells[p2Idx].append((i1,i2)) 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]), lineIndices_)) assert 1 <= len(connections) <= 2 if len(connections) == 1: dangling_lines.append((i1,i2)) 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: directedGraphMap[i1].append(i2) 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) orderedEdges.append(nextEdge) 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])] polylist.append(poly) 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) plt.axis([-0.05,1.05,-0.05,1.05]) polys = polygons(P) for poly in polys: p = matplotlib.patches.Polygon(poly, facecolor=np.random.rand(3,1)) axes.add_patch(p) X,Y = P[:,0],P[:,1] plt.scatter(X, Y, marker='.', zorder=2) plt.axis([-0.05,1.05,-0.05,1.05]) plt.show()

### BrianMacNamee commented May 15, 2015

Hi there. The graphs you are generating rally look great and I'd like to generate some similar ones. I am getting the following error however when I try to run this:

TypeError Traceback (most recent call last)
in ()
161 plt.axis([-0.05,1.05,-0.05,1.05])
162
--> 163 polys = polygons(P)
164
165 for poly in polys:

in polygons(points)
145 '''
146 vertices, lineIndices = voronoi(points)
--> 147 cells = voronoi_cell_lines(points, vertices, lineIndices)
148 polys = voronoi_polygons(cells)
149 polylist = []

in voronoi_cell_lines(points, vertices, lineIndices)
75
76 cells = collections.defaultdict(list)
---> 77 for i1,i2 in lineIndices:
78 v1,v2 = vertices[i1], vertices[i2]
79 mid = (v1+v2)/2

TypeError: 'numpy.int64' object is not iterable

Any suggestions for what I might need to do to fix this?

Thanks.

### letmaik commented May 25, 2015

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

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

vs.

```>>> np.__version__
'1.9.2'
>>> 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.

### 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.

### lukaslueg commented Jul 24, 2015

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

### franciscovargas commented Aug 27, 2015

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

### JoJocoder commented Mar 9, 2018

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 "voronoi_polygons.py", line 169, in
p = matplotlib.patches.Polygon(poly, facecolor=np.random.rand(3,1))
File "G:\Anaconda3\lib\site-packages\matplotlib\patches.py", line 950, in init
Patch.init(self, **kwargs)
File "G:\Anaconda3\lib\site-packages\matplotlib\patches.py", line 126, in init
self.set_facecolor(facecolor)
File "G:\Anaconda3\lib\site-packages\matplotlib\patches.py", line 334, in set_facecolor
self._set_facecolor(color)
File "G:\Anaconda3\lib\site-packages\matplotlib\patches.py", line 324, in _set_facecolor
self._facecolor = colors.to_rgba(color, alpha)
File "G:\Anaconda3\lib\site-packages\matplotlib\colors.py", line 143, in to_rgba
rgba = _to_rgba_no_colorcycle(c, alpha)
File "G:\Anaconda3\lib\site-packages\matplotlib\colors.py", 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.