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