voronoi.py

import numpy as np
from scipy.spatial import Delaunay


#from http://pydec.googlecode.com/svn-history/r10/trunk/pydec/math/circumcenter.py

__all__ = ['is_wellcentered', 'circumcenter', 'circumcenter_barycentric']

from numpy import bmat, hstack, vstack, dot, sqrt, ones, zeros, sum, \
        asarray
from numpy.linalg import solve,norm

def is_wellcentered(pts, tol=1e-8):
    """Determines whether a set of points defines a well-centered simplex.
    """
    barycentric_coordinates = circumcenter_barycentric(pts)    
    return min(barycentric_coordinates) > tol

def circumcenter_barycentric(pts):
    """Barycentric coordinates of the circumcenter of a set of points.
    
    Parameters
    ----------
    pts : array-like
        An N-by-K array of points which define an (N-1)-simplex in K dimensional space.
        N and K must satisfy 1 <= N <= K + 1 and K >= 1.

    Returns
    -------
    coords : ndarray
        Barycentric coordinates of the circumcenter of the simplex defined by pts.
        Stored in an array with shape (K,)
        
    Examples
    --------
    >>> from pydec.math.circumcenter import *
    >>> circumcenter_barycentric([[0],[4]])           # edge in 1D
    array([ 0.5,  0.5])
    >>> circumcenter_barycentric([[0,0],[4,0]])       # edge in 2D
    array([ 0.5,  0.5])
    >>> circumcenter_barycentric([[0,0],[4,0],[0,4]]) # triangle in 2D
    array([ 0. ,  0.5,  0.5])

    See Also
    --------
    circumcenter_barycentric
   
    References
    ----------
    Uses an extension of the method described here:
    http://www.ics.uci.edu/~eppstein/junkyard/circumcenter.html

    """    

    pts = asarray(pts)

    rows,cols = pts.shape

    assert(rows <= cols + 1)    

    A = bmat( [[ 2*dot(pts,pts.T), ones((rows,1)) ],
               [  ones((1,rows)) ,  zeros((1,1))  ]] )

    b = hstack((sum(pts * pts, axis=1),ones((1))))
    x = solve(A,b)
    bary_coords = x[:-1]  

    return bary_coords
    
def circumcenter(pts):
    """Circumcenter and circumradius of a set of points.
    
    Parameters
    ----------
    pts : array-like
        An N-by-K array of points which define an (N-1)-simplex in K dimensional space.
        N and K must satisfy 1 <= N <= K + 1 and K >= 1.

    Returns
    -------
    center : ndarray
        Circumcenter of the simplex defined by pts.  Stored in an array with shape (K,)
    radius : float
        Circumradius of the circumsphere that circumscribes the points defined by pts.
        
    Examples
    --------
    >>> circumcenter([[0],[1]])             # edge in 1D
    (array([ 0.5]), 0.5)
    >>> circumcenter([[0,0],[1,0]])         # edge in 2D
    (array([ 0.5,  0. ]), 0.5)
    >>> circumcenter([[0,0],[1,0],[0,1]])   # triangle in 2D
    (array([ 0.5,  0.5]), 0.70710678118654757)

    See Also
    --------
    circumcenter_barycentric
   
    References
    ----------
    Uses an extension of the method described here:
    http://www.ics.uci.edu/~eppstein/junkyard/circumcenter.html

    """
    pts = asarray(pts)      
    bary_coords = circumcenter_barycentric(pts)
    center = dot(bary_coords,pts)
    radius = norm(pts[0,:] - center)
    return (center,radius)
    
#end

def voronoi(points):
    """
    points is a numpy array n-by-ndim where n is the number of points and ndim is the dimensionality of each point.
    so far this only works with ndim == 2

    based off: http://stackoverflow.com/questions/10650645/python-calculate-voronoi-tesselation-from-scipys-delaunay-triangulation-in-3d
    """
    assert len(points.shape)==2
    n,ndim = points.shape
    
    tri = Delaunay(points)

    #points in the simplex
    p = tri.points[tri.vertices]
    #p[0][0], p[0][1], p[0][2], .. are the vertices of the first simplex
    nsimplices = p.shape[0]
    assert p.shape[1] == ndim+1     #a simplex in ndim-dimensional space has ndim_1 vertices
    assert p.shape[2] == ndim

    #fnd vertices of graph
    ccs = np.array([circumcenter(p[i])[0] for i in range(nsimplices)])
    assert ccs.shape[1] == ndim    

    #tri.neighbors[i] consists of (n_dim+1) indices. these indices are into tri.vertices -- the ID of each simplex that neighbors simplex with ID i.

    #draw Vornoi edges from ccs[i] to vc[i,0] ... vc[i,n_dim+1]
    vc = ccs[tri.neighbors]
    
    vc[tri.neighbors == -1,:] = np.nan # edges at infinity, plotting those would need more work...

    lines = []
    for i in range(ndim+1):
        line = zip(ccs, vc[:,i])
        lines.extend(line)
    
    # Plot it
    import matplotlib.pyplot as plt
    from matplotlib.collections import LineCollection

    lines = LineCollection(lines, edgecolor='k')

    plt.plot(points[:,0], points[:,1], '.')
    plt.plot(ccs[:,0], ccs[:,1], '*')
    plt.gca().add_collection(lines)

    vertices_missing_edge_idx = np.sum(tri.neighbors == -1,axis=1)>0    #using sum as an OR function across all neighbors 
    vertices_missing_edge = ccs[vertices_missing_edge_idx]
    plt.plot(vertices_missing_edge[:,0],vertices_missing_edge[:,1],'x')

    plt.show()

voronoi(np.random.rand(30, 3))