Add lloyd_relaxation arg to Voronoi and Delaunay tellesllations

[jgostick]

https://en.wikipedia.org/wiki/Lloyd%27s_algorithm

[jgostick]

I played around with this today so have worked out the general functions, which I'm putting here instead of on a random branch which I'll probably forget about:

def lloyd_relaxation(vor, mode='rigorous'):
    pts = vor.points
    for i, r in enumerate(vor.point_region):
        if np.all(np.array(vor.regions[r]) > 0):
            pts[i, :pts.shape[1]] = get_centroid(vor.vertices[vor.regions[r]], mode=mode)
    if pts.shape[1] == 2:
        pts = np.c_[pts, np.zeros_like(pts[:, 0])]
    return pts

def get_centroid(pts, mode='rigorous'):
    if mode == 'rigorous':
        tri = sptl.Delaunay(pts)
        centroids = tri.points[tri.simplices].mean(axis=1)
        A = []
        for s in tri.simplices:
            xy = tri.points[s]
            if xy.shape[1] == 2:  # In 2D
                temp = (xy[:, 0].max() - xy[:, 0].min())*(xy[:, 1].max() - xy[:, 1].min())/2
            else:
                d = np.r_[xy.T,[np.ones_like(xy.T[0, :])]]  # In 3D
                temp = np.linalg.det(d)
            A.append(temp)
        A = np.array(A)
        CoM = np.array([(centroids[:, i]*A/A.sum()*len(A)).mean() for i in range(centroids.shape[1])])
    elif mode == 'fast':
        CoM = pts.mean(axis=0)
    return CoM

The lloyd_relaxation function performs a single iteration, so it must be called from within some other for-loop. The get_centroid function returns the centroid of a single voronoi cell, given the xy or xyz coordinates of the cell vertices. The mode controls how the centroid is found. The rigorous mode is quite slow, but finds the actual center of mass of the polyhedra, while the fast method merely takes the average of the vertex coordinates. The speed difference is noticable even on 250 cells, so it's probably quite important for large networks.

The 3D version has not been test yet!

[jgostick]

Also, here is a code snippet showing how to use the lloyd_relaxation function on a network:

import numpy as np
import openpm as op
import matplotlib.pyplot as plt

np.random.seed(3)
domain_size = [1, 1, 0]
n_points = 250
mode = 'rigorous'
pts = op._skgraph.generators.tools.generate_base_points(num_points=n_points, domain_size=domain_size, reflect=True)
pn, vor = op._skgraph.generators.voronoi(points=pts, shape=domain_size)
pts_orig = np.copy(pts)
pts_orig = pts_orig[np.all(pts_orig >= 0, axis=1), :]
pts_orig = pts_orig[np.all(pts_orig <= domain_size, axis=1), :]
for _ in range(5):
    pts = lloyd_relaxation(vor, mode=mode)
    pts = pts[np.all(pts >= 0, axis=1), :]
    pts = pts[np.all(pts <= domain_size, axis=1), :]
    pts = op._skgraph.generators.tools.reflect_base_points(points=pts, domain_size=domain_size)
    pn, vor = op._skgraph.generators.voronoi(points=pts, shape=domain_size)
ax = op._skgraph.visualization.plot_edges(pn)
# ax = op._skgraph.visualization.plot_nodes(pn, ax=ax)
plt.scatter(pts_orig[:, 0], pts_orig[:, 1], marker='x', cmap=plt.cm.rainbow, c=np.arange(n_points), s=10)
pts = pts[np.all(pts >= 0, axis=1), :]
pts = pts[np.all(pts <= domain_size, axis=1), :]
plt.scatter(pts[:, 0], pts[:, 1], marker='x', cmap=plt.cm.rainbow, c=np.arange(n_points), s=10)
for i in range(n_points):
    plt.plot((pts_orig[i, 0], pts[i, 0]), (pts_orig[i, 1], pts[i, 1]), '--', c='grey')
plt.axis(False)
plt.savefig(f'relaxed_voronoi_{mode}.svg');