Skip to content

Instantly share code, notes, and snippets.

@synapticarbors
Created July 14, 2011 21:00

Revisions

  1. synapticarbors created this gist Jul 14, 2011.
    65 changes: 65 additions & 0 deletions voronoi.py
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,65 @@
    #!/usr/bin/env python
    # -----------------------------------------------------------------------------
    # Voronoi diagram from a list of points
    # Copyright (C) 2011 Nicolas P. Rougier
    #
    # Distributed under the terms of the BSD License.
    # -----------------------------------------------------------------------------
    import numpy as np
    import matplotlib
    import matplotlib.pyplot as plt

    def circumcircle(P1,P2,P3):
    '''
    Adapted from:
    http://local.wasp.uwa.edu.au/~pbourke/geometry/circlefrom3/Circle.cpp
    '''
    delta_a = P2 - P1
    delta_b = P3 - P2
    if np.abs(delta_a[0]) <= 0.000000001 and np.abs(delta_b[1]) <= 0.000000001:
    center_x = 0.5*(P2[0] + P3[0])
    center_y = 0.5*(P1[1] + P2[1])
    else:
    aSlope = delta_a[1]/delta_a[0]
    bSlope = delta_b[1]/delta_b[0]
    if np.abs(aSlope-bSlope) <= 0.000000001:
    return None
    center_x= (aSlope*bSlope*(P1[1] - P3[1]) + bSlope*(P1[0] + P2 [0]) \
    - aSlope*(P2[0]+P3[0]) )/(2* (bSlope-aSlope) )
    center_y = -1*(center_x - (P1[0]+P2[0])/2)/aSlope + (P1[1]+P2[1])/2;
    return center_x, center_y

    def voronoi(X,Y):
    P = np.zeros((X.size+4,2))
    P[:X.size,0], P[:Y.size,1] = X, Y
    # We add four points at "infinity"
    m = max(np.abs(X).max(), np.abs(Y).max())*1e5
    P[X.size:,0] = -m, -m, +m, +m
    P[Y.size:,1] = -m, +m, -m, +m
    D = matplotlib.tri.Triangulation(P[:,0],P[:,1])
    T = D.triangles
    n = T.shape[0]
    C = np.zeros((n,2))
    for i in range(n):
    C[i] = circumcircle(P[T[i,0]],P[T[i,1]],P[T[i,2]])
    X,Y = C[:,0], C[:,1]
    segments = []
    for i in range(n):
    for j in range(3):
    k = D.neighbors[i][j]
    if k != -1:
    segments.append( [(X[i],Y[i]), (X[k],Y[k])] )
    return segments

    if __name__ == '__main__':

    X = np.random.random(200)
    Y = np.random.random(200)
    fig = plt.figure(figsize=(10,10))
    axes = plt.subplot(1,1,1)
    plt.scatter(X,Y)
    segments = voronoi(X,Y)
    lines = matplotlib.collections.LineCollection(segments, color='0.75')
    axes.add_collection(lines)
    plt.axis([0,1,0,1])
    plt.show()
    122 changes: 122 additions & 0 deletions voronoi_subdivide.py
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,122 @@
    import numpy as np
    from scipy.spatial.distance import cdist,pdist,squareform
    import matplotlib
    import matplotlib.pyplot as plt
    import matplotlib.tri as tri

    from voronoi import voronoi # http://webloria.loria.fr/~rougier/coding/neural-networks/voronoi.py


    # Code to take a voronoi tesselation of a 2D space and subdivide each voronoi cell by
    # selecting the m furthest neighboring voronoi cell centers and generating n-1 new points between
    # the center and the midpoint to the neighbor, that along with the original center will
    # define a new voronoi parition within the cell.


    k = 10 # number of reference poses
    N = 40000 # number of test points
    m = 2 # number of neighbors to base subdivision off of
    n = 3 # number of slices

    colors = np.arange(m*n+1)

    g = np.linspace(0,1,n+1)
    gdata = np.array([0,1])

    # Define Random reference points
    x_i = np.random.uniform(0,5,size=(k,2))

    # Calculate pairwise distance between references points
    d_ij = squareform(pdist(x_i))

    # Define voronoi
    vorbounds = voronoi(x_i[:,0],x_i[:,1])

    # Create Deluanay Triangulation
    triang = tri.Triangulation(x_i[:,0], x_i[:,1])
    edges = triang.edges

    # Create supplementary poses to subdivide voronoi cells
    y_i = []
    for vori,vor in enumerate(x_i):
    temppoints = []
    # Identify neighboring reference poses
    jj = np.where((edges[:,0] == vori) | (edges[:,1] == vori))
    pj = edges[jj].flatten()
    pj = pj[pj != vori]

    # Get m poses with largest distance
    d_j = d_ij[vori,:]
    neigh_target_j = pj[np.argsort(d_j[pj])[::-1]]

    # Linearly interpolate n points between reference and midpoint with neighbors
    count = 0
    for nj in neigh_target_j:
    if count > m - 1:
    break
    midpoint = 0.5*(vor - x_i[nj,:]) + x_i[nj,:]
    nearest = np.argmin(cdist(np.atleast_2d(midpoint),x_i))
    # exclude centeres connected via the Delaunay triangulation
    # that pass through an intervening voronoi cell
    if not nearest in [nj,vori]:
    continue
    endpoints = np.vstack((vor,midpoint))
    a = np.zeros((g.shape[0]-2,2))
    for dim in xrange(2):
    a[:,dim] = np.interp(g,gdata,endpoints[:,dim])[1:-1]

    temppoints.append(a)
    count += 1

    if len(temppoints):
    y_i.append(np.vstack(temppoints))
    else:
    y_i.append([])

    # Generate a large number of random configs and assign to subdivisions
    x = np.random.uniform(0,5,size=(N,2))
    cassign = np.zeros((N,))
    for ii,p in enumerate(x):
    # Calculate the distance from point to every reference pose
    p = np.atleast_2d(p)
    d_j = cdist(p,x_i,metric='euclidean')

    # Identify voronoi membership
    pi = np.argmin(d_j)

    # Determine membership in subdivisions
    if len(y_i[pi]):
    subpoints = np.vstack((x_i[pi,:],y_i[pi]))
    else:
    subpoints = x_i[pi,:]
    subpoints = np.atleast_2d(subpoints)
    si = np.argmin(cdist(p,subpoints))
    cassign[ii] = colors[si]

    #print cassign
    fig = plt.figure()
    ax = fig.add_subplot(111)

    # Plot reference poses
    ax.plot(x_i[:,0],x_i[:,1],'.',color='yellow',markersize=20)

    # Plot internal points
    y_ii = [a for a in y_i if a != []]
    z = np.vstack(y_ii)
    ax.plot(z[:,0],z[:,1],'.',color='yellow',markersize=10)


    # Plot voronoi edges
    vorlines = matplotlib.collections.LineCollection(vorbounds, color='yellow',linewidths=3)
    ax.add_collection(vorlines)

    # Plot Deluanay triangulation
    ax.triplot(triang, 'wo-')

    # Plot test particles
    ax.scatter(x[:,0],x[:,1],marker='o',c=cassign,s=5)

    plt.axis([-1,6,-1,6])
    plt.show()