Creates a Voronoi diagram with cell polygons using scipy's Delaunay triangulation (scipy >= 0.9)

voronoi_polygons.py

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()

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.

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

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

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

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

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.