Resources Numerical accuracy Calculation policy
Builds a triangulated surface from survey points and breaklines, then its elevations, contours and slopes.
Coordinates are 64-bit and are moved to a local origin first, so state-plane values in the millions keep their digits. The two decisions that shape the mesh (which side of an edge a point is on, and whether a point is inside a triangle's circle) are computed exactly with adaptive-precision arithmetic (robust-predicates 3.0.3). Points closer than 0.0001 project units in plan are one vertex.
6 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.mjs
Whether three points turn counter-clockwise, clockwise, or lie on one line.
sign of (b − a) × (c − a)
orient2d returns a value whose sign is always the true sign, whatever the coordinates. The sign is turned over here because the library counts y downward.
orient(a, b, c) { const { X, Y } = this; return -orient2d(X[a], Y[a], X[b], Y[b], X[c], Y[c]); } // > 0: counter-clockwise (y up)Whether the point across an edge lies inside the circle through the triangle on this side.
sign of the in-circle determinant of a, b, c, d
incircle is exact in the same way. An edge is flipped only when the answer is strictly inside, so four points on one circle never flip back and forth.
illegal(e) {
const o = this.he[e];
if (o < 0) return false;
const { tri, X, Y } = this;
const a = tri[e], b = tri[next(e)], c = tri[prev(e)], d = tri[prev(o)];
return incircle(X[a], Y[a], X[b], Y[b], X[c], Y[c], X[d], Y[d]) > 0 && this.convex(e);
}The triangle a point falls in and its three barycentric weights.
w₁ = ((y_b − y_c)(x − x_c) + (x_c − x_b)(y − y_c)) ÷ d, w₃ = 1 − w₁ − w₂
Ordinary 64-bit arithmetic in local coordinates. A point within 1e-9 of the surface edge is still placed in the nearest triangle, with its weights clamped to the edge.
/** The triangle under (local) x, y and its barycentric weights, or null. */
locate(lx, ly, tol = 1e-9) {
const g = this.grid(), { X, Y, tris } = this;
if (lx < g.b.minX - tol || lx > g.b.maxX + tol || ly < g.b.minY - tol || ly > g.b.maxY + tol) return null;
const i = clampI((lx - g.b.minX) / g.cw, g.nx), j = clampI((ly - g.b.minY) / g.ch, g.ny), k = j * g.nx + i;
let best = null, bestMin = -Infinity;
for (let s = g.start[k]; s < g.start[k + 1]; s++) {
const t = g.items[s], a = tris[3 * t], b = tris[3 * t + 1], c = tris[3 * t + 2];
const d = (Y[b] - Y[c]) * (X[a] - X[c]) + (X[c] - X[b]) * (Y[a] - Y[c]);
if (d === 0) continue;
const w1 = ((Y[b] - Y[c]) * (lx - X[c]) + (X[c] - X[b]) * (ly - Y[c])) / d;
const w2 = ((Y[c] - Y[a]) * (lx - X[c]) + (X[a] - X[c]) * (ly - Y[c])) / d;
const w3 = 1 - w1 - w2, m = Math.min(w1, w2, w3);
if (m >= 0) return { t, w: [w1, w2, w3] };
if (m > bestMin) { bestMin = m; best = { t, w: [w1, w2, w3] }; }
}
if (best && tol > 0) {
// within tol of the triangle in plan (points on the outer edge)
const t = best.t, a = tris[3 * t], b = tris[3 * t + 1], c = tris[3 * t + 2];
const scale = Math.max(Math.hypot(X[b] - X[a], Y[b] - Y[a]), Math.hypot(X[c] - X[b], Y[c] - Y[b]), Math.hypot(X[a] - X[c], Y[a] - Y[c]));
if (-bestMin * scale <= tol) { const w = best.w.map((v) => Math.max(0, v)); const s = w[0] + w[1] + w[2]; return { t, w: w.map((v) => v / s) }; }
}
return null;
}The surface elevation at an easting and northing.
z = w₁·z₁ + w₂·z₂ + w₃·z₃
The local origin is subtracted before the triangle is looked up. The elevation is linear across each triangle.
elevationAt(x, y, tol = 1e-9) {
const r = this.locate(x - this.ox, y - this.oy, tol);
if (!r) return null;
const { tris, Z } = this, t = r.t;
return r.w[0] * Z[tris[3 * t]] + r.w[1] * Z[tris[3 * t + 1]] + r.w[2] * Z[tris[3 * t + 2]];
}Rise over run and the downhill direction of one triangle.
n = u × v; slope = √((n_x ÷ n_z)² + (n_y ÷ n_z)²)
Taken from the plane through the three corners. No finite differences.
/** Triangle t's plane: slope (rise/run), downhill direction (radians), gradient. */
slopeOf(t) {
const { X, Y, Z, tris } = this;
const a = tris[3 * t], b = tris[3 * t + 1], c = tris[3 * t + 2];
const ux = X[b] - X[a], uy = Y[b] - Y[a], uz = Z[b] - Z[a], vx = X[c] - X[a], vy = Y[c] - Y[a], vz = Z[c] - Z[a];
const nx = uy * vz - uz * vy, ny = uz * vx - ux * vz, nz = ux * vy - uy * vx;
const gx = -nx / nz, gy = -ny / nz; // dz/dx, dz/dy
return { slope: Math.hypot(gx, gy), dir: Math.atan2(-gy, -gx), gx, gy };
}The surface's area in plan and along the slope, with its elevation range.
area in plan = Σ |n_z| ÷ 2; surface area = Σ |u × v| ÷ 2
Summed triangle by triangle in local coordinates with ordinary 64-bit addition.
stats() {
if (this._stats) return this._stats;
const { X, Y, Z, tris, ox, oy } = this;
const used = this.usedVertices();
let n = 0, minZ = Infinity, maxZ = -Infinity, sum = 0, minX = Infinity, minY = Infinity, maxX = -Infinity, maxY = -Infinity;
for (let i = 0; i < X.length; i++) {
if (!used[i]) continue;
n++; sum += Z[i];
if (Z[i] < minZ) minZ = Z[i]; if (Z[i] > maxZ) maxZ = Z[i];
if (X[i] < minX) minX = X[i]; if (X[i] > maxX) maxX = X[i]; if (Y[i] < minY) minY = Y[i]; if (Y[i] > maxY) maxY = Y[i];
}
let a2 = 0, a3 = 0;
for (let t = 0; t < tris.length; t += 3) {
const a = tris[t], b = tris[t + 1], c = tris[t + 2];
const ux = X[b] - X[a], uy = Y[b] - Y[a], uz = Z[b] - Z[a], vx = X[c] - X[a], vy = Y[c] - Y[a], vz = Z[c] - Z[a];
const cx = uy * vz - uz * vy, cy = uz * vx - ux * vz, cz = ux * vy - uy * vx;
a2 += Math.abs(cz) / 2;
a3 += Math.hypot(cx, cy, cz) / 2;
}
this._stats = { points: n, triangles: this.nT, minZ, maxZ, meanZ: sum / n, area2d: a2, area3d: a3, bounds: { minX: minX + ox, minY: minY + oy, maxX: maxX + ox, maxY: maxY + oy } };
return this._stats;
}Coordinates are 64-bit. The yes-or-no decisions that shape the geometry (which side of a line, inside or outside a circle) are computed exactly, so rounding cannot flip one.