[974] | 1 | import numpy as np |
---|
| 2 | from scipy.spatial.qhull import _QhullUser |
---|
| 3 | from scipy.spatial.qhull import _Qhull |
---|
| 4 | |
---|
| 5 | from scvt.timer import Timer |
---|
| 6 | from scvt.sphere import unit_vector |
---|
| 7 | from scvt.topology import simplex |
---|
| 8 | |
---|
| 9 | ######################################################## |
---|
| 10 | |
---|
| 11 | # custom Delaunay class, slightly more efficient |
---|
| 12 | |
---|
| 13 | class Delaunay(_QhullUser): |
---|
| 14 | def __init__(self, points): |
---|
| 15 | qhull_options = b"Qbb Qc Qz Q12" |
---|
| 16 | # Run qhull |
---|
| 17 | qhull = _Qhull(b"d", points, qhull_options, required_options=b"Qt") |
---|
| 18 | _QhullUser.__init__(self, qhull) |
---|
| 19 | def _update(self, qhull): |
---|
| 20 | qhull.triangulate() |
---|
| 21 | self.paraboloid_scale, self.paraboloid_shift = \ |
---|
| 22 | qhull.get_paraboloid_shift_scale() |
---|
| 23 | self.simplices, self.neighbors, self.equations, self.coplanar = qhull.get_simplex_facet_array() |
---|
| 24 | _QhullUser._update(self, qhull) |
---|
| 25 | |
---|
| 26 | ######################################################## |
---|
| 27 | |
---|
| 28 | def testn(z): return z>0. |
---|
| 29 | def tests(z): return z<=0. |
---|
| 30 | |
---|
| 31 | def north_south(points, verbose): |
---|
| 32 | with Timer('in stereographic projection', verbose): |
---|
| 33 | north,south = [],[] |
---|
| 34 | for i in range(points.shape[0]): |
---|
| 35 | xyz=unit_vector(points[i,:]) |
---|
| 36 | points[i,:]=xyz |
---|
| 37 | x,y,z=xyz |
---|
| 38 | zn, zs = 1.+z, 1.-z |
---|
| 39 | if(zn>0.): north.append([i,[x/zn,y/zn]]) |
---|
| 40 | if(zs>0.): south.append([i,[x/zs,y/zs]]) |
---|
| 41 | indn, xyn = zip(*north) |
---|
| 42 | inds, xys = zip(*south) |
---|
| 43 | with Timer('in planar Delaunay', verbose): |
---|
| 44 | trin, tris = Delaunay(np.asarray(xyn)), Delaunay(np.asarray(xys)) |
---|
| 45 | |
---|
| 46 | return points[:,2], (indn, trin, testn), (inds, tris, tests) |
---|
| 47 | |
---|
| 48 | def enum_simplices(z,domains): |
---|
| 49 | # domains is a list of tuples (ind, tri, test) |
---|
| 50 | for ind, tri, test in domains: |
---|
| 51 | for i,j,k in tri.simplices: |
---|
| 52 | ii, jj, kk = ind[i], ind[j], ind[k] # global indices |
---|
| 53 | zi, zj, zk = z[ii], z[jj], z[kk] |
---|
| 54 | if test(zi) | test(zj) | test(zk): yield simplex(ii,jj,kk) |
---|
| 55 | |
---|
| 56 | def delaunay(points, verbose=False): |
---|
| 57 | z, north, south = north_south(points, verbose) |
---|
| 58 | with Timer('enumerating simplices', verbose): |
---|
| 59 | simplices = list(set( enum_simplices(z, (north, south) ) )) |
---|
| 60 | return simplices |
---|