source: codes/mesh_generation/scvt/lloyd.py @ 1028

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