import numpy as np from scvt import sphere, sphere_fast from scvt.timer import Timer from scvt.delaunay import delaunay from scvt.topology import all_edges_regions, regions_vertices from scvt.sphere import areas, circum def flatten_regions(regions, edges): degree = [len(region) for region in regions] vert = list(regions_vertices(regions, edges)) return np.asarray(degree), np.asarray(vert) def flatten_simplices(simplices): for i,j,k in simplices: yield i yield j yield k def lloyd(points, weight, verbose): with Timer('in one Lloyd iteration (total)', verbose): N = points.shape[0] simplices = delaunay(points,verbose) # dict_simp=dict(enumerate(simplices)) with Timer('in np.asarray(simplices)', verbose): simp = np.asarray(simplices) with Timer('enumerating edges', verbose): simplex_index={ simplex:n for n,simplex in enumerate(simplices) } edges, regions = all_edges_regions(simplices) # with Timer('computing circumcenters (Python)', verbose): # vertices = [circum(points[a,:],points[b,:],points[c,:]) for a,b,c in simplices] # vertices = np.asarray(vertices) # with Timer('computing circumcenters (jit)', verbose): # vertices = sphere.circum_fast(simp, points) with Timer('computing circumcenters (C)', verbose): vertices=np.zeros((len(simplices),3)) sphere_fast.circum_fast(simp, points,vertices) with Timer('computing barycenters (C)', verbose): regions_size, vert = flatten_regions(regions, edges) # degree, counts = np.unique(regions_size, return_counts=True) # print("Edge counts of Voronoi cells :", dict(zip(degree, counts))) w_i,w_v = weight(points), weight(vertices) bary = points.copy() sphere_fast.barycenters_fast(regions_size, vert, points, w_i, vertices, w_v, bary) # bary = sphere.barycenters_fast(regions_size, vert, points, vertices, weight) # bary = sphere.barycenters(regions, points, edges, vertices) if verbose: area = np.asarray(list(areas(points, simplices))) print('%d cells, %d edges, %d triangles'%(N, len(edges), len(simplices))) print('Total area = 4*pi + ',area.sum()-4*np.pi) if N<50: print(edges) a,b,c=points[simplices[0],:] test(a,b,c) test(a,c,b) test(c,a,b) return bary, vertices, edges