[974] | 1 | import numpy as np |
---|
| 2 | |
---|
| 3 | from scvt import sphere, sphere_fast |
---|
| 4 | from scvt.timer import Timer |
---|
| 5 | from scvt.delaunay import delaunay |
---|
| 6 | from scvt.topology import all_edges_regions, regions_vertices |
---|
| 7 | from scvt.sphere import areas, circum |
---|
| 8 | |
---|
| 9 | def flatten_regions(regions, edges): |
---|
| 10 | degree = [len(region) for region in regions] |
---|
| 11 | vert = list(regions_vertices(regions, edges)) |
---|
| 12 | return np.asarray(degree), np.asarray(vert) |
---|
| 13 | |
---|
| 14 | def flatten_simplices(simplices): |
---|
| 15 | for i,j,k in simplices: |
---|
| 16 | yield i |
---|
| 17 | yield j |
---|
| 18 | yield k |
---|
| 19 | |
---|
| 20 | def lloyd(points, weight, verbose): |
---|
| 21 | with Timer('in one Lloyd iteration (total)', verbose): |
---|
| 22 | N = points.shape[0] |
---|
| 23 | simplices = delaunay(points,verbose) |
---|
| 24 | # dict_simp=dict(enumerate(simplices)) |
---|
| 25 | with Timer('in np.asarray(simplices)', verbose): |
---|
| 26 | simp = np.asarray(simplices) |
---|
| 27 | with Timer('enumerating edges', verbose): |
---|
| 28 | simplex_index={ simplex:n for n,simplex in enumerate(simplices) } |
---|
| 29 | edges, regions = all_edges_regions(simplices) |
---|
| 30 | |
---|
| 31 | # with Timer('computing circumcenters (Python)', verbose): |
---|
| 32 | # vertices = [circum(points[a,:],points[b,:],points[c,:]) for a,b,c in simplices] |
---|
| 33 | # vertices = np.asarray(vertices) |
---|
| 34 | # with Timer('computing circumcenters (jit)', verbose): |
---|
| 35 | # vertices = sphere.circum_fast(simp, points) |
---|
| 36 | with Timer('computing circumcenters (C)', verbose): |
---|
| 37 | vertices=np.zeros((len(simplices),3)) |
---|
| 38 | sphere_fast.circum_fast(simp, points,vertices) |
---|
| 39 | with Timer('computing barycenters (C)', verbose): |
---|
| 40 | regions_size, vert = flatten_regions(regions, edges) |
---|
| 41 | # degree, counts = np.unique(regions_size, return_counts=True) |
---|
| 42 | # print("Edge counts of Voronoi cells :", dict(zip(degree, counts))) |
---|
| 43 | w_i,w_v = weight(points), weight(vertices) |
---|
| 44 | bary = points.copy() |
---|
| 45 | sphere_fast.barycenters_fast(regions_size, vert, points, w_i, vertices, w_v, bary) |
---|
| 46 | # bary = sphere.barycenters_fast(regions_size, vert, points, vertices, weight) |
---|
| 47 | # bary = sphere.barycenters(regions, points, edges, vertices) |
---|
| 48 | |
---|
| 49 | if verbose: |
---|
| 50 | area = np.asarray(list(areas(points, simplices))) |
---|
| 51 | print('%d cells, %d edges, %d triangles'%(N, len(edges), len(simplices))) |
---|
| 52 | print('Total area = 4*pi + ',area.sum()-4*np.pi) |
---|
| 53 | |
---|
| 54 | if N<50: |
---|
| 55 | print(edges) |
---|
| 56 | a,b,c=points[simplices[0],:] |
---|
| 57 | test(a,b,c) |
---|
| 58 | test(a,c,b) |
---|
| 59 | test(c,a,b) |
---|
| 60 | |
---|
| 61 | return bary, vertices, edges |
---|