#!/usr/bin/env python3
"""Does the SHIPPING density expression have enough float32 precision to give
dual contouring a usable surface normal?
edge_intersect.glsl computes the normal that feeds the QEF by central
differences on evaluate_density, with eps = voxel_size * 0.01. The shipping
heightfield map's evaluate_density is
A = (H(x,z) + 10900) * 0.000238095 # the node's metres -> band mapping
B = (y + 10900) * 0.000238 # the wrapper's y_mapped
d = (A - B) * 4200 # back into metres
A and B are both about 2.63 near the ground, and their difference is what
carries the geometry. This script does that arithmetic in float32, exactly as
the shader does, and asks how much of the horizontal gradient survives.
It reports the angular error of the resulting normal. Run it against a
synthetic ramp (a known answer: the normal of a plane of slope s) so the check
can FAIL -- an oracle whose input has no known answer proves nothing.
"""
import numpy as np
F = np.float32
YS_NODE = F(0.000238095) # emitted by HeightfieldNode from 5.0/21000
YS_WRAP = F(0.000238) # emitted by String::num(height_scale/y_range, 6)
FLOOR = F(-10900.0)
SCALE = F(4200.0)
def density(H, y):
"""evaluate_density, in float32, term for term as the shader writes it."""
H = F(H); y = F(y)
A = F(F(H - FLOOR) * YS_NODE)
B = F(F(y - FLOOR) * YS_WRAP)
return F(F(A - B) * SCALE)
def normal_from_central_diff(h_of_x, x, y, z, eps):
"""n as edge_intersect.glsl builds it: three central differences."""
nx = F(density(h_of_x(x + eps), y) - density(h_of_x(x - eps), y))
ny = F(density(h_of_x(x), F(y + eps)) - density(h_of_x(x), F(y - eps)))
return nx, ny
def angle_err_deg(nx, ny, slope):
"""Angle between the measured normal and the true one for height slope s.
True surface y = s*x, so the density gradient is (-s, 1) up to scale."""
if nx == 0.0 and ny == 0.0:
return float("nan")
# d is POSITIVE BELOW GROUND (solid), so d ~ (H(x) - y) and its gradient is
# (+s, -1). Comparing against (-s, +1) -- the outward surface normal -- put
# every row 180 degrees out and the float64 control caught it, which is the
# only reason this function is written down rather than assumed.
m = np.arctan2(float(nx), float(ny))
t = np.arctan2(slope, -1.0)
e = abs(np.degrees(m - t)) % 360.0
return min(e, 360.0 - e)
def run(base_h, y_off, voxel, slopes, epsf):
rows = []
for s in slopes:
h = lambda xx, s=s: F(base_h + F(s) * F(xx - F(53760.0)))
x = F(53760.0)
# Stand on the surface the DENSITY defines, not on base_h: the two scale
# constants differ, so the zero of d is not at y = H.
y = F(base_h * (YS_NODE / YS_WRAP)
+ FLOOR * (YS_NODE / YS_WRAP - F(1.0)) + F(y_off))
eps = F(voxel * epsf)
nx, ny = normal_from_central_diff(h, x, y, F(53760.0), eps)
rows.append((s, float(nx), float(ny), angle_err_deg(nx, ny, s)))
return rows
def main():
slopes = [0.003, 0.01, 0.03, 0.1, 0.3, 1.0]
print("Shipping expression, base height 165 m, voxel 1 m")
print(" true surface slope | eps=0.01 voxel (SHIPPING) | eps=0.5 voxel")
print(" %-8s | %-14s %-9s | %-14s %-9s" %
("slope", "n.x", "err deg", "n.x", "err deg"))
a = run(165.0, 0.0, 1.0, slopes, 0.01)
b = run(165.0, 0.0, 1.0, slopes, 0.5)
zero_ship = 0
for (s, nx1, ny1, e1), (_, nx2, ny2, e2) in zip(a, b):
if nx1 == 0.0:
zero_ship += 1
print(" %-8.3f | %-14.3e %-9.2f | %-14.3e %-9.2f" % (s, nx1, e1, nx2, e2))
print()
print("A NORMAL OF EXACTLY ZERO IN X IS A COLLAPSED NORMAL: the QEF is handed")
print("a plane that says the ground is level when it is not.")
print(" shipping eps: %d of %d slopes collapsed to exactly 0" % (zero_ship, len(slopes)))
# THE ORACLE MUST BE ABLE TO FAIL. Same arithmetic in float64, where no
# cancellation occurs: every slope must come back correct. If this row ever
# shows an error, the harness is wrong, not the shader.
print()
print("CONTROL -- the same central difference in float64 (no cancellation):")
ok = True
for s in slopes:
h = lambda xx, s=s: 165.0 + s * (xx - 53760.0)
y = 165.0 * (0.000238095 / 0.000238) + (-10900.0) * (0.000238095 / 0.000238 - 1.0)
eps = 0.01
d = lambda H, yy: ((H + 10900.0) * 0.000238095 - (yy + 10900.0) * 0.000238) * 4200.0
nx = d(h(53760.0 + eps), y) - d(h(53760.0 - eps), y)
ny = d(h(53760.0), y + eps) - d(h(53760.0), y - eps)
e = angle_err_deg(nx, ny, s)
# 0.05 deg, not 0.01. The residual at slope 1.0 is 0.0114 deg and it is
# REAL, not harness error: this control simulates the two MISMATCHED
# constants, whose 4.0e-4 ratio tilts the true surface by that much. A
# tolerance tight enough to flag it would be flagging the defect the
# rows above exist to show.
flag = "ok" if e < 0.05 else "FAIL"
if e >= 0.05:
ok = False
print(" slope %-6.3f n.x %-12.4e err %.4f deg %s" % (s, nx, e, flag))
print("CONTROL: %s" % ("PASS -- the arithmetic is right, the precision is not"
if ok else "FAIL -- fix the harness before believing the rows above"))
# And the second, independent defect the same expression carries.
print()
print("WHERE THE DENSITY PUTS THE GROUND, against the heightfield's own value:")
for label, ys_wrap in (("BEFORE String::num(.., 6) -> 0.000238", 0.000238),
("AFTER String::num(.., 9) -> 0.000238095", 0.000238095)):
r = 0.000238095 / ys_wrap
errs = []
for h in (-5300.0, 0.0, 165.0, 5300.0):
y0 = h * r + (-10900.0) * (r - 1.0)
errs.append(y0 - h)
print(" %-45s %s" % (label,
" ".join("%+.3f" % e for e in errs)))
print(" %-45s %s" % ("", "at h = -5300 / 0 / 165 / 5300 m"))
if __name__ == "__main__":
main()