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); } /**