PL. 112 · LATTICES / RELAXATION / CENTROIDAL VORONOI
Lloyd Relaxation
Stuart P. Lloyd, 1957 · published 1982
OPEN THE LIVE PLATE ▸DEFINITION
cell(sᵢ) = ⋂ⱼ≠ᵢ { x : (x − sᵢ)·(sⱼ − sᵢ) ≤ ‖sⱼ − sᵢ‖²/2 }
cᵢ = ∫cell x dx / ∫cell dx
sᵢ ← sᵢ + ω(cᵢ − sᵢ), until mean‖cᵢ − sᵢ‖ < tol
NOTES
Scatter points at random, give every point the territory nearer to it than to any other, then move each point to the middle of its own territory and cut the territories again. Repeat and the scatter anneals: the cells even out, and away from the frame they walk into hexagons, which is the tile that minimises mean squared distance to its own centre among all the shapes that tile the plane. Lloyd wrote this down in 1957 at Bell Labs as a way to place the levels of a pulse-code quantizer — where to put the steps of a staircase so the average squared error is least — and showed an optimal quantizer must satisfy two conditions at once: every value goes to the nearest level, and every level sits at the centroid of what it receives. The algorithm is nothing but alternating them, and it is the same iteration later named k-means. The paper stayed an internal memorandum for twenty-five years before it was printed. What is drawn here is the migration rather than the diagram: each bright stub joins a site to the centroid pulling on it, and every stub shrinks to nothing exactly when the tessellation becomes centroidal. Run it far enough and the blue noise is gone — the point set has become a crystal with a few dislocations left in it, which is the outcome every sampling paper reaching for this method warns about.
PROVENANCE
- Origin
- S. P. Lloyd, "Least squares quantization in PCM", IEEE Transactions on Information Theory 28(2), March 1982, 129–137. Written at Bell Telephone Laboratories in 1957 and presented that September to the Institute of Mathematical Statistics, it circulated as a memorandum for twenty-five years and reached print only in the 1982 special issue on quantization
- Found twice
- J. Max, "Quantizing for minimum distortion", IRE Transactions on Information Theory 6(1), March 1960, 7–12, has the same two conditions and the same alternation, arrived at independently while the Lloyd memorandum was still unpublished. That is why the scalar quantizer is called Lloyd-Max, and why 1957 and 1982 are both on this plate
- The planar case
- Lloyd worked in one dimension against a general probability density. The fixed point drawn here — every site at the centroid of its own Voronoi cell — is the centroidal Voronoi tessellation, named and analysed in Q. Du, V. Faber and M. Gunzburger, "Centroidal Voronoi Tessellations: Applications and Algorithms", SIAM Review 41(4), 1999, 637–676. This plate runs the uniform-density case; the general one weights the centroid by the density and crowds sites where it is high
- Why hexagons
- A. Gersho, "Asymptotically optimal block quantization", IEEE Transactions on Information Theory 25(4), 1979, 373–380, conjectured that the optimal quantizer cell is a single polytope tiling the space, which in the plane is the regular hexagon. In two dimensions it is a theorem rather than a conjecture, from the sums-of-moments result of L. Fejes Toth, Acta Mathematica Academiae Scientiarum Hungaricae 10, 1959, 299–304; in three dimensions and above it is still open. The event this plate announces counts how many of its interior cells reached six sides, which is that statement as a measurement rather than a claim
- How the cell is built here
- The cell is defined as an intersection of half planes, so the plate builds it as one: start from the frame and clip by the perpendicular bisector against every other site, after I. E. Sutherland and G. W. Hodgman, "Reentrant polygon clipping", Communications of the ACM 17(1), 1974, 32–42. What survives is the exact polygon, which gives the area centroid in closed form by the shoelace formula and gives real segments to stroke. No lattice is sampled and no pixel is touched, which is why this plate exports as vector where PL. 29 cannot
- Not PL. 29, and not PL. 55
- The Voronoi plate is this partition held still and shaded; the subject here is the fixed-point iteration that moves the sites, and the partition is only what the iteration reads. The overlap with Poisson-disk sampling is real and worth stating: both end in a point set with no clumps and no gaps, and relaxation is a standard way to make blue noise. They differ in what they optimise — Poisson-disk rejects any sample inside a hard exclusion radius and stops, while this minimises quantization energy and keeps going, past blue noise into hexagonal order. That is the failure mode, not a bonus
- Standing
- Public domain — two optimality conditions and the iteration that alternates them, implemented from the papers
- Constants
- ω is the relaxation factor. ω = 1 is the step Lloyd wrote, which converges in well under a second and is over before it can be watched; the default under-relaxes so the migration is the exhibit rather than a transient. The tolerance is a percentage of the mean site spacing, so it means the same thing at any site count and any plate size, and it trades waiting against order: measured over six seeds at 48 sites, 0.9% settles in about 109 steps with 61% of the interior cells hexagonal, and 0.3% takes 273 steps to reach 69%
- Source
- doi:10.1109/TIT.1982.1056489
HOUDINI · VEX
The same published mathematics as a Detail Wrangle body. Paste it into a Wrangle with Run Over set to Detail; every constant is the published value plus a tweak channel, so Create Spare Parameters gives a slider that starts where the paper does.
// FORMA — PL. 112 · LLOYD RELAXATION — Stuart P. Lloyd, 1957 · published 1982
// cell(sᵢ) = ⋂ⱼ≠ᵢ { x : (x − sᵢ)·(sⱼ − sᵢ) ≤ ‖sⱼ − sᵢ‖²/2 }
// cᵢ = ∫cell x dx / ∫cell dx
// sᵢ ← sᵢ + ω(cᵢ − sᵢ), until mean‖cᵢ − sᵢ‖ < tol
// Paste into a Detail Wrangle (Run Over: Detail), no inputs needed.
// Written from the published mathematics, not adapted from any code.
// Constants arrive at their published values. Press the node's Create
// Spare Parameters button and every tweak becomes a slider — starting
// at 0, the published figure, and moving in the constant's own units.
// https://forma-gen.com/#plate=lloyd
float p_sites = 48 + chf('sites_tweak'); // sites · live 20 .. 110
float p_relax = 0.15 + chf('relax_tweak'); // ω — relaxation factor · live 0.1 .. 1
float p_tol = 0.9 + chf('tol_tweak'); // settled below (% of spacing) · live 0.3 .. 3
// The plate's own colour: FORMA's LATTICES accent as a cosine ramp,
// brightest near t = 0 and t = 1, near-black around t = 0.5.
vector forma_ramp(float t){
return set(
0.46 + 0.5 * cos(6.28318530718 * (t + 0)),
0.4185 + 0.4549 * cos(6.28318530718 * (t + 0.05)),
0.1389 + 0.151 * cos(6.28318530718 * (t + 0.1)));
}
// Lloyd relaxation run to its own convergence test in one cook: the cell of a
// site is the intersection of the half planes that put it nearer than each
// other site, so it is built that way — clip the frame by every perpendicular
// bisector in turn (Sutherland and Hodgman) and what survives is the exact
// polygon. Its area centroid closes in the shoelace formula, every site steps
// toward its own centroid by omega, and the sweep stops when the mean step has
// fallen below the tolerance, exactly as the plate stops. The plate does one
// step a frame and shows the migration; a cook has no clock, so this emits the
// settled tessellation the plate holds on screen: cell walls as two-point
// polylines, sites as points. The bright site-to-centroid stubs and the
// migration trails are the plate's motion devices and are zero here by
// construction, so they are not ported.
//
// Domain is the plate's 560-unit reference canvas with its own 5% inset, and
// canvas y is negated onto Houdini's up. Deterministic: the scatter is
// random(counted seed) on the plate's first-epoch seed, so every cook is the
// same tessellation.
float W = 560.0, H = 560.0;
int forma_n = int(rint(p_sites));
int forma_max = 400; // steps this cook will run before it emits
float mg = 0.05;
float x0 = mg * W, x1 = W - mg * W;
float y0 = mg * W, y1 = H - mg * W;
float spacing = sqrt((x1 - x0) * (y1 - y0) / float(forma_n));
float sx[], sy[], cx[], cy[];
resize(sx, forma_n); resize(sy, forma_n);
resize(cx, forma_n); resize(cy, forma_n);
int rc = 1957; // seeded(1957 + epoch*5701) at epoch 0
for (int i = 0; i < forma_n; i++){
sx[i] = x0 + random(rc) * (x1 - x0); rc++;
sy[i] = y0 + random(rc) * (y1 - y0); rc++;
}
// the clipper's buffers, sized once: a convex polygon gains at most one vertex
// per half plane, so four plus n-1 bounds it
int forma_cap = forma_n + 8;
float px[], py[], qx[], qy[], fs[];
int pg[], qg[];
resize(px, forma_cap); resize(py, forma_cap); resize(pg, forma_cap);
resize(qx, forma_cap); resize(qy, forma_cap); resize(qg, forma_cap);
resize(fs, forma_cap);
int lastpass = 0, steps = 0;
for (int pass = 0; pass <= forma_max + 1; pass++){
float drift = 0.0;
for (int i = 0; i < forma_n; i++){
float ax = sx[i], ay = sy[i];
// the frame, tagged -1 so its edges are told apart from cell walls
px[0] = x0; py[0] = y0; pg[0] = -1;
px[1] = x1; py[1] = y0; pg[1] = -1;
px[2] = x1; py[2] = y1; pg[2] = -1;
px[3] = x0; py[3] = y1; pg[3] = -1;
int cnt = 4;
float rad2 = 0.0;
for (int k = 0; k < 4; k++){
float dx = px[k] - ax, dy = py[k] - ay;
float d = dx * dx + dy * dy;
if (d > rad2) rad2 = d;
}
for (int j = 0; j < forma_n; j++){
if (j == i || cnt < 3) continue;
float dx = sx[j] - ax, dy = sy[j] - ay;
float d2 = dx * dx + dy * dy;
// no bisector between coincident sites; and a bisector standing
// further off than the furthest surviving vertex cannot reach what
// is left of the polygon, which is what keeps an n-squared sweep
// cheap without a spatial index
if (d2 < 1e-8 || d2 * 0.25 > rad2) continue;
float lim = d2 * 0.5;
for (int k = 0; k < cnt; k++)
fs[k] = (px[k] - ax) * dx + (py[k] - ay) * dy - lim;
int m = 0;
for (int k = 0; k < cnt; k++){
int k2 = (k + 1 == cnt) ? 0 : k + 1;
float fa = fs[k], fb = fs[k2];
if (fa <= 0.0){
qx[m] = px[k]; qy[m] = py[k]; qg[m] = pg[k]; m++;
if (fb > 0.0){
// leaving: the edge opening here runs along the
// bisector, so it belongs to site j
float u = fa / (fa - fb);
qx[m] = px[k] + (px[k2] - px[k]) * u;
qy[m] = py[k] + (py[k2] - py[k]) * u;
qg[m] = j; m++;
}
} else if (fb <= 0.0){
// re-entering: the edge continues along the one it was cut
// from, so it keeps that tag
float u = fa / (fa - fb);
qx[m] = px[k] + (px[k2] - px[k]) * u;
qy[m] = py[k] + (py[k2] - py[k]) * u;
qg[m] = pg[k]; m++;
}
}
float nr2 = 0.0;
for (int k = 0; k < m; k++){
px[k] = qx[k]; py[k] = qy[k]; pg[k] = qg[k];
float ddx = px[k] - ax, ddy = py[k] - ay;
float dd = ddx * ddx + ddy * ddy;
if (dd > nr2) nr2 = dd;
}
cnt = m; rad2 = nr2;
}
// Shoelace, taken about the site rather than about the origin. The two
// are algebraically identical; on this side they are not, because VEX
// floats carry seven digits and a cell of a few thousand square units
// sitting five hundred units out spends three of them on cancellation.
float a2 = 0.0, gx = 0.0, gy = 0.0;
for (int k = 0; k < cnt; k++){
int k2 = (k + 1 == cnt) ? 0 : k + 1;
float ux = px[k] - ax, uy = py[k] - ay;
float vx = px[k2] - ax, vy = py[k2] - ay;
float cr = ux * vy - vx * uy;
a2 += cr;
gx += (ux + vx) * cr;
gy += (uy + vy) * cr;
}
if (cnt >= 3 && abs(a2) > 1e-6){
gx = ax + gx / (3.0 * a2);
gy = ay + gy / (3.0 * a2);
} else {
gx = ax; gy = ay; // a cell squeezed to nothing stays put
}
cx[i] = gx; cy[i] = gy;
drift += sqrt((gx - ax) * (gx - ax) + (gy - ay) * (gy - ay));
if (lastpass){
// every interior wall belongs to two cells and both compute it;
// the lower-numbered site emits it, so the net is one prim per wall
for (int k = 0; k < cnt; k++){
if (pg[k] >= 0 && pg[k] <= i) continue;
int k2 = (k + 1 == cnt) ? 0 : k + 1;
vector col = forma_ramp(0.9);
int e0 = addpoint(0, set(px[k] - W / 2.0, (H / 2.0) - py[k], 0.0));
int e1 = addpoint(0, set(px[k2] - W / 2.0, (H / 2.0) - py[k2], 0.0));
setpointattrib(0, "Cd", e0, col);
setpointattrib(0, "Cd", e1, col);
setpointattrib(0, "Alpha", e0, 0.66);
setpointattrib(0, "Alpha", e1, 0.66);
addprim(0, "polyline", e0, e1);
}
}
}
if (lastpass) break;
// Settled means the mean distance from site to centroid has fallen below
// the tolerance, measured against the mean site spacing so the test says
// the same thing at any site count. One more pass then emits the cells it
// stopped on. forma_max is the cook's own ceiling: at the slowest corner
// of the declared box the plate itself takes about 690 steps, and a
// wrangle that runs that far is a wrangle nobody waits for.
if (drift / float(forma_n) / spacing < p_tol / 100.0 || steps >= forma_max){
lastpass = 1;
continue;
}
// omega is the relaxation factor: omega = 1 is the step Lloyd wrote, and
// anything below it walks the same descent in smaller moves. In a cook it
// decides the path rather than the destination, which is why the tolerance
// and not the step count is what stops this
for (int i = 0; i < forma_n; i++){
sx[i] += (cx[i] - sx[i]) * p_relax;
sy[i] += (cy[i] - sy[i]) * p_relax;
}
steps++;
}
// the sites themselves, in the ramp's bright lobe like the plate's own dots
for (int i = 0; i < forma_n; i++){
int pt = addpoint(0, set(sx[i] - W / 2.0, (H / 2.0) - sy[i], 0.0));
setpointattrib(0, "Cd", pt, forma_ramp(1.02));
setpointattrib(0, "Alpha", pt, 0.92);
}