Resources Numerical accuracy Calculation policy

Triangulation of solid faces

Splits a face of an imported 3D solid, with its holes, into triangles; the triangles are what the solid's mesh is built from.

How the precision is kept

Points are copied into Float64Array. Orientation and in-circle tests are plain 64-bit determinants, not exact predicates. Side tests use tolerances scaled to the size D of the region: 1e-13 × D for lengths and 1e-13 × D² for areas. The in-circle test is compared with zero. Every loop has an iteration cap.

Calculations under this policy

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.

acis/cdt.js CadCraft 0.8.2 64-bit 6 calculations Reviewed

Tests tools/acis/test-apps.mjs, tools/acis/test.mjs

Region size, enclosing triangle and tolerances

64-bit

Measures the region and sets the tolerances every later test uses.

D = max(Δx, Δy); eps = 1e-13 × D; epsA = 1e-13 × D²

Both tolerances are scaled to the region's size, not absolute. The enclosing triangle's corners are 20 × D to either side of and 10 × D below the center of the box, and 20 × D above it.

  • acis/cdt.js
  • lines 20 to 27
  • CadCraft 0.8.2
const D = Math.max(x1 - x0, y1 - y0) || 1, cx = (x0 + x1) / 2, cy = (y0 + y1) / 2;
// the points, then three far corners of a triangle around all of them
const X = new Float64Array(n + 3), Y = new Float64Array(n + 3);
for (let i = 0; i < n; i++) { X[i] = xy[2 * i]; Y[i] = xy[2 * i + 1]; }
X[n] = cx - 20 * D; Y[n] = cy - 10 * D;
X[n + 1] = cx + 20 * D; Y[n + 1] = cy - 10 * D;
X[n + 2] = cx; Y[n + 2] = cy + 20 * D;
const eps = D * 1e-13, epsA = D * D * 1e-13;   // a length, an area (orient() is twice a triangle's area)

Orientation of three points

64-bit

Twice the signed area of the triangle a, b, c: positive when anticlockwise.

(x_b − x_a) × (y_c − y_a) − (y_b − y_a) × (x_c − x_a)

A 2 × 2 determinant of coordinate differences in 64-bit; the two products and their subtraction are each rounded, and the subtraction cancels for nearly collinear points. Callers compare the result with ±epsA, not with zero.

  • acis/cdt.js
  • line 28
  • CadCraft 0.8.2
const orient = (a, b, c) => (X[b] - X[a]) * (Y[c] - Y[a]) - (Y[b] - Y[a]) * (X[c] - X[a]);

In-circle test

64-bit

Whether point d lies inside the circle through a, b and c.

det of rows (a − d, |a − d|²), (b − d, |b − d|²), (c − d, |c − d|²) > 0

Coordinates are taken relative to d before squaring, then a 3 × 3 determinant is expanded in 64-bit. The caller (legalize) compares the result with exactly 0: there is no tolerance and no exact fallback, so four nearly cocircular points are decided by rounding.

  • acis/cdt.js
  • lines 29 to 32
  • CadCraft 0.8.2
const inCircle = (a, b, c, d) => {
  const adx = X[a] - X[d], ady = Y[a] - Y[d], bdx = X[b] - X[d], bdy = Y[b] - Y[d], cdx = X[c] - X[d], cdy = Y[c] - Y[d];
  return (adx * adx + ady * ady) * (bdx * cdy - cdx * bdy) - (bdx * bdx + bdy * bdy) * (adx * cdy - cdx * ady) + (cdx * cdx + cdy * cdy) * (adx * bdy - bdx * ady);
};

Triangle that holds a point

64-bit

Walks from the last triangle to the one a new point falls in.

The walk crosses an edge when the point is outside it by more than epsA. After 4 × (length of the corner array) steps without an answer it tests every triangle with the same tolerance.

  • acis/cdt.js
  • lines 41 to 57
  • CadCraft 0.8.2
const locate = (p) => {
  let t = last;
  for (let steps = 0; steps < 4 * V.length; steps++) {
    let moved = false;
    for (let k = 0; k < 3; k++) {
      const a = V[t * 3 + k], b = V[t * 3 + (k + 1) % 3];
      if (orient(a, b, p) < -epsA && N[t * 3 + k] >= 0) { t = N[t * 3 + k]; moved = true; break; }
    }
    if (!moved) return t;
  }
  // walking failed (rounding): look at every triangle
  for (let t2 = 0; t2 < V.length / 3; t2++) {
    const a = V[t2 * 3], b = V[t2 * 3 + 1], c = V[t2 * 3 + 2];
    if (orient(a, b, p) >= -epsA && orient(b, c, p) >= -epsA && orient(c, a, p) >= -epsA) return t2;
  }
  return -1;
};

Inserting a point

64-bit

Adds one point: merged into a corner it coincides with, splitting an edge it lies on, or splitting the triangle it is inside.

A point within eps of a triangle corner in both x and y is merged into that corner. A point whose orientation value against an edge is within epsA splits that edge and the triangle across it; otherwise the triangle is split in three.

  • acis/cdt.js
  • lines 98 to 136
  • CadCraft 0.8.2
const add = (p) => {
  const t = locate(p);
  if (t < 0) return false;
  // on top of a corner?
  for (let k = 0; k < 3; k++) { const q = V[t * 3 + k]; if (Math.abs(X[q] - X[p]) <= eps && Math.abs(Y[q] - Y[p]) <= eps) { alias[p] = q < n ? q : -1; return true; } }
  const a = V[t * 3], b = V[t * 3 + 1], c = V[t * 3 + 2];
  const nA = N[t * 3], nB = N[t * 3 + 1], nC = N[t * 3 + 2];
  // on an edge? split it and the triangle across
  let onK = -1;
  for (let k = 0; k < 3; k++) if (Math.abs(orient(V[t * 3 + k], V[t * 3 + (k + 1) % 3], p)) <= epsA) { onK = k; break; }
  if (onK >= 0 && N[t * 3 + onK] >= 0) {
    const k = onK;
    const ea = V[t * 3 + k], eb = V[t * 3 + (k + 1) % 3], ec = V[t * 3 + (k + 2) % 3];
    const u = N[t * 3 + k], ku = edgeOf(u, eb), ed = V[u * 3 + (ku + 2) % 3];
    const t_ca = N[t * 3 + (k + 2) % 3], t_bc = N[t * 3 + (k + 1) % 3];
    const u_ad = N[u * 3 + (ku + 1) % 3], u_db = N[u * 3 + (ku + 2) % 3];
    const t2 = V.length / 3, u2 = t2 + 1;
    // t = (ea, p, ec), t2 = (p, eb, ec), u = (eb, p, ed), u2 = (p, ea, ed)
    V[t * 3] = ea; V[t * 3 + 1] = p; V[t * 3 + 2] = ec; N[t * 3] = u2; N[t * 3 + 1] = t2; N[t * 3 + 2] = t_ca;
    V.push(p, eb, ec); N.push(u, t_bc, t);
    V[u * 3] = eb; V[u * 3 + 1] = p; V[u * 3 + 2] = ed; N[u * 3] = t2; N[u * 3 + 1] = u2; N[u * 3 + 2] = u_db;
    V.push(p, ea, ed); N.push(t, u_ad, u);
    setN(t_bc, ec, eb, t2); setN(u_ad, ed, ea, u2);
    // a boundary edge split in two stays boundary in its halves
    if (fixed.has(ek(ea, eb))) { fixed.delete(ek(ea, eb)); fixed.add(ek(ea, p)); fixed.add(ek(p, eb)); }
    last = t;
    legalize([[t, p], [t2, p], [u, p], [u2, p]]);
    return true;
  }
  const t2 = V.length / 3, t3 = t2 + 1;
  // t = (a, b, p), t2 = (b, c, p), t3 = (c, a, p)
  V[t * 3] = a; V[t * 3 + 1] = b; V[t * 3 + 2] = p; N[t * 3] = nA; N[t * 3 + 1] = t2; N[t * 3 + 2] = t3;
  V.push(b, c, p); N.push(nB, t3, t);
  V.push(c, a, p); N.push(nC, t, t2);
  setN(nB, c, b, t2); setN(nC, a, c, t3);
  last = t;
  legalize([[t, p], [t2, p], [t3, p]]);
  return true;
};

Do two segments cross

64-bit

Whether edge c–d properly crosses boundary segment a–b.

All four orientation values must be beyond ±epsA with opposite signs in each pair. Segments that touch at an end or are nearly collinear do not count as crossing.

  • acis/cdt.js
  • lines 143 to 146
  • CadCraft 0.8.2
const crosses = (a, b, c, d) => {
  const o1 = orient(a, b, c), o2 = orient(a, b, d), o3 = orient(c, d, a), o4 = orient(c, d, b);
  return ((o1 > epsA && o2 < -epsA) || (o1 < -epsA && o2 > epsA)) && ((o3 > epsA && o4 < -epsA) || (o3 < -epsA && o4 > epsA));
};

In detail

Method
Constrained Delaunay triangulation: points inserted one by one with edge flips (Lawson), boundary segments forced in by flipping crossing edges (Sloan, 1993, as the file's header cites), inside triangles kept by parity across boundary segments.
Number types
Float64Array coordinates; plain arrays of integer indexes.
Tolerances
eps = 1e-13 × D (coincident points); epsA = 1e-13 × D² (orientation, twice a triangle's area). D is the larger side of the bounding box, or 1 when that is 0. The in-circle test has no tolerance.
Limits
Not exact arithmetic: for nearly collinear or nearly cocircular points the sign of a determinant can be decided by rounding. The point walk falls back to testing every triangle; flip loops stop at 1e7 steps (legalize) and 100000 steps per boundary segment.

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.

Every policy that takes this approach