| [e4c61dd] | 1 | import Path from "./path.js";
|
|---|
| 2 | import Polygon from "./polygon.js";
|
|---|
| 3 |
|
|---|
| 4 | export default class Voronoi {
|
|---|
| 5 | constructor(delaunay, [xmin, ymin, xmax, ymax] = [0, 0, 960, 500]) {
|
|---|
| 6 | if (!((xmax = +xmax) >= (xmin = +xmin)) || !((ymax = +ymax) >= (ymin = +ymin))) throw new Error("invalid bounds");
|
|---|
| 7 | this.delaunay = delaunay;
|
|---|
| 8 | this._circumcenters = new Float64Array(delaunay.points.length * 2);
|
|---|
| 9 | this.vectors = new Float64Array(delaunay.points.length * 2);
|
|---|
| 10 | this.xmax = xmax, this.xmin = xmin;
|
|---|
| 11 | this.ymax = ymax, this.ymin = ymin;
|
|---|
| 12 | this._init();
|
|---|
| 13 | }
|
|---|
| 14 | update() {
|
|---|
| 15 | this.delaunay.update();
|
|---|
| 16 | this._init();
|
|---|
| 17 | return this;
|
|---|
| 18 | }
|
|---|
| 19 | _init() {
|
|---|
| 20 | const {delaunay: {points, hull, triangles}, vectors} = this;
|
|---|
| 21 | let bx, by; // lazily computed barycenter of the hull
|
|---|
| 22 |
|
|---|
| 23 | // Compute circumcenters.
|
|---|
| 24 | const circumcenters = this.circumcenters = this._circumcenters.subarray(0, triangles.length / 3 * 2);
|
|---|
| 25 | for (let i = 0, j = 0, n = triangles.length, x, y; i < n; i += 3, j += 2) {
|
|---|
| 26 | const t1 = triangles[i] * 2;
|
|---|
| 27 | const t2 = triangles[i + 1] * 2;
|
|---|
| 28 | const t3 = triangles[i + 2] * 2;
|
|---|
| 29 | const x1 = points[t1];
|
|---|
| 30 | const y1 = points[t1 + 1];
|
|---|
| 31 | const x2 = points[t2];
|
|---|
| 32 | const y2 = points[t2 + 1];
|
|---|
| 33 | const x3 = points[t3];
|
|---|
| 34 | const y3 = points[t3 + 1];
|
|---|
| 35 |
|
|---|
| 36 | const dx = x2 - x1;
|
|---|
| 37 | const dy = y2 - y1;
|
|---|
| 38 | const ex = x3 - x1;
|
|---|
| 39 | const ey = y3 - y1;
|
|---|
| 40 | const ab = (dx * ey - dy * ex) * 2;
|
|---|
| 41 |
|
|---|
| 42 | if (Math.abs(ab) < 1e-9) {
|
|---|
| 43 | // For a degenerate triangle, the circumcenter is at the infinity, in a
|
|---|
| 44 | // direction orthogonal to the halfedge and away from the “center” of
|
|---|
| 45 | // the diagram <bx, by>, defined as the hull’s barycenter.
|
|---|
| 46 | if (bx === undefined) {
|
|---|
| 47 | bx = by = 0;
|
|---|
| 48 | for (const i of hull) bx += points[i * 2], by += points[i * 2 + 1];
|
|---|
| 49 | bx /= hull.length, by /= hull.length;
|
|---|
| 50 | }
|
|---|
| 51 | const a = 1e9 * Math.sign((bx - x1) * ey - (by - y1) * ex);
|
|---|
| 52 | x = (x1 + x3) / 2 - a * ey;
|
|---|
| 53 | y = (y1 + y3) / 2 + a * ex;
|
|---|
| 54 | } else {
|
|---|
| 55 | const d = 1 / ab;
|
|---|
| 56 | const bl = dx * dx + dy * dy;
|
|---|
| 57 | const cl = ex * ex + ey * ey;
|
|---|
| 58 | x = x1 + (ey * bl - dy * cl) * d;
|
|---|
| 59 | y = y1 + (dx * cl - ex * bl) * d;
|
|---|
| 60 | }
|
|---|
| 61 | circumcenters[j] = x;
|
|---|
| 62 | circumcenters[j + 1] = y;
|
|---|
| 63 | }
|
|---|
| 64 |
|
|---|
| 65 | // Compute exterior cell rays.
|
|---|
| 66 | let h = hull[hull.length - 1];
|
|---|
| 67 | let p0, p1 = h * 4;
|
|---|
| 68 | let x0, x1 = points[2 * h];
|
|---|
| 69 | let y0, y1 = points[2 * h + 1];
|
|---|
| 70 | vectors.fill(0);
|
|---|
| 71 | for (let i = 0; i < hull.length; ++i) {
|
|---|
| 72 | h = hull[i];
|
|---|
| 73 | p0 = p1, x0 = x1, y0 = y1;
|
|---|
| 74 | p1 = h * 4, x1 = points[2 * h], y1 = points[2 * h + 1];
|
|---|
| 75 | vectors[p0 + 2] = vectors[p1] = y0 - y1;
|
|---|
| 76 | vectors[p0 + 3] = vectors[p1 + 1] = x1 - x0;
|
|---|
| 77 | }
|
|---|
| 78 | }
|
|---|
| 79 | render(context) {
|
|---|
| 80 | const buffer = context == null ? context = new Path : undefined;
|
|---|
| 81 | const {delaunay: {halfedges, inedges, hull}, circumcenters, vectors} = this;
|
|---|
| 82 | if (hull.length <= 1) return null;
|
|---|
| 83 | for (let i = 0, n = halfedges.length; i < n; ++i) {
|
|---|
| 84 | const j = halfedges[i];
|
|---|
| 85 | if (j < i) continue;
|
|---|
| 86 | const ti = Math.floor(i / 3) * 2;
|
|---|
| 87 | const tj = Math.floor(j / 3) * 2;
|
|---|
| 88 | const xi = circumcenters[ti];
|
|---|
| 89 | const yi = circumcenters[ti + 1];
|
|---|
| 90 | const xj = circumcenters[tj];
|
|---|
| 91 | const yj = circumcenters[tj + 1];
|
|---|
| 92 | this._renderSegment(xi, yi, xj, yj, context);
|
|---|
| 93 | }
|
|---|
| 94 | let h0, h1 = hull[hull.length - 1];
|
|---|
| 95 | for (let i = 0; i < hull.length; ++i) {
|
|---|
| 96 | h0 = h1, h1 = hull[i];
|
|---|
| 97 | const t = Math.floor(inedges[h1] / 3) * 2;
|
|---|
| 98 | const x = circumcenters[t];
|
|---|
| 99 | const y = circumcenters[t + 1];
|
|---|
| 100 | const v = h0 * 4;
|
|---|
| 101 | const p = this._project(x, y, vectors[v + 2], vectors[v + 3]);
|
|---|
| 102 | if (p) this._renderSegment(x, y, p[0], p[1], context);
|
|---|
| 103 | }
|
|---|
| 104 | return buffer && buffer.value();
|
|---|
| 105 | }
|
|---|
| 106 | renderBounds(context) {
|
|---|
| 107 | const buffer = context == null ? context = new Path : undefined;
|
|---|
| 108 | context.rect(this.xmin, this.ymin, this.xmax - this.xmin, this.ymax - this.ymin);
|
|---|
| 109 | return buffer && buffer.value();
|
|---|
| 110 | }
|
|---|
| 111 | renderCell(i, context) {
|
|---|
| 112 | const buffer = context == null ? context = new Path : undefined;
|
|---|
| 113 | const points = this._clip(i);
|
|---|
| 114 | if (points === null || !points.length) return;
|
|---|
| 115 | context.moveTo(points[0], points[1]);
|
|---|
| 116 | let n = points.length;
|
|---|
| 117 | while (points[0] === points[n-2] && points[1] === points[n-1] && n > 1) n -= 2;
|
|---|
| 118 | for (let i = 2; i < n; i += 2) {
|
|---|
| 119 | if (points[i] !== points[i-2] || points[i+1] !== points[i-1])
|
|---|
| 120 | context.lineTo(points[i], points[i + 1]);
|
|---|
| 121 | }
|
|---|
| 122 | context.closePath();
|
|---|
| 123 | return buffer && buffer.value();
|
|---|
| 124 | }
|
|---|
| 125 | *cellPolygons() {
|
|---|
| 126 | const {delaunay: {points}} = this;
|
|---|
| 127 | for (let i = 0, n = points.length / 2; i < n; ++i) {
|
|---|
| 128 | const cell = this.cellPolygon(i);
|
|---|
| 129 | if (cell) cell.index = i, yield cell;
|
|---|
| 130 | }
|
|---|
| 131 | }
|
|---|
| 132 | cellPolygon(i) {
|
|---|
| 133 | const polygon = new Polygon;
|
|---|
| 134 | this.renderCell(i, polygon);
|
|---|
| 135 | return polygon.value();
|
|---|
| 136 | }
|
|---|
| 137 | _renderSegment(x0, y0, x1, y1, context) {
|
|---|
| 138 | let S;
|
|---|
| 139 | const c0 = this._regioncode(x0, y0);
|
|---|
| 140 | const c1 = this._regioncode(x1, y1);
|
|---|
| 141 | if (c0 === 0 && c1 === 0) {
|
|---|
| 142 | context.moveTo(x0, y0);
|
|---|
| 143 | context.lineTo(x1, y1);
|
|---|
| 144 | } else if (S = this._clipSegment(x0, y0, x1, y1, c0, c1)) {
|
|---|
| 145 | context.moveTo(S[0], S[1]);
|
|---|
| 146 | context.lineTo(S[2], S[3]);
|
|---|
| 147 | }
|
|---|
| 148 | }
|
|---|
| 149 | contains(i, x, y) {
|
|---|
| 150 | if ((x = +x, x !== x) || (y = +y, y !== y)) return false;
|
|---|
| 151 | return this.delaunay._step(i, x, y) === i;
|
|---|
| 152 | }
|
|---|
| 153 | *neighbors(i) {
|
|---|
| 154 | const ci = this._clip(i);
|
|---|
| 155 | if (ci) for (const j of this.delaunay.neighbors(i)) {
|
|---|
| 156 | const cj = this._clip(j);
|
|---|
| 157 | // find the common edge
|
|---|
| 158 | if (cj) loop: for (let ai = 0, li = ci.length; ai < li; ai += 2) {
|
|---|
| 159 | for (let aj = 0, lj = cj.length; aj < lj; aj += 2) {
|
|---|
| 160 | if (ci[ai] === cj[aj]
|
|---|
| 161 | && ci[ai + 1] === cj[aj + 1]
|
|---|
| 162 | && ci[(ai + 2) % li] === cj[(aj + lj - 2) % lj]
|
|---|
| 163 | && ci[(ai + 3) % li] === cj[(aj + lj - 1) % lj]) {
|
|---|
| 164 | yield j;
|
|---|
| 165 | break loop;
|
|---|
| 166 | }
|
|---|
| 167 | }
|
|---|
| 168 | }
|
|---|
| 169 | }
|
|---|
| 170 | }
|
|---|
| 171 | _cell(i) {
|
|---|
| 172 | const {circumcenters, delaunay: {inedges, halfedges, triangles}} = this;
|
|---|
| 173 | const e0 = inedges[i];
|
|---|
| 174 | if (e0 === -1) return null; // coincident point
|
|---|
| 175 | const points = [];
|
|---|
| 176 | let e = e0;
|
|---|
| 177 | do {
|
|---|
| 178 | const t = Math.floor(e / 3);
|
|---|
| 179 | points.push(circumcenters[t * 2], circumcenters[t * 2 + 1]);
|
|---|
| 180 | e = e % 3 === 2 ? e - 2 : e + 1;
|
|---|
| 181 | if (triangles[e] !== i) break; // bad triangulation
|
|---|
| 182 | e = halfedges[e];
|
|---|
| 183 | } while (e !== e0 && e !== -1);
|
|---|
| 184 | return points;
|
|---|
| 185 | }
|
|---|
| 186 | _clip(i) {
|
|---|
| 187 | // degenerate case (1 valid point: return the box)
|
|---|
| 188 | if (i === 0 && this.delaunay.hull.length === 1) {
|
|---|
| 189 | return [this.xmax, this.ymin, this.xmax, this.ymax, this.xmin, this.ymax, this.xmin, this.ymin];
|
|---|
| 190 | }
|
|---|
| 191 | const points = this._cell(i);
|
|---|
| 192 | if (points === null) return null;
|
|---|
| 193 | const {vectors: V} = this;
|
|---|
| 194 | const v = i * 4;
|
|---|
| 195 | return this._simplify(V[v] || V[v + 1]
|
|---|
| 196 | ? this._clipInfinite(i, points, V[v], V[v + 1], V[v + 2], V[v + 3])
|
|---|
| 197 | : this._clipFinite(i, points));
|
|---|
| 198 | }
|
|---|
| 199 | _clipFinite(i, points) {
|
|---|
| 200 | const n = points.length;
|
|---|
| 201 | let P = null;
|
|---|
| 202 | let x0, y0, x1 = points[n - 2], y1 = points[n - 1];
|
|---|
| 203 | let c0, c1 = this._regioncode(x1, y1);
|
|---|
| 204 | let e0, e1 = 0;
|
|---|
| 205 | for (let j = 0; j < n; j += 2) {
|
|---|
| 206 | x0 = x1, y0 = y1, x1 = points[j], y1 = points[j + 1];
|
|---|
| 207 | c0 = c1, c1 = this._regioncode(x1, y1);
|
|---|
| 208 | if (c0 === 0 && c1 === 0) {
|
|---|
| 209 | e0 = e1, e1 = 0;
|
|---|
| 210 | if (P) P.push(x1, y1);
|
|---|
| 211 | else P = [x1, y1];
|
|---|
| 212 | } else {
|
|---|
| 213 | let S, sx0, sy0, sx1, sy1;
|
|---|
| 214 | if (c0 === 0) {
|
|---|
| 215 | if ((S = this._clipSegment(x0, y0, x1, y1, c0, c1)) === null) continue;
|
|---|
| 216 | [sx0, sy0, sx1, sy1] = S;
|
|---|
| 217 | } else {
|
|---|
| 218 | if ((S = this._clipSegment(x1, y1, x0, y0, c1, c0)) === null) continue;
|
|---|
| 219 | [sx1, sy1, sx0, sy0] = S;
|
|---|
| 220 | e0 = e1, e1 = this._edgecode(sx0, sy0);
|
|---|
| 221 | if (e0 && e1) this._edge(i, e0, e1, P, P.length);
|
|---|
| 222 | if (P) P.push(sx0, sy0);
|
|---|
| 223 | else P = [sx0, sy0];
|
|---|
| 224 | }
|
|---|
| 225 | e0 = e1, e1 = this._edgecode(sx1, sy1);
|
|---|
| 226 | if (e0 && e1) this._edge(i, e0, e1, P, P.length);
|
|---|
| 227 | if (P) P.push(sx1, sy1);
|
|---|
| 228 | else P = [sx1, sy1];
|
|---|
| 229 | }
|
|---|
| 230 | }
|
|---|
| 231 | if (P) {
|
|---|
| 232 | e0 = e1, e1 = this._edgecode(P[0], P[1]);
|
|---|
| 233 | if (e0 && e1) this._edge(i, e0, e1, P, P.length);
|
|---|
| 234 | } else if (this.contains(i, (this.xmin + this.xmax) / 2, (this.ymin + this.ymax) / 2)) {
|
|---|
| 235 | return [this.xmax, this.ymin, this.xmax, this.ymax, this.xmin, this.ymax, this.xmin, this.ymin];
|
|---|
| 236 | }
|
|---|
| 237 | return P;
|
|---|
| 238 | }
|
|---|
| 239 | _clipSegment(x0, y0, x1, y1, c0, c1) {
|
|---|
| 240 | // for more robustness, always consider the segment in the same order
|
|---|
| 241 | const flip = c0 < c1;
|
|---|
| 242 | if (flip) [x0, y0, x1, y1, c0, c1] = [x1, y1, x0, y0, c1, c0];
|
|---|
| 243 | while (true) {
|
|---|
| 244 | if (c0 === 0 && c1 === 0) return flip ? [x1, y1, x0, y0] : [x0, y0, x1, y1];
|
|---|
| 245 | if (c0 & c1) return null;
|
|---|
| 246 | let x, y, c = c0 || c1;
|
|---|
| 247 | if (c & 0b1000) x = x0 + (x1 - x0) * (this.ymax - y0) / (y1 - y0), y = this.ymax;
|
|---|
| 248 | else if (c & 0b0100) x = x0 + (x1 - x0) * (this.ymin - y0) / (y1 - y0), y = this.ymin;
|
|---|
| 249 | else if (c & 0b0010) y = y0 + (y1 - y0) * (this.xmax - x0) / (x1 - x0), x = this.xmax;
|
|---|
| 250 | else y = y0 + (y1 - y0) * (this.xmin - x0) / (x1 - x0), x = this.xmin;
|
|---|
| 251 | if (c0) x0 = x, y0 = y, c0 = this._regioncode(x0, y0);
|
|---|
| 252 | else x1 = x, y1 = y, c1 = this._regioncode(x1, y1);
|
|---|
| 253 | }
|
|---|
| 254 | }
|
|---|
| 255 | _clipInfinite(i, points, vx0, vy0, vxn, vyn) {
|
|---|
| 256 | let P = Array.from(points), p;
|
|---|
| 257 | if (p = this._project(P[0], P[1], vx0, vy0)) P.unshift(p[0], p[1]);
|
|---|
| 258 | if (p = this._project(P[P.length - 2], P[P.length - 1], vxn, vyn)) P.push(p[0], p[1]);
|
|---|
| 259 | if (P = this._clipFinite(i, P)) {
|
|---|
| 260 | for (let j = 0, n = P.length, c0, c1 = this._edgecode(P[n - 2], P[n - 1]); j < n; j += 2) {
|
|---|
| 261 | c0 = c1, c1 = this._edgecode(P[j], P[j + 1]);
|
|---|
| 262 | if (c0 && c1) j = this._edge(i, c0, c1, P, j), n = P.length;
|
|---|
| 263 | }
|
|---|
| 264 | } else if (this.contains(i, (this.xmin + this.xmax) / 2, (this.ymin + this.ymax) / 2)) {
|
|---|
| 265 | P = [this.xmin, this.ymin, this.xmax, this.ymin, this.xmax, this.ymax, this.xmin, this.ymax];
|
|---|
| 266 | }
|
|---|
| 267 | return P;
|
|---|
| 268 | }
|
|---|
| 269 | _edge(i, e0, e1, P, j) {
|
|---|
| 270 | while (e0 !== e1) {
|
|---|
| 271 | let x, y;
|
|---|
| 272 | switch (e0) {
|
|---|
| 273 | case 0b0101: e0 = 0b0100; continue; // top-left
|
|---|
| 274 | case 0b0100: e0 = 0b0110, x = this.xmax, y = this.ymin; break; // top
|
|---|
| 275 | case 0b0110: e0 = 0b0010; continue; // top-right
|
|---|
| 276 | case 0b0010: e0 = 0b1010, x = this.xmax, y = this.ymax; break; // right
|
|---|
| 277 | case 0b1010: e0 = 0b1000; continue; // bottom-right
|
|---|
| 278 | case 0b1000: e0 = 0b1001, x = this.xmin, y = this.ymax; break; // bottom
|
|---|
| 279 | case 0b1001: e0 = 0b0001; continue; // bottom-left
|
|---|
| 280 | case 0b0001: e0 = 0b0101, x = this.xmin, y = this.ymin; break; // left
|
|---|
| 281 | }
|
|---|
| 282 | // Note: this implicitly checks for out of bounds: if P[j] or P[j+1] are
|
|---|
| 283 | // undefined, the conditional statement will be executed.
|
|---|
| 284 | if ((P[j] !== x || P[j + 1] !== y) && this.contains(i, x, y)) {
|
|---|
| 285 | P.splice(j, 0, x, y), j += 2;
|
|---|
| 286 | }
|
|---|
| 287 | }
|
|---|
| 288 | return j;
|
|---|
| 289 | }
|
|---|
| 290 | _project(x0, y0, vx, vy) {
|
|---|
| 291 | let t = Infinity, c, x, y;
|
|---|
| 292 | if (vy < 0) { // top
|
|---|
| 293 | if (y0 <= this.ymin) return null;
|
|---|
| 294 | if ((c = (this.ymin - y0) / vy) < t) y = this.ymin, x = x0 + (t = c) * vx;
|
|---|
| 295 | } else if (vy > 0) { // bottom
|
|---|
| 296 | if (y0 >= this.ymax) return null;
|
|---|
| 297 | if ((c = (this.ymax - y0) / vy) < t) y = this.ymax, x = x0 + (t = c) * vx;
|
|---|
| 298 | }
|
|---|
| 299 | if (vx > 0) { // right
|
|---|
| 300 | if (x0 >= this.xmax) return null;
|
|---|
| 301 | if ((c = (this.xmax - x0) / vx) < t) x = this.xmax, y = y0 + (t = c) * vy;
|
|---|
| 302 | } else if (vx < 0) { // left
|
|---|
| 303 | if (x0 <= this.xmin) return null;
|
|---|
| 304 | if ((c = (this.xmin - x0) / vx) < t) x = this.xmin, y = y0 + (t = c) * vy;
|
|---|
| 305 | }
|
|---|
| 306 | return [x, y];
|
|---|
| 307 | }
|
|---|
| 308 | _edgecode(x, y) {
|
|---|
| 309 | return (x === this.xmin ? 0b0001
|
|---|
| 310 | : x === this.xmax ? 0b0010 : 0b0000)
|
|---|
| 311 | | (y === this.ymin ? 0b0100
|
|---|
| 312 | : y === this.ymax ? 0b1000 : 0b0000);
|
|---|
| 313 | }
|
|---|
| 314 | _regioncode(x, y) {
|
|---|
| 315 | return (x < this.xmin ? 0b0001
|
|---|
| 316 | : x > this.xmax ? 0b0010 : 0b0000)
|
|---|
| 317 | | (y < this.ymin ? 0b0100
|
|---|
| 318 | : y > this.ymax ? 0b1000 : 0b0000);
|
|---|
| 319 | }
|
|---|
| 320 | _simplify(P) {
|
|---|
| 321 | if (P && P.length > 4) {
|
|---|
| 322 | for (let i = 0; i < P.length; i+= 2) {
|
|---|
| 323 | const j = (i + 2) % P.length, k = (i + 4) % P.length;
|
|---|
| 324 | if (P[i] === P[j] && P[j] === P[k] || P[i + 1] === P[j + 1] && P[j + 1] === P[k + 1]) {
|
|---|
| 325 | P.splice(j, 2), i -= 2;
|
|---|
| 326 | }
|
|---|
| 327 | }
|
|---|
| 328 | if (!P.length) P = null;
|
|---|
| 329 | }
|
|---|
| 330 | return P;
|
|---|
| 331 | }
|
|---|
| 332 | }
|
|---|
