#!/usr/bin/env python3
"""
THE REFERENCE vs OUR PATH, one window, one statistic.
erode_heightmap.cpp is the tool from #50 that ran the filter on the diffusion
field and produced a result worth keeping. Our engine path runs the SAME
erosion_filter.inc through erosion_wrapper_glsl.h's heightmap_gn. They differ
only in the unit conversion handed to the filter:
erode_heightmap heightmap_gn
vertical H 0.5*(hi-lo), the TILE y_range*EN_HAMP*0.5, a CONSTANT
horizontal L gully/cell = 85.7 m M_TH = 60000 m
slope factor L/H ~0.025 ~22.9
If ~900x more slope saturates the filter's slope-gated masks, the reference
should FOLLOW the field's existing valleys and ours should not.
THE STATISTIC IS THE SAME FOR BOTH and is computed against the SOURCE field, so
neither implementation can share a bug with its own judge:
corr(|delta|, concavity) concavity = -laplacian; positive = already a valley
THE SHUFFLE CONTROL must collapse to ~0 for both, or the number is measuring the
window rather than the filter.
python compare_reference.py
"""
import argparse, os, subprocess
import numpy as np
from sweep_assumed_slope import render, MAP_PARAMS
from stitch_block import read_hblk
HERE = os.path.dirname(os.path.abspath(__file__))
DC = os.path.abspath(os.path.join(HERE, "..", "dc-terrain"))
HBLK = os.path.abspath(os.path.join(HERE, "..", "hyper-tps", "assets",
"terrain", "diffusion_gully_30km.hblk"))
MAP = dict(MAP_PARAMS)
MEAS_RMS, MEAS_STRENGTH = 0.227, MAP["en_strength"]
def args_with(**over):
d = dict(MAP_PARAMS); d.update(over)
out = []
for k, v in d.items():
out += ["--set", k, str(v)]
return out
def stats(delta, concav, slope, rng):
ae = np.abs(delta)
return (float(np.sqrt((delta**2).mean())),
float(np.corrcoef(ae.ravel(), concav.ravel())[0, 1]),
float(np.corrcoef(ae.ravel(), slope.ravel())[0, 1]),
float(np.corrcoef(rng.permutation(ae.ravel()), concav.ravel())[0, 1]))
def main():
ap = argparse.ArgumentParser()
ap.add_argument("--ox", type=float, default=50000.0)
ap.add_argument("--oz", type=float, default=34000.0)
ap.add_argument("--n", type=int, default=512)
ap.add_argument("--mpp", type=float, default=30.0)
a = ap.parse_args()
rng = np.random.default_rng(0)
# ---- the source window, at the field's own 30 m pitch ------------------
arr, meta = read_hblk(HBLK)
i0 = int(round((a.ox - meta["origin_x"]) / meta["mpp"]))
j0 = int(round((a.oz - meta["origin_z"]) / meta["mpp"]))
step = int(round(a.mpp / meta["mpp"])) or 1
src = arr[j0:j0+a.n*step:step, i0:i0+a.n*step:step].astype(np.float32)
assert src.shape == (a.n, a.n), src.shape
src.tofile("/tmp/cr_in.f32")
h = src.astype(np.float64)
sz, sx = np.gradient(h, a.mpp)
slope = np.sqrt(sx*sx + sz*sz)
concav = -(np.gradient(np.gradient(h, a.mpp, axis=0), a.mpp, axis=0) +
np.gradient(np.gradient(h, a.mpp, axis=1), a.mpp, axis=1))
print(f"window {a.n*a.mpp/1000:.2f} km at ({a.ox:.0f},{a.oz:.0f}) @ {a.mpp} m/px")
print(f" relief {h.min():.0f} .. {h.max():.0f} m\n")
rows = []
# ---- A. the REFERENCE, at its own defaults ----------------------------
for gully, strength, oct_ in [(900.0, 0.22, 5), (60.0, 0.22, 7)]:
r = subprocess.run(["/tmp/erode_heightmap", "--in", "/tmp/cr_in.f32",
"--n", str(a.n), "--mpp", str(a.mpp),
"--out", "/tmp/cr_ref.f32",
"--gully-metres", str(gully),
"--strength", str(strength),
"--octaves", str(oct_)],
cwd=DC, capture_output=True, text=True)
if r.returncode != 0:
print(r.stdout, r.stderr); raise SystemExit("erode_heightmap failed")
out = np.fromfile("/tmp/cr_ref.f32", dtype="<f4",
count=a.n*a.n).astype(np.float64).reshape(a.n, a.n)
rows.append((f"reference gully {gully:.0f} m, str {strength}",) +
stats(out - h, concav, slope, rng))
# ---- B. OUR path, at the map's settings and at 10 m depth --------------
off = render("/tmp/erosion_render", a.ox, a.oz, a.n, a.mpp,
"/tmp/cr_off.f32", args_with(en_strength=0))
for label, st in [("ours map settings", MEAS_STRENGTH),
("ours at 10 m depth", MEAS_STRENGTH*(10.0/MEAS_RMS))]:
g = render("/tmp/erosion_render", a.ox, a.oz, a.n, a.mpp,
f"/tmp/cr_our_{st:.4g}.f32", args_with(en_strength=st))
rows.append((label,) + stats(g - off, concav, slope, rng))
print(f"{'':>34} {'rms m':>8} {'corr concav':>12} {'corr slope':>11} {'shuffle':>9}")
for name, rms, cc, cs, sh in rows:
print(f"{name:>34} {rms:>8.2f} {cc:>12.3f} {cs:>11.3f} {sh:>9.3f}")
# THE VERDICT KEYS ON SLOPE, NOT CONCAVITY.
#
# The first version of this keyed on concavity and printed "the unit
# conversion is NOT the explanation" -- while the slope column in the SAME
# table showed 0.559/0.368 for the reference against 0.020 for ours. The
# conclusion was false and the evidence refuting it was already on screen.
#
# Concavity was the wrong question: NEITHER implementation deepens ground
# that is already a valley, so "enhances existing lines" is not what the
# reference does either. What it does is cut HARDER ON STEEPER GROUND, and
# that is what reads as following the terrain.
ref_s = max(abs(r[3]) for r in rows[:2])
our_s = max(abs(r[3]) for r in rows[2:])
ref_c = max(abs(r[2]) for r in rows[:2])
our_c = max(abs(r[2]) for r in rows[2:])
print(f"\ncorr SLOPE reference {ref_s:.3f} ours {our_s:.3f}"
f" ({ref_s/our_s if our_s else float('inf'):.1f}x)")
print(f"corr CONCAV reference {ref_c:.3f} ours {our_c:.3f}"
" (both ~0 -- neither follows existing valleys)")
if our_s < 0.5 * ref_s:
print("\nREADS AS: ours cuts far more uniformly than the reference. The\n"
" erosion is placed without regard to the terrain under it.")
else:
print("\nREADS AS: ours now follows the terrain's slope comparably to the\n"
" reference.")
if __name__ == "__main__":
main()