Builds the triangulated ground surface from points, breaklines and a boundary, and works out its contours, elevation bands and outline.
How the precision is kept
Coordinates are moved to a whole-number local origin before triangulating and held in Float64Arrays; the file's header gives the reason, that state-plane sized values keep their precision. The triangulation's orientation and in-circle decisions use the robust-predicates library; the steps recorded here (merging points, splitting crossing breaklines, clipping, contours) are plain 64-bit arithmetic with fixed tolerances in the caller's units.
Calculations under this policy
17 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.
The origin is the whole number at or below the smallest coordinate, and every vertex is stored as its difference from it. Points are merged when they round to the same cell of a 1e-4 grid, so two points closer than 1e-4 that straddle a cell edge stay separate. The first point's position is kept; a breakline's elevation replaces a point's when they differ by more than 1e-3.
tin.js
lines 67 to 81
BimCraft 0.8.1
ox = Math.floor(ox); oy = Math.floor(oy);
const X = [], Y = [], Z = [], key = new Map();
const kOf = (x, y) => Math.round(x / SNAP) + ',' + Math.round(y / SNAP);
const addV = (x, y, z, hard) => {
const lx = x - ox, ly = y - oy, k = kOf(lx, ly);
const at = key.get(k);
if (at !== undefined) {
if (Math.abs(Z[at] - z) > 1e-3) { report.zConflicts++; if (hard) Z[at] = z; }
else if (!hard) report.merged++;
return at;
}
key.set(k, X.length); X.push(lx); Y.push(ly); Z.push(z);
return X.length - 1;
};
Coordinates handed to the triangulator
64-bit
Packs the local coordinates for Delaunator and runs the Delaunay triangulation.
coords = (x, −y) per vertex
The coordinates go into a Float64Array unchanged except for the sign of y, which is exact. Delaunator's own arithmetic is in that library.
tin.js
lines 103 to 105
BimCraft 0.8.1
const coords = new Float64Array(X.length * 2);
for (let i = 0; i < X.length; i++) { coords[2 * i] = X[i]; coords[2 * i + 1] = -Y[i]; } // mirrored: Delaunator's triangles come out counter-clockwise in y-up
const del = new Delaunator(coords);
Boundary clipped to the data and given elevations
64-bit
Clips the outer boundary to the convex hull of the data and takes each corner's elevation from a first surface built without it.
z_corner = elevationAt(x, y) with a 1e-6 edge tolerance
The clipped corners are computed by line intersection in project coordinates, not local ones. A corner whose elevation cannot be found stops the build with an error.
Drops the triangles whose center lies outside the boundary.
keep when the centroid (Σx ÷ 3, Σy ÷ 3) is inside the boundary polygon
One point-in-polygon test per triangle at its centroid, with no tolerance. The boundary was forced in as edges first, so no triangle straddles it unless a boundary edge could not be forced.
tin.js
lines 115 to 119
BimCraft 0.8.1
const poly = boundary.map((p) => ({ x: p.x - ox, y: p.y - oy }));
for (let t = 0; t < keep.length; t++) {
const a = mesh.tri[3 * t], b = mesh.tri[3 * t + 1], c = mesh.tri[3 * t + 2];
if (!pointInPoly((X[a] + X[b] + X[c]) / 3, (Y[a] + Y[b] + Y[c]) / 3, poly)) keep[t] = 0;
}
Surface from given triangles
64-bit
Builds a surface from vertices and faces read from a file, keeping every triangle counter-clockwise.
o = (x_b − x_a)(y_c − y_a) − (y_b − y_a)(x_c − x_a); swap b, c when o < 0
The orientation here is a plain 64-bit cross product on local coordinates, not the exact predicate. A face is dropped only when the product is exactly zero, so a nearly flat sliver can be kept with either orientation.
tin.js
lines 135 to 157
BimCraft 0.8.1
/** A TIN given as vertices and faces (from LandXML). Faces with a missing vertex are dropped. */
export function surfaceFromFaces({ pts, faces }) {
const report = { notes: [], errors: [] };
let ox = Infinity, oy = Infinity;
for (const p of pts) if (p) { ox = Math.min(ox, p.x); oy = Math.min(oy, p.y); }
if (!Number.isFinite(ox)) { report.errors.push(t('tin.err.nopts')); return { surface: null, report }; }
ox = Math.floor(ox); oy = Math.floor(oy);
const X = new Float64Array(pts.length), Y = new Float64Array(pts.length), Z = new Float64Array(pts.length);
pts.forEach((p, i) => { if (p) { X[i] = p.x - ox; Y[i] = p.y - oy; Z[i] = p.z; } });
const tris = [];
let bad = 0;
for (const f of faces) {
if (f.some((i) => !Number.isInteger(i) || i < 0 || i >= pts.length || !pts[i]) || f[0] === f[1] || f[1] === f[2] || f[0] === f[2]) { bad++; continue; }
let [a, b, c] = f;
const o = (X[b] - X[a]) * (Y[c] - Y[a]) - (Y[b] - Y[a]) * (X[c] - X[a]);
if (o === 0) { bad++; continue; }
if (o < 0) [b, c] = [c, b]; // keep every triangle counter-clockwise
tris.push(a, b, c);
}
if (bad) report.notes.push(t('tin.badfaces', { n: bad }));
if (!tris.length) { report.errors.push(t('tin.err.nofaces')); return { surface: null, report }; }
return { surface: new Surface(X, Y, Z, Int32Array.from(tris), ox, oy), report };
}
Crossing of two breakline segments
64-bit
Finds where two segments properly cross, as a fraction along each.
t = cross(c − a, s) ÷ cross(r, s); u = cross(c − a, r) ÷ cross(r, s)
Segments are parallel when |cross(r, s)| < 1e-12, an absolute test on a product of two lengths. A crossing within 1e-9 of either end of either segment, as a fraction, is not a crossing.
tin.js
lines 201 to 210
BimCraft 0.8.1
/** Proper crossing of two segments (not at their ends): {t along the first, u along the second}. */
function segHit(ax, ay, bx, by, cx, cy, dx, dy) {
const rx = bx - ax, ry = by - ay, sx = dx - cx, sy = dy - cy;
const den = rx * sy - ry * sx;
if (Math.abs(den) < 1e-12) return null;
const t = ((cx - ax) * sy - (cy - ay) * sx) / den, u = ((cx - ax) * ry - (cy - ay) * rx) / den;
const eps = 1e-9;
if (t <= eps || t >= 1 - eps || u <= eps || u >= 1 - eps) return null;
return { t, u };
}
Vertex at a breakline crossing
64-bit
Places a new vertex where two breaklines cross and gives it an elevation.
x = x_a + (x_b − x_a) t; z from the first breakline
The position is interpolated along the first segment only. When both are breaklines the one given first sets the elevation, and a difference over 0.01 between the two is counted as a conflict. The new vertex is merged with an existing one on the same 1e-4 cell.
tin.js
lines 177 to 185
BimCraft 0.8.1
const x = X[a] + (X[b] - X[a]) * hit.t, y = Y[a] + (Y[b] - Y[a]) * hit.t;
const zi = Z[a] + (Z[b] - Z[a]) * hit.t, zj = Z[c] + (Z[d] - Z[c]) * hit.u;
const k = kOf(x, y);
let v = key.get(k);
// a breakline gives the crossing its elevation; a boundary edge does not
// between two breaklines the one given first wins
const z = segs[i][2] && segs[j][2] ? (i < j ? zi : zj) : segs[i][2] ? zi : segs[j][2] ? zj : zi;
if (v === undefined) { v = X.length; key.set(k, v); X.push(x); Y.push(y); Z.push(z); }
if (segs[i][2] && segs[j][2] && Math.abs(zi - zj) > 0.01) report.zConflicts++;
Side of a breakline a vertex is on
64-bit
Decides whether a vertex lies left of, right of, or on the line of a breakline being forced in.
on the line when |orient| ÷ L < 1e-7 × (1 + L)
The sign comes from the exact orientation predicate, but a vertex closer to the line than 1e-7 × (1 + length) is treated as on it. That tolerance grows with the breakline's length and is in the caller's unit.
tin.js
lines 291 to 299
BimCraft 0.8.1
/** Which side of p→q v is on (1 left, -1 right), 0 when it is within a hair of the line. */
side(p, q, v) {
const { X, Y } = this;
const o = this.orient(p, q, v);
if (o === 0) return 0;
const L = Math.hypot(X[q] - X[p], Y[q] - Y[p]);
if (Math.abs(o) / L < 1e-7 * (1 + L)) return 0;
return o > 0 ? 1 : -1;
}
Whether a diagonal can be flipped
Exact
Tests that the two triangles on an edge form a strictly convex quadrilateral.
orient(c, d, a) and orient(c, d, b) have opposite signs
Both signs come from robust-predicates' orient2d through the mesh's orient method, so the decision is exact for the stored coordinates.
tin.js
lines 382 to 390
BimCraft 0.8.1
/** The two triangles on e make a strictly convex quadrilateral (so e can be flipped). */
convex(e) {
const o = this.he[e];
if (o < 0) return false;
const { tri } = this;
const a = tri[e], b = tri[next(e)], c = tri[prev(e)], d = tri[prev(o)];
const s1 = this.orient(c, d, a), s2 = this.orient(c, d, b);
return (s1 > 0 && s2 < 0) || (s1 < 0 && s2 > 0);
}
Whether an edge crosses a breakline
Exact
Tests whether two segments properly cross during the flips that force a breakline in.
each segment's ends lie on opposite sides of the other
Four exact orientation signs; a touching end counts as not crossing.
tin.js
lines 378 to 381
BimCraft 0.8.1
crosses(p, q, u, v) {
const a = this.orient(p, q, u), b = this.orient(p, q, v), c = this.orient(u, v, p), d = this.orient(u, v, q);
return ((a > 0 && b < 0) || (a < 0 && b > 0)) && ((c > 0 && d < 0) || (c < 0 && d > 0));
}
Long outer triangles removed
64-bit
Removes triangles on the outside whose outer edge is longer than a limit, working inward.
remove when outer edge length > maxEdge
Edge lengths use Math.hypot and are compared with the limit with no tolerance. An edge that is a breakline is never peeled. The bitwise operations in this method act on triangle indexes only.
tin.js
lines 418 to 431
BimCraft 0.8.1
/** Removes outside triangles whose outer edge is longer than maxEdge, working inward. */
peel(keep, maxEdge) {
const { tri, he, X, Y, con } = this;
const len = (e) => Math.hypot(X[tri[next(e)]] - X[tri[e]], Y[tri[next(e)]] - Y[tri[e]]);
const outer = (e) => { const o = he[e]; return o < 0 || !keep[(o / 3) | 0]; };
const stack = [];
for (let e = 0; e < tri.length; e++) if (keep[(e / 3) | 0] && outer(e)) stack.push(e);
while (stack.length) {
const e = stack.pop(), t = (e / 3) | 0;
if (!keep[t] || con[e] || !outer(e) || len(e) <= maxEdge) continue;
keep[t] = 0;
for (const f of [3 * t, 3 * t + 1, 3 * t + 2]) { const o = he[f]; if (o >= 0 && keep[(o / 3) | 0]) stack.push(o); }
}
}
Contour lines
64-bit
The contour lines of the surface at every multiple of an interval.
level = base + k × interval; t = (L − zᵢ) ÷ (zⱼ − zᵢ); point = pᵢ + (pⱼ − pᵢ) t
Each level is base + k × interval, a product, not a running sum. A vertex within interval × 1e-7 of a level is treated as that much above it, which moves the contour there by at most that height. A crossing on an edge shared by two triangles is computed in each of them and matched by the edge's two vertex numbers, not by coordinates. Loops smaller than 1e-6 of the surface's diagonal are dropped.
tin.js
lines 563 to 592
BimCraft 0.8.1
/** Contour polylines every `interval` (levels at base + k·interval), joined across triangles. */
contours(interval, { base = 0, maxLevels = 2000 } = {}) {
if (!(interval > 0)) return [];
const { X, Y, Z, tris, ox, oy } = this, s = this.stats();
const k0 = Math.ceil((s.minZ - base) / interval - 1e-9), k1 = Math.floor((s.maxZ - base) / interval + 1e-9);
if (k1 - k0 > maxLevels) return null;
const out = [];
const eps = interval * 1e-7;
const tiny = 1e-6 * Math.hypot(s.bounds.maxX - s.bounds.minX, s.bounds.maxY - s.bounds.minY);
for (let k = k0; k <= k1; k++) {
const L = base + k * interval;
const zv = (i) => (Math.abs(Z[i] - L) < eps ? L + eps : Z[i]); // a vertex exactly on the level counts as just above it
const segs = []; // [keyA, keyB, ax, ay, bx, by]
const cross = (i, j) => {
const zi = zv(i), zj = zv(j);
if ((zi > L) === (zj > L)) return null;
const t = (L - zi) / (zj - zi);
return { k: i < j ? i + '_' + j : j + '_' + i, x: X[i] + (X[j] - X[i]) * t, y: Y[i] + (Y[j] - Y[i]) * t };
};
for (let t = 0; t < tris.length; t += 3) {
const a = tris[t], b = tris[t + 1], c = tris[t + 2];
const za = zv(a), zb = zv(b), zc = zv(c);
if ((za > L && zb > L && zc > L) || (za < L && zb < L && zc < L)) continue;
const hits = [cross(a, b), cross(b, c), cross(c, a)].filter(Boolean);
if (hits.length === 2) segs.push(hits);
}
for (const line of chain(segs)) if (extent(line.pts) > tiny) out.push({ z: L, pts: line.pts.map((p) => ({ x: p.x + ox, y: p.y + oy })), closed: line.closed });
}
return out;
}
Triangle clipped at an elevation
64-bit
Keeps the part of a triangle above or below a level, for the elevation bands.
t = f_p ÷ (f_p − f_q); new point = p + (q − p) t at z = L
Linear interpolation on edges whose ends are on opposite sides; the new point is given the level itself as its height. A vertex exactly on the level counts as kept.
tin.js
lines 665 to 675
BimCraft 0.8.1
/** Keeps the part of a triangle polygon [x,y,z] with sign·(z − L) ≥ 0. */
function clipZ(poly, L, sign) {
const out = [];
for (let i = 0; i < poly.length; i++) {
const p = poly[i], q = poly[(i + 1) % poly.length];
const fp = sign * (p[2] - L), fq = sign * (q[2] - L);
if (fp >= 0) out.push(p);
if ((fp >= 0) !== (fq >= 0)) { const t = fp / (fp - fq); out.push([p[0] + (q[0] - p[0]) * t, p[1] + (q[1] - p[1]) * t, L]); }
}
return out;
}
Convex hull
64-bit
The convex hull of a set of points, counter-clockwise.
monotone chain; pop while cross(o, a, b) ≤ 0
The turn test is a plain 64-bit cross product in project coordinates, not the exact predicate, so for nearly collinear points far from the origin the sign can be decided by rounding. Collinear points are left off the hull.
tin.js
lines 715 to 724
BimCraft 0.8.1
export function convexHull(pts) {
const p = pts.slice().sort((a, b) => a.x - b.x || a.y - b.y);
if (p.length < 3) return p;
const cr = (o, a, b) => (a.x - o.x) * (b.y - o.y) - (a.y - o.y) * (b.x - o.x);
const lo = [], hi = [];
for (const q of p) { while (lo.length >= 2 && cr(lo[lo.length - 2], lo[lo.length - 1], q) <= 0) lo.pop(); lo.push(q); }
for (let i = p.length - 1; i >= 0; i--) { const q = p[i]; while (hi.length >= 2 && cr(hi[hi.length - 2], hi[hi.length - 1], q) <= 0) hi.pop(); hi.push(q); }
lo.pop(); hi.pop();
return lo.concat(hi);
}
Polygon area
64-bit
The signed area of a polygon.
A = ½ Σ (xᵢ yᵢ₊₁ − xᵢ₊₁ yᵢ)
The shoelace sum on the coordinates as given, added in order. With large project coordinates each product is far larger than the area, so leading digits cancel.
tin.js
line 725
BimCraft 0.8.1
export function polyArea(pts) { let a = 0; for (let i = 0; i < pts.length; i++) { const p = pts[i], q = pts[(i + 1) % pts.length]; a += p.x * q.y - q.x * p.y; } return a / 2; }
Point inside a polygon
64-bit
Tests whether a point is inside a polygon.
count crossings of a ray to the right; odd = inside
The crossing's x is one division and one multiplication per edge, with no tolerance, so a point on or within rounding of an edge can fall either way.
tin.js
lines 727 to 734
BimCraft 0.8.1
export function pointInPoly(x, y, pts) {
let inside = false;
for (let i = 0, j = pts.length - 1; i < pts.length; j = i++) {
const a = pts[i], b = pts[j];
if ((a.y > y) !== (b.y > y) && x < ((b.x - a.x) * (y - a.y)) / (b.y - a.y) + a.x) inside = !inside;
}
return inside;
}
Polygon clipped by a convex polygon
64-bit
Clips a polygon to a convex polygon, used to keep a boundary inside the data's hull.
inside when cross(b − a, p − a) ≥ −1e-9; t = cross(a − p, e) ÷ cross(d, e)
The inside test allows −1e-9 on a cross product, an absolute value in squared length units that does not scale with the polygon's size. Points closer than 1e-9 to the next one are dropped from the result.
tin.js
lines 735 to 755
BimCraft 0.8.1
/** Sutherland–Hodgman: subject polygon clipped by a convex counter-clockwise polygon. */
export function clipConvex(subject, clip) {
let out = subject;
for (let i = 0; i < clip.length && out.length; i++) {
const a = clip[i], b = clip[(i + 1) % clip.length];
const inside = (p) => (b.x - a.x) * (p.y - a.y) - (b.y - a.y) * (p.x - a.x) >= -1e-9;
const inp = out; out = [];
for (let k = 0; k < inp.length; k++) {
const p = inp[k], q = inp[(k + 1) % inp.length];
const ip = inside(p), iq = inside(q);
if (ip) out.push(p);
if (ip !== iq) {
const dx = q.x - p.x, dy = q.y - p.y, ex = b.x - a.x, ey = b.y - a.y;
const t = ((a.x - p.x) * ey - (a.y - p.y) * ex) / (dx * ey - dy * ex);
out.push({ x: p.x + dx * t, y: p.y + dy * t });
}
}
}
// drop repeats
return out.filter((p, i) => { const q = out[(i + 1) % out.length]; return Math.hypot(p.x - q.x, p.y - q.y) > 1e-9; });
}
In detail
Method
Delaunator for the Delaunay triangulation; breaklines forced in by edge flipping (the file cites Sloan 1993) and Delaunay restored round them by Lawson flips; contours by linear interpolation along triangle edges; hull by the monotone chain.
Number types
Float64Array coordinates, Int32Array triangle and half-edge indexes. Libraries: delaunator and robust-predicates (orient2d, incircle).
Tolerances
Points closer than 1e-4 in plan are one vertex. An elevation difference over 1e-3 at a merged vertex, or over 0.01 where two breaklines cross, is counted as a conflict. A vertex is on a breakline when its distance from it is under 1e-7 × (1 + length). Contour levels treat a vertex within interval × 1e-7 of the level as just above it. All are in the caller's length unit; BimCraft passes millimeters.
Limits
The fixed tolerances do not scale with the unit, so the same file read in meters and in millimeters merges points at different real distances. The convex hull and the polygon helpers use plain floating-point cross products, not the exact predicates.
Approach: 64-bit, direct
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.