OverviewFilesCommitsIssuesDocsReleasesPeers

dc-terrain commits / 01fdec3

qef/polish: the Newton step was missing a factor of voxel_size
gradient() returns the raw DIFFERENCE across 2*eps = voxel_size world units and never divides by the spacing, so it is grad_d * voxel_size, not grad_d. Newton to the isosurface wants dx = -d * grad_d / |grad_d|^2 which in terms of g = voxel_size * grad_d is dx = -d * g * voxel_size / |g|^2 and the voxel_size was absent. The step was therefore 1/voxel_size times too large. Measured on an analytic sphere (tests/gate_polish_units.py), mean distance from the true surface after the 3 iterations: voxel shipping fixed 0.25 0.0143 0.0000 0.50 0.0304 0.0000 1.00 0.0000 0.0000 2.00 0.0608 0.0000 4.00 0.4091 0.0000 8.00 1.2759 0.0005 Exact at exactly one voxel size, 1.0. That is the signature of a missing linear factor and of nothing else. The two regimes look like different bugs and are one. At the 0.25 m LOD the step overshoots 4x, so Newton oscillates, saturates the +/-0.5*voxel clamp, and adjacent cells settle half a voxel apart -- a regular lattice of tents, one per cell, which is the "upside-down pyramid" floor. On a flat plane the measured vertical spread between neighbouring cells is 0.1875 at voxel 0.25 and 1.2656 at voxel 4, against 0.0 fixed. At coarse LOD the step is instead far too small and the vertex is left up to 1.3 world units off the surface, which is detail going blocky with distance. Only ORGANIC cells reach this code -- hard (convex-brush) cells skip the polish entirely. That is why a Quake CONVEX brush (edit_type 8/9, exact ray-vs-plane Hermite plus the plane-set snap) was unaffected while a SPHERE edit terraced: dc_compute.cpp:937 maps spheres to type 0/1/6/7, which is the organic path. BOX (2/3) and MESH SDF (4/5) are on that path too and should have shown the same artifact. max_step was also clamping the wrong quantity. It is documented as bounding the DISPLACEMENT but bounded the coefficient, which is then multiplied by grad, so the actual bound was scaled by |grad|. The displacement itself is now clamped. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
01fdec300ecf3c499fc5b02a399f3486a263ad10
by plankatron <fleeky@gmail.com> · 2026-09-08T01:38:37.000Z

Changed files (3) +167 -6

shaders/qef_solve.glsl +29 -2
……@@ -464,8 +464,35 @@ void main() {
464464 float grad_sq = dot(grad, grad);
465465 if (grad_sq < 1e-12) break;
466466
467 float step = clamp(d / grad_sq, -max_step, max_step);
468 local -= step * grad;
467 // NEWTON STEP TO THE ISOSURFACE.
468 //
469 // gradient() returns the raw DIFFERENCE across 2*eps = voxel_size
470 // world units -- it is never divided by the spacing -- so it is
471 // grad_d * voxel_size, NOT grad_d. Newton wants
472 // dx = -d * grad_d / |grad_d|^2
473 // which in terms of g = voxel_size * grad_d is
474 // dx = -d * g * voxel_size / |g|^2.
475 // The voxel_size was missing, so the step was 1/voxel_size times
476 // too large: 4x overshoot at the 0.25 m LOD (Newton oscillates,
477 // saturates the clamp, and adjacent cells land half a voxel apart
478 // -- one tent per cell, the "upside-down pyramid" floor), and badly
479 // UNDER-stepped at coarse LOD (1.3 world units off at voxel 8 --
480 // detail going blocky with distance).
481 //
482 // It was exact at exactly one voxel size, 1.0, which is the
483 // signature of a missing linear factor and nothing else.
484 //
485 // Only ORGANIC cells reach here -- hard (convex-brush) cells skip
486 // the polish -- which is why a Quake brush was unaffected while a
487 // SPHERE edit, which is type 0/1/6/7 and organic, terraced.
488 //
489 // max_step also clamped the wrong quantity: it is documented as
490 // bounding the DISPLACEMENT but bounded the coefficient, which is
491 // then multiplied by grad. Clamp the displacement itself.
492 vec3 dx = -(d * voxel_size / grad_sq) * grad;
493 float dxl = length(dx);
494 if (dxl > max_step) dx *= max_step / dxl;
495 local += dx;
469496 }
470497
471498 // Re-clamp to cell bounds after polish
shaders/vertex_polish.glsl +29 -4
……@@ -121,10 +121,35 @@ void main() {
121121 float grad_sq = dot(grad, grad);
122122 if (grad_sq < 1e-12) break; // flat region, no improvement possible
123123
124 // Newton step: d / dot(grad, grad), clamped for safety
125 float step = clamp(d / grad_sq, -max_step, max_step);
126
127 local -= step * grad;
124 // NEWTON STEP TO THE ISOSURFACE.
125 //
126 // gradient() returns the raw DIFFERENCE across 2*eps = voxel_size
127 // world units -- it is never divided by the spacing -- so it is
128 // grad_d * voxel_size, NOT grad_d. Newton wants
129 // dx = -d * grad_d / |grad_d|^2
130 // which in terms of g = voxel_size * grad_d is
131 // dx = -d * g * voxel_size / |g|^2.
132 // The voxel_size was missing, so the step was 1/voxel_size times
133 // too large: 4x overshoot at the 0.25 m LOD (Newton oscillates,
134 // saturates the clamp, and adjacent cells land half a voxel apart
135 // -- one tent per cell, the "upside-down pyramid" floor), and badly
136 // UNDER-stepped at coarse LOD (1.3 world units off at voxel 8 --
137 // detail going blocky with distance).
138 //
139 // It was exact at exactly one voxel size, 1.0, which is the
140 // signature of a missing linear factor and nothing else.
141 //
142 // Only ORGANIC cells reach here -- hard (convex-brush) cells skip
143 // the polish -- which is why a Quake brush was unaffected while a
144 // SPHERE edit, which is type 0/1/6/7 and organic, terraced.
145 //
146 // max_step also clamped the wrong quantity: it is documented as
147 // bounding the DISPLACEMENT but bounded the coefficient, which is
148 // then multiplied by grad. Clamp the displacement itself.
149 vec3 dx = -(d * voxel_size / grad_sq) * grad;
150 float dxl = length(dx);
151 if (dxl > max_step) dx *= max_step / dxl;
152 local += dx;
128153 }
129154
130155 // Clamp result to voxel cell bounds
tests/gate_polish_units.py +109 -0
……@@ -0,0 +1,109 @@
1"""Verify the CORRECTED polish step, including the clamp fix.
2
3Two units errors in the shipping loop:
4 1. gradient() returns the raw DIFFERENCE across voxel_size world units, so it
5 is grad_d * voxel_size. Newton needs dx = -d*grad_d/|grad_d|^2, which in
6 terms of g is -d*g*voxel_size/|g|^2. The voxel_size is missing.
7 2. max_step is documented as bounding the DISPLACEMENT, but it clamps `step`,
8 which is then multiplied by grad. So the bound is scaled by |grad|.
9"""
10import numpy as np
11
12R, ITERS = 8.0, 3
13
14
15def density(p):
16 return R - np.linalg.norm(p)
17
18
19def gradient(p, vs):
20 eps = vs * 0.5
21 g = np.empty(3)
22 for a in range(3):
23 o = np.zeros(3); o[a] = eps
24 g[a] = density(p + o) - density(p - o)
25 return g
26
27
28def polish(start, vs, mode):
29 local = np.array(start, dtype=float)
30 max_step = 0.5 * vs
31 for _ in range(ITERS):
32 d = density(local)
33 g = gradient(local, vs)
34 gsq = g @ g
35 if gsq < 1e-12:
36 break
37 if mode == "shipping":
38 local -= np.clip(d / gsq, -max_step, max_step) * g
39 elif mode == "scale_only":
40 local -= np.clip(d * vs / gsq, -max_step, max_step) * g
41 else: # "fixed": scale AND clamp the displacement, not the coefficient
42 dx = -(d * vs / gsq) * g
43 dl = np.linalg.norm(dx)
44 if dl > max_step:
45 dx *= max_step / dl
46 local += dx
47 return local
48
49
50rng = np.random.default_rng(3)
51print("%-8s %-14s %-14s %-14s" % ("voxel", "shipping", "scale only", "scale+clamp"))
52for vs in (0.25, 0.5, 1.0, 2.0, 4.0, 8.0):
53 out = {}
54 for mode in ("shipping", "scale_only", "fixed"):
55 errs = []
56 r2 = np.random.default_rng(3)
57 for _ in range(400):
58 dv = r2.normal(size=3); dv /= np.linalg.norm(dv)
59 start = dv * (R + r2.uniform(-0.5, 0.5) * vs)
60 p = polish(start, vs, mode)
61 errs.append(abs(np.linalg.norm(p) - R))
62 out[mode] = np.mean(errs)
63 print("%-8.2f %-14.5f %-14.5f %-14.5f"
64 % (vs, out["shipping"], out["scale_only"], out["fixed"]))
65
66# The tent signature: on a FLAT floor, how far apart do neighbouring cells end
67# up vertically? That difference IS the pyramid height.
68print()
69print("FLAT FLOOR — vertical spread between neighbouring cells (= tent height)")
70
71
72def flat_density(p):
73 return -p[2] + 0.37 # a plane that does NOT sit on a grid line
74
75
76def flat_polish(start, vs, mode):
77 local = np.array(start, dtype=float)
78 max_step = 0.5 * vs
79 for _ in range(ITERS):
80 d = flat_density(local)
81 eps = vs * 0.5
82 g = np.array([0.0, 0.0,
83 flat_density(local + np.array([0, 0, eps]))
84 - flat_density(local - np.array([0, 0, eps]))])
85 gsq = g @ g
86 if gsq < 1e-12:
87 break
88 if mode == "shipping":
89 local -= np.clip(d / gsq, -max_step, max_step) * g
90 else:
91 dx = -(d * vs / gsq) * g
92 dl = np.linalg.norm(dx)
93 if dl > max_step:
94 dx *= max_step / dl
95 local += dx
96 return local
97
98
99print("%-8s %-22s %-22s" % ("voxel", "shipping spread", "fixed spread"))
100for vs in (0.25, 0.5, 1.0, 2.0, 4.0):
101 res = {}
102 for mode in ("shipping", "fixed"):
103 zs = []
104 for k in range(24):
105 z0 = (k % 4) * vs * 0.25 # different start height per cell
106 p = flat_polish(np.array([0.0, 0.0, z0]), vs, mode)
107 zs.append(p[2])
108 res[mode] = max(zs) - min(zs)
109 print("%-8.2f %-22.5f %-22.5f" % (vs, res["shipping"], res["fixed"]))