#!/usr/bin/env python3
"""
SWEEP assumed_slope_mix over the diffusion heightfield.
WHY THIS KNOB. erosion_filter.inc:547 is
gully_slope = mix(slope2, slope2/slope_length * assumed_slope.x, assumed_slope.y)
At mix = 1.0 -- what we ship -- that collapses to slope2/slope_length * 0.7:
the heightmap's slope vector keeps its DIRECTION and loses its MAGNITUDE to a
constant. The filter's own comment (:228) says the real slope then "survives
only through slope_length in the initial mask". So incision is the same depth
on a valley wall as on a plain, which is what "the erosion doesn't follow the
heightmap's large-scale structures" looks like.
At mix = 0 the real slope magnitude drives the gullies.
WHAT IS MEASURED, so the picture is not the only evidence:
relief sd of the surface -- does the knob change the terrain at all
|erosion| rms of (this render - the erosion-OFF render), in metres:
how much the filter is actually moving the surface
corr(|ero|,slope) correlation between per-pixel erosion magnitude and the
UNDERLYING heightmap slope. THIS is the question: at mix 1.0
it should be near zero (erosion ignores steepness), and it
should climb as mix falls. If it does not move, the knob is
not the cause and this sweep says so.
python sweep_assumed_slope.py [--mpp 15] [--n 512]
"""
import argparse, os, subprocess, sys
import numpy as np
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"))
MIXES = [1.0, 0.75, 0.5, 0.25, 0.0]
# THE MAP'S CONSTANTS, not the renderer's defaults.
#
# erosion_render's ErosionWrapperParams defaults are the filter's own, and they
# put 633 m rms of incision on the field. map_dc_diffusion.gd asks for 30 m.
# Sweeping a knob at the wrong operating point measures a terrain nobody plays,
# so every value below is transcribed from scenes/maps/map_dc_diffusion.gd and
# the two DERIVED ones are derived here the same way the map derives them.
FLOOR_Y, CEILING_Y, TERRAIN_HEIGHT = -10900.0, 10100.0, 60000.0
GULLY_METRES, DEPTH_METRES = 60.0, 30.0
CELL_SCALE, OCTAVES, GAIN, HEIGHT_SCALE = 0.7, 7, 0.5, 5.0
# The FIELD's half-relief, the SLOPE's vertical unit. Measured over all 841
# cached tiles: -4829..5110 m. NOT the playable y_range, which is 2.1x larger.
FIELD_HALF_RANGE = 4970.0
EROSION_H = (CEILING_Y - FLOOR_Y) * 0.25 * 0.5 # map :313
MAP_PARAMS = [
("floor_y", FLOOR_Y),
("ceiling_y", CEILING_Y),
("terrain_height", TERRAIN_HEIGHT),
("en_height_scale", HEIGHT_SCALE),
("en_strength", DEPTH_METRES / (2.0 * EROSION_H)), # map :602
("en_erosion_scale", GULLY_METRES / (TERRAIN_HEIGHT * CELL_SCALE)),
("en_cell_scale", CELL_SCALE),
("en_octaves", OCTAVES),
("en_gain", GAIN),
("field_half_range", FIELD_HALF_RANGE),
]
def map_args():
out = []
for k, v in MAP_PARAMS:
out += ["--set", k, str(v)]
return out
def render(binpath, ox, oz, n, mpp, out, extra=()):
cmd = [binpath, "--heightmap", HBLK, "--ox", str(ox), "--oz", str(oz),
"--n", str(n), "--mpp", str(mpp), "--out", out] + list(extra)
r = subprocess.run(cmd, cwd=DC, capture_output=True, text=True)
if r.returncode != 0:
print(r.stdout, r.stderr); raise SystemExit("renderer failed: " + " ".join(cmd))
return np.fromfile(out, dtype="<f4", count=n*n).astype(np.float64).reshape(n, n)
def tint(h, lo, hi, step):
t = np.clip((h - lo) / max(hi - lo, 1e-6), 0, 1)
stops = np.array([[0.05,0.12,0.28],[0.20,0.35,0.45],[0.30,0.42,0.22],
[0.55,0.50,0.30],[0.60,0.55,0.52],[0.92,0.92,0.94]])
idx = t * (len(stops) - 1)
i0 = np.floor(idx).astype(int); i1 = np.minimum(i0 + 1, len(stops) - 1)
f = (idx - i0)[..., None]
rgb = stops[i0] * (1 - f) + stops[i1] * f
gy, gx = np.gradient(h, step)
nrm = 1.0 / np.sqrt(gx*gx + gy*gy + 1.0)
shade = np.clip(0.45 + 0.75 * ((-gx - gy) * nrm * 0.5 + 0.5), 0, 1.35)
return np.clip(rgb * shade[..., None], 0, 1)
def main():
ap = argparse.ArgumentParser()
ap.add_argument("--bin", default="/tmp/erosion_render")
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=15.0)
ap.add_argument("--out", default="/tmp/sweep_assumed_slope.png")
ap.add_argument("--defaults", action="store_true",
help="use erosion_render's own defaults instead of the map's")
a = ap.parse_args()
base_args = [] if a.defaults else map_args()
print("parameters: " + ("erosion_render DEFAULTS" if a.defaults
else "map_dc_diffusion.gd (gully %g m, depth %g m, strength %.6g)"
% (GULLY_METRES, DEPTH_METRES, DEPTH_METRES/(2.0*EROSION_H))))
span = a.n * a.mpp
print(f"window {span/1000:.2f} km at ({a.ox:.0f}, {a.oz:.0f}), "
f"{a.n}x{a.n} @ {a.mpp} m/px\n")
# The UNDERLYING heightmap slope, from the source field itself -- not from
# any render -- so it is independent of what the filter did.
arr, meta = read_hblk(HBLK)
ii = ((np.arange(a.n) * a.mpp + a.ox) - meta["origin_x"]) / meta["mpp"]
jj = ((np.arange(a.n) * a.mpp + a.oz) - meta["origin_z"]) / meta["mpp"]
gi = np.clip(np.round(ii).astype(int), 0, meta["side"]-1)
gj = np.clip(np.round(jj).astype(int), 0, meta["side"]-1)
hsrc = arr[np.ix_(gj, gi)].astype(np.float64)
sz, sx = np.gradient(hsrc, a.mpp)
src_slope = np.sqrt(sx*sx + sz*sz)
# erosion OFF: the baseline the filter is displacing from.
base = render(a.bin, a.ox, a.oz, a.n, a.mpp, "/tmp/sw_off.f32",
base_args + ["--set", "en_strength", "0"])
print(f"{'mix':>6} {'relief sd':>11} {'|erosion| rms':>15} "
f"{'max |ero|':>11} {'corr(|ero|,slope)':>19}")
rows, imgs = [], []
for m in MIXES:
h = render(a.bin, a.ox, a.oz, a.n, a.mpp, f"/tmp/sw_{m}.f32",
base_args + ["--set", "assumed_slope_mix", str(m)])
ero = h - base
aero = np.abs(ero)
c = np.corrcoef(aero.ravel(), src_slope.ravel())[0, 1]
print(f"{m:>6.2f} {h.std():>11.1f} {np.sqrt((ero**2).mean()):>15.2f} "
f"{aero.max():>11.2f} {c:>19.3f}")
rows.append((m, c)); imgs.append(h)
lo = min(x.min() for x in imgs); hi = max(x.max() for x in imgs)
strip = np.concatenate([tint(x, lo, hi, a.mpp) for x in imgs], axis=1)
try:
from PIL import Image
Image.fromarray((strip*255).astype(np.uint8)).save(a.out)
print(f"\nwrote {a.out} (mix {', '.join(str(m) for m in MIXES)}, left to right)")
except ImportError:
print("\n(PIL missing — numbers only)")
d = rows[-1][1] - rows[0][1]
print(f"\ncorr moves {rows[0][1]:+.3f} (mix 1.0) -> {rows[-1][1]:+.3f} (mix 0.0), "
f"delta {d:+.3f}")
print("READS AS: the knob " + ("DOES" if abs(d) > 0.05 else "DOES NOT") +
" change how far the erosion follows the heightmap's slope.")
if __name__ == "__main__":
main()