1 | import numpy as np |
---|
2 | import matplotlib.pyplot as plt |
---|
3 | from numba import jit |
---|
4 | |
---|
5 | from scvt import topology |
---|
6 | from scvt import sphere |
---|
7 | from scvt import delaunay |
---|
8 | from scvt import lloyd |
---|
9 | |
---|
10 | def test(a,b,c): |
---|
11 | p=circum(a,b,c) |
---|
12 | print(norm(a-p),norm(b-p),norm(c-p)) # should be equal |
---|
13 | |
---|
14 | @jit |
---|
15 | def weight(xyz): |
---|
16 | y = xyz[:,1]; |
---|
17 | return 1.+15.*np.exp(-16.*y*y) |
---|
18 | |
---|
19 | def main(nbp,inner,outer,filename=None): |
---|
20 | N=10*nbp*nbp+2 |
---|
21 | points=np.random.randn(N,3) |
---|
22 | for i in range(outer): |
---|
23 | verbose = True |
---|
24 | print('Outer Lloyd iteration %d/%d'%(i,outer)) |
---|
25 | for j in range(inner): |
---|
26 | points, vertices, edges = lloyd.lloyd(points, weight, verbose) |
---|
27 | verbose = False |
---|
28 | plt.figure(figsize=(10,10)) |
---|
29 | sphere.plot_mesh(vertices, topology.dual_edges(edges), orient, "k") |
---|
30 | if N<1000: |
---|
31 | sphere.plot_mesh(points,edges,orient,"r") |
---|
32 | if filename is None: |
---|
33 | plt.show() |
---|
34 | else: |
---|
35 | plt.savefig(filename) |
---|
36 | |
---|
37 | |
---|
38 | |
---|
39 | orient=np.asarray([0,0,1.]) |
---|
40 | |
---|
41 | ################## Manual tests with octahedron #################### |
---|
42 | |
---|
43 | points=np.asarray([[1.,0.,0.],[0.,1.,0.],[0.,0.,1.], |
---|
44 | [-1.,0.,0.],[0.,-1.,0.],[0.,0.,-1.]]) |
---|
45 | simplices = delaunay.delaunay(points) |
---|
46 | print("== Simplices ==", simplices) |
---|
47 | edges, regions = topology.all_edges_regions(simplices) |
---|
48 | print('%d cells, %d edges, %d triangles'%(points.shape[0], len(edges), len(simplices))) |
---|
49 | print(edges) |
---|
50 | for n,region in enumerate(regions): |
---|
51 | print("Voronoi cell %d :"%n, region) |
---|
52 | for n,region in enumerate(regions): |
---|
53 | print("Voronoi cell %d :"%n, [edges[edge] for edge in region]) |
---|
54 | |
---|
55 | vertices = [sphere.circum(points[a,:],points[b,:],points[c,:]) for a,b,c in simplices] |
---|
56 | vertices = np.asarray(vertices) |
---|
57 | print("== Vertices ==", vertices) |
---|
58 | bary = sphere.barycenters(regions, points, edges, vertices) |
---|
59 | |
---|
60 | regions_size, vert = lloyd.flatten_regions(regions, edges) |
---|
61 | degree, counts = np.unique(regions_size, return_counts=True) |
---|
62 | print("Edge counts of Voronoi cells :", dict(zip(degree, counts))) |
---|
63 | |
---|
64 | ####################### Larger grid ####################### |
---|
65 | #main(8, 1,1) |
---|
66 | #main(8, 100,5) |
---|
67 | main(16, 100,5, 'nbp16.png') |
---|
68 | main(40, 100,10, 'nbp40.png') |
---|