Resources Numerical accuracy Calculation policy
Cuts the ground surface along a polyline and gives the distance, offset and elevation at every triangle edge it crosses; also measures along, locates on and crosses polylines.
The cut has a point wherever the line crosses a triangle edge, so it follows the triangulated surface with no sampling. Coordinates are shifted by the surface's origin before the crossings are computed. Distances along the line are a running 64-bit sum of segment lengths plus a fraction of the current segment. Points closer than 1e-9 of the line's length are merged.
Open item. Scale the 1e-15 parallel test in segCross to the lengths of the two segments.
8 calculations in 1 source file. The code is read from the released app when this page loads. It is the calculation people run, not a copy kept for this page.
Tests bin/test-civilcraft-river.mjs
The plan length of a polyline.
L = Σ √(Δx² + Δy²)
A running 64-bit sum of Math.hypot segment lengths in order, with no compensation.
export function lineLength(line) {
let L = 0;
for (let i = 1; i < (line?.length || 0); i++) L += Math.hypot(line[i].x - line[i - 1].x, line[i].y - line[i - 1].y);
return L;
}The point and direction at a distance along a polyline.
u = (d − length before the segment) ÷ segment length, clamped to [0, 1]; p = a + (b − a)·u
The distance before the segment is a running sum, then one division and one interpolation. Zero-length segments are skipped before the division. The direction is Math.atan2 of the segment.
export function pointAt(line, d) {
let acc = 0;
for (let i = 1; i < line.length; i++) {
const p = line[i - 1], q = line[i], l = Math.hypot(q.x - p.x, q.y - p.y);
if (l <= 0) continue;
if (d <= acc + l || i === line.length - 1) {
const u = Math.max(0, Math.min(1, (d - acc) / l));
return { x: p.x + (q.x - p.x) * u, y: p.y + (q.y - p.y) * u, dir: Math.atan2(q.y - p.y, q.x - p.x) };
}
acc += l;
}
const p = line[0];
return { x: p.x, y: p.y, dir: 0 };
}The distance along a polyline and the offset to the side of the point nearest to a given point.
u = ((p − a)·(b − a)) ÷ |b − a|², clamped; d = length before + u·|b − a|; offset = ± distance, sign from the cross product
Dot and cross products of coordinate differences in 64-bit; the differences are taken from the segment's start, so they stay small. A later segment replaces the best only when it is closer by more than 1e-12.
export function locateOn(line, p) {
let best = null, acc = 0;
for (let i = 1; i < line.length; i++) {
const a = line[i - 1], b = line[i], dx = b.x - a.x, dy = b.y - a.y, l2 = dx * dx + dy * dy, l = Math.sqrt(l2);
if (l2 <= 0) continue;
const u = Math.max(0, Math.min(1, ((p.x - a.x) * dx + (p.y - a.y) * dy) / l2));
const x = a.x + dx * u, y = a.y + dy * u, dist = Math.hypot(p.x - x, p.y - y);
if (!best || dist < best.dist - 1e-12) {
const cross = dx * (p.y - a.y) - dy * (p.x - a.x); // > 0: p is to the left
best = { d: acc + u * l, off: cross > 0 ? -dist : dist, x, y, dir: Math.atan2(dy, dx), dist };
}
acc += l;
}
if (!best) return null;
delete best.dist;
return best;
}Where two segments cross, as fractions along each.
den = r × s; t = ((q − p) × s) ÷ den; u = ((q − p) × r) ÷ den
Plain 64-bit cross products. Segments count as parallel when |den| < 1e-15, an absolute value not scaled to their lengths. Fractions within 1e-12 outside [0, 1] are accepted and clamped.
function segCross(px, py, rx, ry, qx, qy, sx, sy) { // p + t·r vs q + u·s → [t, u] in [0,1]² or null
const den = rx * sy - ry * sx;
if (Math.abs(den) < 1e-15) return null;
const t = ((qx - px) * sy - (qy - py) * sx) / den, u = ((qx - px) * ry - (qy - py) * rx) / den;
return t >= -EPS && t <= 1 + EPS && u >= -EPS && u <= 1 + EPS ? [Math.max(0, Math.min(1, t)), Math.max(0, Math.min(1, u))] : null;
}The first place one polyline crosses another, with the distance along each.
d_a = length before + t·|segment a|; d_b = length before + u·|segment b|
Distances are running sums of segment lengths plus the crossing fraction times the segment length. The crossing point is interpolated on the first polyline.
export function crossing(a, b) {
let accA = 0;
for (let i = 1; i < a.length; i++) {
const p = a[i - 1], q = a[i], la = Math.hypot(q.x - p.x, q.y - p.y);
let accB = 0, hit = null;
for (let j = 1; j < b.length; j++) {
const r = b[j - 1], s = b[j], lb = Math.hypot(s.x - r.x, s.y - r.y);
const c = segCross(p.x, p.y, q.x - p.x, q.y - p.y, r.x, r.y, s.x - r.x, s.y - r.y);
if (c && (!hit || c[0] < hit.t)) hit = { t: c[0], da: accA + c[0] * la, db: accB + c[1] * lb, x: p.x + (q.x - p.x) * c[0], y: p.y + (q.y - p.y) * c[0] };
accB += lb;
}
if (hit) { delete hit.t; return hit; }
accA += la;
}
return null;
}A straight line square to a polyline at a distance along it, a set length to each side.
right = (sin dir, −cos dir); p = c − right × left; q = c + right × right length
64-bit sine and cosine of the segment direction, one multiply and add per coordinate.
export function squareAcross(line, d, left, right) {
const p = pointAt(line, d);
const rx = Math.sin(p.dir), ry = -Math.cos(p.dir); // right of the direction (y up)
return [{ x: p.x - rx * left, y: p.y - ry * left }, { x: p.x + rx * right, y: p.y + ry * right }];
}The points of the ground along a line: distance, offset, coordinates and elevation at every triangle edge crossing and line vertex.
d = length before + t·|segment|; o = d − origin; z from the surface at (x, y)
Crossings are found in coordinates local to the surface's origin, while the point itself is interpolated in project coordinates. Points closer than max(1e-9 × line length, 1e-12) are merged, a tolerance that scales with the line. Elevations come from the surface's lookup with tolerance 1e-7.
export function cutSection(surface, line, { origin = 0 } = {}) {
const length = lineLength(line);
const res = { length, origin, runs: [], breaks: 0 };
if (!surface || !Array.isArray(line) || line.length < 2 || !(length > 0)) return res;
const { ox, oy } = surface;
const cand = []; // { d, x, y, k }
let acc = 0;
for (let i = 1; i < line.length; i++) {
const a = line[i - 1], b = line[i], l = Math.hypot(b.x - a.x, b.y - a.y);
if (!(l > 0)) continue;
cand.push({ d: acc, x: a.x, y: a.y, k: i === 1 ? 'end' : 'vertex' });
for (const t of edgeParams(surface, a.x - ox, a.y - oy, b.x - ox, b.y - oy)) {
if (t <= EPS || t >= 1 - EPS) continue;
cand.push({ d: acc + t * l, x: a.x + (b.x - a.x) * t, y: a.y + (b.y - a.y) * t, k: 'edge' });
}
acc += l;
}
const last = line[line.length - 1];
cand.push({ d: length, x: last.x, y: last.y, k: 'end' });
cand.sort((u, v) => u.d - v.d);
const tol = Math.max(1e-9 * length, 1e-12);
let run = null, prev = -Infinity;
for (const c of cand) {
if (c.d - prev < tol) { // merge: keep the stronger kind (end > vertex > edge)
const r = run?.[run.length - 1];
if (r && c.k !== 'edge') r.k = r.k === 'end' || c.k === 'end' ? 'end' : 'vertex';
continue;
}
prev = c.d;
const z = surface.elevationAt(c.x, c.y, 1e-7);
if (z === null) { run = null; continue; }
if (!run) res.runs.push((run = []));
run.push({ d: c.d, o: c.d - origin, x: c.x, y: c.y, z, k: c.k });
if (c.k === 'edge') res.breaks++;
}
res.runs = res.runs.filter((r) => r.length >= 2);
return res;
}The elevation of a cut section at any offset, between its points.
z = z₀ + (z₁ − z₀)·(o − o₀) ÷ (o₁ − o₀)
Linear interpolation in 64-bit; a zero-width interval returns the second point's elevation instead of dividing. Offsets are matched with an absolute 1e-9.
export function zAt(runs, o) {
for (const r of runs) {
if (o < r[0].o - 1e-9 || o > r[r.length - 1].o + 1e-9) continue;
for (let i = 1; i < r.length; i++) {
if (o <= r[i].o + 1e-9) { const p = r[i - 1], q = r[i], d = q.o - p.o; return d > 0 ? p.z + ((q.z - p.z) * (o - p.o)) / d : q.z; }
}
}
return null;
}Every value is a 64-bit number from input to result. Formulas are closed form or a direct solve, and every comparison uses a stated tolerance.