Fix chamfered crisp edges: anchor the QEM quadric solve at the edge midpoint

src/worker/sdf/qef.ts

@@ -101,26 +101,22 @@ function eigenSymmetric3(
 }
 
 /**
- * Solve the QEF for a set of Hermite samples (surface points + normals),
- * anchored at the mass point.
- *
- * Returns the position minimizing sum_i (n_i . (x - p_i))^2, where
- * rank-deficient directions of A^T A stay at the mass point.
+ * Solve the 3x3 symmetric system A x = b by truncated pseudo-inverse,
+ * anchored at `anchor`: directions of A whose eigenvalue falls below
+ * SINGULAR_TOLERANCE x the largest keep the anchor's coordinate instead
+ * of being dragged by noise.  This is the rank-safe replacement for a
+ * determinant solve wherever A can be legitimately singular — flat
+ * regions (rank 1) and straight creases (rank 2), where an unregularized
+ * solution slides arbitrarily far along the null direction.
  */
-export function solveQEF(points: Vec3[], normals: Vec3[], massPoint: Vec3): Vec3 {
-  let a00 = 0, a01 = 0, a02 = 0, a11 = 0, a12 = 0, a22 = 0;
-  let b0 = 0, b1 = 0, b2 = 0;
-
-  for (let i = 0; i < points.length; i++) {
-    const [nx, ny, nz] = normals[i];
-    const d = nx * points[i][0] + ny * points[i][1] + nz * points[i][2];
-    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;
-  }
-
-  // Shift to mass-point-relative coordinates: solve ATA * dx = ATb - ATA * mp
-  const [mx, my, mz] = massPoint;
+export function solveSymmetric3Anchored(
+  a00: number, a01: number, a02: number,
+  a11: number, a12: number, a22: number,
+  b0: number, b1: number, b2: number,
+  anchor: Vec3,
+): Vec3 {
+  const [mx, my, mz] = anchor;
+  // Shift to anchor-relative coordinates: solve A * dx = b - A * anchor
   const r0 = b0 - (a00 * mx + a01 * my + a02 * mz);
   const r1 = b1 - (a01 * mx + a11 * my + a12 * mz);
   const r2 = b2 - (a02 * mx + a12 * my + a22 * mz);
@@ -141,3 +137,25 @@ export function solveQEF(points: Vec3[], normals: Vec3[], massPoint: Vec3): Vec3
 
   return [mx + dx, my + dy, mz + dz];
 }
+
+/**
+ * Solve the QEF for a set of Hermite samples (surface points + normals),
+ * anchored at the mass point.
+ *
+ * Returns the position minimizing sum_i (n_i . (x - p_i))^2, where
+ * rank-deficient directions of A^T A stay at the mass point.
+ */
+export function solveQEF(points: Vec3[], normals: Vec3[], massPoint: Vec3): Vec3 {
+  let a00 = 0, a01 = 0, a02 = 0, a11 = 0, a12 = 0, a22 = 0;
+  let b0 = 0, b1 = 0, b2 = 0;
+
+  for (let i = 0; i < points.length; i++) {
+    const [nx, ny, nz] = normals[i];
+    const d = nx * points[i][0] + ny * points[i][1] + nz * points[i][2];
+    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;
+  }
+
+  return solveSymmetric3Anchored(a00, a01, a02, a11, a12, a22, b0, b1, b2, massPoint);
+}

src/worker/sdf/simplify.test.ts

@@ -1,6 +1,7 @@
 import { describe, it, expect } from 'vitest';
 import { simplifyMesh, splitCreaseEdges } from './simplify';
 import { marchingCubes } from './marchingCubes';
+import { dualContour } from './dualContour';
 import { evaluateSDF } from './evaluate';
 import type { SDFNode, BBox } from './types';
 import type { MeshResult } from './marchingCubes';
@@ -26,6 +27,51 @@ function makeGrid(node: SDFNode, resolution: number, bbox: BBox): Float32Array {
   return grid;
 }
 
+/** Exact point-to-triangle distance (Ericson, Real-Time Collision Detection) */
+function pointTriangleDistance(
+  p: [number, number, number], positions: Float32Array,
+  ia: number, ib: number, ic: number,
+): number {
+  const a = [positions[ia * 3], positions[ia * 3 + 1], positions[ia * 3 + 2]];
+  const b = [positions[ib * 3], positions[ib * 3 + 1], positions[ib * 3 + 2]];
+  const c = [positions[ic * 3], positions[ic * 3 + 1], positions[ic * 3 + 2]];
+  const sub = (u: number[], v: number[]) => [u[0] - v[0], u[1] - v[1], u[2] - v[2]];
+  const dot = (u: number[], v: number[]) => u[0] * v[0] + u[1] * v[1] + u[2] * v[2];
+  const len = (u: number[]) => Math.hypot(u[0], u[1], u[2]);
+
+  const ab = sub(b, a), ac = sub(c, a), ap = sub(p, a);
+  const d1 = dot(ab, ap), d2 = dot(ac, ap);
+  if (d1 <= 0 && d2 <= 0) return len(ap);
+  const bp = sub(p, b);
+  const d3 = dot(ab, bp), d4 = dot(ac, bp);
+  if (d3 >= 0 && d4 <= d3) return len(bp);
+  const vc = d1 * d4 - d3 * d2;
+  if (vc <= 0 && d1 >= 0 && d3 <= 0) {
+    const v = d1 / (d1 - d3);
+    return len(sub(p, [a[0] + v * ab[0], a[1] + v * ab[1], a[2] + v * ab[2]]));
+  }
+  const cp = sub(p, c);
+  const d5 = dot(ab, cp), d6 = dot(ac, cp);
+  if (d6 >= 0 && d5 <= d6) return len(cp);
+  const vb = d5 * d2 - d1 * d6;
+  if (vb <= 0 && d2 >= 0 && d6 <= 0) {
+    const w = d2 / (d2 - d6);
+    return len(sub(p, [a[0] + w * ac[0], a[1] + w * ac[1], a[2] + w * ac[2]]));
+  }
+  const va = d3 * d6 - d5 * d4;
+  if (va <= 0 && d4 - d3 >= 0 && d5 - d6 >= 0) {
+    const w = (d4 - d3) / (d4 - d3 + (d5 - d6));
+    return len(sub(p, [b[0] + w * (c[0] - b[0]), b[1] + w * (c[1] - b[1]), b[2] + w * (c[2] - b[2])]));
+  }
+  const denom = 1 / (va + vb + vc);
+  const v = vb * denom, w = vc * denom;
+  return len(sub(p, [
+    a[0] + ab[0] * v + ac[0] * w,
+    a[1] + ab[1] * v + ac[1] * w,
+    a[2] + ab[2] * v + ac[2] * w,
+  ]));
+}
+
 function meshFromSDF(node: SDFNode, res: number): ReturnType<typeof marchingCubes> {
   const bbox: BBox = { min: [-8, -8, -8], max: [8, 8, 8] };
   const grid = makeGrid(node, res, bbox);
@@ -171,6 +217,59 @@ describe('simplifyMesh error-bounded mode', () => {
     expect(boxKept).toBeLessThan(sphereKept * 0.5);
   });
 
+  it('keeps the mesh covering sharp creases (no chamfered edges)', () => {
+    // Regression: the quadric cost along a straight crease is zero in the
+    // direction of the crease line, so an unregularized optimal-point solve
+    // let collapse targets slide arbitrarily far along the edge.  Vertices
+    // ended up out of order along the crease and triangles chamfered across
+    // the corner — vertices all exactly on the surface, mesh watertight,
+    // but the crisp edge visibly cut off.  Measure actual coverage: every
+    // point of the true crease line must be close to the simplified mesh.
+    const bbox: BBox = { min: [-8, -8, -8], max: [8, 8, 8] };
+    const res = 48;
+    const voxel = 16 / res;
+    const box: SDFNode = { kind: 'box', size: [10, 10, 10] };
+    const grid = makeGrid(box, res, bbox);
+    const mesh = dualContour(grid, res, bbox, box);
+    const simplified = simplifyMesh(mesh, { maxError: voxel * 0.05 });
+
+    // All 12 crease edges of the box
+    const creases: [number[], number[]][] = [];
+    for (const axis of [0, 1, 2]) {
+      for (const s1 of [-5, 5]) {
+        for (const s2 of [-5, 5]) {
+          const p0 = [0, 0, 0], p1 = [0, 0, 0];
+          p0[axis] = -5; p1[axis] = 5;
+          p0[(axis + 1) % 3] = p1[(axis + 1) % 3] = s1;
+          p0[(axis + 2) % 3] = p1[(axis + 2) % 3] = s2;
+          creases.push([p0, p1]);
+        }
+      }
+    }
+
+    let worst = 0;
+    for (const [p0, p1] of creases) {
+      for (let s = 0; s <= 40; s++) {
+        const t = s / 40;
+        const p: [number, number, number] = [
+          p0[0] + (p1[0] - p0[0]) * t,
+          p0[1] + (p1[1] - p0[1]) * t,
+          p0[2] + (p1[2] - p0[2]) * t,
+        ];
+        let best = Infinity;
+        const { positions, indices } = simplified;
+        for (let f = 0; f < indices.length; f += 3) {
+          best = Math.min(best, pointTriangleDistance(p,
+            positions, indices[f], indices[f + 1], indices[f + 2]));
+        }
+        worst = Math.max(worst, best);
+      }
+    }
+    // Before the anchored solve this was ~1.6 (half a face chamfered off);
+    // grid-accurate creases sit essentially at 0.
+    expect(worst).toBeLessThan(voxel * 0.5);
+  });
+
   it('respects targetRatio as a floor when both constraints are given', () => {
     const mesh = meshFromSDF({ kind: 'box', size: [10, 10, 10] }, 24);
     const simplified = simplifyMesh(mesh, { maxError: 10, targetRatio: 0.5 });

src/worker/sdf/simplify.ts

@@ -12,6 +12,7 @@
  */
 
 import type { MeshResult } from './marchingCubes';
+import { solveSymmetric3Anchored } from './qef';
 
 /** Symmetric 4x4 matrix stored as 10 floats (upper triangle) */
 type Quadric = [number, number, number, number, number, number, number, number, number, number];
@@ -45,25 +46,26 @@ function evalQ(q: Quadric, x: number, y: number, z: number): number {
        + q[9];
 }
 
-/** Find the point minimizing the quadric error. Returns null if the 3x3 system is singular. */
-function optimalQ(q: Quadric): [number, number, number] | null {
-  // Solve the 3x3 linear system:
-  // [a00 a01 a02] [x]   [-a03]
-  // [a01 a11 a12] [y] = [-a13]
-  // [a02 a12 a22] [z]   [-a23]
-  const a = q[0], b = q[1], c = q[2], d = q[3];
-  const e = q[4], f = q[5], g = q[6];
-  const h = q[7], k = q[8];
-
-  const det = a*(e*h - f*f) - b*(b*h - f*c) + c*(b*f - e*c);
-  if (Math.abs(det) < 1e-12) return null;
-
-  const inv = 1 / det;
-  return [
-    ((e*h - f*f)*(-d) + (c*f - b*h)*(-g) + (b*f - c*e)*(-k)) * inv,
-    ((c*f - b*h)*(-d) + (a*h - c*c)*(-g) + (b*c - a*f)*(-k)) * inv,
-    ((b*f - c*e)*(-d) + (b*c - a*f)*(-g) + (a*e - b*b)*(-k)) * inv,
-  ];
+/**
+ * Find the point minimizing the quadric error, anchored at `anchor`.
+ *
+ * The 3x3 system is legitimately singular exactly at the features that
+ * matter most: a flat region is rank 1 and a straight crease is rank 2,
+ * where any point on the plane / crease line has identical cost.  A
+ * determinant solve there is ill-conditioned and lets the vertex slide
+ * arbitrarily far along the null direction — on a box this shuffled crease
+ * vertices out of order along the edge line and the triangulation chamfered
+ * across the corner.  The truncated eigen solve keeps null-direction
+ * coordinates at the anchor (the edge midpoint), which also preserves the
+ * ordering of vertices along a crease.
+ */
+function optimalQ(q: Quadric, anchor: [number, number, number]): [number, number, number] {
+  const [x, y, z] = solveSymmetric3Anchored(
+    q[0], q[1], q[2], q[4], q[5], q[7],
+    -q[3], -q[6], -q[8],
+    anchor,
+  );
+  return [x, y, z];
 }
 
 /** Binary min-heap keyed by cost */
@@ -204,19 +206,33 @@ export function simplifyMesh(mesh: MeshResult, target: number | SimplifyOptions,
   // Compute cost and optimal position for collapsing edge (ra, rb)
   function computeEdgeCost(ra: number, rb: number): { cost: number; pos: [number, number, number] } {
     const q = addQ(quadrics[ra], quadrics[rb]);
-    const opt = optimalQ(q);
-    if (opt) {
-      return { cost: evalQ(q, opt[0], opt[1], opt[2]), pos: opt };
-    }
     const mid: [number, number, number] = [
       (vx[ra] + vx[rb]) * 0.5, (vy[ra] + vy[rb]) * 0.5, (vz[ra] + vz[rb]) * 0.5,
     ];
-    const ea = evalQ(q, vx[ra], vy[ra], vz[ra]);
-    const eb = evalQ(q, vx[rb], vy[rb], vz[rb]);
-    const em = evalQ(q, mid[0], mid[1], mid[2]);
-    if (ea <= eb && ea <= em) return { cost: ea, pos: [vx[ra], vy[ra], vz[ra]] };
-    if (eb <= ea && eb <= em) return { cost: eb, pos: [vx[rb], vy[rb], vz[rb]] };
-    return { cost: em, pos: mid };
+    const opt = optimalQ(q, mid);
+
+    // Even the anchored solve can move a long way when the strong directions
+    // genuinely demand it; a collapse target far outside the edge's own
+    // neighborhood is never geometrically sensible, so fall back to the
+    // endpoint/midpoint candidates in that case.
+    const ex = vx[rb] - vx[ra], ey = vy[rb] - vy[ra], ez = vz[rb] - vz[ra];
+    const edgeLenSq = ex * ex + ey * ey + ez * ez;
+    const ox = opt[0] - mid[0], oy = opt[1] - mid[1], oz = opt[2] - mid[2];
+    const optDistSq = ox * ox + oy * oy + oz * oz;
+
+    let best: [number, number, number] = opt;
+    let bestCost = optDistSq <= Math.max(edgeLenSq * 4, 1e-12)
+      ? evalQ(q, opt[0], opt[1], opt[2])
+      : Infinity;
+    for (const cand of [
+      mid,
+      [vx[ra], vy[ra], vz[ra]] as [number, number, number],
+      [vx[rb], vy[rb], vz[rb]] as [number, number, number],
+    ]) {
+      const c = evalQ(q, cand[0], cand[1], cand[2]);
+      if (c < bestCost) { bestCost = c; best = cand; }
+    }
+    return { cost: bestCost, pos: best };
   }
 
   // Priority queue — entries carry the original vertex pair and a generation