"""Verify the CORRECTED polish step, including the clamp fix.
Two units errors in the shipping loop:
1. gradient() returns the raw DIFFERENCE across voxel_size world units, so it
is grad_d * voxel_size. Newton needs dx = -d*grad_d/|grad_d|^2, which in
terms of g is -d*g*voxel_size/|g|^2. The voxel_size is missing.
2. max_step is documented as bounding the DISPLACEMENT, but it clamps `step`,
which is then multiplied by grad. So the bound is scaled by |grad|.
"""
import numpy as np
R, ITERS = 8.0, 3
def density(p):
return R - np.linalg.norm(p)
def gradient(p, vs):
eps = vs * 0.5
g = np.empty(3)
for a in range(3):
o = np.zeros(3); o[a] = eps
g[a] = density(p + o) - density(p - o)
return g
def polish(start, vs, mode):
local = np.array(start, dtype=float)
max_step = 0.5 * vs
for _ in range(ITERS):
d = density(local)
g = gradient(local, vs)
gsq = g @ g
if gsq < 1e-12:
break
if mode == "shipping":
local -= np.clip(d / gsq, -max_step, max_step) * g
elif mode == "scale_only":
local -= np.clip(d * vs / gsq, -max_step, max_step) * g
else: # "fixed": scale AND clamp the displacement, not the coefficient
dx = -(d * vs / gsq) * g
dl = np.linalg.norm(dx)
if dl > max_step:
dx *= max_step / dl
local += dx
return local
rng = np.random.default_rng(3)
print("%-8s %-14s %-14s %-14s" % ("voxel", "shipping", "scale only", "scale+clamp"))
for vs in (0.25, 0.5, 1.0, 2.0, 4.0, 8.0):
out = {}
for mode in ("shipping", "scale_only", "fixed"):
errs = []
r2 = np.random.default_rng(3)
for _ in range(400):
dv = r2.normal(size=3); dv /= np.linalg.norm(dv)
start = dv * (R + r2.uniform(-0.5, 0.5) * vs)
p = polish(start, vs, mode)
errs.append(abs(np.linalg.norm(p) - R))
out[mode] = np.mean(errs)
print("%-8.2f %-14.5f %-14.5f %-14.5f"
% (vs, out["shipping"], out["scale_only"], out["fixed"]))
# The tent signature: on a FLAT floor, how far apart do neighbouring cells end
# up vertically? That difference IS the pyramid height.
print()
print("FLAT FLOOR — vertical spread between neighbouring cells (= tent height)")
def flat_density(p):
return -p[2] + 0.37 # a plane that does NOT sit on a grid line
def flat_polish(start, vs, mode):
local = np.array(start, dtype=float)
max_step = 0.5 * vs
for _ in range(ITERS):
d = flat_density(local)
eps = vs * 0.5
g = np.array([0.0, 0.0,
flat_density(local + np.array([0, 0, eps]))
- flat_density(local - np.array([0, 0, eps]))])
gsq = g @ g
if gsq < 1e-12:
break
if mode == "shipping":
local -= np.clip(d / gsq, -max_step, max_step) * g
else:
dx = -(d * vs / gsq) * g
dl = np.linalg.norm(dx)
if dl > max_step:
dx *= max_step / dl
local += dx
return local
print("%-8s %-22s %-22s" % ("voxel", "shipping spread", "fixed spread"))
for vs in (0.25, 0.5, 1.0, 2.0, 4.0):
res = {}
for mode in ("shipping", "fixed"):
zs = []
for k in range(24):
z0 = (k % 4) * vs * 0.25 # different start height per cell
p = flat_polish(np.array([0.0, 0.0, z0]), vs, mode)
zs.append(p[2])
res[mode] = max(zs) - min(zs)
print("%-8.2f %-22.5f %-22.5f" % (vs, res["shipping"], res["fixed"]))