Instantly share code, notes, and snippets.

Embed
What would you like to do?
Unconstrained Lloyd Iteration
from scipy.spatial import Voronoi, voronoi_plot_2d
import matplotlib.pyplot as plt
import numpy as np
import umap, os
def find_centroid(verts):
'''Return the centroid of a polygon described by `verts`'''
area = 0
x = 0
y = 0
for i in range(len(verts)-1):
step = (verts[i, 0] * verts[i+1, 1]) - (verts[i+1, 0] * verts[i, 1])
area += step
x += (verts[i, 0] + verts[i+1, 0]) * step
y += (verts[i, 1] + verts[i+1, 1]) * step
if area == 0: area += 0.01
return np.array([ (1/(3*area))*x, (1/(3*area))*y ])
def lloyd_iterate(X):
voronoi = Voronoi(X, qhull_options='Qbb Qc Qx')
centroids = []
for i in voronoi.regions:
region = voronoi.vertices[i + [i[0]]]
centroids.append( find_centroid( region ) )
return np.array(centroids)
def plot(X, name):
'''Plot the Voronoi map of 2D numpy array X'''
v = Voronoi(X, qhull_options='Qbb Qc Qx')
plot = voronoi_plot_2d(v, show_vertices=False, line_colors='y', line_alpha=0.5, point_size=5)
plot.set_figheight(14)
plot.set_figwidth(20)
plt.axis([-10, 10, -10, 10])
if not os.path.exists('plots'): os.makedirs('plots')
plot.savefig( 'plots/' + str(name) + '.png' )
# get 1000 observations in two dimensions and plot their Voronoi map
X = np.random.rand(1000, 4)
X = umap.UMAP().fit_transform(X)
plot(X, 0)
# run 4 iterations, plotting each result
for i in range(100):
X = lloyd_iterate(X)
plot(X, i)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment