Resources Numerical accuracy Calculation policy
Adds points inside a triangulated face and keeps the triangles well shaped, for curved faces read from a .3dm file.
Orientation and in-circle tests are evaluated in 64-bit floating point, as written, without exact arithmetic. The in-circle test subtracts the tested point from the other three first, which keeps its products small. A point too close to a side or corner is left out by a tolerance relative to the triangle's area.
Open item. Compute the orientation and in-circle signs with exact arithmetic, so nearly collinear or cocircular points cannot flip a decision.
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 None found by the review
Twice the signed area of a triangle: which way it turns.
(b − a) × (c − a) = (bₓ − aₓ)(c_y − a_y) − (b_y − a_y)(cₓ − aₓ)
Differences from point a are taken first, then two products are subtracted. For nearly collinear points the two products are nearly equal and the sign can be wrong; it is not computed exactly.
const orient = (X, Y, a, b, c) => (X[b] - X[a]) * (Y[c] - Y[a]) - (Y[b] - Y[a]) * (X[c] - X[a]);Whether a fourth point lies inside the circle through three others.
determinant of the rows (pₓ − dₓ, p_y − d_y, |p − d|²) for p = a, b, c; inside when it is positive
All three points are taken relative to the tested point before squaring, so the size of the coordinates does not matter, only the size of the triangle. The determinant's sign is a floating-point result with no error bound or exact fallback.
function inCircle(X, Y, a, b, c, d) {
const ax = X[a] - X[d], ay = Y[a] - Y[d], bx = X[b] - X[d], by = Y[b] - Y[d], cx = X[c] - X[d], cy = Y[c] - Y[d];
return (ax * ax + ay * ay) * (bx * cy - cx * by) - (bx * bx + by * by) * (ax * cy - cx * ay) + (cx * cx + cy * cy) * (ax * by - bx * ay);
}Walks from the last triangle toward the point, then searches all triangles if the walk fails.
A point is inside when all three edge tests are zero or positive, a plain comparison with zero and no tolerance. The walk is capped at 4096 steps before the full search.
/** The triangle holding (x, y): a walk from the last one, then every triangle when the walk gets lost. */
locate(x, y) {
const X = this.X, Y = this.Y;
const inside = (t) => {
const T = this.T[t];
if (!T) return -1;
const [a, b, c] = T;
const o = (p, q) => (X[q] - X[p]) * (y - Y[p]) - (Y[q] - Y[p]) * (x - X[p]);
const e0 = o(a, b), e1 = o(b, c), e2 = o(c, a);
if (e0 >= 0 && e1 >= 0 && e2 >= 0) return -2;
// step across the side it lies most outside of
const m = Math.min(e0, e1, e2);
const [p, q] = m === e0 ? [a, b] : m === e1 ? [b, c] : [c, a];
const n = this.he.get(this.key(q, p));
return n === undefined ? -1 : n;
};
let t = this.T[this.last] ? this.last : this.T.findIndex(Boolean);
for (let step = 0; step < 4096 && t >= 0; step++) {
const r = inside(t);
if (r === -2) return t;
t = r;
}
for (let k = 0; k < this.T.length; k++) if (inside(k) === -2) return k;
return -1;
}Splits the triangle a new point falls in, or the side it falls on.
ε = 10⁻⁶ × orient(a, b, c); on a side when a sub-area ≤ ε
The tolerance scales with the triangle's own area. A point within it of exactly one side splits that side and the triangle beyond; a point near a corner, on a boundary side or outside by more than the tolerance is left out, and the function returns false.
insert(x, y) {
const t = this.locate(x, y);
if (t < 0) return false;
const X = this.X, Y = this.Y;
const [a, b, c] = this.T[t];
// too close to a side or corner: leave it out (the refinement after this fills in where needed)
const area = orient(X, Y, a, b, c);
const p = X.length;
X.push(x); Y.push(y);
const w0 = orient(X, Y, a, b, p), w1 = orient(X, Y, b, c, p), w2 = orient(X, Y, c, a, p);
const eps = area * 1e-6;
if (Math.min(w0, w1, w2) <= eps) {
// on a side (the other two well inside): split it and the triangle beyond; on a corner or a polygon side: leave it out
const on = [w0, w1, w2].filter((w) => w <= eps).length;
const [q, r, o] = w0 <= eps ? [a, b, c] : w1 <= eps ? [b, c, a] : [c, a, b];
const n = this.he.get(this.key(r, q));
if (on !== 1 || n === undefined || Math.min(w0, w1, w2) < -eps) { X.pop(); Y.pop(); return false; }
const N = this.T[n];
const d = N[0] === r ? N[2] : N[1] === r ? N[0] : N[1];
if (orient(X, Y, r, p, d) <= 0 || orient(X, Y, p, q, d) <= 0) { X.pop(); Y.pop(); return false; }
this.drop(n);
this.put(t, q, p, o); this.put(n, p, r, o);
const t3 = this.T.length; this.put(t3, r, p, d);
const t4 = this.T.length; this.put(t4, p, q, d);
this.last = t;
this.legalize(r, o, p); this.legalize(o, q, p); this.legalize(q, d, p); this.legalize(d, r, p);
return true;
}
this.put(t, a, b, p);
const t1 = this.T.length; this.put(t1, b, c, p);
const t2 = this.T.length; this.put(t2, c, a, p);
this.last = t;
this.legalize(a, b, p); this.legalize(b, c, p); this.legalize(c, a, p);
return true;
}Flips an edge while the opposite corner lies inside the circumcircle.
A flip needs a strictly positive in-circle value and two strictly positive orientation tests for the new triangles, so a flip never creates a turned-over triangle as far as the floating-point signs are right. The recursion stops at depth 64.
/** Edge a→b of triangle (a, b, p): flip it while the far corner lies in the circle. */
legalize(a, b, p, depth = 0) {
if (depth > 64) return;
const t = this.he.get(this.key(a, b));
const n = this.he.get(this.key(b, a));
if (t === undefined || n === undefined) return;
const N = this.T[n];
const d = N[0] === b ? N[2] : N[1] === b ? N[0] : N[1]; // the corner of (b, a, d) not on the edge
if (inCircle(this.X, this.Y, a, b, p, d) <= 0) return;
// the new edge p–d must leave both triangles turning the same way
if (orient(this.X, this.Y, a, d, p) <= 0 || orient(this.X, this.Y, d, b, p) <= 0) return;
this.drop(n);
this.put(t, a, d, p);
this.put(n, d, b, p);
this.legalize(a, d, p, depth + 1);
this.legalize(d, b, p, depth + 1);
}Makes an existing triangulation Delaunay by flipping edges.
The same three tests as the edge flip, over every triangle, for at most 8 passes. It stops early when a pass makes no flip; after 8 passes edges that still fail the test are left as they are.
export function flipAll(X, Y, tris) {
const m = new Mesh2(X, Y, tris);
for (let pass = 0; pass < 8; pass++) {
let flips = 0;
for (let t = 0; t < m.T.length; t++) {
const T = m.T[t];
if (!T) continue;
for (let k = 0; k < 3; k++) {
const a = T[k], b = T[(k + 1) % 3], p = T[(k + 2) % 3];
const n = m.he.get(m.key(b, a));
if (n === undefined) continue;
const N = m.T[n];
const d = N[0] === b ? N[2] : N[1] === b ? N[0] : N[1];
if (inCircle(X, Y, a, b, p, d) > 0 && orient(X, Y, a, d, p) > 0 && orient(X, Y, d, b, p) > 0) {
m.drop(n); m.put(t, a, d, p); m.put(n, d, b, p); flips++;
break;
}
}
}
if (!flips) break;
}
return m.out();
}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.