[974] | 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') |
---|