31:T67d,description = md`# Lloyd’s Algorithm With Velocities A fork showing a variant on over-relaxation that adds velocities. The basic formula is essentially: > ${tex`\mathbf{v}_{n+1} = k_1\mathbf{v}_{n} + k_2\times`}*displacement from current position to centroid* Where ${tex`k_1`} is some value between 0 and 1 (fully damped to undamped respectively), and ${tex`k_2`} a value close to two. [See the appendix below](#appendix) for a discussion of choosing good values, and other details of the implementation (we also add clamping, jittering, a speed limit). The result is that our system converges much faster: ![](https://i.imgur.com/WmTe8oj.png) In this particular screenshot, the speed is over 10 times faster. This matches the amount of blue cells covered by a single red cell, which also seem to be close to a ratio of 10:1. This difference appears to increase with the number of points - this is mostly due to normal relaxation becoming slower and slower the more points there are. In fact, given enough points, the standard algorithm never reaches equilibrium *at all* in the set-up below. **Original Description** > [Lloyd’s algorithm](https://en.wikipedia.org/wiki/Lloyd's_algorithm) iteratively relaxes a Voronoi diagram: at each step, a new Voronoi diagram is constructed from the centroids of the cells from the previous step. This implementation uses [d3-delaunay](https://github.com/d3/d3-delaunay) and [d3-polygon](https://github.com/d3/d3-polygon); see also [Jason Davies’](https://www.jasondavies.com/lloyd/). The red lines in the animation below represent the Voronoi cells, while the black lines show the movement of the centroids.`32:T19f9,appendix = md`## Choosing ${tex`k_1`} and ${tex`k_2`} On the one hand, higher values of ${tex`k_1`} and ${tex`k_2`} lead to a faster convergence on an even distribution of points. On the other, it appears that there is a certain tipping point after which the speed of the points in the system only increases. The reason is that we are essentially adding momentum to our system: a positive value of ${tex`k_1`} means any displacement towards the centroid represents an acceleration of our point. This would be equivalent to adding kinetic energy faster than the system can lose. Interestingly, for ${tex`k_2 = 1`}, the tipping point for ${tex`k_1`} appears to be any value greater than ${tex`\frac{1}{\sqrt{2}}`}. In other words, the ideal value in that case is: > ${tex`\mathbf{v}_{n+1} = \frac{1}{\sqrt{2}}\mathbf{v}_{n} + `} *displacement from current position to centroid* Why ${tex`\frac{1}{\sqrt{2}}`}? No idea! It's a very suspicious value though, the kind that makes me think that there is some geometric explanation for it. Regardless, it turns out we can add other systems to push ${tex`k`} closer to 1. ## Boundaries and Overshoot With standard relaxation points can only converge to their own center of mass. That means they can never overshoot, which here means being displaced outside their own cell. Therefore they can never escape the bounding box of our Voronoi diagram. With velocities this changes however: after a displacement step the generator might end up outside its previous cell, and therefore the diagram's bounding box. Points being out of bounds can crash the algorithms in \`d3-delaunay\`. To fix that, we must clamp points our bounding box. This leads to another issue: with a simple clamp the points might end up overlapping - especially in the corners this is likely to happen. To prevent that we should also jitter their position. If we finally reset the velocities to zero as well, we reduce the overal speed of the system. This effectively is an [inelastic collision](https://en.wikipedia.org/wiki/Inelastic_collision), which loses some kinetic energy. This turns out to help in practice! So putting all of that together: // clamp to boundaries. If a point is over the edge, jitter position // to reduce chance of delaunaytor breaking due to having many points // in, say, (0,0). On top of that, set velocity to zero to lose some // kinetic energy and prevent bounce-back. let x1 = x0 + vx; let y1 = y0 + vy; if (x1 <= 0) { x1 = Math.random(); vx = 0; vy = 0; } else if (x1 >= width - 1) { x1 = width - 1 - Math.random(); vx = 0; vy = 0; } if (y1 <= 0) { y1 = Math.random(); vx = 0; vy = 0; } else if (y1 >= height - 1) { y1 = height - 1 - Math.random(); vx = 0; vy = 0; } Another option to both limit overshoot is adding a speed limit. But how to determine a good value for that limit? Well, lets think backwards from equilibrium. At this point, the average area of each cell in the Voronoi diagram will be: > ${tex`\frac{area}{cells} = \frac{width \times height}{cells}`} The Voronoi diagram will be approximately hexagonal, which is means the cells are approximately circular. So a decent enough approximation of the radius ${tex`r`} of the cell will be: > ${tex`r \approx \sqrt{\frac{width \times height}{cells \times \Pi}}`} If we limit the speed to this value, we ensure that there is an upper bound to the total energy of the system (namely, all points moving at this speed). This actually lets us reduce the damping factor even further! The reason this works is that when a distribution of points is close to an equilibrium, their displacement to the centroid will be smaller. This means the addition of kinetic energy to the sysem is reduced. As a result, as long as there is *some* damping factor, and *some* speed limit to avoid the energy from building up too high before it can go down again, eventually our system should lose its the kinetic energy. ## More thoughts on ${tex`\mathbf{v}_{n+1} = \frac{1}{\sqrt{2}}\mathbf{v}_{n} +`} *displacement* Why *that* constant value? I _suspect_ it might be a close enough approximation of [a critically damped oscillatory system](https://en.wikipedia.org/wiki/Damping_ratio). It most certainly isn't, though: when rendering very few points it visibly overshoots (perhaps the friction should be a product of speed instead? Or perhaps we should include the area of the cell in the equation somehow, like as a scaling factor?) I would love to know, but I am not sure how to determine what. Perhaps by reasoning backwards from a fully relaxed equilibrium, treating the generators as a system of connected simple springs with equal spring constant, and then reasoning about what would happen if we pulled all generators to the center? I dunno, I dropped out of physics over a decade ago. I did a quick cursory search on Google Scholar. These papers may be worth looking into (am still looking for more): Li, Hongwei, et al. [*"Fast capacity constrained Voronoi tessellation."*](https://dl.acm.org/citation.cfm?id=1730985) Proceedings of the 2010 ACM SIGGRAPH symposium on Interactive 3D Graphics and Games. ACM, 2010. - \`pdf mirrors:\` [microsoft.com](https://www.microsoft.com/en-us/research/publication/fast-capacity-constrained-voronoi-tessellation/), [researchgate.net](https://www.researchgate.net/profile/Chi_Wing_Fu/publication/220792129_Fast_Capacity_Constrained_Voronoi_Tessellation/links/547f2c190cf2de80e7cc77e1/Fast-Capacity-Constrained-Voronoi-Tessellation.pdf), [psu.edu](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.167.5579&rep=rep1&type=pdf)) Hateley, James C., Huayi Wei, and Long Chen. [*"Fast methods for computing centroidal Voronoi tessellations."*](http://web.math.ucsb.edu/~jhateley/paper/fmccvt.pdf) Journal of Scientific Computing 63.1 (2015): 185-212. - \`pdf mirrors:\` [researchgate.net](https://www.researchgate.net/profile/Huayi_Wei/publication/271952077_Fast_Methods_for_Computing_Centroidal_Voronoi_Tessellations/links/55b22d0508ae092e9650344a/Fast-Methods-for-Computing-Centroidal-Voronoi-Tessellations.pdf), [academia.edu](https://s3.amazonaws.com/academia.edu.documents/46742763/Fast_Methods_for_Computing_Centroidal_Vo20160623-1028-1392fm.pdf?AWSAccessKeyId=AKIAIWOWYYGZ2Y53UL3A&Expires=1540387269&Signature=CUFnqT2Zcpd7K75acNId7P%2Fo4wM%3D&response-content-disposition=inline%3B%20filename%3DFast_Methods_for_Computing_Centroidal_Vo.pdf) [↩back](#description) --- `33:T958,voronoiGenerator1 = { replay; const n = totalPoints; const points = new Float64Array(n*2); for(let i = 0; i < points.length;){ const t = Math.random() * Math.PI * 2; const a = Math.random(); points[i++] = a * Math.sin(t) + width * 0.5; points[i++] = a * Math.cos(t) + height * 0.5; } context1.font = '20pt sans'; context1.fillStyle = "white" context1.fillRect(0, 0, width, height); context1.lineCap = "square"; let delaunay, voronoi; let vsum = 0; let minSpeed = Math.sqrt(width*height/n)*0.005; let avg = minSpeed*2; let maxSpeed = minSpeed*2; let strIterations = ''; let iterations = 0; let rendering = renderDefault; delaunay = new Delaunay(points); while(rendering) { context1.fillStyle = clearView ? "white" : "#FFF4"; context1.fillRect(0, 0, width, height); voronoi = delaunay.voronoi([0, 0, width, height]); voronoi.render(context1); context1.lineWidth = 1; context1.strokeStyle = clearView ? "#D00" : "#8002"; context1.stroke(); context1.beginPath(); vsum = 0; maxSpeed = 0; for (let i = 0; i < n; ++i) { const x0 = points[i * 2]; const y0 = points[i * 2 + 1]; let [x1, y1] = d3.polygonCentroid(voronoi.cellPolygon(i)); points[i * 2] = x1; points[i * 2 + 1] = y1; const dx = x1 - x0; const dy = y1 - y0; let speed = Math.sqrt(dx*dx + dy*dy); vsum += speed || 0; if (speed > maxSpeed) maxSpeed = speed; speed = speed < 3 ? 3/speed : 1; context1.moveTo(x0, y0); context1.lineTo(dx*speed + x0, dy*speed + y0); } context1.lineWidth = 1; context1.stroke(); avg = (avg*3 + vsum/n)*0.25; rendering = maxSpeed > minSpeed; strIterations = ' - ' + (rendering ? '' : 'finished in ') + iterations + ' iterations'; const avgSpeedStr = `Average speed: ${avg.toPrecision(4)} (default${strIterations})`; const maxSpeedStr = `Max speed: ${maxSpeed.toPrecision(4)}`; const ctx = renderBlend && renderVelocities ? statsContext1 : context1; statsContext1.clearRect(0, 0, width, height); ctx.lineWidth = 8; ctx.strokeStyle = "#FFF8"; ctx.fillStyle = "black"; ctx.strokeText(avgSpeedStr, 10, 40); ctx.fillText(avgSpeedStr, 10, 40); ctx.strokeText(maxSpeedStr, 10, 80); ctx.fillText(maxSpeedStr, 10, 80); yield ++iterations; delaunay.update(); } }34:Te21,voronoiGenerator2 = { replay; const n = totalPoints; const points = new Float64Array(n*2); for(let i = 0; i < points.length;){ const t = Math.random() * Math.PI * 2; const a = Math.random(); points[i++] = a * Math.sin(t) + width * 0.5; points[i++] = a * Math.cos(t) + height * 0.5; } const velocities = new Float64Array(n*2); for(let i = 0, x0 = width*0.5, y0 = height*0.5; i < points.length; i+=2) { velocities[i] = (points[i] - x0)*2 velocities[i + 1] = (points[i + 1] - y0)*2; } context2.font = '20pt sans'; context2.fillStyle = "white" context2.fillRect(0, 0, width, height); context2.lineCap = "square"; let delaunay, voronoi; let vsum = 0; let minSpeed = Math.sqrt(width*height/n)*0.005; let avg = minSpeed*2; let maxSpeed = minSpeed*2; let strIterations = ''; let iterations = 0; let rendering = renderVelocities delaunay = new Delaunay(points); while(rendering) { context2.fillStyle = clearView ? "white" : "#FFF4"; context2.fillRect(0, 0, width, height); voronoi = delaunay.voronoi([0, 0, width, height]); context2.beginPath(); voronoi.render(context2); context2.lineWidth = 1; context2.strokeStyle = clearView ? "#44F" : "#44F2"; context2.stroke(); context2.beginPath(); vsum = 0; maxSpeed = 0; for (let i = 0; i < n; ++i) { const [x0, y0] = [points[i * 2], points[i * 2 + 1]]; let [cx, cy] = d3.polygonCentroid(voronoi.cellPolygon(i)); let [vx, vy] = [velocities[i * 2], velocities[i * 2 + 1]]; const dx = cx - x0; const dy = cy - y0; vx = (vx*damping + 1.925*dx); vy = (vy*damping + 1.925*dy); // clamp to boundaries. Jitter edges to reduce chances // of delaunaytor breaking due to having many points in, say, (0,0). let x1 = x0 + vx; if (x1 <= 1) { x1 = 1 + Math.random(); vx = 0; vy = 0; } else if (x1 >= width - 2) { x1 = width - 2 - Math.random(); vx = 0; vy = 0; } let y1 = y0 + vy; if (y1 <= 1) { y1 = 1 + Math.random(); vx = 0; vy = 0; } else if (y1 >= height - 2) { y1 = height - 2 - Math.random(); vx = 0; vy = 0; } let speed = Math.sqrt(vx*vx + vy*vy) || 0; vsum += speed; if (speed > speedLimit) { vx *= speedLimit/speed; vy *= speedLimit/speed; speed = speedLimit; } if (speed > maxSpeed) maxSpeed = speed; velocities[i * 2] = vx; velocities[i * 2 + 1] = vy; points[i * 2] = x1; points[i * 2 + 1] = y1; // emphasize movement if les than a pixel speed = speed < 3 ? 3 / speed : 1; context2.moveTo(x0, y0), context2.lineTo(speed*(x1-x0) + x0, speed*(y1-y0) + y0); } context2.lineWidth = 1; context2.stroke(); avg = (avg*3 + vsum/n)*0.25; rendering = maxSpeed > minSpeed; strIterations = ' - ' + (rendering ? '' : 'finished in ') + iterations + ' iterations '; const avgSpeedStr = `Average speed: ${avg.toPrecision(4)} (momentum${strIterations})`; const maxSpeedStr = `Max speed: ${maxSpeed.toPrecision(4)}`; const ctx = renderBlend && renderDefault ? statsContext2 : context2; statsContext2.clearRect(0, 0, width, height); ctx.lineWidth = 8; ctx.strokeStyle = "#FFF8"; ctx.fillStyle = "black"; ctx.strokeText(avgSpeedStr, 10, 120); ctx.fillText(avgSpeedStr, 10, 120); ctx.strokeText(maxSpeedStr, 10, 160); ctx.fillText(maxSpeedStr, 10, 160); yield ++iterations; delaunay.update(); } }35:T3908,Delaunator = { const EPSILON = Math.pow(2, -52); const EDGE_STACK = new Uint32Array(512); class Delaunator { static from(points, getX = defaultGetX, getY = defaultGetY) { const n = points.length; const coords = new Float64Array(n * 2); for (let i = 0; i < n; i++) { const p = points[i]; coords[2 * i] = getX(p); coords[2 * i + 1] = getY(p); } return new Delaunator(coords); } constructor(coords) { const n = coords.length >> 1; if (n > 0 && typeof coords[0] !== 'number') throw new Error('Expected coords to contain numbers.'); this.coords = coords; // arrays that will store the triangulation graph const maxTriangles = Math.max(2 * n - 5, 0); this._triangles = new Uint32Array(maxTriangles * 3); this._halfedges = new Int32Array(maxTriangles * 3); // temporary arrays for tracking the edges of the advancing convex hull this._hashSize = Math.ceil(Math.sqrt(n)); this._hullPrev = new Uint32Array(n); // edge to prev edge this._hullNext = new Uint32Array(n); // edge to next edge this._hullTri = new Uint32Array(n); // edge to adjacent triangle this._hullHash = new Int32Array(this._hashSize).fill(-1); // angular edge hash // temporary arrays for sorting points this._ids = new Uint32Array(n); this._dists = new Float64Array(n); this.update(); } update() { const {coords, _hullPrev: hullPrev, _hullNext: hullNext, _hullTri: hullTri, _hullHash: hullHash} = this; const n = coords.length >> 1; // populate an array of point indices; calculate input data bbox let minX = Infinity; let minY = Infinity; let maxX = -Infinity; let maxY = -Infinity; for (let i = 0; i < n; i++) { const x = coords[2 * i]; const y = coords[2 * i + 1]; if (x < minX) minX = x; if (y < minY) minY = y; if (x > maxX) maxX = x; if (y > maxY) maxY = y; this._ids[i] = i; } const cx = (minX + maxX) / 2; const cy = (minY + maxY) / 2; let minDist = Infinity; let i0, i1, i2; // pick a seed point close to the center for (let i = 0; i < n; i++) { const d = dist(cx, cy, coords[2 * i], coords[2 * i + 1]); if (d < minDist) { i0 = i; minDist = d; } } const i0x = coords[2 * i0]; const i0y = coords[2 * i0 + 1]; minDist = Infinity; // find the point closest to the seed for (let i = 0; i < n; i++) { if (i === i0) continue; const d = dist(i0x, i0y, coords[2 * i], coords[2 * i + 1]); if (d < minDist && d > 0) { i1 = i; minDist = d; } } let i1x = coords[2 * i1]; let i1y = coords[2 * i1 + 1]; let minRadius = Infinity; // find the third point which forms the smallest circumcircle with the first two for (let i = 0; i < n; i++) { if (i === i0 || i === i1) continue; const r = circumradius(i0x, i0y, i1x, i1y, coords[2 * i], coords[2 * i + 1]); if (r < minRadius) { i2 = i; minRadius = r; } } let i2x = coords[2 * i2]; let i2y = coords[2 * i2 + 1]; if (minRadius === Infinity) { // order collinear points by dx (or dy if all x are identical) // and return the list as a hull for (let i = 0; i < n; i++) { this._dists[i] = (coords[2 * i] - coords[0]) || (coords[2 * i + 1] - coords[1]); } quicksort(this._ids, this._dists, 0, n - 1); const hull = new Uint32Array(n); let j = 0; for (let i = 0, d0 = -Infinity; i < n; i++) { const id = this._ids[i]; if (this._dists[id] > d0) { hull[j++] = id; d0 = this._dists[id]; } } this.hull = hull.subarray(0, j); this.triangles = new Uint32Array(0); this.halfedges = new Uint32Array(0); return; } // swap the order of the seed points for counter-clockwise orientation if (orient(i0x, i0y, i1x, i1y, i2x, i2y)) { const i = i1; const x = i1x; const y = i1y; i1 = i2; i1x = i2x; i1y = i2y; i2 = i; i2x = x; i2y = y; } const center = circumcenter(i0x, i0y, i1x, i1y, i2x, i2y); this._cx = center.x; this._cy = center.y; for (let i = 0; i < n; i++) { this._dists[i] = dist(coords[2 * i], coords[2 * i + 1], center.x, center.y); } // sort the points by distance from the seed triangle circumcenter quicksort(this._ids, this._dists, 0, n - 1); // set up the seed triangle as the starting hull this._hullStart = i0; let hullSize = 3; hullNext[i0] = hullPrev[i2] = i1; hullNext[i1] = hullPrev[i0] = i2; hullNext[i2] = hullPrev[i1] = i0; hullTri[i0] = 0; hullTri[i1] = 1; hullTri[i2] = 2; hullHash.fill(-1); hullHash[this._hashKey(i0x, i0y)] = i0; hullHash[this._hashKey(i1x, i1y)] = i1; hullHash[this._hashKey(i2x, i2y)] = i2; this.trianglesLen = 0; this._addTriangle(i0, i1, i2, -1, -1, -1); for (let k = 0, xp, yp; k < this._ids.length; k++) { const i = this._ids[k]; const x = coords[2 * i]; const y = coords[2 * i + 1]; // skip near-duplicate points if (k > 0 && Math.abs(x - xp) <= EPSILON && Math.abs(y - yp) <= EPSILON) continue; xp = x; yp = y; // skip seed triangle points if (i === i0 || i === i1 || i === i2) continue; // find a visible edge on the convex hull using edge hash let start = 0; for (let j = 0, key = this._hashKey(x, y); j < this._hashSize; j++) { start = hullHash[(key + j) % this._hashSize]; if (start !== -1 && start !== hullNext[start]) break; } start = hullPrev[start]; let e = start, q; while (q = hullNext[e], !orient(x, y, coords[2 * e], coords[2 * e + 1], coords[2 * q], coords[2 * q + 1])) { e = q; if (e === start) { e = -1; break; } } if (e === -1) continue; // likely a near-duplicate point; skip it // add the first triangle from the point let t = this._addTriangle(e, i, hullNext[e], -1, -1, hullTri[e]); // recursively flip triangles from the point until they satisfy the Delaunay condition hullTri[i] = this._legalize(t + 2); hullTri[e] = t; // keep track of boundary triangles on the hull hullSize++; // walk forward through the hull, adding more triangles and flipping recursively let n = hullNext[e]; while (q = hullNext[n], orient(x, y, coords[2 * n], coords[2 * n + 1], coords[2 * q], coords[2 * q + 1])) { t = this._addTriangle(n, i, q, hullTri[i], -1, hullTri[n]); hullTri[i] = this._legalize(t + 2); hullNext[n] = n; // mark as removed hullSize--; n = q; } // walk backward from the other side, adding more triangles and flipping if (e === start) { while (q = hullPrev[e], orient(x, y, coords[2 * q], coords[2 * q + 1], coords[2 * e], coords[2 * e + 1])) { t = this._addTriangle(q, i, e, -1, hullTri[e], hullTri[q]); this._legalize(t + 2); hullTri[q] = t; hullNext[e] = e; // mark as removed hullSize--; e = q; } } // update the hull indices this._hullStart = hullPrev[i] = e; hullNext[e] = hullPrev[n] = i; hullNext[i] = n; // save the two new edges in the hash table hullHash[this._hashKey(x, y)] = i; hullHash[this._hashKey(coords[2 * e], coords[2 * e + 1])] = e; } this.hull = new Uint32Array(hullSize); for (let i = 0, e = this._hullStart; i < hullSize; i++) { this.hull[i] = e; e = hullNext[e]; } // trim typed triangle mesh arrays this.triangles = this._triangles.subarray(0, this.trianglesLen); this.halfedges = this._halfedges.subarray(0, this.trianglesLen); } _hashKey(x, y) { return Math.floor(pseudoAngle(x - this._cx, y - this._cy) * this._hashSize) % this._hashSize; } _legalize(a) { const {_triangles: triangles, _halfedges: halfedges, coords} = this; let i = 0; let ar = 0; // recursion eliminated with a fixed-size stack while (true) { const b = halfedges[a]; /* if the pair of triangles doesn't satisfy the Delaunay condition * (p1 is inside the circumcircle of [p0, pl, pr]), flip them, * then do the same check/flip recursively for the new pair of triangles * * pl pl * /||\ / \ * al/ || \bl al/ \a * / || \ / \ * / a||b \ flip /___ar___\ * p0\ || /p1 => p0\---bl---/p1 * \ || / \ / * ar\ || /br b\ /br * \||/ \ / * pr pr */ const a0 = a - a % 3; ar = a0 + (a + 2) % 3; if (b === -1) { // convex hull edge if (i === 0) break; a = EDGE_STACK[--i]; continue; } const b0 = b - b % 3; const al = a0 + (a + 1) % 3; const bl = b0 + (b + 2) % 3; const p0 = triangles[ar]; const pr = triangles[a]; const pl = triangles[al]; const p1 = triangles[bl]; const illegal = inCircle( coords[2 * p0], coords[2 * p0 + 1], coords[2 * pr], coords[2 * pr + 1], coords[2 * pl], coords[2 * pl + 1], coords[2 * p1], coords[2 * p1 + 1]); if (illegal) { triangles[a] = p1; triangles[b] = p0; const hbl = halfedges[bl]; // edge swapped on the other side of the hull (rare); fix the halfedge reference if (hbl === -1) { let e = this._hullStart; do { if (this._hullTri[e] === bl) { this._hullTri[e] = a; break; } e = this._hullPrev[e]; } while (e !== this._hullStart); } this._link(a, hbl); this._link(b, halfedges[ar]); this._link(ar, bl); const br = b0 + (b + 1) % 3; // don't worry about hitting the cap: it can only happen on extremely degenerate input if (i < EDGE_STACK.length) { EDGE_STACK[i++] = br; } } else { if (i === 0) break; a = EDGE_STACK[--i]; } } return ar; } _link(a, b) { this._halfedges[a] = b; if (b !== -1) this._halfedges[b] = a; } // add a new triangle given vertex indices and adjacent half-edge ids _addTriangle(i0, i1, i2, a, b, c) { const t = this.trianglesLen; this._triangles[t] = i0; this._triangles[t + 1] = i1; this._triangles[t + 2] = i2; this._link(t, a); this._link(t + 1, b); this._link(t + 2, c); this.trianglesLen += 3; return t; } } // monotonically increases with real angle, but doesn't need expensive trigonometry function pseudoAngle(dx, dy) { const p = dx / (Math.abs(dx) + Math.abs(dy)); return (dy > 0 ? 3 - p : 1 + p) / 4; // [0..1] } function dist(ax, ay, bx, by) { const dx = ax - bx; const dy = ay - by; return dx * dx + dy * dy; } function orient(px, py, qx, qy, rx, ry) { return (qy - py) * (rx - qx) - (qx - px) * (ry - qy) < 0; } function inCircle(ax, ay, bx, by, cx, cy, px, py) { const dx = ax - px; const dy = ay - py; const ex = bx - px; const ey = by - py; const fx = cx - px; const fy = cy - py; const ap = dx * dx + dy * dy; const bp = ex * ex + ey * ey; const cp = fx * fx + fy * fy; return dx * (ey * cp - bp * fy) - dy * (ex * cp - bp * fx) + ap * (ex * fy - ey * fx) < 0; } function circumradius(ax, ay, bx, by, cx, cy) { const dx = bx - ax; const dy = by - ay; const ex = cx - ax; const ey = cy - ay; const bl = dx * dx + dy * dy; const cl = ex * ex + ey * ey; const d = 0.5 / (dx * ey - dy * ex); const x = (ey * bl - dy * cl) * d; const y = (dx * cl - ex * bl) * d; return x * x + y * y; } function circumcenter(ax, ay, bx, by, cx, cy) { const dx = bx - ax; const dy = by - ay; const ex = cx - ax; const ey = cy - ay; const bl = dx * dx + dy * dy; const cl = ex * ex + ey * ey; const d = 0.5 / (dx * ey - dy * ex); const x = ax + (ey * bl - dy * cl) * d; const y = ay + (dx * cl - ex * bl) * d; return {x, y}; } function quicksort(ids, dists, left, right) { if (right - left <= 20) { for (let i = left + 1; i <= right; i++) { const temp = ids[i]; const tempDist = dists[temp]; let j = i - 1; while (j >= left && dists[ids[j]] > tempDist) ids[j + 1] = ids[j--]; ids[j + 1] = temp; } } else { const median = (left + right) >> 1; let i = left + 1; let j = right; swap(ids, median, i); if (dists[ids[left]] > dists[ids[right]]) swap(ids, left, right); if (dists[ids[i]] > dists[ids[right]]) swap(ids, i, right); if (dists[ids[left]] > dists[ids[i]]) swap(ids, left, i); const temp = ids[i]; const tempDist = dists[temp]; while (true) { do i++; while (dists[ids[i]] < tempDist); do j--; while (dists[ids[j]] > tempDist); if (j < i) break; swap(ids, i, j); } ids[left + 1] = ids[j]; ids[j] = temp; if (right - i + 1 >= j - left) { quicksort(ids, dists, i, right); quicksort(ids, dists, left, j - 1); } else { quicksort(ids, dists, left, j - 1); quicksort(ids, dists, i, right); } } } function swap(arr, i, j) { const tmp = arr[i]; arr[i] = arr[j]; arr[j] = tmp; } function defaultGetX(p) { return p[0]; } function defaultGetY(p) { return p[1]; } return Delaunator; }36:T277b, class Voronoi { constructor(delaunay, [xmin, ymin, xmax, ymax] = [0, 0, 960, 500]) { if (!((xmax = +xmax) >= (xmin = +xmin)) || !((ymax = +ymax) >= (ymin = +ymin))) throw new Error("invalid bounds"); const {points, hull, triangles} = this.delaunay = delaunay; const circumcenters = this.circumcenters = new Float64Array(triangles.length / 3 * 2); const vectors = this.vectors = new Float64Array(points.length * 2); this.xmax = xmax, this.xmin = xmin; this.ymax = ymax, this.ymin = ymin; // Compute circumcenters. for (let i = 0, j = 0, n = triangles.length; i < n; i += 3, j += 2) { const t1 = triangles[i] * 2; const t2 = triangles[i + 1] * 2; const t3 = triangles[i + 2] * 2; const x1 = points[t1]; const y1 = points[t1 + 1]; const x2 = points[t2]; const y2 = points[t2 + 1]; const x3 = points[t3]; const y3 = points[t3 + 1]; const a2 = x1 - x2; const a3 = x1 - x3; const b2 = y1 - y2; const b3 = y1 - y3; const d1 = x1 * x1 + y1 * y1; const d2 = d1 - x2 * x2 - y2 * y2; const d3 = d1 - x3 * x3 - y3 * y3; const ab = (a3 * b2 - a2 * b3) * 2; // degenerate case (2 points) if (!ab) { circumcenters[j] = (x1 + x3) / 2 + 1e8 * b3; circumcenters[j + 1] = (y1 + y3) / 2 - 1e8 * a3; } else { circumcenters[j] = (b2 * d3 - b3 * d2) / ab; circumcenters[j + 1] = (a3 * d2 - a2 * d3) / ab; } } // Compute exterior cell rays. let h = hull[hull.length - 1]; let p0, p1 = h * 4; let x0, x1 = points[2 * h]; let y0, y1 = points[2 * h + 1]; for (let i = 0; i < hull.length; ++i) { h = hull[i]; p0 = p1, x0 = x1, y0 = y1; p1 = h * 4, x1 = points[2 * h], y1 = points[2 * h + 1]; vectors[p0 + 2] = vectors[p1] = y0 - y1; vectors[p0 + 3] = vectors[p1 + 1] = x1 - x0; } } render(context) { const buffer = context == null ? context = new Path : undefined; const {delaunay: {halfedges, inedges, hull}, circumcenters, vectors} = this; if (hull.length <= 1) return null; for (let i = 0, n = halfedges.length; i < n; ++i) { const j = halfedges[i]; if (j < i) continue; const ti = Math.floor(i / 3) * 2; const tj = Math.floor(j / 3) * 2; const xi = circumcenters[ti]; const yi = circumcenters[ti + 1]; const xj = circumcenters[tj]; const yj = circumcenters[tj + 1]; this._renderSegment(xi, yi, xj, yj, context); } let h0, h1 = hull[hull.length - 1]; for (let i = 0; i < hull.length; ++i) { h0 = h1, h1 = hull[i]; const t = Math.floor(inedges[h1] / 3) * 2; const x = circumcenters[t]; const y = circumcenters[t + 1]; const v = h0 * 4; const p = this._project(x, y, vectors[v + 2], vectors[v + 3]); if (p) this._renderSegment(x, y, p[0], p[1], context); } return buffer && buffer.value(); } renderBounds(context) { const buffer = context == null ? context = new Path : undefined; context.rect(this.xmin, this.ymin, this.xmax - this.xmin, this.ymax - this.ymin); return buffer && buffer.value(); } renderCell(i, context) { const buffer = context == null ? context = new Path : undefined; const points = this._clip(i); if (points === null) return; context.moveTo(points[0], points[1]); for (let i = 2, n = points.length; i < n; i += 2) { context.lineTo(points[i], points[i + 1]); } context.closePath(); return buffer && buffer.value(); } *cellPolygons() { const {delaunay: {points}} = this; for (let i = 0, n = points.length / 2; i < n; ++i) { const cell = this.cellPolygon(i); if (cell) yield cell; } } cellPolygon(i) { const polygon = new Polygon; this.renderCell(i, polygon); return polygon.value(); } _renderSegment(x0, y0, x1, y1, context) { let S; const c0 = this._regioncode(x0, y0); const c1 = this._regioncode(x1, y1); if (c0 === 0 && c1 === 0) { context.moveTo(x0, y0); context.lineTo(x1, y1); } else if (S = this._clipSegment(x0, y0, x1, y1, c0, c1)) { context.moveTo(S[0], S[1]); context.lineTo(S[2], S[3]); } } contains(i, x, y) { if ((x = +x, x !== x) || (y = +y, y !== y)) return false; return this.delaunay._step(i, x, y) === i; } _cell(i) { const {circumcenters, delaunay: {inedges, halfedges, triangles}} = this; const e0 = inedges[i]; if (e0 === -1) return null; // coincident point const points = []; let e = e0; do { const t = Math.floor(e / 3); points.push(circumcenters[t * 2], circumcenters[t * 2 + 1]); e = e % 3 === 2 ? e - 2 : e + 1; if (triangles[e] !== i) break; // bad triangulation e = halfedges[e]; } while (e !== e0 && e !== -1); return points; } _clip(i) { // degenerate case (1 valid point: return the box) if (i === 0 && this.delaunay.hull.length === 1) { return [this.xmax, this.ymin, this.xmax, this.ymax, this.xmin, this.ymax, this.xmin, this.ymin]; } const points = this._cell(i); if (points === null) return null; const {vectors: V} = this; const v = i * 4; return V[v] || V[v + 1] ? this._clipInfinite(i, points, V[v], V[v + 1], V[v + 2], V[v + 3]) : this._clipFinite(i, points); } _clipFinite(i, points) { const n = points.length; let P = null; let x0, y0, x1 = points[n - 2], y1 = points[n - 1]; let c0, c1 = this._regioncode(x1, y1); let e0, e1; for (let j = 0; j < n; j += 2) { x0 = x1, y0 = y1, x1 = points[j], y1 = points[j + 1]; c0 = c1, c1 = this._regioncode(x1, y1); if (c0 === 0 && c1 === 0) { e0 = e1, e1 = 0; if (P) P.push(x1, y1); else P = [x1, y1]; } else { let S, sx0, sy0, sx1, sy1; if (c0 === 0) { if ((S = this._clipSegment(x0, y0, x1, y1, c0, c1)) === null) continue; [sx0, sy0, sx1, sy1] = S; } else { if ((S = this._clipSegment(x1, y1, x0, y0, c1, c0)) === null) continue; [sx1, sy1, sx0, sy0] = S; e0 = e1, e1 = this._edgecode(sx0, sy0); if (e0 && e1) this._edge(i, e0, e1, P, P.length); if (P) P.push(sx0, sy0); else P = [sx0, sy0]; } e0 = e1, e1 = this._edgecode(sx1, sy1); if (e0 && e1) this._edge(i, e0, e1, P, P.length); if (P) P.push(sx1, sy1); else P = [sx1, sy1]; } } if (P) { e0 = e1, e1 = this._edgecode(P[0], P[1]); if (e0 && e1) this._edge(i, e0, e1, P, P.length); } else if (this.contains(i, (this.xmin + this.xmax) / 2, (this.ymin + this.ymax) / 2)) { return [this.xmax, this.ymin, this.xmax, this.ymax, this.xmin, this.ymax, this.xmin, this.ymin]; } return P; } _clipSegment(x0, y0, x1, y1, c0, c1) { while (true) { if (c0 === 0 && c1 === 0) return [x0, y0, x1, y1]; if (c0 & c1) return null; let x, y, c = c0 || c1; if (c & 0b1000) x = x0 + (x1 - x0) * (this.ymax - y0) / (y1 - y0), y = this.ymax; else if (c & 0b0100) x = x0 + (x1 - x0) * (this.ymin - y0) / (y1 - y0), y = this.ymin; else if (c & 0b0010) y = y0 + (y1 - y0) * (this.xmax - x0) / (x1 - x0), x = this.xmax; else y = y0 + (y1 - y0) * (this.xmin - x0) / (x1 - x0), x = this.xmin; if (c0) x0 = x, y0 = y, c0 = this._regioncode(x0, y0); else x1 = x, y1 = y, c1 = this._regioncode(x1, y1); } } _clipInfinite(i, points, vx0, vy0, vxn, vyn) { let P = Array.from(points), p; if (p = this._project(P[0], P[1], vx0, vy0)) P.unshift(p[0], p[1]); if (p = this._project(P[P.length - 2], P[P.length - 1], vxn, vyn)) P.push(p[0], p[1]); if (P = this._clipFinite(i, P)) { for (let j = 0, n = P.length, c0, c1 = this._edgecode(P[n - 2], P[n - 1]); j < n; j += 2) { c0 = c1, c1 = this._edgecode(P[j], P[j + 1]); if (c0 && c1) j = this._edge(i, c0, c1, P, j), n = P.length; } } else if (this.contains(i, (this.xmin + this.xmax) / 2, (this.ymin + this.ymax) / 2)) { P = [this.xmin, this.ymin, this.xmax, this.ymin, this.xmax, this.ymax, this.xmin, this.ymax]; } return P; } _edge(i, e0, e1, P, j) { while (e0 !== e1) { let x, y; switch (e0) { case 0b0101: e0 = 0b0100; continue; // top-left case 0b0100: e0 = 0b0110, x = this.xmax, y = this.ymin; break; // top case 0b0110: e0 = 0b0010; continue; // top-right case 0b0010: e0 = 0b1010, x = this.xmax, y = this.ymax; break; // right case 0b1010: e0 = 0b1000; continue; // bottom-right case 0b1000: e0 = 0b1001, x = this.xmin, y = this.ymax; break; // bottom case 0b1001: e0 = 0b0001; continue; // bottom-left case 0b0001: e0 = 0b0101, x = this.xmin, y = this.ymin; break; // left } if ((P[j] !== x || P[j + 1] !== y) && this.contains(i, x, y)) { P.splice(j, 0, x, y), j += 2; } } return j; } _project(x0, y0, vx, vy) { let t = Infinity, c, x, y; if (vy < 0) { // top if (y0 <= this.ymin) return null; if ((c = (this.ymin - y0) / vy) < t) y = this.ymin, x = x0 + (t = c) * vx; } else if (vy > 0) { // bottom if (y0 >= this.ymax) return null; if ((c = (this.ymax - y0) / vy) < t) y = this.ymax, x = x0 + (t = c) * vx; } if (vx > 0) { // right if (x0 >= this.xmax) return null; if ((c = (this.xmax - x0) / vx) < t) x = this.xmax, y = y0 + (t = c) * vy; } else if (vx < 0) { // left if (x0 <= this.xmin) return null; if ((c = (this.xmin - x0) / vx) < t) x = this.xmin, y = y0 + (t = c) * vy; } return [x, y]; } _edgecode(x, y) { return (x === this.xmin ? 0b0001 : x === this.xmax ? 0b0010 : 0b0000) | (y === this.ymin ? 0b0100 : y === this.ymax ? 0b1000 : 0b0000); } _regioncode(x, y) { return (x < this.xmin ? 0b0001 : x > this.xmax ? 0b0010 : 0b0000) | (y < this.ymin ? 0b0100 : y > this.ymax ? 0b1000 : 0b0000); } }37:T1d68,Delaunay = { const tau = 2 * Math.PI; function pointX(p) { return p[0]; } function pointY(p) { return p[1]; } function area(hull, points) { let n = hull.length, x0, y0, x1 = points[2 * hull[n - 1]], y1 = points[2 * hull[n - 1] + 1], area = 0; for (let i = 0; i < n; i ++) { x0 = x1, y0 = y1; x1 = points[2 * hull[i]]; y1 = points[2 * hull[i] + 1]; area += y0 * x1 - x0 * y1; } return area / 2; } function jitter(x, y, r) { return [x + Math.sin(x + y) * r, y + Math.cos(x - y) * r]; } class Delaunay { constructor(points) { let d = new Delaunator(points), a = area(d.hull, points); if (a < 1e-10 && d.hull && d.hull.length > 2) { this.collinear = Int32Array.from({length: points.length/2}, (_,i) => i) .sort((i, j) => points[2 * i] - points[2 * j] || points[2 * i + 1] - points[2 * j + 1]); // for exact neighbors const e = this.collinear[0], f = this.collinear[this.collinear.length - 1], bounds = [ points[2 * e], points[2 * e + 1], points[2 * f], points[2 * f + 1] ], r = 1e-8 * Math.sqrt((bounds[3] - bounds[1])**2 + (bounds[2] - bounds[0])**2); for (let i = 0, n = points.length / 2; i < n; ++i) { const p = jitter(points[2 * i], points[2 * i + 1], r); points[2 * i] = p[0]; points[2 * i + 1] = p[1]; } d = new Delaunator(points); } this._delaunator = d; this.inedges = new Int32Array(points.length / 2); this._hullIndex = new Int32Array(points.length / 2); this.points = this._delaunator.coords; this._init(); } update() { this._delaunator.update(); this._init(); } _init() { const halfedges = this.halfedges = this._delaunator.halfedges; const hull = this.hull = this._delaunator.hull; const triangles = this.triangles = this._delaunator.triangles; const inedges = this.inedges.fill(-1); const hullIndex = this._hullIndex.fill(-1); // Compute an index from each point to an (arbitrary) incoming halfedge // Used to give the first neighbor of each point; for this reason, // on the hull we give priority to exterior halfedges for (let e = 0, n = halfedges.length; e < n; ++e) { const p = triangles[e % 3 === 2 ? e - 2 : e + 1]; if (halfedges[e] === -1 || inedges[p] === -1) inedges[p] = e; } for (let i = 0, n = hull.length; i < n; ++i) { hullIndex[hull[i]] = i; } // degenerate case: 1 or 2 (distinct) points if (hull.length <= 2 && hull.length > 0) { this.triangles = new Int32Array(3).fill(-1); this.halfedges = new Int32Array(3).fill(-1); this.triangles[0] = hull[0]; this.triangles[1] = hull[1]; this.triangles[2] = hull[1]; inedges[hull[0]] = 1; if (hull.length === 2) inedges[hull[1]] = 0; } } voronoi(bounds) { return new Voronoi(this, bounds); } *neighbors(i) { const {inedges, hull, _hullIndex, halfedges, triangles} = this; // degenerate case with several collinear points if (this.collinear) { const l = this.collinear.indexOf(i); if (l > 0) yield this.collinear[l - 1]; if (l < this.collinear.length - 1) yield this.collinear[l + 1]; return; } const e0 = inedges[i]; if (e0 === -1) return; // coincident point let e = e0, p0 = -1; do { yield p0 = triangles[e]; e = e % 3 === 2 ? e - 2 : e + 1; if (triangles[e] !== i) return; // bad triangulation e = halfedges[e]; if (e === -1) { const p = hull[(_hullIndex[i] + 1) % hull.length]; if (p !== p0) yield p; return; } } while (e !== e0); } find(x, y, i = 0) { if ((x = +x, x !== x) || (y = +y, y !== y)) return -1; const i0 = i; let c; while ((c = this._step(i, x, y)) >= 0 && c !== i && c !== i0) i = c; return c; } _step(i, x, y) { const {inedges, points} = this; if (inedges[i] === -1 || !points.length) return (i + 1) % (points.length >> 1); let c = i; let dc = (x - points[i * 2]) ** 2 + (y - points[i * 2 + 1]) ** 2; for (const t of this.neighbors(i)) { const dt = (x - points[t * 2]) ** 2 + (y - points[t * 2 + 1]) ** 2; if (dt < dc) dc = dt, c = t; } return c; } render(context) { const buffer = context == null ? context = new Path : undefined; const {points, halfedges, triangles} = this; for (let i = 0, n = halfedges.length; i < n; ++i) { const j = halfedges[i]; if (j < i) continue; const ti = triangles[i] * 2; const tj = triangles[j] * 2; context.moveTo(points[ti], points[ti + 1]); context.lineTo(points[tj], points[tj + 1]); } this.renderHull(context); return buffer && buffer.value(); } renderPoints(context, r = 2) { const buffer = context == null ? context = new Path : undefined; const {points} = this; for (let i = 0, n = points.length; i < n; i += 2) { const x = points[i], y = points[i + 1]; context.moveTo(x + r, y); context.arc(x, y, r, 0, tau); } return buffer && buffer.value(); } renderHull(context) { const buffer = context == null ? context = new Path : undefined; const {hull, points} = this; const h = hull[0] * 2, n = hull.length; context.moveTo(points[h], points[h + 1]); for (let i = 1; i < n; ++i) { const h = 2 * hull[i]; context.lineTo(points[h], points[h + 1]); } context.closePath(); return buffer && buffer.value(); } hullPolygon() { const polygon = new Polygon; this.renderHull(polygon); return polygon.value(); } renderTriangle(i, context) { const buffer = context == null ? context = new Path : undefined; const {points, triangles} = this; const t0 = triangles[i *= 3] * 2; const t1 = triangles[i + 1] * 2; const t2 = triangles[i + 2] * 2; context.moveTo(points[t0], points[t0 + 1]); context.lineTo(points[t1], points[t1 + 1]); context.lineTo(points[t2], points[t2 + 1]); context.closePath(); return buffer && buffer.value(); } *trianglePolygons() { const {triangles} = this; for (let i = 0, n = triangles.length / 3; i < n; ++i) { yield this.trianglePolygon(i); } } trianglePolygon(i) { const polygon = new Polygon; this.renderTriangle(i, polygon); return polygon.value(); } } Delaunay.from = function(points, fx = pointX, fy = pointY, that) { return new Delaunay("length" in points ? flatArray(points, fx, fy, that) : Float64Array.from(flatIterable(points, fx, fy, that))); }; function flatArray(points, fx, fy, that) { const n = points.length; const array = new Float64Array(n * 2); for (let i = 0; i < n; ++i) { const p = points[i]; array[i * 2] = fx.call(that, p, i, points); array[i * 2 + 1] = fy.call(that, p, i, points); } return array; } function* flatIterable(points, fx, fy, that) { let i = 0; for (const p of points) { yield fx.call(that, p, i, points); yield fy.call(that, p, i, points); ++i; } } return Delaunay; }38:T492,Path = { const epsilon = 1e-6; class Path { constructor() { this._x0 = this._y0 = // start of current subpath this._x1 = this._y1 = null; // end of current subpath this._ = ""; } moveTo(x, y) { this._ += `M${this._x0 = this._x1 = +x},${this._y0 = this._y1 = +y}`; } closePath() { if (this._x1 !== null) { this._x1 = this._x0, this._y1 = this._y0; this._ += "Z"; } } lineTo(x, y) { this._ += `L${this._x1 = +x},${this._y1 = +y}`; } arc(x, y, r) { x = +x, y = +y, r = +r; const x0 = x + r; const y0 = y; if (r < 0) throw new Error("negative radius"); if (this._x1 === null) this._ += `M${x0},${y0}`; else if (Math.abs(this._x1 - x0) > epsilon || Math.abs(this._y1 - y0) > epsilon) this._ += "L" + x0 + "," + y0; if (!r) return; this._ += `A${r},${r},0,1,1,${x - r},${y}A${r},${r},0,1,1,${this._x1 = x0},${this._y1 = y0}`; } rect(x, y, w, h) { this._ += `M${this._x0 = this._x1 = +x},${this._y0 = this._y1 = +y}h${+w}v${+h}h${-w}Z`; } value() { return this._ || null; } } return Path; } 39:[] 2f:{"title":"Lloyd’s Algorithm With Velocities","cells":[{"id":0,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":195,"value":"viewof totalPoints = slider({\n title: \"Number of generators\",\n min: 10,\n max: 20000,\n step: 1,\n value: 2000,\n submit: true,\n format: \",\",\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":559,"value":"viewof damping = slider({\n title: \"Velocity damping (0 = fully damped, 1 = undamped)\",\n min: 0,\n max: 1,\n precision: 4,\n value: 0.925,\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":682,"value":"viewof speedLimitFactor = slider({\n title: \"Max Speed Limit multiplier (up to 10x)\",\n min: 0,\n max: 10,\n precision: 4,\n value: 7.5,\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":697,"value":"viewof renderBlend = {\n await visibility();\n return checkbox({\n options: [{ value: \"blended\", label: \"Render overlay of both methods\" }],\n value: \"blended\"\n });\n}\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":84,"value":"viewof renderDefault = checkbox({\n options: [{ value: \"normal\", label: \"Render Normal Relaxation\" }],\n value: \"normal\"\n})\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":100,"value":"viewof renderVelocities = checkbox({\n options: [{ value: \"velocities\", label: \"Render Relaxation With Velocities\" }],\n value: \"velocities\"\n})\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":130,"value":"viewof clearView = checkbox({\n options: [{ value: \"clear\", label: \"Clear view each frame (uncheck to emphasize movement)\" }],\n value: \"clear\"\n})\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":28,"value":"viewof replay = html`