CadCraft™ v0.8.1
Use it now
Menu

Resources Numerical accuracy Calculation policy

Distance grid: baking, sampling and volume

Stores a shape as signed distances on a sparse voxel lattice, built from a formula or from a closed triangle mesh, and reads values, gradients and the enclosed volume back from it.

How the precision is kept

Distances are computed in 64-bit and then stored in Float32Array bricks of 8 × 8 × 8 nodes, clamped to ± band × h. Everything read back (samples, gradients, the volume, the contour mesh and the exports) starts from those 32-bit values. A stored value is at most band × h in size, so its 32-bit rounding is at most about 6e-8 of band × h; the method's own error is set by the voxel size h.

Open item. Store brick values in a Float64Array, or state the 32-bit grid next to the voxel size as the method's resolution, so sampled distances, the volume and the contour do not start from 32-bit numbers.

Calculations under this policy

14 calculations in 4 source files. 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/bake.js PartCraft 0.8.1 32-bit in part 2 calculations Reviewed

Tests tools/field-test/volumes.mjs

Bricks to compute

64-bit

Which bricks of the lattice can hold the surface.

brick range = floor((box ∓ band × h) ÷ h ÷ 8); ball radius = 3.5√3 × h

Brick indices are floors of the box limits divided by h and by 8, in 64-bit. Each brick is classified with a ball around its center at (index × 8 + 3.5) × h; only bricks that may hold the surface are computed.

  • field/bake.js
  • lines 16 to 30
  • PartCraft 0.8.1
export function plan(spec, { h, bbox = null, band = 3 } = {}) {
  const node = compile(spec);
  const box = bbox || node.bounds;
  if (!box.every(Number.isFinite)) throw new FieldError('nobox', 'field: this field has no finite box; pass bbox or intersect it with a box');
  const pad = band * h;
  const lo = [0, 1, 2].map((a) => Math.floor((box[a] - pad) / h / 8)), hi = [0, 1, 2].map((a) => Math.floor((box[a + 3] + pad) / h / 8));
  for (let a = 0; a < 3; a++) if (!inRange(lo[a]) || !inRange(hi[a])) throw new FieldError('toolarge', 'field: part too large for the grid (more than 16 000 bricks along an axis)');
  const keys = [], inside = [], r = R_BRICK * h;
  for (let bx = lo[0]; bx <= hi[0]; bx++) for (let by = lo[1]; by <= hi[1]; by++) for (let bz = lo[2]; bz <= hi[2]; bz++) {
    const cls = node.cls((bx * 8 + 3.5) * h, (by * 8 + 3.5) * h, (bz * 8 + 3.5) * h, r, band * h);
    if (cls === OUT) continue;
    if (cls === IN) inside.push(keyOf(bx, by, bz)); else keys.push(keyOf(bx, by, bz));
  }
  return { h, band, keys, inside, bricks: (hi[0] - lo[0] + 1) * (hi[1] - lo[1] + 1) * (hi[2] - lo[2] + 1) };
}

Node values of a brick

32-bit in part

The signed distance at each of a brick's 512 nodes.

value = clamp(f(i·h, j·h, k·h), −band × h, band × h)

The field is evaluated at whole-number multiples of h in 64-bit, clamped to ± band × h, then stored in a Float32Array. Two grids with the same h share exactly the same sample points.

  • field/bake.js
  • lines 32 to 46
  • PartCraft 0.8.1
export function bakeKeys(spec, h, band, keys) {
  const f = compile(spec).f, bw = band * h, out = [];
  for (const key of keys) {
    const bz = Math.floor(key % 131072), rest = (key - bz) / 131072, by = rest % 131072, bx = (rest - by) / 131072;
    const x0 = (bx - 65536) * 8, y0 = (by - 65536) * 8, z0 = (bz - 65536) * 8;
    const arr = new Float32Array(BN);
    let n = 0;
    for (let k = 0; k < 8; k++) for (let j = 0; j < 8; j++) for (let i = 0; i < 8; i++, n++) {
      let v = f((x0 + i) * h, (y0 + j) * h, (z0 + k) * h);
      arr[n] = v > bw ? bw : v < -bw ? -bw : v;
    }
    out.push(arr);
  }
  return out;
}
field/grid.js PartCraft 0.8.1 32-bit in part 5 calculations Reviewed

Tests tools/field-test/volumes.mjs

Brick key

Exact

One whole number that names a brick.

key = ((bx + 65536) × 131072 + (by + 65536)) × 131072 + (bz + 65536)

A whole number below 2⁵¹, held exactly in a 64-bit number; brickOf undoes it with remainders and exact divisions.

  • field/grid.js
  • line 16
  • PartCraft 0.8.1
export const keyOf = (bx, by, bz) => ((bx + OFF) * S + (by + OFF)) * S + (bz + OFF);

Value at a point

32-bit in part

The signed distance at any point, between the lattice nodes.

trilinear: a + (b − a) × t along x, then y, then z

The eight corner values are read from Float32Array bricks, or are ± band × h for absent and inside bricks. The interpolation itself is 64-bit.

  • field/grid.js
  • lines 37 to 47
  • PartCraft 0.8.1
/** Trilinear sample at a world point. */
sample(x, y, z) {
  const h = this.h, fx = x / h, fy = y / h, fz = z / h;
  const i = Math.floor(fx), j = Math.floor(fy), k = Math.floor(fz);
  const u = fx - i, v = fy - j, w = fz - k;
  const c000 = this.at(i, j, k), c100 = this.at(i + 1, j, k), c010 = this.at(i, j + 1, k), c110 = this.at(i + 1, j + 1, k);
  const c001 = this.at(i, j, k + 1), c101 = this.at(i + 1, j, k + 1), c011 = this.at(i, j + 1, k + 1), c111 = this.at(i + 1, j + 1, k + 1);
  const a0 = c000 + (c100 - c000) * u, a1 = c010 + (c110 - c010) * u, b0 = c001 + (c101 - c001) * u, b1 = c011 + (c111 - c011) * u;
  const d0 = a0 + (a1 - a0) * v, d1 = b0 + (b1 - b0) * v;
  return d0 + (d1 - d0) * w;
}

Gradient at a point

32-bit in part

The direction and rate of change of the distance.

g = (f(p + h/2) − f(p − h/2)) ÷ h, per axis

The difference of two samples one voxel apart, divided by h. Both samples come from 32-bit node values, so the difference of two nearly equal numbers carries their 32-bit rounding.

  • field/grid.js
  • lines 48 to 52
  • PartCraft 0.8.1
/** Gradient (central differences over one voxel) → [gx, gy, gz]. */
gradient(x, y, z) {
  const e = this.h * 0.5;
  return [this.sample(x + e, y, z) - this.sample(x - e, y, z), this.sample(x, y + e, z) - this.sample(x, y - e, z), this.sample(x, y, z + e) - this.sample(x, y, z - e)].map((g) => g / (2 * e));
}

Box of the stored bricks

64-bit

The world box that holds every stored node.

lo = brick × 8 × h; hi = (brick × 8 + 7) × h

Whole numbers times h in 64-bit.

  • field/grid.js
  • lines 64 to 70
  • PartCraft 0.8.1
/** World box of the stored bricks (node positions). */
worldBox() {
  const e = this.extent();
  if (!e) return null;
  const h = this.h;
  return [e.lo[0] * 8 * h, e.lo[1] * 8 * h, e.lo[2] * 8 * h, (e.hi[0] * 8 + 7) * h, (e.hi[1] * 8 + 7) * h, (e.hi[2] * 8 + 7) * h];
}

Enclosed volume

32-bit in part

The volume inside the surface.

V = h³ × Σ clamp(½ − d ÷ h, 0, 1)

Node values are read from Float32Array bricks; the fractions are added in brick-key order in 64-bit and multiplied by h³ once at the end. A wholly inside brick adds 512. The file's comment: exact for a flat surface, second-order for a curved one.

  • field/grid.js
  • lines 80 to 93
  • PartCraft 0.8.1
/**
 * Volume enclosed (world units³): every node counts its share of a voxel, clamp(½ − d/h, 0, 1) — exact for a flat surface,
 * second-order for a curved one. INSIDE bricks count whole.
 */
volume() {
  const h = this.h, h3 = h * h * h;
  let v = 0;
  for (const k of this.sortedKeys()) {
    const b = this.bricks.get(k);
    if (b === INSIDE) { v += BN; continue; }
    for (let n = 0; n < BN; n++) { const f = 0.5 - b[n] / h; v += f <= 0 ? 0 : f >= 1 ? 1 : f; }
  }
  return v * h3;
}
field/mesh-sdf.js PartCraft 0.8.1 32-bit in part 5 calculations Reviewed

Tests tools/field-test/volumes.mjs

Welding corners

64-bit

Merges mesh corners that sit at the same position.

Two corners merge only when their three coordinates are identical; there is no distance tolerance. The welded positions are copied into a Float64Array.

  • field/mesh-sdf.js
  • lines 14 to 23
  • PartCraft 0.8.1
function weld(mesh) {
  const P = mesh.positions, I = mesh.indices, map = new Map(), pos = [], idx = new Int32Array(I.length);
  for (let n = 0; n < I.length; n++) {
    const v = I[n], x = P[v * 3], y = P[v * 3 + 1], z = P[v * 3 + 2], key = x + ',' + y + ',' + z;
    let w = map.get(key);
    if (w === undefined) { w = pos.length / 3; map.set(key, w); pos.push(x, y, z); }
    idx[n] = w;
  }
  return { P: Float64Array.from(pos), I: idx };
}

Angle-weighted normals

64-bit

The pseudo-normal at each corner and edge, used for the sign.

angle = acos(e₁·e₂ ÷ (|e₁||e₂|)); N_vertex = Σ angle × n_face

The cosine is clamped to [−1, 1] before acos. Triangles of zero area are skipped first, so the two edge lengths in the divisor are not zero. The comment notes that any monotone weight gives the same sign.

  • field/mesh-sdf.js
  • lines 46 to 48
  • PartCraft 0.8.1
const cs = (e1x * e2x + e1y * e2y + e1z * e2z) / (l1 * l2);
const ang = Math.acos(cs < -1 ? -1 : cs > 1 ? 1 : cs);        // only the weight of the vertex normal; any monotone weight is fine for the sign
VN[p0 * 3] += ang * nx; VN[p0 * 3 + 1] += ang * ny; VN[p0 * 3 + 2] += ang * nz;

Signed distance at a node

32-bit in part

The distance from a lattice node to the mesh, negative inside.

d = √(least squared distance); sign = sign((p − closest)·N)

The least squared distance over the brick's triangles is found in 64-bit and one square root is taken. The sign is the dot product of the vector from the closest point with the pseudo-normal of the closest face, edge or corner. The clamped value is stored in a Float32Array.

  • field/mesh-sdf.js
  • lines 83 to 101
  • PartCraft 0.8.1
const d = Math.sqrt(best);
// pseudo-normal of the feature the closest point lies on
const a = I[bt * 3], b = I[bt * 3 + 1], c = I[bt * 3 + 2];
let nx, ny, nz;
switch (bq) {
  case 1: nx = VN[a * 3]; ny = VN[a * 3 + 1]; nz = VN[a * 3 + 2]; break;
  case 2: nx = VN[b * 3]; ny = VN[b * 3 + 1]; nz = VN[b * 3 + 2]; break;
  case 3: nx = VN[c * 3]; ny = VN[c * 3 + 1]; nz = VN[c * 3 + 2]; break;
  case 4: { const e = EN.get(ekey(a, b)); nx = e[0]; ny = e[1]; nz = e[2]; break; }
  case 5: { const e = EN.get(ekey(b, c)); nx = e[0]; ny = e[1]; nz = e[2]; break; }
  case 6: { const e = EN.get(ekey(c, a)); nx = e[0]; ny = e[1]; nz = e[2]; break; }
  default: nx = FN[bt * 3]; ny = FN[bt * 3 + 1]; nz = FN[bt * 3 + 2];
}
// vector from the closest point to the node = p − closest; recompute the closest point once for the winner
closest(px, py, pz, P, a * 3, b * 3, c * 3, q);
const dx = px - q[4], dy = py - q[5], dz = pz - q[6];
const sgn = dx * nx + dy * ny + dz * nz;
const v = sgn < 0 ? -d : d;
arr[n] = v > bw ? bw : v < -bw ? -bw : v;

Closest point on a triangle

64-bit

The nearest point of a triangle to a point, and which face, edge or corner it lies on.

Region tests on dot products of edge vectors; edge and face cases divide by differences of those dot products. The face case divides by va + vb + vc without a check; triangles of zero area are skipped before they get here.

  • field/mesh-sdf.js
  • lines 144 to 175
  • PartCraft 0.8.1
/** out = [d², region, 0, 0, cx, cy, cz]; region 0 face, 1/2/3 vertex a/b/c, 4/5/6 edge ab/bc/ca. */
function closest(px, py, pz, P, a, b, c, out) {
  const ax = P[a], ay = P[a + 1], az = P[a + 2], bx = P[b], by = P[b + 1], bz = P[b + 2], cx = P[c], cy = P[c + 1], cz = P[c + 2];
  const abx = bx - ax, aby = by - ay, abz = bz - az, acx = cx - ax, acy = cy - ay, acz = cz - az;
  const apx = px - ax, apy = py - ay, apz = pz - az;
  const d1 = abx * apx + aby * apy + abz * apz, d2 = acx * apx + acy * apy + acz * apz;
  let reg = 0, qx, qy, qz;
  if (d1 <= 0 && d2 <= 0) { reg = 1; qx = ax; qy = ay; qz = az; }
  else {
    const bpx = px - bx, bpy = py - by, bpz = pz - bz, d3 = abx * bpx + aby * bpy + abz * bpz, d4 = acx * bpx + acy * bpy + acz * bpz;
    if (d3 >= 0 && d4 <= d3) { reg = 2; qx = bx; qy = by; qz = bz; }
    else {
      const vc = d1 * d4 - d3 * d2;
      if (vc <= 0 && d1 >= 0 && d3 <= 0) { reg = 4; const v = d1 / (d1 - d3); qx = ax + v * abx; qy = ay + v * aby; qz = az + v * abz; }
      else {
        const cpx = px - cx, cpy = py - cy, cpz = pz - cz, d5 = abx * cpx + aby * cpy + abz * cpz, d6 = acx * cpx + acy * cpy + acz * cpz;
        if (d6 >= 0 && d5 <= d6) { reg = 3; qx = cx; qy = cy; qz = cz; }
        else {
          const vb = d5 * d2 - d1 * d6;
          if (vb <= 0 && d2 >= 0 && d6 <= 0) { reg = 6; const w = d2 / (d2 - d6); qx = ax + w * acx; qy = ay + w * acy; qz = az + w * acz; }
          else {
            const va = d3 * d6 - d5 * d4;
            if (va <= 0 && d4 - d3 >= 0 && d5 - d6 >= 0) { reg = 5; const w = (d4 - d3) / ((d4 - d3) + (d5 - d6)); qx = bx + w * (cx - bx); qy = by + w * (cy - by); qz = bz + w * (cz - bz); }
            else { const den = 1 / (va + vb + vc), v = vb * den, w = vc * den; qx = ax + abx * v + acx * w; qy = ay + aby * v + acy * w; qz = az + abz * v + acz * w; }
          }
        }
      }
    }
  }
  const dx = px - qx, dy = py - qy, dz = pz - qz;
  out[0] = dx * dx + dy * dy + dz * dz; out[1] = reg; out[4] = qx; out[5] = qy; out[6] = qz;
}

Inside or outside far from the surface

64-bit

Marks bricks no triangle comes near as inside or outside.

x_cross = (w₀a_x + w₁b_x + w₂c_x) ÷ (w₀ + w₁ + w₂); inside when an odd number of crossings lie before the brick

One ray along +x per row of bricks. A triangle counts when its three edge functions have the same strict sign, so the divisor is not zero. The ray is shifted by a few thousandths of a voxel, up to three times, until the number of crossings is even; a row that never gets an even count is left unmarked.

  • field/mesh-sdf.js
  • lines 121 to 129
  • PartCraft 0.8.1
const y = (by * 8 + 3.5 + [1.2345e-3, -3.14159e-3, 5.7721e-3][attempt]) * h, z = (bz * 8 + 3.5 + [2.71828e-3, 1.41421e-3, -4.6692e-3][attempt]) * h;
xs.length = 0;
for (const t of l) {
  const a = I[t * 3] * 3, b = I[t * 3 + 1] * 3, c = I[t * 3 + 2] * 3;
  const ay = P[a + 1] - y, az = P[a + 2] - z, byy = P[b + 1] - y, bzz = P[b + 2] - z, cy = P[c + 1] - y, cz = P[c + 2] - z;
  const w0 = byy * cz - bzz * cy, w1 = cy * az - cz * ay, w2 = ay * bzz - az * byy;
  if ((w0 > 0 && w1 > 0 && w2 > 0) || (w0 < 0 && w1 < 0 && w2 < 0)) xs.push((w0 * P[a] + w1 * P[b] + w2 * P[c]) / (w0 + w1 + w2));
}
ok = xs.length % 2 === 0;
field/field.js PartCraft 0.8.1 32-bit in part 2 calculations Reviewed

Tests None found by the review

Density map as stored

32-bit in part

Keeps a density map that steers a wall thickness or cell size.

The map's values are copied into a Float32Array when the shape is described.

  • field/field.js
  • lines 80 to 84
  • PartCraft 0.8.1
/**
 * A number that follows a density map: `data` (nx·ny·nz values in 0…1 at cell centres, x fastest) over the box [x0,y0,z0,x1,y1,z1], linearly interpolated
 * between centres and held at the edge value outside; the profile is v0 + (v1 − v0)·density. For a wall that thickens where a topology result is solid.
 */
export const density = (data, n, box, v0, v1) => ({ t: 'dens', data: Float32Array.from(data), n: [...n], box: [...box], v0, v1 });

Number that follows a density map

32-bit in part

The wall thickness or cell size at a point, from the density map.

v = v₀ + (v₁ − v₀) × density(p), density by trilinear interpolation between cell centers

The eight map values are read from the Float32Array and interpolated in 64-bit; outside the box the edge value is held. Then one linear blend between v₀ and v₁.

  • field/field.js
  • lines 93 to 103
  • PartCraft 0.8.1
if (s.t === 'dens') {
  const [nx, ny, nz] = s.n, D = s.data, [x0, y0, z0, x1, y1, z1] = s.box, v0 = s.v0, v1 = s.v1;
  const ax = (x, lo, hi, n) => { const u = n > 1 ? (x - lo) / (hi - lo) * n - 0.5 : 0; const c = u < 0 ? 0 : u > n - 1 ? n - 1 : u; const i = Math.min(n - 2 < 0 ? 0 : n - 2, Math.floor(c)); return [i, n > 1 ? c - i : 0]; };
  return (x, y, z) => {
    const [i, fx] = ax(x, x0, x1, nx), [j, fy] = ax(y, y0, y1, ny), [k, fz] = ax(z, z0, z1, nz);
    const i1 = nx > 1 ? i + 1 : 0, j1 = ny > 1 ? j + 1 : 0, k1 = nz > 1 ? k + 1 : 0, at = (a, b, c) => D[a + nx * (b + ny * c)];
    const d00 = at(i, j, k) * (1 - fx) + at(i1, j, k) * fx, d10 = at(i, j1, k) * (1 - fx) + at(i1, j1, k) * fx, d01 = at(i, j, k1) * (1 - fx) + at(i1, j, k1) * fx, d11 = at(i, j1, k1) * (1 - fx) + at(i1, j1, k1) * fx;
    const d = (d00 * (1 - fy) + d10 * fy) * (1 - fz) + (d01 * (1 - fy) + d11 * fy) * fz;
    return v0 + (v1 - v0) * d;
  };
}

In detail

Method
Nodes on the world lattice x = i × h. Baking samples the compiled field at each node. From a mesh: point to triangle distance by Voronoi regions (the file cites Ericson, Real-Time Collision Detection 5.1.5) with the sign from angle-weighted pseudo-normals (Bærentzen and Aanæs 2005). Sampling is trilinear; the gradient is a central difference over one voxel; the volume counts clamp(½ − d ÷ h, 0, 1) per node.
Number types
64-bit arithmetic; Float32Array bricks; whole-number brick keys of 51 bits held exactly in a 64-bit number.
Tolerances
Band of 3 voxels by default. A brick is skipped when a ball of 3.5√3 voxels around its center is wholly outside or inside the band. Ray parity for far bricks tries up to three slightly shifted rays.
Limits
grid.js says the volume rule is exact for a flat surface and second-order for a curved one. A density map is also kept in a Float32Array. Edge keys in mesh-sdf.js (a × 2³² + b) are exact only below 2²¹ welded corners.

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