source: codes/mesh_generation/scvt/delaunay.py @ 1035

import numpy as np
from scipy.spatial.qhull import _QhullUser
from scipy.spatial.qhull import _Qhull

from scvt.timer import Timer
from scvt.sphere import unit_vector
from scvt.topology import simplex

########################################################

# custom Delaunay class, slightly more efficient

class Delaunay(_QhullUser):
    def __init__(self, points):
        qhull_options = b"Qbb Qc Qz Q12"
        # Run qhull
        qhull = _Qhull(b"d", points, qhull_options, required_options=b"Qt")
        _QhullUser.__init__(self, qhull)
    def _update(self, qhull):
        qhull.triangulate()
        self.paraboloid_scale, self.paraboloid_shift = \
                               qhull.get_paraboloid_shift_scale()
        self.simplices, self.neighbors, self.equations, self.coplanar = qhull.get_simplex_facet_array()
        _QhullUser._update(self, qhull)

########################################################        

def testn(z): return z>0.
def tests(z): return z<=0.

def north_south(points, verbose):
    with Timer('in stereographic projection', verbose):
        north,south = [],[]
        for i in range(points.shape[0]):
            xyz=unit_vector(points[i,:])
            points[i,:]=xyz
            x,y,z=xyz
            zn, zs = 1.+z, 1.-z
            if(zn>0.): north.append([i,[x/zn,y/zn]])
            if(zs>0.):  south.append([i,[x/zs,y/zs]])
        indn, xyn = zip(*north)
        inds, xys = zip(*south)
    with Timer('in planar Delaunay', verbose):
        trin, tris = Delaunay(np.asarray(xyn)), Delaunay(np.asarray(xys))

    return points[:,2], (indn, trin, testn), (inds, tris, tests)

def enum_simplices(z,domains):
    # domains is a list of tuples (ind, tri, test)
    for ind, tri, test in domains:
        for i,j,k in tri.simplices:
            ii, jj, kk = ind[i], ind[j], ind[k] # global indices
            zi, zj, zk = z[ii], z[jj], z[kk]
            if test(zi) | test(zj) | test(zk): yield simplex(ii,jj,kk)

def delaunay(points, verbose=False):
    z, north, south = north_south(points, verbose)
    with Timer('enumerating simplices', verbose):
        simplices = list(set( enum_simplices(z, (north, south) ) ))
    return simplices