Resources Numerical accuracy Calculation policy

Surface mesh from a distance grid (dual contouring)

Turns a distance grid into a closed triangle mesh, one vertex per surface cell.

How the precision is kept

Edge crossings, normals and the least-squares vertex are computed in 64-bit from the grid's 32-bit node values and collected in a Float64Array. The finished positions are then converted to a Float32Array, and the bounding box is taken from the converted values. Connectivity is decided by signs and whole-number tables only.

Open item. Return the positions as the Float64Array they were computed in.

Calculations under this policy

5 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.

field/contour.js PartCraft 0.8.1 32-bit in part 5 calculations Reviewed Nov 6, 2026

Tests tools/field-test/watertight.mjs

Edge crossing and its normal

32-bit in part

Where the surface crosses a cell edge, and the surface direction there.

t = vᵢ ÷ (vᵢ − vₒ); p = p_in + (p_out − p_in) × t

The corner values come from the grid's Float32Array bricks. vᵢ is below 0 and vₒ is not, so the divisor is never zero. The normal is the trilinear gradient divided by its length; below 1e-12 the edge direction is used.

field/contour.js
lines 190 to 196
PartCraft 0.8.1

const vi = cv[inC], vo = cv[outC], t = vi / (vi - vo);
const x0 = (i + (inC & 1)) * h, y0 = (j + (inC >> 1 & 1)) * h, z0 = (k + (inC >> 2 & 1)) * h;
const x1 = (i + (outC & 1)) * h, y1 = (j + (outC >> 1 & 1)) * h, z1 = (k + (outC >> 2 & 1)) * h;
const px = x0 + (x1 - x0) * t, py = y0 + (y1 - y0) * t, pz = z0 + (z1 - z0) * t;
let g = cellGradient(cv, (px - i * h) / h, (py - j * h) / h, (pz - k * h) / h); let gl = Math.sqrt(g[0] * g[0] + g[1] * g[1] + g[2] * g[2]);
if (!(gl > 1e-12)) { g = [x1 - x0, y1 - y0, z1 - z0]; gl = h; }
P.push(px, py, pz); N.push(g[0] / gl, g[1] / gl, g[2] / gl);

Vertex of a surface cell

64-bit

The point that best fits the crossing planes of one cell.

(AᵀA + λI)y = Aᵀb, x = m + y, λ = 0.1

The normal equations are summed in 64-bit relative to the mean crossing m and solved by Cramer's rule. λ on the diagonal keeps the system well posed; if the determinant is not above 1e-12 the mean is used. The vertex is clamped to its cell and kept in a Float64Array at this stage.

field/contour.js
lines 202 to 216
PartCraft 0.8.1

let a00 = LAMBDA, a01 = 0, a02 = 0, a11 = LAMBDA, a12 = 0, a22 = LAMBDA, b0 = 0, b1 = 0, b2 = 0;
for (let e = 0; e < n; e++) {
  const nx = N[e * 3], ny = N[e * 3 + 1], nz = N[e * 3 + 2], d = nx * (P[e * 3] - mx) + ny * (P[e * 3 + 1] - my) + nz * (P[e * 3 + 2] - mz);
  a00 += nx * nx; a01 += nx * ny; a02 += nx * nz; a11 += ny * ny; a12 += ny * nz; a22 += nz * nz;
  b0 += nx * d; b1 += ny * d; b2 += nz * d;
}
const det = a00 * (a11 * a22 - a12 * a12) - a01 * (a01 * a22 - a12 * a02) + a02 * (a01 * a12 - a11 * a02);
let x = mx, y = my, z = mz;
if (det > 1e-12) {
  x += (b0 * (a11 * a22 - a12 * a12) - a01 * (b1 * a22 - a12 * b2) + a02 * (b1 * a12 - a11 * b2)) / det;
  y += (a00 * (b1 * a22 - a12 * b2) - b0 * (a01 * a22 - a12 * a02) + a02 * (a01 * b2 - b1 * a02)) / det;
  z += (a00 * (a11 * b2 - b1 * a12) - a01 * (a01 * b2 - b1 * a02) + b0 * (a01 * a12 - a11 * a02)) / det;
}
const clamp = (v, lo) => (v < lo ? lo : v > lo + h ? lo + h : v);
pushPos(clamp(x, i * h), clamp(y, j * h), clamp(z, k * h));

Gradient inside a cell

64-bit

The gradient of the interpolated distance inside one cell.

Differences of the eight corner values weighted by the other two coordinates: the gradient of the trilinear interpolant, per voxel.

field/contour.js
lines 264 to 272
PartCraft 0.8.1

/** Gradient of the trilinear interpolant of a cell's eight corner values (corner c = x | y<<1 | z<<2) at (u, v, w) in [0,1]³, per voxel. */
function cellGradient(c, u, v, w) {
  const iu = 1 - u, iv = 1 - v, iw = 1 - w;
  return [
    (c[1] - c[0]) * iv * iw + (c[3] - c[2]) * v * iw + (c[5] - c[4]) * iv * w + (c[7] - c[6]) * v * w,
    (c[2] - c[0]) * iu * iw + (c[3] - c[1]) * u * iw + (c[6] - c[4]) * iu * w + (c[7] - c[5]) * u * w,
    (c[4] - c[0]) * iu * iv + (c[5] - c[1]) * u * iv + (c[6] - c[2]) * iu * v + (c[7] - c[3]) * u * v,
  ];
}

Which diagonal splits a quad

64-bit

Chooses how each four-sided face becomes two triangles.

use diagonal 0–2 when |p₀ − p₂|² ≤ |p₁ − p₃|²

Squared lengths of the two diagonals are compared in 64-bit, with ties going to the first. A diagonal that is already an edge is avoided so that no edge gets three triangles.

field/contour.js
lines 155 to 161
PartCraft 0.8.1

const D2 = (p, q) => (pa[p * 3] - pa[q * 3]) ** 2 + (pa[p * 3 + 1] - pa[q * 3 + 1]) ** 2 + (pa[p * 3 + 2] - pa[q * 3 + 2]) ** 2;
const finish = (Q) => {                                  // choose the diagonals of a slab's quads, in order, and emit the triangles
  const { qa, qn, slab } = Q, near = [slab - 1, slab, slab + 1], dnear = [slab - 1, slab], dset = setOf(diags, slab);
  for (let q = 0; q < qn; q += 4) {
    const o0 = qa[q], o1 = qa[q + 1], o2 = qa[q + 2], o3 = qa[q + 3];
    const bad02 = hasIn(perims, near, o0, o2) || hasIn(diags, dnear, o0, o2), bad13 = hasIn(perims, near, o1, o3) || hasIn(diags, dnear, o1, o3);
    const use02 = bad02 === bad13 ? D2(o0, o2) <= D2(o1, o3) : !bad02;

Positions and box as returned

32-bit in part

The mesh's corner positions and its bounding box.

The 64-bit positions are converted to a Float32Array here, and the bounding box is read from the converted array. A coordinate near 100 mm is then held to about 4e-6 mm.

field/contour.js
lines 258 to 261
PartCraft 0.8.1

const positions = Float32Array.from(pa.subarray(0, pn)), indices = ta.slice(0, tn);
const bbox = { min: [Infinity, Infinity, Infinity], max: [-Infinity, -Infinity, -Infinity] };
for (let n = 0; n < positions.length; n += 3) for (let a = 0; a < 3; a++) { if (positions[n + a] < bbox.min[a]) bbox.min[a] = positions[n + a]; if (positions[n + a] > bbox.max[a]) bbox.max[a] = positions[n + a]; }
return { positions, indices, bbox };

In detail

Method: Dual contouring (the file cites Ju, Losasso, Schaefer, Warren 2002): the vertex minimizes Σ(nᵢ·(x − pᵢ))² + λ|x − m|², solved by Cramer's rule and clamped to its cell. One quad per sign-changing edge, split along the shorter diagonal unless that edge already exists.
Number types: Float32Array node values in; Float64Array work buffer; Float32Array positions and Uint32Array indices out.
Tolerances: λ = 0.1 pulls the vertex toward the mean crossing. A gradient shorter than 1e-12 falls back to the edge direction. The 3 × 3 solve is skipped when the determinant is not above 1e-12.
Limits: The mesh follows the grid, so its accuracy is set by the voxel size h; positions are stored in 32-bit.

Approach: 32-bit step, flagged

One step stores or exchanges 32-bit values that reach a number a person reads. It is listed here with its open item, and the step is either moved to 64-bit or the app will warn where the number is read.

Every policy that takes this approach