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