basant307/AI_Governance_Project

1 2const EPSILON = Math.pow(2, -52);3const EDGE_STACK = new Uint32Array(512);4 5import {orient2d} from 'robust-predicates';6 7/** @template {ArrayLike<number>} T */8export default class Delaunator {9 10    /**11     * Constructs a delaunay triangulation object given an array of points (`[x, y]` by default).12     * `getX` and `getY` are optional functions of the form `(point) => value` for custom point formats.13     *14     * @template P15     * @param {P[]} points16     * @param {(p: P) => number} [getX]17     * @param {(p: P) => number} [getY]18     */19    // @ts-expect-error TS232220    static from(points, getX = defaultGetX, getY = defaultGetY) {21        const n = points.length;22        const coords = new Float64Array(n * 2);23 24        for (let i = 0; i < n; i++) {25            const p = points[i];26            coords[2 * i] = getX(p);27            coords[2 * i + 1] = getY(p);28        }29 30        return new Delaunator(coords);31    }32 33    /**34     * Constructs a delaunay triangulation object given an array of point coordinates of the form:35     * `[x0, y0, x1, y1, ...]` (use a typed array for best performance). Duplicate points are skipped.36     *37     * @param {T} coords38     */39    constructor(coords) {40        const n = coords.length >> 1;41        if (n > 0 && typeof coords[0] !== 'number') throw new Error('Expected coords to contain numbers.');42 43        this.coords = coords;44 45        // arrays that will store the triangulation graph46        const maxTriangles = Math.max(2 * n - 5, 0);47        /** @private */ this._triangles = new Uint32Array(maxTriangles * 3);48        /** @private */ this._halfedges = new Int32Array(maxTriangles * 3);49 50        // temporary arrays for tracking the edges of the advancing convex hull51        /** @private */ this._hashSize = Math.ceil(Math.sqrt(n));52        /** @private */ this._hullPrev = new Uint32Array(n); // edge to prev edge53        /** @private */ this._hullNext = new Uint32Array(n); // edge to next edge54        /** @private */ this._hullTri = new Uint32Array(n); // edge to adjacent triangle55        /** @private */ this._hullHash = new Int32Array(this._hashSize); // angular edge hash56 57        // temporary arrays for sorting points58        /** @private */ this._ids = new Uint32Array(n);59        /** @private */ this._dists = new Float64Array(n);60 61        /** @private */ this.trianglesLen = 0;62        /** @private */ this._cx = 0;63        /** @private */ this._cy = 0;64        /** @private */ this._hullStart = 0;65 66 67        /** A `Uint32Array` array of indices that reference points on the convex hull of the input data, counter-clockwise. */68        this.hull = this._triangles;69        /** A `Uint32Array` array of triangle vertex indices (each group of three numbers forms a triangle). All triangles are directed counterclockwise. */70        this.triangles = this._triangles;71        /**72         * A `Int32Array` array of triangle half-edge indices that allows you to traverse the triangulation.73         * `i`-th half-edge in the array corresponds to vertex `triangles[i]` the half-edge is coming from.74         * `halfedges[i]` is the index of a twin half-edge in an adjacent triangle (or `-1` for outer half-edges on the convex hull).75         */76        this.halfedges = this._halfedges;77 78        this.update();79    }80 81    /**82     * Updates the triangulation if you modified `delaunay.coords` values in place, avoiding expensive memory allocations.83     * Useful for iterative relaxation algorithms such as Lloyd's.84     */85    update() {86        const {coords, _hullPrev: hullPrev, _hullNext: hullNext, _hullTri: hullTri, _hullHash: hullHash} =  this;87        const n = coords.length >> 1;88 89        // populate an array of point indices; calculate input data bbox90        let minX = Infinity;91        let minY = Infinity;92        let maxX = -Infinity;93        let maxY = -Infinity;94 95        for (let i = 0; i < n; i++) {96            const x = coords[2 * i];97            const y = coords[2 * i + 1];98            if (x < minX) minX = x;99            if (y < minY) minY = y;100            if (x > maxX) maxX = x;101            if (y > maxY) maxY = y;102            this._ids[i] = i;103        }104        const cx = (minX + maxX) / 2;105        const cy = (minY + maxY) / 2;106 107        let i0 = 0, i1 = 0, i2 = 0;108 109        // pick a seed point close to the center110        for (let i = 0, minDist = Infinity; i < n; i++) {111            const d = dist(cx, cy, coords[2 * i], coords[2 * i + 1]);112            if (d < minDist) {113                i0 = i;114                minDist = d;115            }116        }117        const i0x = coords[2 * i0];118        const i0y = coords[2 * i0 + 1];119 120        // find the point closest to the seed121        for (let i = 0, minDist = Infinity; i < n; i++) {122            if (i === i0) continue;123            const d = dist(i0x, i0y, coords[2 * i], coords[2 * i + 1]);124            if (d < minDist && d > 0) {125                i1 = i;126                minDist = d;127            }128        }129        let i1x = coords[2 * i1];130        let i1y = coords[2 * i1 + 1];131 132        let minRadius = Infinity;133 134        // find the third point which forms the smallest circumcircle with the first two135        for (let i = 0; i < n; i++) {136            if (i === i0 || i === i1) continue;137            const r = circumradius(i0x, i0y, i1x, i1y, coords[2 * i], coords[2 * i + 1]);138            if (r < minRadius) {139                i2 = i;140                minRadius = r;141            }142        }143        let i2x = coords[2 * i2];144        let i2y = coords[2 * i2 + 1];145 146        if (minRadius === Infinity) {147            // order collinear points by dx (or dy if all x are identical)148            // and return the list as a hull149            for (let i = 0; i < n; i++) {150                this._dists[i] = (coords[2 * i] - coords[0]) || (coords[2 * i + 1] - coords[1]);151            }152            quicksort(this._ids, this._dists, 0, n - 1);153            const hull = new Uint32Array(n);154            let j = 0;155            for (let i = 0, d0 = -Infinity; i < n; i++) {156                const id = this._ids[i];157                const d = this._dists[id];158                if (d > d0) {159                    hull[j++] = id;160                    d0 = d;161                }162            }163            this.hull = hull.subarray(0, j);164            this.triangles = new Uint32Array(0);165            this.halfedges = new Int32Array(0);166            return;167        }168 169        // swap the order of the seed points for counter-clockwise orientation170        if (orient2d(i0x, i0y, i1x, i1y, i2x, i2y) < 0) {171            const i = i1;172            const x = i1x;173            const y = i1y;174            i1 = i2;175            i1x = i2x;176            i1y = i2y;177            i2 = i;178            i2x = x;179            i2y = y;180        }181 182        const center = circumcenter(i0x, i0y, i1x, i1y, i2x, i2y);183        this._cx = center.x;184        this._cy = center.y;185 186        for (let i = 0; i < n; i++) {187            this._dists[i] = dist(coords[2 * i], coords[2 * i + 1], center.x, center.y);188        }189 190        // sort the points by distance from the seed triangle circumcenter191        quicksort(this._ids, this._dists, 0, n - 1);192 193        // set up the seed triangle as the starting hull194        this._hullStart = i0;195        let hullSize = 3;196 197        hullNext[i0] = hullPrev[i2] = i1;198        hullNext[i1] = hullPrev[i0] = i2;199        hullNext[i2] = hullPrev[i1] = i0;200 201        hullTri[i0] = 0;202        hullTri[i1] = 1;203        hullTri[i2] = 2;204 205        hullHash.fill(-1);206        hullHash[this._hashKey(i0x, i0y)] = i0;207        hullHash[this._hashKey(i1x, i1y)] = i1;208        hullHash[this._hashKey(i2x, i2y)] = i2;209 210        this.trianglesLen = 0;211        this._addTriangle(i0, i1, i2, -1, -1, -1);212 213        for (let k = 0, xp = 0, yp = 0; k < this._ids.length; k++) {214            const i = this._ids[k];215            const x = coords[2 * i];216            const y = coords[2 * i + 1];217 218            // skip near-duplicate points219            if (k > 0 && Math.abs(x - xp) <= EPSILON && Math.abs(y - yp) <= EPSILON) continue;220            xp = x;221            yp = y;222 223            // skip seed triangle points224            if (i === i0 || i === i1 || i === i2) continue;225 226            // find a visible edge on the convex hull using edge hash227            let start = 0;228            for (let j = 0, key = this._hashKey(x, y); j < this._hashSize; j++) {229                start = hullHash[(key + j) % this._hashSize];230                if (start !== -1 && start !== hullNext[start]) break;231            }232 233            start = hullPrev[start];234            let e = start, q;235            while (q = hullNext[e], orient2d(x, y, coords[2 * e], coords[2 * e + 1], coords[2 * q], coords[2 * q + 1]) >= 0) {236                e = q;237                if (e === start) {238                    e = -1;239                    break;240                }241            }242            if (e === -1) continue; // likely a near-duplicate point; skip it243 244            // add the first triangle from the point245            let t = this._addTriangle(e, i, hullNext[e], -1, -1, hullTri[e]);246 247            // recursively flip triangles from the point until they satisfy the Delaunay condition248            hullTri[i] = this._legalize(t + 2);249            hullTri[e] = t; // keep track of boundary triangles on the hull250            hullSize++;251 252            // walk forward through the hull, adding more triangles and flipping recursively253            let n = hullNext[e];254            while (q = hullNext[n], orient2d(x, y, coords[2 * n], coords[2 * n + 1], coords[2 * q], coords[2 * q + 1]) < 0) {255                t = this._addTriangle(n, i, q, hullTri[i], -1, hullTri[n]);256                hullTri[i] = this._legalize(t + 2);257                hullNext[n] = n; // mark as removed258                hullSize--;259                n = q;260            }261 262            // walk backward from the other side, adding more triangles and flipping263            if (e === start) {264                while (q = hullPrev[e], orient2d(x, y, coords[2 * q], coords[2 * q + 1], coords[2 * e], coords[2 * e + 1]) < 0) {265                    t = this._addTriangle(q, i, e, -1, hullTri[e], hullTri[q]);266                    this._legalize(t + 2);267                    hullTri[q] = t;268                    hullNext[e] = e; // mark as removed269                    hullSize--;270                    e = q;271                }272            }273 274            // update the hull indices275            this._hullStart = hullPrev[i] = e;276            hullNext[e] = hullPrev[n] = i;277            hullNext[i] = n;278 279            // save the two new edges in the hash table280            hullHash[this._hashKey(x, y)] = i;281            hullHash[this._hashKey(coords[2 * e], coords[2 * e + 1])] = e;282        }283 284        this.hull = new Uint32Array(hullSize);285        for (let i = 0, e = this._hullStart; i < hullSize; i++) {286            this.hull[i] = e;287            e = hullNext[e];288        }289 290        // trim typed triangle mesh arrays291        this.triangles = this._triangles.subarray(0, this.trianglesLen);292        this.halfedges = this._halfedges.subarray(0, this.trianglesLen);293    }294 295    /**296     * Calculate an angle-based key for the edge hash used for advancing convex hull.297     *298     * @param {number} x299     * @param {number} y300     * @private301     */302    _hashKey(x, y) {303        return Math.floor(pseudoAngle(x - this._cx, y - this._cy) * this._hashSize) % this._hashSize;304    }305 306    /**307     * Flip an edge in a pair of triangles if it doesn't satisfy the Delaunay condition.308     *309     * @param {number} a310     * @private311     */312    _legalize(a) {313        const {_triangles: triangles, _halfedges: halfedges, coords} = this;314 315        let i = 0;316        let ar = 0;317 318        // recursion eliminated with a fixed-size stack319        while (true) {320            const b = halfedges[a];321 322            /* if the pair of triangles doesn't satisfy the Delaunay condition323             * (p1 is inside the circumcircle of [p0, pl, pr]), flip them,324             * then do the same check/flip recursively for the new pair of triangles325             *326             *           pl                    pl327             *          /||\                  /  \328             *       al/ || \bl            al/    \a329             *        /  ||  \              /      \330             *       /  a||b  \    flip    /___ar___\331             *     p0\   ||   /p1   =>   p0\---bl---/p1332             *        \  ||  /              \      /333             *       ar\ || /br             b\    /br334             *          \||/                  \  /335             *           pr                    pr336             */337            const a0 = a - a % 3;338            ar = a0 + (a + 2) % 3;339 340            if (b === -1) { // convex hull edge341                if (i === 0) break;342                a = EDGE_STACK[--i];343                continue;344            }345 346            const b0 = b - b % 3;347            const al = a0 + (a + 1) % 3;348            const bl = b0 + (b + 2) % 3;349 350            const p0 = triangles[ar];351            const pr = triangles[a];352            const pl = triangles[al];353            const p1 = triangles[bl];354 355            const illegal = inCircle(356                coords[2 * p0], coords[2 * p0 + 1],357                coords[2 * pr], coords[2 * pr + 1],358                coords[2 * pl], coords[2 * pl + 1],359                coords[2 * p1], coords[2 * p1 + 1]);360 361            if (illegal) {362                triangles[a] = p1;363                triangles[b] = p0;364 365                const hbl = halfedges[bl];366 367                // edge swapped on the other side of the hull (rare); fix the half-edge reference368                if (hbl === -1) {369                    let e = this._hullStart;370                    do {371                        if (this._hullTri[e] === bl) {372                            this._hullTri[e] = a;373                            break;374                        }375                        e = this._hullPrev[e];376                    } while (e !== this._hullStart);377                }378                this._link(a, hbl);379                this._link(b, halfedges[ar]);380                this._link(ar, bl);381 382                const br = b0 + (b + 1) % 3;383 384                // don't worry about hitting the cap: it can only happen on extremely degenerate input385                if (i < EDGE_STACK.length) {386                    EDGE_STACK[i++] = br;387                }388            } else {389                if (i === 0) break;390                a = EDGE_STACK[--i];391            }392        }393 394        return ar;395    }396 397    /**398     * Link two half-edges to each other.399     * @param {number} a400     * @param {number} b401     * @private402     */403    _link(a, b) {404        this._halfedges[a] = b;405        if (b !== -1) this._halfedges[b] = a;406    }407 408    /**409     * Add a new triangle given vertex indices and adjacent half-edge ids.410     *411     * @param {number} i0412     * @param {number} i1413     * @param {number} i2414     * @param {number} a415     * @param {number} b416     * @param {number} c417     * @private418     */419    _addTriangle(i0, i1, i2, a, b, c) {420        const t = this.trianglesLen;421 422        this._triangles[t] = i0;423        this._triangles[t + 1] = i1;424        this._triangles[t + 2] = i2;425 426        this._link(t, a);427        this._link(t + 1, b);428        this._link(t + 2, c);429 430        this.trianglesLen += 3;431 432        return t;433    }434}435 436/**437 * Monotonically increases with real angle, but doesn't need expensive trigonometry.438 *439 * @param {number} dx440 * @param {number} dy441 */442function pseudoAngle(dx, dy) {443    const p = dx / (Math.abs(dx) + Math.abs(dy));444    return (dy > 0 ? 3 - p : 1 + p) / 4; // [0..1]445}446 447/**448 * Squared distance between two points.449 *450 * @param {number} ax451 * @param {number} ay452 * @param {number} bx453 * @param {number} by454 */455function dist(ax, ay, bx, by) {456    const dx = ax - bx;457    const dy = ay - by;458    return dx * dx + dy * dy;459}460 461/**462 * Check whether point P is inside a circle formed by points A, B, C.463 *464 * @param {number} ax465 * @param {number} ay466 * @param {number} bx467 * @param {number} by468 * @param {number} cx469 * @param {number} cy470 * @param {number} px471 * @param {number} py472 */473function inCircle(ax, ay, bx, by, cx, cy, px, py) {474    const dx = ax - px;475    const dy = ay - py;476    const ex = bx - px;477    const ey = by - py;478    const fx = cx - px;479    const fy = cy - py;480 481    const ap = dx * dx + dy * dy;482    const bp = ex * ex + ey * ey;483    const cp = fx * fx + fy * fy;484 485    return dx * (ey * cp - bp * fy) -486           dy * (ex * cp - bp * fx) +487           ap * (ex * fy - ey * fx) < 0;488}489 490/**491 * Squared radius of the circle formed by points A, B, C.492 *493 * @param {number} ax494 * @param {number} ay495 * @param {number} bx496 * @param {number} by497 * @param {number} cx498 * @param {number} cy499 */500function circumradius(ax, ay, bx, by, cx, cy) {501    const dx = bx - ax;502    const dy = by - ay;503    const ex = cx - ax;504    const ey = cy - ay;505 506    const bl = dx * dx + dy * dy;507    const cl = ex * ex + ey * ey;508    const d = 0.5 / (dx * ey - dy * ex);509 510    const x = (ey * bl - dy * cl) * d;511    const y = (dx * cl - ex * bl) * d;512 513    return x * x + y * y;514}515 516/**517 * Get coordinates of a circumcenter for points A, B, C.518 *519 * @param {number} ax520 * @param {number} ay521 * @param {number} bx522 * @param {number} by523 * @param {number} cx524 * @param {number} cy525 */526function circumcenter(ax, ay, bx, by, cx, cy) {527    const dx = bx - ax;528    const dy = by - ay;529    const ex = cx - ax;530    const ey = cy - ay;531 532    const bl = dx * dx + dy * dy;533    const cl = ex * ex + ey * ey;534    const d = 0.5 / (dx * ey - dy * ex);535 536    const x = ax + (ey * bl - dy * cl) * d;537    const y = ay + (dx * cl - ex * bl) * d;538 539    return {x, y};540}541 542/**543 * Sort points by distance via an array of point indices and an array of calculated distances.544 *545 * @param {Uint32Array} ids546 * @param {Float64Array} dists547 * @param {number} left548 * @param {number} right549 */550function quicksort(ids, dists, left, right) {551    if (right - left <= 20) {552        for (let i = left + 1; i <= right; i++) {553            const temp = ids[i];554            const tempDist = dists[temp];555            let j = i - 1;556            while (j >= left && dists[ids[j]] > tempDist) ids[j + 1] = ids[j--];557            ids[j + 1] = temp;558        }559    } else {560        const median = (left + right) >> 1;561        let i = left + 1;562        let j = right;563        swap(ids, median, i);564        if (dists[ids[left]] > dists[ids[right]]) swap(ids, left, right);565        if (dists[ids[i]] > dists[ids[right]]) swap(ids, i, right);566        if (dists[ids[left]] > dists[ids[i]]) swap(ids, left, i);567 568        const temp = ids[i];569        const tempDist = dists[temp];570        while (true) {571            do i++; while (dists[ids[i]] < tempDist);572            do j--; while (dists[ids[j]] > tempDist);573            if (j < i) break;574            swap(ids, i, j);575        }576        ids[left + 1] = ids[j];577        ids[j] = temp;578 579        if (right - i + 1 >= j - left) {580            quicksort(ids, dists, i, right);581            quicksort(ids, dists, left, j - 1);582        } else {583            quicksort(ids, dists, left, j - 1);584            quicksort(ids, dists, i, right);585        }586    }587}588 589/**590 * @param {Uint32Array} arr591 * @param {number} i592 * @param {number} j593 */594function swap(arr, i, j) {595    const tmp = arr[i];596    arr[i] = arr[j];597    arr[j] = tmp;598}599 600/** @param {[number, number]} p */601function defaultGetX(p) {602    return p[0];603}604/** @param {[number, number]} p */605function defaultGetY(p) {606    return p[1];607}608