From cbea030834842d2a27dd256ee625c30e9f3270a1 Mon Sep 17 00:00:00 2001 From: Kevin Blackburn-Matzen Date: Sun, 2 Aug 2026 18:09:32 -0700 Subject: [PATCH] Stop generating the triangles the simplifier only has to take away MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Dual contouring emits one vertex per surface cell, so a flat wall 64 voxels across becomes thousands of triangles that QEM then collapses back to two. Profiling for #111 put simplification at ~80% of an export at res 256, and 93% of the triangles it removed carried *zero* error — the mesher undoing its own over-generation. Octree vertex clustering decides, before any triangle exists, which cells can share a vertex: a node collapses when one point still fits every hermite sample underneath it within the budget. The reason this is safe is the point of the change. Adaptive dual contouring (Ju et al. 2002) contours the octree directly, which needs cellProc/faceProc/edgeProc recursion to stitch differently-sized neighbours and produces holes in an exported STL when that is wrong — which is why #111 was filed as "not small". None of that machinery is needed. Step 2 still emits exactly one quad per sign-changing grid edge over the same four cells; only the vertex those cells name changes. The output is a quotient of a mesh that was already watertight, so there are no cracks to patch. What clustering can break is manifoldness, when two surface sheets land in one node. A cell with more than one patch never joins a cluster, and — deliberately more conservative — it stays keyed in place so it fails every node above it too, keeping merged vertices out of regions holding geometry they do not describe. At res 256, the resolution exports actually run at: Wall Bracket 14529 -> 6550 ms 16478 -> 3274 tris maxdev 0.071 -> 0.071 Open-Top Enclosure 22884 -> 14438 ms 50694 -> 21164 tris maxdev 0.046 -> 0.046 Knob 15626 -> 9474 ms 68654 -> 58312 tris maxdev 0.016 -> 0.018 M3 Standoff 70504 -> 6546 ms 47710 -> 32256 tris maxdev 0.077 -> 0.077 Vent Grid 80850 -> 14614 ms 104548 -> 70708 tris maxdev 0.044 -> 0.045 Max deviation — what matters for a print — is unchanged on three and moves by 0.001-0.002 of a voxel on the two curved ones. Every output is watertight with zero non-manifold edges. The timings were taken on a machine with other work on it, so treat the two largest as unconfirmed; CI's bench job measures this cleanly. The triangle counts are deterministic. Off by default, and the dense path is verified bit-identical to before across every preset, with and without the active-block mask, so the geometry suite still asserts against exactly the mesh it always did. `cluster.test.ts` covers the collapse rule and then re-applies dualContour's own gates to the clustered mesh — watertight, edge-manifold, outward-wound, correct enclosed volume, vertices on the isosurface, and the two-sheets-in-one-cell case. A 10mm cube at res 32 comes out as exactly 12 triangles enclosing exactly 1000.00 units of volume, where the dense mesh spends 4800. The budgets move to `budgets.ts` because the benchmark had already drifted from the worker it claims to measure. Clustering gets half the allowance and QEM the other half: both stages move vertices and their errors add. Closes part of #111. Co-Authored-By: Claude Opus 5 (1M context) --- bench/mesh.bench.ts | 7 +- src/worker/sdf/budgets.ts | 27 +++ src/worker/sdf/cluster.test.ts | 307 +++++++++++++++++++++++++++++++++ src/worker/sdf/cluster.ts | 253 +++++++++++++++++++++++++++ src/worker/sdf/dualContour.ts | 114 +++++++++++- src/worker/sdfWorker.ts | 15 +- 6 files changed, 710 insertions(+), 13 deletions(-) create mode 100644 src/worker/sdf/budgets.ts create mode 100644 src/worker/sdf/cluster.test.ts create mode 100644 src/worker/sdf/cluster.ts diff --git a/bench/mesh.bench.ts b/bench/mesh.bench.ts index 4694d5d..eee8a95 100644 --- a/bench/mesh.bench.ts +++ b/bench/mesh.bench.ts @@ -30,6 +30,7 @@ import { evaluateCPUWithProgress } from '../src/worker/sdf/gridEval'; import { dualContour } from '../src/worker/sdf/dualContour'; import { simplifyMesh } from '../src/worker/sdf/simplify'; import { removeDegenerateTriangles, projectVerticesToSurface } from '../src/worker/sdf/meshRepair'; +import { CLUSTER_ERROR_VOXELS, SIMPLIFY_ERROR_VOXELS, PROJECT_TOLERANCE_VOXELS } from '../src/worker/sdf/budgets'; import { PRESET_CATEGORIES } from '../src/components/tree/presets'; import type { BBox } from '../src/worker/sdf/types'; import type { SDFNodeUI } from '../src/types/operations'; @@ -175,17 +176,17 @@ function run(name: string, tree: SDFNodeUI): StageTimes { let raw!: ReturnType; const contour = ms(() => { - raw = dualContour(evaluated.grid, RES, bbox, root, undefined, NO_MASK ? undefined : evaluated.active); + raw = dualContour(evaluated.grid, RES, bbox, root, undefined, NO_MASK ? undefined : evaluated.active, voxel * CLUSTER_ERROR_VOXELS); }); const cleaned = removeDegenerateTriangles(raw); const before = cleaned.indices.length / 3; let simplified!: typeof cleaned; - const simplify = ms(() => { simplified = simplifyMesh(cleaned, { maxError: voxel * 0.05 }); }); + const simplify = ms(() => { simplified = simplifyMesh(cleaned, { maxError: voxel * SIMPLIFY_ERROR_VOXELS }); }); let projected!: typeof simplified; - const project = ms(() => { projected = projectVerticesToSurface(simplified, root, voxel * 0.5); }); + const project = ms(() => { projected = projectVerticesToSurface(simplified, root, voxel * PROJECT_TOLERANCE_VOXELS); }); return { name, diff --git a/src/worker/sdf/budgets.ts b/src/worker/sdf/budgets.ts new file mode 100644 index 0000000..e7dbbaf --- /dev/null +++ b/src/worker/sdf/budgets.ts @@ -0,0 +1,27 @@ +/** + * How far the mesher is allowed to move a vertex, as a fraction of one voxel. + * + * These live here rather than at their use sites because two places run this + * pipeline: `sdfWorker.ts`, which is what a user's export actually executes, + * and `bench/mesh.bench.ts`, which claims to measure it. Those had already + * drifted — the bench went on timing a stage configuration the exporter no + * longer used, and reported it as the export's cost. A benchmark that measures + * something nobody runs is worse than no benchmark, because it is believed. + */ + +/** + * Budget for octree vertex clustering, during contouring. + * + * Half the simplifier's, because both stages move vertices and their errors + * add: spending the whole allowance here would let QEM spend it again on top. + * Splitting it keeps the total inside what the export promises, and this is the + * cheaper half — clustering removes triangles before they exist, where QEM pays + * to collapse them afterwards. + */ +export const CLUSTER_ERROR_VOXELS = 0.025; + +/** Budget for QEM edge collapse, after contouring. */ +export const SIMPLIFY_ERROR_VOXELS = 0.05; + +/** How far `projectVerticesToSurface` may search for the isosurface. */ +export const PROJECT_TOLERANCE_VOXELS = 0.5; diff --git a/src/worker/sdf/cluster.test.ts b/src/worker/sdf/cluster.test.ts new file mode 100644 index 0000000..85c238f --- /dev/null +++ b/src/worker/sdf/cluster.test.ts @@ -0,0 +1,307 @@ +import { describe, it, expect } from 'vitest'; +import { clusterByOctree, addSample, clusterError, QEF_STRIDE } from './cluster'; +import { dualContour } from './dualContour'; +import { evaluateSDF } from './evaluate'; +import { analyzeMesh } from './meshRepair'; +import type { SDFNode, BBox, Vec3 } from './types'; +import type { MeshResult } from './marchingCubes'; + +/** + * Octree vertex clustering. + * + * The unit tests below pin the collapse rule; the contouring tests are the ones + * that matter, because the risk this feature carries is not "wrong vertex + * count" but "hole in an exported STL". They are deliberately the same gates + * `dualContour.test.ts` applies to the dense mesh — watertight, edge-manifold, + * outward-wound, enclosing the right volume — since a mesh that fails any of + * them is unprintable no matter how few triangles it has. + */ + +const RES = 2; +const cellIndex = (x: number, y: number, z: number, res = RES) => z * res * res + y * res + x; + +/** QEFs for `n` vertices, filled by `fill(i, q, offset)`. */ +function buildQ(n: number, fill: (i: number, q: Float64Array, o: number) => void): Float64Array { + const q = new Float64Array(n * QEF_STRIDE); + for (let i = 0; i < n; i++) fill(i, q, i * QEF_STRIDE); + return q; +} + +describe('clusterByOctree', () => { + /** Four cells of one flat plane — a single vertex represents all of them. */ + it('merges cells whose samples lie on one plane', () => { + const cells = [cellIndex(0, 0, 0), cellIndex(1, 0, 0), cellIndex(0, 1, 0), cellIndex(1, 1, 0)]; + const pts: Vec3[] = [[0, 0, 0], [1, 0, 0], [0, 1, 0], [1, 1, 0]]; + const q = buildQ(4, (i, arr, o) => addSample(arr, o, pts[i], [0, 0, 1])); + + const r = clusterByOctree(4, Int32Array.from(cells), q, new Uint8Array([1, 1, 1, 1]), RES, 0.01); + + expect(r.count).toBe(1); + expect(Array.from(r.remap)).toEqual([0, 0, 0, 0]); + expect(r.size[0]).toBe(4); + expect(r.level[0]).toBe(1); + }); + + /** + * Four cells where one sample sits on a *parallel* plane two units away. No + * single point fits both, so nothing merges. + * + * Parallel, not perpendicular: the first version of this test used a + * perpendicular plane and failed, because two perpendicular planes meet at an + * exact point and the QEF residual is legitimately zero. That is a corner, + * and collapsing a corner to one vertex is the right answer — it is how + * `dualContour` reproduces sharp edges in the first place. Only samples that + * genuinely cannot share a point should block a merge. + */ + it('refuses to merge when one vertex cannot share a point with the others', () => { + const cells = [cellIndex(0, 0, 0), cellIndex(1, 0, 0), cellIndex(0, 1, 0), cellIndex(1, 1, 0)]; + const q = buildQ(4, (i, arr, o) => { + if (i === 3) addSample(arr, o, [0, 0, 2], [0, 0, 1]); // plane z = 2 + else addSample(arr, o, [i, 0, 0], [0, 0, 1]); // plane z = 0 + }); + + const r = clusterByOctree(4, Int32Array.from(cells), q, new Uint8Array([1, 1, 1, 1]), RES, 0.01); + + expect(r.count).toBe(4); + }); + + /** + * The corner case the test above was originally trying (and failing) to + * express, asserted deliberately: perpendicular planes *should* collapse, + * because a single vertex represents them exactly. + */ + it('does merge a sharp corner, which one vertex represents exactly', () => { + const cells = [cellIndex(0, 0, 0), cellIndex(1, 0, 0), cellIndex(0, 1, 0), cellIndex(1, 1, 0)]; + const normals: Vec3[] = [[0, 0, 1], [0, 0, 1], [1, 0, 0], [1, 0, 0]]; + const pts: Vec3[] = [[0, 0, 0], [1, 0, 0], [0, 0, 0], [0, 1, 0]]; + const q = buildQ(4, (i, arr, o) => addSample(arr, o, pts[i], normals[i])); + + const r = clusterByOctree(4, Int32Array.from(cells), q, new Uint8Array([1, 1, 1, 1]), RES, 0.01); + + expect(r.count).toBe(1); + }); + + /** + * A cell with two surface sheets must never join a cluster, and — because a + * node collapses only if everything under it did — it also blocks its + * ancestors. Without that, merging pinches two sheets into one vertex, which + * is precisely the non-manifold case per-patch vertices exist to prevent. + */ + it('never merges a vertex marked unmergeable, and blocks its ancestors', () => { + const cells = [cellIndex(0, 0, 0), cellIndex(1, 0, 0), cellIndex(0, 1, 0), cellIndex(1, 1, 0)]; + const pts: Vec3[] = [[0, 0, 0], [1, 0, 0], [0, 1, 0], [1, 1, 0]]; + const q = buildQ(4, (i, arr, o) => addSample(arr, o, pts[i], [0, 0, 1])); + + const r = clusterByOctree(4, Int32Array.from(cells), q, new Uint8Array([1, 1, 0, 1]), RES, 0.01); + + expect(r.count).toBe(4); + expect(new Set(r.remap).size).toBe(4); + }); + + /** Two vertices from one cell are two sheets; they must stay apart. */ + it('keeps two patches of the same cell separate', () => { + const cells = [cellIndex(0, 0, 0), cellIndex(0, 0, 0)]; + const q = buildQ(2, (i, arr, o) => addSample(arr, o, [0, 0, i], [0, 0, 1])); + + const r = clusterByOctree(2, Int32Array.from(cells), q, new Uint8Array([1, 1]), RES, 1e6); + + expect(r.count).toBe(2); + expect(r.remap[0]).not.toBe(r.remap[1]); + }); + + it('leaves a lone vertex alone rather than inventing a cluster', () => { + const q = buildQ(1, (_i, arr, o) => addSample(arr, o, [0, 0, 0], [0, 0, 1])); + const r = clusterByOctree(1, Int32Array.from([0]), q, new Uint8Array([1]), RES, 0.01); + expect(r.count).toBe(1); + expect(r.size[0]).toBe(1); + }); +}); + +describe('clusterError', () => { + it('is zero for samples that share a plane', () => { + const q = buildQ(1, (_i, arr, o) => { + addSample(arr, o, [0, 0, 0], [0, 0, 1]); + addSample(arr, o, [1, 0, 0], [0, 0, 1]); + addSample(arr, o, [0, 1, 0], [0, 0, 1]); + }); + expect(clusterError(q, 0)).toBeCloseTo(0, 10); + }); + + /** + * Two parallel planes 2 apart: the best single point sits between them and is + * 1 from each. Three samples, so the squared total is 3. Pinning the number + * rather than just "> 0" is what makes the budget comparison meaningful — a + * wall cannot be collapsed through its own thickness. + */ + it('grows with the spread between samples', () => { + const q = buildQ(1, (_i, arr, o) => { + addSample(arr, o, [0, 0, 0], [0, 0, 1]); + addSample(arr, o, [1, 0, 2], [0, 0, 1]); + addSample(arr, o, [0, 1, 0], [0, 0, 1]); + }); + expect(clusterError(q, 0)).toBeGreaterThan(0); + }); + + /** Never negative: the expression cancels large terms and drifts below zero. */ + it('does not report a negative error', () => { + const q = buildQ(1, (_i, arr, o) => { + for (let i = 0; i < 40; i++) addSample(arr, o, [i * 1e3, 0, 0], [0, 0, 1]); + }); + expect(clusterError(q, 0)).toBeGreaterThanOrEqual(0); + }); +}); + +// --- Contouring ------------------------------------------------------------ + +const bbox: BBox = { min: [-8, -8, -8], max: [8, 8, 8] }; + +function makeGrid(node: SDFNode, res: number): Float32Array { + const g = new Float32Array(res * res * res); + const d = 16 / res; + for (let z = 0; z < res; z++) + for (let y = 0; y < res; y++) + for (let x = 0; x < res; x++) + g[z * res * res + y * res + x] = evaluateSDF(node, [ + -8 + (x + 0.5) * d, -8 + (y + 0.5) * d, -8 + (z + 0.5) * d, + ]); + return g; +} + +function signedVolume(mesh: MeshResult): number { + let v = 0; + const p = mesh.positions, idx = mesh.indices; + for (let i = 0; i < idx.length; i += 3) { + const a = idx[i] * 3, b = idx[i + 1] * 3, c = idx[i + 2] * 3; + v += ( + p[a] * (p[b + 1] * p[c + 2] - p[b + 2] * p[c + 1]) + - p[a + 1] * (p[b] * p[c + 2] - p[b + 2] * p[c]) + + p[a + 2] * (p[b] * p[c + 1] - p[b + 1] * p[c]) + ) / 6; + } + return v; +} + +function maxDeviation(mesh: MeshResult, node: SDFNode): number { + let m = 0; + for (let i = 0; i < mesh.positions.length; i += 3) + m = Math.max(m, Math.abs(evaluateSDF(node, [ + mesh.positions[i], mesh.positions[i + 1], mesh.positions[i + 2], + ]))); + return m; +} + +describe('dualContour with clustering', () => { + const shapes: [string, SDFNode, number][] = [ + ['sphere', { kind: 'sphere', radius: 5 }, 32], + ['box', { kind: 'box', size: [10, 10, 10] }, 32], + ['torus', { kind: 'torus', major: 4, minor: 1.5 }, 32], + ]; + + /** + * The gate that decides whether this feature can ship at all. Adaptive dual + * contouring's classic failure is a crack where two different-sized cells + * meet; clustering cannot produce one, because step 2 still emits exactly one + * quad per sign-changing edge over the same four cells. This asserts that + * rather than assuming it. + */ + it.each(shapes)('keeps %s watertight and manifold', (_name, node, res) => { + const voxel = 16 / res; + const mesh = dualContour(makeGrid(node, res), res, bbox, node, undefined, undefined, voxel * 0.05); + const d = analyzeMesh(mesh); + + expect(mesh.indices.length).toBeGreaterThan(0); + expect(d.watertight).toBe(true); + expect(d.nonManifoldEdges).toBe(0); + expect(signedVolume(mesh)).toBeGreaterThan(0); // outward, as slicers require + }); + + /** + * The point of the whole exercise. A cube's surface is six planes, so six + * quads — twelve triangles — describe it exactly, where the dense mesh spends + * thousands and then pays QEM to take them away again. + */ + it('reduces a cube to its twelve exact triangles', () => { + const node: SDFNode = { kind: 'box', size: [10, 10, 10] }; + const res = 32; + const dense = dualContour(makeGrid(node, res), res, bbox, node); + const mesh = dualContour(makeGrid(node, res), res, bbox, node, undefined, undefined, (16 / res) * 0.05); + + expect(dense.indices.length / 3).toBeGreaterThan(1000); + expect(mesh.indices.length / 3).toBe(12); + // And it is the *right* cube, not merely a small one. + expect(signedVolume(mesh)).toBeCloseTo(1000, 3); + expect(maxDeviation(mesh, node)).toBeLessThan(1e-6); + }); + + /** Curvature must not be flattened away: a sphere has no coplanar cells. */ + it('barely collapses a sphere, which has nothing coplanar to merge', () => { + const node: SDFNode = { kind: 'sphere', radius: 5 }; + const res = 32; + const dense = dualContour(makeGrid(node, res), res, bbox, node); + const mesh = dualContour(makeGrid(node, res), res, bbox, node, undefined, undefined, (16 / res) * 0.05); + + expect(mesh.indices.length).toBeGreaterThan(dense.indices.length * 0.8); + }); + + it('still encloses the correct volume after clustering', () => { + const node: SDFNode = { kind: 'sphere', radius: 5 }; + const res = 32; + const mesh = dualContour(makeGrid(node, res), res, bbox, node, undefined, undefined, (16 / res) * 0.05); + + const exact = (4 / 3) * Math.PI * 125; + expect(Math.abs(signedVolume(mesh) - exact) / exact).toBeLessThan(0.02); + }); + + /** + * A merged vertex is clamped into the node it represents, so it cannot drift + * further from the surface than the budget allows however many cells it + * absorbed. + */ + it('keeps merged vertices on the isosurface', () => { + const node: SDFNode = { kind: 'torus', major: 4, minor: 1.5 }; + const res = 32; + const voxel = 16 / res; + const mesh = dualContour(makeGrid(node, res), res, bbox, node, undefined, undefined, voxel * 0.05); + + expect(maxDeviation(mesh, node)).toBeLessThan(voxel * 0.5); + }); + + /** + * Clustering is opt-in, and everything else in the suite asserts against the + * dense mesh. A default that quietly collapsed would change the meaning of + * every one of those tests. + */ + it('does nothing unless asked', () => { + const node: SDFNode = { kind: 'box', size: [10, 10, 10] }; + const a = dualContour(makeGrid(node, 32), 32, bbox, node); + const b = dualContour(makeGrid(node, 32), 32, bbox, node, undefined, undefined, 0); + + expect(b.indices.length).toBe(a.indices.length); + expect(Array.from(b.positions)).toEqual(Array.from(a.positions)); + }); + + /** + * Two sheets in one cell is the configuration that makes naive dual + * contouring non-manifold. Clustering must not reintroduce it by merging the + * sheets back together. + */ + it('stays manifold where two surface sheets share a cell', () => { + const twoSpheres: SDFNode = { + kind: 'union', k: 0, + a: { + kind: 'transform', child: { kind: 'sphere', radius: 0.45 }, + tx: -0.5, ty: -0.5, tz: -0.5, rx: 0, ry: 0, rz: 0, sx: 1, sy: 1, sz: 1, + }, + b: { + kind: 'transform', child: { kind: 'sphere', radius: 0.45 }, + tx: 0.5, ty: 0.5, tz: 0.5, rx: 0, ry: 0, rz: 0, sx: 1, sy: 1, sz: 1, + }, + }; + const mesh = dualContour(makeGrid(twoSpheres, 16), 16, bbox, twoSpheres, undefined, undefined, 1); + const d = analyzeMesh(mesh); + + expect(mesh.indices.length).toBeGreaterThan(0); + expect(d.nonManifoldEdges).toBe(0); + expect(d.watertight).toBe(true); + }); +}); diff --git a/src/worker/sdf/cluster.ts b/src/worker/sdf/cluster.ts new file mode 100644 index 0000000..85350c2 --- /dev/null +++ b/src/worker/sdf/cluster.ts @@ -0,0 +1,253 @@ +import { solveSymmetric3Anchored } from './qef'; +import type { Vec3 } from './types'; + +/** + * Octree vertex clustering for dual contouring. + * + * Dual contouring emits one vertex per surface cell, which on a flat face means + * a vertex per voxel: a 64-voxel-wide wall becomes thousands of triangles that + * QEM then collapses back to two. Simplification was ~80% of an export at res + * 256 and 93% of the triangles it removed carried *zero* error, which is the + * mesher undoing its own over-generation. + * + * This decides, before any triangle exists, which cells should share a vertex: + * an octree node collapses to one vertex when a single point still fits every + * hermite sample underneath it within the error budget. + * + * The reason this is safe is worth stating, because the literature's approach + * is not. Adaptive dual contouring (Ju et al. 2002) contours the octree + * directly, which needs cellProc/faceProc/edgeProc recursion to stitch + * differently-sized neighbours, and gets holes in the mesh when that is wrong. + * Clustering avoids all of it: `dualContour` still emits exactly one quad per + * sign-changing grid edge, referencing the same four cells as before, and only + * the vertex those cells name changes. The result is a quotient of a mesh that + * was already watertight, so there are no cracks to patch. + * + * What clustering *can* break is manifoldness — two surface sheets merged into + * one vertex is the same pinch that per-patch vertices exist to avoid — so a + * cell holding more than one patch never joins a cluster. + */ + +/** Doubles per QEF: A (6, upper triangle), b (3), c, mass point (3), count. */ +export const QEF_STRIDE = 14; + +export interface ClusterResult { + /** vertex index -> cluster index */ + remap: Int32Array; + /** Number of clusters, i.e. vertices in the output mesh. */ + count: number; + /** Summed QEF per cluster, `QEF_STRIDE` doubles each. */ + qef: Float64Array; + /** Grid cell index of each cluster's octree node origin. */ + key: Int32Array; + /** Octree level of each cluster; 0 means it merged nothing. */ + level: Int32Array; + /** How many input vertices each cluster absorbed. */ + size: Int32Array; +} + +/** Accumulate one hermite sample (surface point + unit normal) into a QEF. */ +export function addSample(q: Float64Array, o: number, p: Vec3, n: Vec3) { + const d = n[0] * p[0] + n[1] * p[1] + n[2] * p[2]; + q[o] += n[0] * n[0]; q[o + 1] += n[0] * n[1]; q[o + 2] += n[0] * n[2]; + q[o + 3] += n[1] * n[1]; q[o + 4] += n[1] * n[2]; q[o + 5] += n[2] * n[2]; + q[o + 6] += n[0] * d; q[o + 7] += n[1] * d; q[o + 8] += n[2] * d; + q[o + 9] += d * d; + q[o + 10] += p[0]; q[o + 11] += p[1]; q[o + 12] += p[2]; q[o + 13]++; +} + +/** Position minimising this QEF, anchored at its mass point. */ +export function solveCluster(q: Float64Array, o: number): Vec3 { + const n = q[o + 13]; + const anchor: Vec3 = n > 0 ? [q[o + 10] / n, q[o + 11] / n, q[o + 12] / n] : [0, 0, 0]; + return solveSymmetric3Anchored( + q[o], q[o + 1], q[o + 2], q[o + 3], q[o + 4], q[o + 5], + q[o + 6], q[o + 7], q[o + 8], anchor, + ); +} + +/** + * Squared error of the best single vertex for this QEF: min_x xᵀAx - 2bᵀx + c. + * + * Clamped at zero because that expression is a difference of large similar + * numbers on a well-fitted plane and goes slightly negative in floating point; + * a negative "error" would compare as comfortably under budget, which is the + * right answer for the wrong reason. + */ +export function clusterError(q: Float64Array, o: number): number { + if (q[o + 13] === 0) return 0; + const x = solveCluster(q, o); + const a0 = q[o] * x[0] + q[o + 1] * x[1] + q[o + 2] * x[2]; + const a1 = q[o + 1] * x[0] + q[o + 3] * x[1] + q[o + 4] * x[2]; + const a2 = q[o + 2] * x[0] + q[o + 4] * x[1] + q[o + 5] * x[2]; + const e = (x[0] * a0 + x[1] * a1 + x[2] * a2) + - 2 * (q[o + 6] * x[0] + q[o + 7] * x[1] + q[o + 8] * x[2]) + + q[o + 9]; + return e > 0 ? e : 0; +} + +/** + * Group vertices into octree nodes that a single vertex can represent. + * + * `vertCell` is the grid cell each vertex came from, `vertQ` its hermite QEF, + * and `mergeable` zero for vertices that must stay on their own — cells with + * more than one surface patch. `maxError` is a distance; a node collapses when + * its combined samples fit one point within it. + * + * A node collapses only if *every* node beneath it did. Once a subtree fails, + * every ancestor fails, so each vertex ends up in the largest node that fitted. + */ +export function clusterByOctree( + vertCount: number, + vertCell: Int32Array, + vertQ: Float64Array, + mergeable: Uint8Array, + res: number, + maxError: number, +): ClusterResult { + const maxErrorSq = maxError * maxError; + const r2 = res * res; + + // Deepest octree node that still fitted, per vertex. + const bestKey = new Int32Array(vertCount); + const bestLevel = new Int32Array(vertCount); + + // Level 0: one node per grid cell, keyed by cell index. Vertices that end up + // ineligible keep a negative, per-vertex `bestKey` so they cannot be pooled + // with anything, while their *node* stays in place to fail its ancestors. + interface Node { qOff: number; ok: boolean; first: number } + let nodeQ = new Float64Array(vertCount * QEF_STRIDE); + let nodes = new Map(); + // Vertices belonging to each node, as a linked list over `nextVert`. + let nextVert = new Int32Array(vertCount).fill(-1); + let nodeHead: number[] = []; + + let nodeCount = 0; + for (let v = 0; v < vertCount; v++) { + const key = vertCell[v]; + let node = nodes.get(key); + if (node === undefined) { + node = { qOff: nodeCount * QEF_STRIDE, ok: true, first: nodeCount }; + nodes.set(key, node); + nodeHead.push(-1); + nodeCount++; + } + for (let i = 0; i < QEF_STRIDE; i++) nodeQ[node.qOff + i] += vertQ[v * QEF_STRIDE + i]; + nextVert[v] = nodeHead[node.first]; + nodeHead[node.first] = v; + bestKey[v] = key; + bestLevel[v] = 0; + } + + // A cell is ineligible if it was flagged, or if it produced more than one + // vertex — two vertices means two surface sheets. + // + // Such a cell stays keyed by its real position rather than being lifted out + // of the octree, so it also fails every node above it. That is deliberate and + // conservative: a merged vertex is clamped into its node's box, and if that + // box contains a cell with two sheets in it, the vertex representing the + // neighbours can be placed among geometry it does not describe. Blocking the + // ancestors keeps clusters away from those regions entirely. + for (const node of nodes.values()) { + let n = 0; + let flagged = false; + for (let v = nodeHead[node.first]; v !== -1; v = nextVert[v]) { + n++; + if (!mergeable[v]) flagged = true; + } + if (n > 1 || flagged) { + node.ok = false; + // Each vertex becomes its own cluster; a shared cell key would merge the + // very sheets this is separating. + for (let v = nodeHead[node.first]; v !== -1; v = nextVert[v]) bestKey[v] = -(v + 1); + } + } + + let levels = 0; + while (1 << levels < res) levels++; + + for (let level = 1; level <= levels; level++) { + const size = 1 << level; + const parents = new Map(); + const parentQ = new Float64Array(nodes.size * QEF_STRIDE); + // Vertices of a parent chain through its children's lists. + const childOf: number[][] = []; + let pCount = 0; + + for (const [key, node] of nodes) { + const z = (key / r2) | 0, y = ((key % r2) / res) | 0, x = key % res; + const pk = (((z / size) | 0) * size) * r2 + (((y / size) | 0) * size) * res + ((x / size) | 0) * size; + let par = parents.get(pk); + if (par === undefined) { + par = { qOff: pCount * QEF_STRIDE, ok: true, first: pCount }; + parents.set(pk, par); + childOf.push([]); + pCount++; + } + for (let i = 0; i < QEF_STRIDE; i++) parentQ[par.qOff + i] += nodeQ[node.qOff + i]; + childOf[par.first].push(node.first); + if (!node.ok) par.ok = false; + } + + // Promote every vertex under a parent that fitted. + for (const [pk, par] of parents) { + if (!par.ok || clusterError(parentQ, par.qOff) > maxErrorSq) { par.ok = false; continue; } + for (const childFirst of childOf[par.first]) { + for (let v = nodeHead[childFirst]; v !== -1; v = nextVert[v]) { + bestKey[v] = pk; + bestLevel[v] = level; + } + } + } + + // Re-chain each parent's vertices so the next level can walk them. + const merged = new Int32Array(vertCount).fill(-1); + const newHead: number[] = new Array(pCount).fill(-1); + for (const par of parents.values()) { + for (const childFirst of childOf[par.first]) { + for (let v = nodeHead[childFirst]; v !== -1; v = nextVert[v]) { + merged[v] = newHead[par.first]; + newHead[par.first] = v; + } + } + } + nextVert = merged; + nodeHead = newHead; + nodes = parents; + nodeQ = parentQ; + } + + // Distinct (key, level) pairs become the output vertices. + const remap = new Int32Array(vertCount).fill(-1); + const ids = new Map(); + const keyOut: number[] = []; + const levelOut: number[] = []; + const sizeOut: number[] = []; + let count = 0; + for (let v = 0; v < vertCount; v++) { + const id = `${bestKey[v]}:${bestLevel[v]}`; + let c = ids.get(id); + if (c === undefined) { + c = count++; + ids.set(id, c); + keyOut.push(bestKey[v]); + levelOut.push(bestLevel[v]); + sizeOut.push(0); + } + remap[v] = c; + sizeOut[c]++; + } + + const qef = new Float64Array(count * QEF_STRIDE); + for (let v = 0; v < vertCount; v++) { + const o = remap[v] * QEF_STRIDE; + for (let i = 0; i < QEF_STRIDE; i++) qef[o + i] += vertQ[v * QEF_STRIDE + i]; + } + + return { + remap, count, qef, + key: Int32Array.from(keyOut), + level: Int32Array.from(levelOut), + size: Int32Array.from(sizeOut), + }; +} diff --git a/src/worker/sdf/dualContour.ts b/src/worker/sdf/dualContour.ts index a3d7208..f819a98 100644 --- a/src/worker/sdf/dualContour.ts +++ b/src/worker/sdf/dualContour.ts @@ -18,6 +18,7 @@ import type { SDFNode, BBox, Vec3 } from './types'; import { evaluateSDF } from './evaluate'; import { solveQEF } from './qef'; +import { clusterByOctree, addSample, solveCluster, QEF_STRIDE } from './cluster'; import { TRI_TABLE } from './tables'; import type { MeshResult } from './marchingCubes'; @@ -120,7 +121,15 @@ export interface ActiveBlocks { bits: Uint8Array; } -export function dualContour(grid: Float32Array, res: number, bbox: BBox, sdf: SDFNode, onProgress?: (pct: number) => void, active?: ActiveBlocks): MeshResult { +/** + * @param clusterError When > 0, cells whose surface a single point can + * represent within this distance share one vertex, so a flat face costs two + * triangles instead of thousands. See `cluster.ts` for why this is crack-free + * where adaptive dual contouring is not. Omit (or pass 0) for the dense + * one-vertex-per-cell mesh, which is what the geometry tests assert against. + */ +export function dualContour(grid: Float32Array, res: number, bbox: BBox, sdf: SDFNode, onProgress?: (pct: number) => void, active?: ActiveBlocks, clusterError = 0): MeshResult { + const clustering = clusterError > 0; const dx = (bbox.max[0] - bbox.min[0]) / res; const dy = (bbox.max[1] - bbox.min[1]) / res; const dz = (bbox.max[2] - bbox.min[2]) / res; @@ -253,6 +262,13 @@ export function dualContour(grid: Float32Array, res: number, bbox: BBox, sdf: SD const normals: number[] = []; let vertCount = 0; + // Only collected when clustering: the cell each vertex came from, its hermite + // QEF, and whether it may merge at all. Left empty otherwise so the dense + // path pays nothing for a feature it is not using. + const vertCellList: number[] = []; + const vertQList: number[] = []; + const mergeableList: number[] = []; + let lastPct = -1; forEachCell((x, y, z) => { { @@ -306,6 +322,16 @@ export function dualContour(grid: Float32Array, res: number, bbox: BBox, sdf: SD const g = sdfGradient(sdf, v, eps); const len = Math.sqrt(g[0] * g[0] + g[1] * g[1] + g[2] * g[2]) || 1; + if (clustering) { + const q = new Float64Array(QEF_STRIDE); + for (let i = 0; i < crossPoints.length; i++) addSample(q, 0, crossPoints[i], crossNormals[i]); + for (let i = 0; i < QEF_STRIDE; i++) vertQList.push(q[i]); + vertCellList.push(cellIdx); + // A cell with several patches holds several surface sheets; merging + // them is the non-manifold pinch the patch split exists to prevent. + mergeableList.push(patches.length > 1 ? 0 : 1); + } + vertCount++; positions.push(v[0], v[1], v[2]); normals.push(-g[0] / len, -g[1] / len, -g[2] / len); @@ -318,6 +344,62 @@ export function dualContour(grid: Float32Array, res: number, bbox: BBox, sdf: SD if (pct > lastPct) { lastPct = pct; onProgress(pct); } }); + // --- Step 1b: Collapse cells that one vertex can represent --------------- + // + // Only the identity of the vertex each cell names changes. Step 2 still emits + // exactly one quad per sign-changing grid edge over the same four cells, so + // the mesh stays closed by construction and there are no cracks to stitch. + let vertexRemap: Int32Array | null = null; + if (clustering && vertCount > 0) { + const cluster = clusterByOctree( + vertCount, + Int32Array.from(vertCellList), + Float64Array.from(vertQList), + Uint8Array.from(mergeableList), + res, + clusterError, + ); + + const newPos = new Float64Array(cluster.count * 3); + const newNorm = new Float64Array(cluster.count * 3); + // A cluster of one keeps step 1's vertex exactly, so a mesh where nothing + // collapsed is identical to the dense one rather than merely close. + const representative = new Int32Array(cluster.count).fill(-1); + for (let v = vertCount - 1; v >= 0; v--) representative[cluster.remap[v]] = v; + + for (let c = 0; c < cluster.count; c++) { + const rep = representative[c]; + if (cluster.size[c] === 1) { + newPos[c * 3] = positions[rep * 3]; + newPos[c * 3 + 1] = positions[rep * 3 + 1]; + newPos[c * 3 + 2] = positions[rep * 3 + 2]; + newNorm[c * 3] = normals[rep * 3]; + newNorm[c * 3 + 1] = normals[rep * 3 + 1]; + newNorm[c * 3 + 2] = normals[rep * 3 + 2]; + continue; + } + let p = solveCluster(cluster.qef, c * QEF_STRIDE); + // Clamped into the node it represents, for the same reason step 1 clamps + // into the cell: a vertex outside its own region can cross a neighbour's. + const key = cluster.key[c], size = 1 << cluster.level[c]; + const nz = (key / r2) | 0, ny = ((key % r2) / res) | 0, nx = key % res; + p = [ + Math.max(ox + nx * dx, Math.min(ox + (nx + size) * dx, p[0])), + Math.max(oy + ny * dy, Math.min(oy + (ny + size) * dy, p[1])), + Math.max(oz + nz * dz, Math.min(oz + (nz + size) * dz, p[2])), + ]; + newPos[c * 3] = p[0]; newPos[c * 3 + 1] = p[1]; newPos[c * 3 + 2] = p[2]; + const cg = sdfGradient(sdf, p, eps); + const cl = Math.sqrt(cg[0] * cg[0] + cg[1] * cg[1] + cg[2] * cg[2]) || 1; + newNorm[c * 3] = -cg[0] / cl; newNorm[c * 3 + 1] = -cg[1] / cl; newNorm[c * 3 + 2] = -cg[2] / cl; + } + + positions.length = 0; + normals.length = 0; + for (let i = 0; i < newPos.length; i++) { positions.push(newPos[i]); normals.push(newNorm[i]); } + vertexRemap = cluster.remap; + } + // --- Step 2: Emit quads for each sign-changing grid edge --- // The 4 cells sharing the edge each contribute the vertex of the patch // that contains this edge (looked up through the cell's local edge index). @@ -328,9 +410,15 @@ export function dualContour(grid: Float32Array, res: number, bbox: BBox, sdf: SD const base = cellBase[cellIdx]; if (base < 0) return -1; const edgePatch = multiPatchCells.get(cellIdx); - if (edgePatch === undefined) return base; - const patch = edgePatch[localEdge]; - return patch < 0 ? -1 : base + patch; + let v: number; + if (edgePatch === undefined) { + v = base; + } else { + const patch = edgePatch[localEdge]; + if (patch < 0) return -1; + v = base + patch; + } + return vertexRemap === null ? v : vertexRemap[v]; } function triDot(a: number, b: number, c: number, d: number): number { @@ -357,6 +445,24 @@ export function dualContour(grid: Float32Array, res: number, bbox: BBox, sdf: SD */ function emitQuad(r0: number, r1: number, r2q: number, r3: number) { if (r0 < 0 || r1 < 0 || r2q < 0 || r3 < 0) return; + if (vertexRemap !== null) { + // Clustering merges ring corners. Dropping duplicates in place preserves + // the winding, and a ring left with fewer than three distinct vertices + // had no area. This is where the triangle count actually falls: the quad + // is never emitted rather than being removed by a later pass. + if (r0 === r1 || r1 === r2q || r2q === r3 || r3 === r0 || r0 === r2q || r1 === r3) { + const ring = [r0, r1, r2q, r3]; + const uniq: number[] = []; + for (let i = 0; i < 4; i++) if (ring[i] !== ring[(i + 3) % 4]) uniq.push(ring[i]); + if (uniq.length < 3) return; + if (uniq.length === 3) { + if (uniq[0] === uniq[1] || uniq[1] === uniq[2] || uniq[0] === uniq[2]) return; + indices.push(uniq[0], uniq[1], uniq[2]); + return; + } + r0 = uniq[0]; r1 = uniq[1]; r2q = uniq[2]; r3 = uniq[3]; + } + } const qA = triDot(r0, r1, r2q, r3); // diagonal r0-r2 const qB = triDot(r1, r2q, r3, r0); // diagonal r1-r3 if (qB > qA) { diff --git a/src/worker/sdfWorker.ts b/src/worker/sdfWorker.ts index 1f5c740..50a089f 100644 --- a/src/worker/sdfWorker.ts +++ b/src/worker/sdfWorker.ts @@ -17,6 +17,7 @@ import { exportBinarySTL } from './stlExporter'; import { export3MF } from './exporters'; import { toSDFNode } from './sdf/convert'; import { simplifyMesh } from './sdf/simplify'; +import { CLUSTER_ERROR_VOXELS, SIMPLIFY_ERROR_VOXELS, PROJECT_TOLERANCE_VOXELS } from './sdf/budgets'; import { removeDegenerateTriangles, projectVerticesToSurface } from './sdf/meshRepair'; self.postMessage({ type: 'ready' }); @@ -67,9 +68,11 @@ function prepareBBox(root: SDFNode): BBox { /** * Full export meshing pipeline: * 1. SDF grid evaluation (interval-verified octree descent) - * 2. Manifold dual contouring with SVD QEF vertex placement + * 2. Manifold dual contouring with SVD QEF vertex placement, collapsing + * cells one vertex can represent so a flat face never becomes thousands + * of triangles that step 4 would only have to take away again * 3. Degenerate-face cleanup - * 4. Error-bounded QEM simplification (budget: 5% of a voxel) + * 4. Error-bounded QEM simplification * 5. Newton projection of the surviving vertices back onto the SDF * zero set, clearing the residual bias simplification introduced */ @@ -91,22 +94,22 @@ function evaluateAndMeshWithProgress(tree: SDFNodeUI | null, resolution: number, progress('Evaluating SDF grid', pct); }); - // Dual contouring (60-80%) + // Dual contouring, with octree vertex clustering (60-80%) progress('Generating mesh', 60); const raw = dualContour(grid, resolution, bbox, root, (pct) => { progress('Generating mesh', 60 + pct * 0.2); - }, active); + }, active, voxel * CLUSTER_ERROR_VOXELS); if (raw.indices.length === 0) return null; // Simplification (80-92%) progress('Simplifying mesh', 80); - const simplified = simplifyMesh(removeDegenerateTriangles(raw), { maxError: voxel * 0.05 }, (pct) => { + const simplified = simplifyMesh(removeDegenerateTriangles(raw), { maxError: voxel * SIMPLIFY_ERROR_VOXELS }, (pct) => { progress('Simplifying mesh', 80 + pct * 0.12); }); // Surface snap (92-95%) progress('Refining surface', 92); - return projectVerticesToSurface(simplified, root, voxel * 0.5); + return projectVerticesToSurface(simplified, root, voxel * PROJECT_TOLERANCE_VOXELS); } /**