Voronoi Tesselation of Random Points on the Unit Square

1-README.md

Voronoi Tesselation

10000 pixels partitions in to 30 regions as follows:

415,  91, 596, 420, 347, 239, 313, 298, 447, 259, 
414, 229, 137, 347, 353, 449, 744, 453, 479, 313, 
362, 363, 494,  87, 343, 395, 327, 140,  62, 285

Right now, this uses LLoyd's algorithm or k-means clustering on N=10⁴ points, where we have decided the centers in advance.

• In an earlier version, the mapping of the points to the clusters was not parallel. Fixed
• Naive implementations tend to run in N² time while Fortune's Sweep Line runs in N log N.
• There may be other time saving tricks as well: scicomp.stackexchange.com.

We can apply this to get Picasso-like images from normal photographs.

1-Voronoi.py

import numpy as np
import matplotlib.pyplot as plt

# N random points in the unit squre
vt = np.array([np.random.random() + 1j*np.random.random() for k in range(N)])

# grid of points to cluster
dz = 0.01
pts = np.arange(0,1.01, dz)[..., None] + 1j*np.arange(0,1.01,dz)[None, ...]

# LLOYD's ALGORITHM
d = np.abs(pts[...,...,None] - vt[None,None, ...])
col = d.argmin(axis=2)
 
# http://stackoverflow.com/questions/8218032/how-to-turn-a-boolean-array-into-index-array-in-numpy
# http://stackoverflow.com/questions/7179532/using-multiple-levels-of-boolean-index-mask-in-numpy

for n in range(N):
    mask = np.where(col == n)
    plt.plot(pts[mask].real, pts[mask].imag, '.', markersize=5, alpha = 0.75, color=colors[n])

plt.axis("Equal")
plt.axis("off")
plt.show()

2-fortune.py

N = 15
vt = np.array([np.random.random() + 1j*np.random.random() for k in range(N)])
plt.plot(vt.real, vt.imag, 'ro', markersize=3)


y = np.arange(0,1,0.005)

for x in np.arange(-1,1,0.005):
    # http://stackoverflow.com/questions/8218032/how-to-turn-a-boolean-array-into-index-array-in-numpy
    mask = np.where(vt.real - x >  0)    

    if len(mask[0]) > 0:      
        w = 0.5*((x + vt[mask].real) +  (y[..., None] - vt[mask].imag)**2/(vt[mask].real - x))
        plt.plot(np.amin(w, axis=1), y,  'k-')

plt.xlim([0,1])
plt.ylim([0,1])

plt.axis("Equal")
plt.axis("off")
plt.grid(True)

plt.show()

2-README.md

Fortune's Sweepline

You can get the Voronoi tessellation in quadratic time using naive methods or n log n time using Fortune's sweep line.

Amazing! You are an artist.