entity/geometry/triangulateOnGrid.tstypescript
import earcut from "earcut";
import Delaunator from "delaunator";
import Constrainautor from "@kninnug/constrainautor";
import type { CoordinatesType } from "./Geometry";
type Triangulation = { points: CoordinatesType[]; indexes: number[] };
/**
* Grid points closer to a ring than this part of the step are left out, a point on a ring edge would stop the edge
* from being constrained.
*/
const GRID_POINT_MARGIN = 0.05;
const triangleArea = (a: CoordinatesType, b: CoordinatesType, c: CoordinatesType) =>
((b[0] - a[0]) * (c[1] - a[1]) - (c[0] - a[0]) * (b[1] - a[1])) * 0.5;
function ringArea(ring: CoordinatesType[]) {
let s = 0;
for (let i = 0, j = ring.length - 1; i < ring.length; j = i++) {
s += (ring[j][0] - ring[i][0]) * (ring[j][1] + ring[i][1]);
}
return Math.abs(s * 0.5);
}
function distanceToSegment(x: number, y: number, a: CoordinatesType, b: CoordinatesType) {
const dx = b[0] - a[0],
dy = b[1] - a[1],
l = dx * dx + dy * dy;
let t = l > 0 ? ((x - a[0]) * dx + (y - a[1]) * dy) / l : 0;
t = Math.max(0, Math.min(1, t));
return Math.hypot(a[0] + t * dx - x, a[1] + t * dy - y);
}
/**
* The nodes of the grid inside the polygon, not too close to its rings: a scan along each grid row with the ring
* segments that reach it.
*/
function gridPointsInside(rings: CoordinatesType[][], step: number): CoordinatesType[] {
const rows = new Map<number, [CoordinatesType, CoordinatesType][]>(),
margin = GRID_POINT_MARGIN * step,
out: CoordinatesType[] = [];
for (const ring of rings) {
for (let i = 0, j = ring.length - 1; i < ring.length; j = i++) {
const a = ring[j],
b = ring[i],
from = Math.ceil((Math.min(a[1], b[1]) - margin) / step),
to = Math.floor((Math.max(a[1], b[1]) + margin) / step);
for (let r = from; r <= to; r++) {
const row = rows.get(r);
row ? row.push([a, b]) : rows.set(r, [[a, b]]);
}
}
}
for (const [r, segments] of rows) {
const y = r * step,
xs: number[] = [];
for (const [a, b] of segments) {
if (a[1] > y !== b[1] > y) {
xs.push(a[0] + ((b[0] - a[0]) * (y - a[1])) / (b[1] - a[1]));
}
}
xs.sort((p, q) => p - q);
for (let k = 0; k + 1 < xs.length; k += 2) {
for (let c = Math.ceil((xs[k] + margin) / step); c * step < xs[k + 1] - margin; c++) {
const x = c * step;
let near = false;
for (const [a, b] of segments) {
if (distanceToSegment(x, y, a, b) < margin) {
near = true;
break;
}
}
if (!near) out.push([x, y]);
}
}
}
return out;
}
/**
* Splits the triangles with an edge longer than the step, but for the rounding, at the middle of their longest edge until there is none.
* Every triangle of a long edge splits it at the same point, so the triangles stay joined.
*/
function splitLongEdges(points: CoordinatesType[], indexes: number[], step: number): number[] {
const out: number[] = [],
middles = new Map<number, number>(),
stack = indexes.slice(),
longest = step * (1 + 1e-9);
const length = (a: number, b: number) =>
Math.max(Math.abs(points[b][0] - points[a][0]), Math.abs(points[b][1] - points[a][1]));
const middle = (a: number, b: number) => {
const lo = Math.min(a, b),
hi = Math.max(a, b),
key = (hi * (hi + 1)) / 2 + lo;
let m = middles.get(key);
if (m === undefined) {
m = points.length;
points.push([(points[a][0] + points[b][0]) * 0.5, (points[a][1] + points[b][1]) * 0.5]);
middles.set(key, m);
}
return m;
};
while (stack.length) {
const c = stack.pop()!,
b = stack.pop()!,
a = stack.pop()!,
ab = length(a, b),
bc = length(b, c),
ca = length(c, a);
if (ab <= longest && bc <= longest && ca <= longest) {
out.push(a, b, c);
} else if (ab >= bc && ab >= ca) {
const m = middle(a, b);
stack.push(a, m, c, m, b, c);
} else if (bc >= ca) {
const m = middle(b, c);
stack.push(b, m, a, m, c, a);
} else {
const m = middle(c, a);
stack.push(c, m, b, m, a, b);
}
}
return out;
}
/**
* The constrained Delaunay triangulation of the ring points and of the grid points inside, the triangles out of the
* polygon and in its holes are left out. Throws when the rings are not a valid polygon.
*/
function triangulateDelaunay(rings: CoordinatesType[][], step: number): Triangulation {
const points: CoordinatesType[] = rings.flat(),
first = new Map<string, number>(),
// The point of each Delaunator vertex and the vertex of each point
vertexPoint: number[] = [],
same: number[] = [];
for (let i = 0; i < points.length; i++) {
const k = points[i][0] + "," + points[i][1];
let v = first.get(k);
if (v === undefined) {
v = vertexPoint.length;
vertexPoint.push(i);
first.set(k, v);
}
same.push(v);
}
const edges: [number, number][] = [];
let offset = 0;
for (const ring of rings) {
for (let i = 0; i < ring.length; i++) {
const a = same[offset + i],
b = same[offset + ((i + 1) % ring.length)];
if (a !== b) edges.push([a, b]);
}
offset += ring.length;
}
for (const p of gridPointsInside(rings, step)) {
vertexPoint.push(points.length);
points.push(p);
}
const coords = new Float64Array(vertexPoint.length * 2);
for (let v = 0; v < vertexPoint.length; v++) {
coords[2 * v] = points[vertexPoint[v]][0];
coords[2 * v + 1] = points[vertexPoint[v]][1];
}
const del = new Delaunator(coords),
con = new Constrainautor(del, edges),
{ triangles, halfedges } = del,
count = triangles.length / 3,
inside = new Int8Array(count).fill(-1),
queue: number[] = [];
// Over a hull edge is out of the polygon, over a constrained edge the side changes
for (let e = 0; e < halfedges.length; e++) {
if (halfedges[e] === -1) {
const t = Math.floor(e / 3);
if (inside[t] === -1) {
inside[t] = con.isConstrained(e) ? 1 : 0;
queue.push(t);
}
}
}
while (queue.length) {
const t = queue.pop()!;
for (let e = 3 * t; e < 3 * t + 3; e++) {
const twin = halfedges[e];
if (twin === -1) continue;
const n = Math.floor(twin / 3),
side = con.isConstrained(e) ? 1 - inside[t] : inside[t];
if (inside[n] === -1) {
inside[n] = side;
queue.push(n);
} else if (inside[n] !== side) {
throw new Error("The rings are not a valid polygon");
}
}
}
let indexes: number[] = [],
area = 0;
for (let t = 0; t < count; t++) {
if (inside[t] === 1) {
const a = vertexPoint[triangles[3 * t]],
b = vertexPoint[triangles[3 * t + 1]],
c = vertexPoint[triangles[3 * t + 2]];
indexes.push(a, b, c);
area += triangleArea(points[a], points[b], points[c]);
}
}
let expected = ringArea(rings[0]);
for (let j = 1; j < rings.length; j++) expected -= ringArea(rings[j]);
if (!(Math.abs(Math.abs(area) - expected) <= 1e-9 * expected)) {
throw new Error("The rings are not a valid polygon");
}
indexes = splitLongEdges(points, indexes, step);
return { points, indexes };
}
/**
* Splits a convex polygon by the line where the coordinate on the axis is the value, into the parts below and above.
* The point on the line is the same for both parts and for every polygon with that edge, whichever way it goes.
*/
function splitConvex(poly: CoordinatesType[], axis: number, value: number): [CoordinatesType[], CoordinatesType[]] {
const below: CoordinatesType[] = [],
above: CoordinatesType[] = [],
other = 1 - axis;
for (let i = 0; i < poly.length; i++) {
const a = poly[i],
b = poly[(i + 1) % poly.length],
da = a[axis] - value,
db = b[axis] - value;
if (da <= 0) below.push(a);
if (da >= 0) above.push(a);
if ((da < 0 && db > 0) || (da > 0 && db < 0)) {
const first = a[0] < b[0] || (a[0] === b[0] && a[1] < b[1]),
p = first ? a : b,
q = first ? b : a,
t = (value - p[axis]) / (q[axis] - p[axis]),
c: CoordinatesType = [0, 0];
c[axis] = value;
c[other] = p[other] + t * (q[other] - p[other]);
below.push(c);
above.push(c);
}
}
return [below, above];
}
/**
* Earcut triangles cut by the grid. Earcut takes the rings that are not a valid polygon too.
*/
function triangulateEarcut(rings: CoordinatesType[][], step: number): Triangulation {
const points: CoordinatesType[] = rings.flat(),
vertices: number[] = [],
holes: number[] = [];
let offset = 0;
for (let j = 0; j < rings.length; j++) {
if (j > 0) holes.push(offset);
offset += rings[j].length;
}
for (const p of points) vertices.push(p[0], p[1]);
const tri = earcut(vertices, holes, 2);
const index = new Map<string, number>();
for (let i = 0; i < points.length; i++) {
const k = points[i][0] + "," + points[i][1];
if (!index.has(k)) index.set(k, i);
}
const indexOf = (p: CoordinatesType) => {
const k = p[0] + "," + p[1];
let i = index.get(k);
if (i === undefined) {
i = points.length;
points.push([p[0], p[1]]);
index.set(k, i);
}
return i;
};
const indexes: number[] = [];
const addCell = (cell: CoordinatesType[]) => {
if (cell.length < 3) return;
const i0 = indexOf(cell[0]);
let i1 = indexOf(cell[1]);
for (let k = 2; k < cell.length; k++) {
const i2 = indexOf(cell[k]);
if (i0 !== i1 && i1 !== i2 && i2 !== i0) indexes.push(i0, i1, i2);
i1 = i2;
}
};
for (let t = 0; t < tri.length; t += 3) {
const triangle = [points[tri[t]], points[tri[t + 1]], points[tri[t + 2]]];
let minX = Infinity,
maxX = -Infinity,
minY = Infinity,
maxY = -Infinity;
for (const p of triangle) {
minX = Math.min(minX, p[0]);
maxX = Math.max(maxX, p[0]);
minY = Math.min(minY, p[1]);
maxY = Math.max(maxY, p[1]);
}
let rest = triangle;
for (let i = Math.floor(minX / step) + 1; rest.length >= 3; i++) {
let column: CoordinatesType[];
if (i * step < maxX) {
[column, rest] = splitConvex(rest, 0, i * step);
} else {
column = rest;
rest = [];
}
let up = column;
for (let j = Math.floor(minY / step) + 1; up.length >= 3; j++) {
let cell: CoordinatesType[];
if (j * step < maxY) {
[cell, up] = splitConvex(up, 1, j * step);
} else {
cell = up;
up = [];
}
addCell(cell);
}
}
}
return { points, indexes };
}
/**
* Triangulates a polygon with no edge longer than the step, for the cube tiles.
* @param {Array.<Array.<CoordinatesType>>} rings - Geodetic rings, the outer one and the holes.
* @param {number} step - Longest edge in degrees.
* @returns {{points: Array.<CoordinatesType>, indexes: Array.<number>}} - The ring points, then the new ones.
*/
export function triangulateOnGrid(rings: CoordinatesType[][], step: number): Triangulation {
try {
return triangulateDelaunay(rings, step);
} catch {
return triangulateEarcut(rings, step);
}
}