refine: convert any triangulation into constrained Delaunay

README.md

@@ -1,19 +1,19 @@
 # Earcut
 
-The fastest and smallest JavaScript **polygon triangulation** library. 3.5KB gzipped.
+The fastest and smallest JavaScript **polygon triangulation** library. 4KB gzipped.
+[![Node](https://github.com/mapbox/earcut/actions/workflows/node.yml/badge.svg)](https://github.com/mapbox/earcut/actions/workflows/node.yml)
+[![](https://img.shields.io/badge/simply-awesome-brightgreen.svg)](https://github.com/mourner/projects)
+
+![Earcut triangulation example](earcut.png)
 
-It is designed to be fast enough for real-time triangulation in the browser,
-sacrificing triangulation quality for raw speed and simplicity,
-while being robust enough to handle most practical datasets without crashing or producing garbage.
+Earcut is designed to be fast enough for real-time [triangulation](https://en.wikipedia.org/wiki/Polygon_triangulation) in the browser,
+favoring raw speed and simplicity over triangulation quality,
+while being robust enough to handle most practical datasets without crashing or producing garbage,
+with an option to refine the result to [Delaunay](https://en.wikipedia.org/wiki/Delaunay_triangulation) quality at a small cost.
 Originally built for [Mapbox GL JS](https://www.mapbox.com/mapbox-gljs) (WebGL-based interactive maps),
 it's now also used by [Three.js](https://threejs.org/) and many other projects.
 
-![Earcut triangulation example](earcut.png)
-
-[![Node](https://github.com/mapbox/earcut/actions/workflows/node.yml/badge.svg)](https://github.com/mapbox/earcut/actions/workflows/node.yml)
-[![Average time to resolve an issue](http://isitmaintained.com/badge/resolution/mapbox/earcut.svg)](http://isitmaintained.com/project/mapbox/earcut "Average time to resolve an issue")
-[![Percentage of issues still open](http://isitmaintained.com/badge/open/mapbox/earcut.svg)](http://isitmaintained.com/project/mapbox/earcut "Percentage of issues still open")
-[![](https://img.shields.io/badge/simply-awesome-brightgreen.svg)](https://github.com/mourner/projects)
+## [Demo](https://mapbox.github.io/earcut/viz/)
 
 ## Usage
 
@@ -39,7 +39,8 @@ and [Triangulation by Ear Clipping](http://www.geometrictools.com/Documentation/
 Earcut is heavily optimized for its primary workload — triangulating polygons from
 [Mapbox Vector Tiles](https://github.com/mapbox/vector-tile-spec). On a representative
 benchmark of **119,680 real-world polygons** (1.9M vertices) drawn from a window of map tiles
-through zooms 4–16, it triangulates the whole set in **~480 ms** on a Macbook Pro M1 Pro (2021).
+through zooms 4–16, it triangulates the whole set in **~445 ms** on a Macbook Pro M1 Pro (2021),
+with optional Delaunay refinement taking additional **~184 ms**.
 You can run the benchmark yourself with `npm run bench`.
 
 ## Robustness
@@ -62,7 +63,7 @@ conforming mesh, [remove T-junctions in a post-process](https://github.com/mapbo
 Install with NPM: `npm install earcut`, then import as a module:
 
 ```js
-import earcut, {flatten, deviation} from 'earcut';
+import earcut, {flatten, deviation, refine} from 'earcut';
 ```
 
 Or use as a module directly in the browser with [jsDelivr](https://www.jsdelivr.com/esm):
@@ -104,6 +105,23 @@ If you pass a single vertex as a hole, Earcut treats it as a Steiner point.
 
 Output triangles always have a consistent **winding order** regardless of the input polygon's winding — counter-clockwise in a y-up coordinate system (clockwise in y-down/screen space). If you need the opposite orientation (e.g. for back-face culling or normals in 3D), call `.reverse()` on the result.
 
+### `refine(triangles, vertices[, dimensions = 2])`
+
+If triangle quality matters, you can run an optional refinement pass after triangulation:
+
+```js
+const triangles = earcut(vertices, holes, dimensions);
+refine(triangles, vertices, dimensions);
+```
+
+This mutates `triangles` in place, legalizing interior edges with Delaunay flips while preserving
+the polygon boundary and holes. It keeps the same number of triangles and the same index format,
+but usually removes many skinny triangles and reduces total triangle edge length.
+
+Refinement is a post-process, so it doesn't affect the speed of normal `earcut` calls unless you
+explicitly call it. It assumes a valid manifold triangulation, such as the output of `earcut`, and
+doesn't repair invalid polygon input or make the mesh conforming.
+
 ### `flatten(data)`
 
 If your input is a multi-dimensional array (e.g. [GeoJSON Polygon](http://geojson.org/geojson-spec.html#polygon)),

bench/bench-refine.js

@@ -0,0 +1,75 @@
+// Benchmark the optional Delaunay refinement post-pass (refine() in src/earcut.js)
+// over the same realistic MVT fixture as bench-tiles.js. Measures the added cost of
+// refine() relative to earcut alone, and reports the resulting triangle quality.
+import earcut, {refine} from '../src/earcut.js';
+import {readTilesFixture} from './tiles-fixture.js';
+
+const polys = readTilesFixture();
+const totalVerts = polys.reduce((s, d) => s + d.vertices.length / d.dimensions, 0);
+
+// pre-triangulate once per poly so the timed paths only measure their own work;
+// refine mutates the index array, so each timed pass re-triangulates a fresh copy.
+function triOnly(d)    { return earcut(d.vertices, d.holes, d.dimensions).length; }
+function triRefine(d)  { const t = earcut(d.vertices, d.holes, d.dimensions); refine(t, d.vertices, d.dimensions); return t.length; }
+
+function run(fn) { let s = 0; for (const d of polys) s += fn(d); return s; }
+
+// interleaved A/B (ABAB): cancels thermal drift instead of attributing it to one side.
+function timeAB(a, b) {
+    let ba = Infinity, bb = Infinity;
+    for (let i = 0; i < 7; i++) {
+        let s = performance.now(); run(a); ba = Math.min(ba, performance.now() - s);
+        s = performance.now(); run(b); bb = Math.min(bb, performance.now() - s);
+    }
+    return {a: ba, b: bb};
+}
+
+// --- triangle quality (normalized q = 4√3·area / Σ edge², in [0,1]) ---
+function quality(getIndices) {
+    let slivers = 0, tinyAngle = 0, perim = 0, sumQ = 0, tris = 0;
+    for (const d of polys) {
+        const c = d.vertices, dim = d.dimensions;
+        const idx = getIndices(d);
+        for (let i = 0; i < idx.length; i += 3) {
+            const ax = c[idx[i] * dim], ay = c[idx[i] * dim + 1];
+            const bx = c[idx[i + 1] * dim], by = c[idx[i + 1] * dim + 1];
+            const cx = c[idx[i + 2] * dim], cy = c[idx[i + 2] * dim + 1];
+            const area = Math.abs((bx - ax) * (cy - ay) - (by - ay) * (cx - ax)) / 2;
+            const e1 = (bx - ax) ** 2 + (by - ay) ** 2;
+            const e2 = (cx - bx) ** 2 + (cy - by) ** 2;
+            const e3 = (ax - cx) ** 2 + (ay - cy) ** 2;
+            const sumSq = e1 + e2 + e3;
+            const q = sumSq > 0 ? 4 * Math.sqrt(3) * area / sumSq : 0;
+            if (q < 0.1) slivers++;
+            // smallest angle via law of cosines on the shortest-opposite side
+            const la = Math.sqrt(e1), lb = Math.sqrt(e2), lc = Math.sqrt(e3);
+            const angA = Math.acos(Math.min(1, Math.max(-1, (lb * lb + lc * lc - la * la) / (2 * lb * lc || 1))));
+            const angB = Math.acos(Math.min(1, Math.max(-1, (la * la + lc * lc - lb * lb) / (2 * la * lc || 1))));
+            const minAng = Math.min(angA, angB, Math.PI - angA - angB) * 180 / Math.PI;
+            if (minAng < 1) tinyAngle++;
+            perim += la + lb + lc;
+            sumQ += q; tris++;
+        }
+    }
+    return {slivers, tinyAngle, perim, meanQ: sumQ / tris, tris};
+}
+
+// warm both paths
+run(triOnly); run(triRefine);
+
+const {a: tEarcut, b: tRefine} = timeAB(triOnly, triRefine);
+const overhead = (tRefine - tEarcut) / tEarcut * 100;
+
+console.log(`fixture:  ${polys.length.toLocaleString()} polygons, ${totalVerts.toLocaleString()} vertices`);
+console.log(`\nearcut:          ${tEarcut.toFixed(1)} ms`);
+console.log(`earcut+refine:   ${tRefine.toFixed(1)} ms   (+${overhead.toFixed(0)}%, ${(tRefine - tEarcut).toFixed(1)} ms refine)`);
+
+const qa = quality(d => earcut(d.vertices, d.holes, d.dimensions));
+const qb = quality((d) => { const t = earcut(d.vertices, d.holes, d.dimensions); refine(t, d.vertices, d.dimensions); return t; });
+const pct = (x, y) => `${((y - x) / x * 100).toFixed(0)}%`;
+console.log('\nquality:              earcut     refined    delta');
+console.log(`  slivers (q<0.1)   ${String(qa.slivers).padStart(8)}  ${String(qb.slivers).padStart(8)}    ${pct(qa.slivers, qb.slivers)}`);
+console.log(`  min-angle < 1deg  ${String(qa.tinyAngle).padStart(8)}  ${String(qb.tinyAngle).padStart(8)}    ${pct(qa.tinyAngle, qb.tinyAngle)}`);
+console.log(`  total perimeter   ${qa.perim.toExponential(2).padStart(8)}  ${qb.perim.toExponential(2).padStart(8)}    ${pct(qa.perim, qb.perim)}`);
+console.log(`  mean q            ${qa.meanQ.toFixed(3).padStart(8)}  ${qb.meanQ.toFixed(3).padStart(8)}    +${pct(qa.meanQ, qb.meanQ)}`);
+console.log(`  triangles         ${qa.tris.toLocaleString()} (invariant: ${qa.tris === qb.tris})`);

package.json

@@ -1,9 +1,10 @@
 {
   "name": "earcut",
-  "version": "3.1.1",
+  "version": "3.2.0",
   "description": "The fastest and smallest JavaScript polygon triangulation library for your WebGL apps",
   "main": "src/earcut.js",
   "type": "module",
+  "sideEffects": false,
   "exports": "./src/earcut.js",
   "types": "src/earcut.d.ts",
   "files": [

src/earcut.js

@@ -873,3 +873,124 @@ export function flatten(data) {
     }
     return {vertices, holes, dimensions};
 }
+
+// Reusable module-level scratch for refine():
+//   he      = twin half-edge of each edge, or -1 on the polygon boundary
+//   hTable  = open-addressing hash, slot -> half-edge index, valid iff hStamp[slot] === gen
+/** @type {Int32Array} */ let edgeStack;
+/** @type {Int32Array} */ let he;
+/** @type {Int32Array} */ let hTable;
+/** @type {Uint32Array} */ let hStamp;
+let hMask = 0, gen = 0;
+
+/**
+ * Refine a triangulation toward the constrained Delaunay triangulation by legalizing every
+ * interior edge in place with Lawson flips — maximizing the minimum angle and removing most
+ * slivers. An optional post-pass for {@link earcut} output, or any manifold triangle-index array
+ * indexing into `coords`. Adapted from delaunator's edge legalization.
+ *
+ * Uses non-robust predicates: float input is fine, and the worst case is a not-quite-Delaunay
+ * edge, never an invalid mesh.
+ *
+ * @param {number[]} triangles triangle indices, as returned by {@link earcut}; mutated in place
+ * @param {ArrayLike<number>} coords the flat vertex coordinates passed to {@link earcut}
+ * @param {number} [dim=2] number of coordinates per vertex in `coords`
+ * @example refine(earcut(data), data);
+ */
+export function refine(triangles, coords, dim = 2) {
+    const t = triangles;
+    const n = t.length;
+    const V = coords.length / dim;
+    ensureScratch(n);
+    gen++;              // bumping the generation logically empties the hash (no clearing)
+    he.fill(-1, 0, n);
+
+    // Build half-edge twins with an undirected-edge hash; consumed slots mark linked pairs.
+    for (let e = 0; e < n; e++) {
+        const a = t[e], b = t[nextHE(e)];
+        const lo = a < b ? a : b, hi = a < b ? b : a;
+        let h = ((hi * V + lo) * 0x9e3779b1 >>> 0) & hMask;
+        while (hStamp[h] === gen) {
+            const s = hTable[h];
+            // s === -1 marks a consumed slot (a pair already linked) — skip past it
+            if (s !== -1) {
+                const sa = t[s], sb = t[nextHE(s)];
+                if ((sa === lo && sb === hi) || (sa === hi && sb === lo)) {
+                    he[e] = s; he[s] = e; hTable[h] = -1; // link, then consume the slot
+                    break;
+                }
+            }
+            h = (h + 1) & hMask;
+        }
+        if (hStamp[h] !== gen) { hTable[h] = e; hStamp[h] = gen; } // first occurrence: insert
+    }
+
+    for (let e0 = 0; e0 < n; e0++) {
+        if (he[e0] === -1) continue;
+        let i = 0, a = e0;
+        while (true) {
+            const b = he[a];
+            const a0 = a - a % 3;
+            const ar = a0 + (a + 2) % 3;
+            if (b === -1) { if (i === 0) break; a = edgeStack[--i]; continue; }
+
+            const b0 = b - b % 3;
+            const al = a0 + (a + 1) % 3;
+            const bl = b0 + (b + 2) % 3;
+            const p0 = t[ar], pr = t[a], pl = t[al], p1 = t[bl];
+
+            const x0 = coords[p0 * dim], y0 = coords[p0 * dim + 1];
+            const xr = coords[pr * dim], yr = coords[pr * dim + 1];
+            const xl = coords[pl * dim], yl = coords[pl * dim + 1];
+            const x1 = coords[p1 * dim], y1 = coords[p1 * dim + 1];
+
+            // Both triangles of the flipped diagonal p0-p1 must be CCW (i.e. the quad is convex).
+            // Flipping a reflex quad would push a triangle outside the polygon; this guards against
+            // it. Boundary/hole edges need no guard — they self-protect via he === -1.
+            const convex = orient(x0, y0, xr, yr, x1, y1) > 0 && orient(x0, y0, x1, y1, xl, yl) > 0;
+
+            if (convex && !inCircle(x0, y0, xr, yr, xl, yl, x1, y1)) {
+                t[a] = p1; t[b] = p0;
+                const hbl = he[bl], har = he[ar];
+                he[a] = hbl; if (hbl !== -1) he[hbl] = a;
+                he[b] = har; if (har !== -1) he[har] = b;
+                he[ar] = bl; he[bl] = ar;
+                if (i < edgeStack.length) edgeStack[i++] = b0 + (b + 1) % 3;
+            } else {
+                if (i === 0) break;
+                a = edgeStack[--i];
+            }
+        }
+    }
+}
+
+/** @param {number} ax @param {number} ay @param {number} bx @param {number} by @param {number} cx @param {number} cy */
+function orient(ax, ay, bx, by, cx, cy) {
+    return (bx - ax) * (cy - ay) - (by - ay) * (cx - ax);
+}
+
+// Whether p is inside the circumcircle of triangle (a, b, c). Sign is negated vs the usual
+// predicate to match earcut's CCW winding — the standard sign would build the anti-Delaunay mesh.
+/** @param {number} ax @param {number} ay @param {number} bx @param {number} by @param {number} cx @param {number} cy @param {number} px @param {number} py */
+function inCircle(ax, ay, bx, by, cx, cy, px, py) {
+    const dx = ax - px, dy = ay - py, ex = bx - px, ey = by - py, fx = cx - px, fy = cy - py;
+    const ap = dx * dx + dy * dy, bp = ex * ex + ey * ey, cp = fx * fx + fy * fy;
+    return dx * (ey * cp - bp * fy) - dy * (ex * cp - bp * fx) + ap * (ex * fy - ey * fx) < 0;
+}
+
+/** @param {number} e */
+function nextHE(e) { // next half-edge within the same triangle
+    return e - e % 3 + (e + 1) % 3;
+}
+
+// Grow the scratch arrays on demand (like earcut's z-order arrays). Allocating lazily here rather
+// than at module load lets the whole refine() block tree-shake away for callers who don't use it.
+/** @param {number} n */
+function ensureScratch(n) {
+    if (!edgeStack) edgeStack = new Int32Array(512);
+    if (!he || he.length < n) he = new Int32Array(n);
+    let size = 1;
+    while (size < n * 4) size <<= 1; // power-of-two table, load factor <= 0.25
+    if (!hTable || hTable.length < size) { hTable = new Int32Array(size); hStamp = new Uint32Array(size); }
+    hMask = size - 1;
+}

test/expected.json

@@ -1,5 +1,6 @@
 {
     "triangles": {
+        "earcut": 548,
         "building": 13,
         "dude": 106,
         "water": 2482,