d72cb94 · terrain-lookdev

Changed files (6) +734 -0

amplitude_law.py +112 -0

…	…	@@ -0,0 +1,112 @@
	1	#!/usr/bin/env python3
	2	"""
	3	WHAT ACTUALLY SETS THE EROSION'S DEPTH?
	4
	5	map_dc_diffusion.gd:602 derives strength from a depth in metres by
	6
	7	strength = DEPTH_METRES / (2 * EROSION_H), EROSION_H = 2625 m
	8
	9	i.e. it believes depth = 2 * strength * H. At the map's settings that
	10	predicts 30 m. sweep_assumed_slope.py measured 0.23 m rms / 0.39 m peak on the
	11	same field with the same parameters -- about 100x short. So the formula is
	12	missing something.
	13
	14	THE SUSPECT is en_erosion_scale, the gully BAND. It shrank by 105x when the
	15	band moved from 6300 m to 60 m (#55 section 3), and 105 is close to the 130x
	16	shortfall. If depth tracks the band as well as the strength, the map's algebra
	17	is wrong by exactly that factor and nobody would see erosion at a 60 m band.
	18
	19	THIS DOES NOT ASSUME THAT. It varies each one alone, over a wide range, and
	20	reports the measured slope in log-log -- an exponent near 1 means proportional,
	21	near 0 means the amplitude does not care about that parameter at all. The two
	22	sweeps cross at the map's own operating point, which is the control: both must
	23	report the same amplitude there or the harness is inconsistent.
	24
	25	python amplitude_law.py
	26	"""
	27	import argparse, os, subprocess
	28	import numpy as np
	29	from sweep_assumed_slope import (render, MAP_PARAMS, EROSION_H, GULLY_METRES,
	30	DEPTH_METRES, TERRAIN_HEIGHT, CELL_SCALE)
	31
	32	MAP = dict(MAP_PARAMS)
	33	MAP_STRENGTH = MAP["en_strength"]
	34	MAP_ESCALE = MAP["en_erosion_scale"]
	35	FACTORS = [1, 2, 4, 8, 16, 32, 64, 128]
	36
	37
	38	def args_with(**over):
	39	d = dict(MAP_PARAMS); d.update(over)
	40	out = []
	41	for k, v in d.items():
	42	out += ["--set", k, str(v)]
	43	return out
	44
	45
	46	def amplitude(binpath, ox, oz, n, mpp, base, **over):
	47	h = render(binpath, ox, oz, n, mpp, "/tmp/amp.f32", args_with(**over))
	48	e = h - base
	49	return float(np.sqrt((e**2).mean())), float(np.abs(e).max())
	50
	51
	52	def loglog_slope(xs, ys):
	53	xs = np.asarray(xs, float); ys = np.asarray(ys, float)
	54	ok = (xs > 0) & (ys > 0)
	55	if ok.sum() < 2: return float("nan")
	56	return float(np.polyfit(np.log(xs[ok]), np.log(ys[ok]), 1)[0])
	57
	58
	59	def main():
	60	ap = argparse.ArgumentParser()
	61	ap.add_argument("--bin", default="/tmp/erosion_render")
	62	ap.add_argument("--ox", type=float, default=50000.0)
	63	ap.add_argument("--oz", type=float, default=34000.0)
	64	ap.add_argument("--n", type=int, default=512)
	65	ap.add_argument("--mpp", type=float, default=15.0)
	66	a = ap.parse_args()
	67
	68	print(f"window {a.n*a.mpp/1000:.2f} km at ({a.ox:.0f}, {a.oz:.0f}) @ {a.mpp} m/px")
	69	print(f"map operating point: strength {MAP_STRENGTH:.6g}, "
	70	f"erosion_scale {MAP_ESCALE:.6g} (gully {GULLY_METRES:g} m)")
	71	print(f"map PREDICTS depth = 2strengthH = {2MAP_STRENGTHEROSION_H:.2f} m\n")
	72
	73	base = render(a.bin, a.ox, a.oz, a.n, a.mpp, "/tmp/amp_off.f32",
	74	args_with(en_strength=0))
	75
	76	print("A. STRENGTH alone (erosion_scale pinned at the map's)")
	77	print(f"{'xmap':>6} {'strength':>12} {'rms m':>10} {'peak m':>10} {'2sH':>10}")
	78	xs, ys = [], []
	79	for f in FACTORS:
	80	s = MAP_STRENGTH * f
	81	rms, pk = amplitude(a.bin, a.ox, a.oz, a.n, a.mpp, base, en_strength=s)
	82	print(f"{f:>6} {s:>12.6g} {rms:>10.3f} {pk:>10.3f} {2sEROSION_H:>10.2f}")
	83	xs.append(f); ys.append(rms)
	84	slope_s = loglog_slope(xs, ys)
	85
	86	print("\nB. EROSION_SCALE alone (strength pinned at the map's)")
	87	print(f"{'xmap':>6} {'gully m':>10} {'esc':>12} {'rms m':>10} {'peak m':>10}")
	88	xs2, ys2 = [], []
	89	for f in FACTORS:
	90	esc = MAP_ESCALE * f
	91	rms, pk = amplitude(a.bin, a.ox, a.oz, a.n, a.mpp, base, en_erosion_scale=esc)
	92	print(f"{f:>6} {GULLY_METRES*f:>10.0f} {esc:>12.6g} {rms:>10.3f} {pk:>10.3f}")
	93	xs2.append(f); ys2.append(rms)
	94	slope_e = loglog_slope(xs2, ys2)
	95
	96	print(f"\nlog-log slope vs strength : {slope_s:+.3f}"
	97	" (1 = proportional, 0 = ignores it)")
	98	print(f"log-log slope vs erosion_scale : {slope_e:+.3f}")
	99	print(f"CONTROL, both sweeps at x1 : {ys[0]:.3f} vs {ys2[0]:.3f} m"
	100	f" {'consistent' if abs(ys[0]-ys2[0]) < 1e-6 else 'INCONSISTENT — harness bug'}")
	101
	102	if ys[0] > 0:
	103	need = DEPTH_METRES / ys[0]
	104	print(f"\nTo reach the map's requested {DEPTH_METRES:g} m at the CURRENT band,"
	105	f" amplitude must rise {need:.1f}x")
	106	if abs(slope_s) > 0.5:
	107	print(f" -> via strength: {MAP_STRENGTH:.6g} * {need**(1/slope_s):.1f}"
	108	f" = {MAP_STRENGTHneed*(1/slope_s):.6g}")
	109
	110
	111	if __name__ == "__main__":
	112	main()

band_strip.py +114 -0

…	…	@@ -0,0 +1,114 @@
	1	#!/usr/bin/env python3
	2	"""
	3	THE GULLY BAND, swept, with DEPTH HELD CONSTANT.
	4
	5	The depth strip settled depth at 10 m. What it could not show is that all seven
	6	octaves sit BELOW the band: at 60 m they run 60, 30, 15, 7.5, 3.8, 1.9, 0.94 m,
	7	so nothing the filter draws is larger than a 60 m gully and the result reads as
	8	one frequency of texture laid over the diffusion field rather than as a valley
	9	system. That is the "not multifractal enough" complaint, and it is a property of
	10	where the band sits, not of how deep the cuts are.
	11
	12	DEPTH IS HELD, WHICH IS THE WHOLE POINT. amplitude_law.py measured incision
	13	proportional to strength AND to erosion_scale, so widening the band at fixed
	14	strength deepens everything and the panels would differ in two ways at once.
	15	Strength is therefore scaled by the inverse of the band, and the achieved rms
	16	is printed per panel so you can see it held.
	17
	18	THE OCTAVE LADDER is printed for each band -- that is the thing being chosen.
	19
	20	python band_strip.py [--depth 10]
	21	"""
	22	import argparse, os
	23	import numpy as np
	24	from sweep_assumed_slope import (render, MAP_PARAMS, GULLY_METRES,
	25	TERRAIN_HEIGHT, CELL_SCALE)
	26
	27	MAP = dict(MAP_PARAMS)
	28	MEAS_RMS, MEAS_STRENGTH = 0.227, MAP["en_strength"] # amplitude_law.py, at 60 m
	29	BANDS = [60.0, 120.0, 180.0, 360.0]
	30
	31
	32	def args_with(**over):
	33	d = dict(MAP_PARAMS); d.update(over)
	34	out = []
	35	for k, v in d.items():
	36	out += ["--set", k, str(v)]
	37	return out
	38
	39
	40	def tint(h, lo, hi, step):
	41	t = np.clip((h - lo) / max(hi - lo, 1e-6), 0, 1)
	42	stops = np.array([[0.05,0.12,0.28],[0.20,0.35,0.45],[0.30,0.42,0.22],
	43	[0.55,0.50,0.30],[0.60,0.55,0.52],[0.92,0.92,0.94]])
	44	idx = t * (len(stops) - 1)
	45	i0 = np.floor(idx).astype(int); i1 = np.minimum(i0 + 1, len(stops) - 1)
	46	f = (idx - i0)[..., None]
	47	rgb = stops[i0] * (1 - f) + stops[i1] * f
	48	gy, gx = np.gradient(h, step)
	49	n = 1.0 / np.sqrt(gxgx + gygy + 1.0)
	50	shade = np.clip(0.40 + 0.85 * ((-gx - gy) * n * 0.5 + 0.5), 0, 1.4)
	51	return np.clip(rgb * shade[..., None], 0, 1)
	52
	53
	54	def main():
	55	ap = argparse.ArgumentParser()
	56	ap.add_argument("--bin", default="/tmp/erosion_render")
	57	ap.add_argument("--ox", type=float, default=50000.0)
	58	ap.add_argument("--oz", type=float, default=34000.0)
	59	ap.add_argument("--n", type=int, default=768)
	60	ap.add_argument("--mpp", type=float, default=8.0)
	61	ap.add_argument("--depth", type=float, default=10.0, help="target rms, metres")
	62	ap.add_argument("--octaves", type=int, default=int(MAP["en_octaves"]))
	63	ap.add_argument("--outdir", default="/tmp/bandstrip")
	64	a = ap.parse_args()
	65	os.makedirs(a.outdir, exist_ok=True)
	66
	67	print(f"window {a.n*a.mpp/1000:.2f} km at ({a.ox:.0f}, {a.oz:.0f}) @ {a.mpp} m/px")
	68	print(f"depth held at {a.depth:g} m rms, {a.octaves} octaves\n")
	69
	70	off = render(a.bin, a.ox, a.oz, a.n, a.mpp, "/tmp/bs_off.f32",
	71	args_with(en_strength=0))
	72	panels = [("erosion OFF", off, 0.0, 0.0, 0.0, "")]
	73
	74	for band in BANDS:
	75	esc = band / (TERRAIN_HEIGHT * CELL_SCALE)
	76	# depth ~ strength * erosion_scale, so hold the product.
	77	s = MEAS_STRENGTH * (a.depth / MEAS_RMS) * (GULLY_METRES / band)
	78	h = render(a.bin, a.ox, a.oz, a.n, a.mpp, f"/tmp/bs_{band:g}.f32",
	79	args_with(en_strength=s, en_erosion_scale=esc,
	80	en_octaves=a.octaves))
	81	e = h - off
	82	ladder = " ".join(f"{band/(2**k):.0f}" for k in range(a.octaves))
	83	panels.append((f"band {band:g} m", h, s,
	84	float(np.sqrt((e**2).mean())), float(np.abs(e).max()),
	85	ladder))
	86
	87	print(f"{'panel':>14} {'strength':>11} {'rms m':>8} {'peak m':>8} {'d:w':>7} octaves (m)")
	88	for name, _, s, rms, pk, lad in panels:
	89	band = float(name.split()[1]) if "band" in name else 1.0
	90	print(f"{name:>14} {s:>11.5g} {rms:>8.2f} {pk:>8.2f} "
	91	f"{(pk/band if 'band' in name else 0):>7.3f} {lad}")
	92
	93	lo = min(p[1].min() for p in panels); hi = max(p[1].max() for p in panels)
	94	imgs = [tint(p[1], lo, hi, a.mpp) for p in panels]
	95	strip = np.concatenate(imgs, axis=1)
	96
	97	from PIL import Image, ImageDraw
	98	im = Image.fromarray((strip*255).astype(np.uint8))
	99	d = ImageDraw.Draw(im)
	100	for k, (name, _, s, rms, pk, lad) in enumerate(panels):
	101	x = k * a.n + 8
	102	d.rectangle([x-6, 4, x + 300, 50], fill=(0, 0, 0))
	103	d.text((x, 8), name, fill=(255, 255, 255))
	104	d.text((x, 22), f"rms {rms:.2f} m peak {pk:.2f} m", fill=(210, 210, 210))
	105	d.text((x, 36), f"oct: {lad}", fill=(150, 150, 150))
	106	Image.fromarray((imgs[k]*255).astype(np.uint8)).save(
	107	os.path.join(a.outdir, f"band{k}.png"))
	108	out = os.path.join(a.outdir, "strip.png")
	109	im.save(out)
	110	print(f"\nwrote {out} ({im.width}x{im.height}, shared tint {lo:.0f}..{hi:.0f} m)")
	111
	112
	113	if __name__ == "__main__":
	114	main()

compare_reference.py +126 -0

…	…	@@ -0,0 +1,126 @@
	1	#!/usr/bin/env python3
	2	"""
	3	THE REFERENCE vs OUR PATH, one window, one statistic.
	4
	5	erode_heightmap.cpp is the tool from #50 that ran the filter on the diffusion
	6	field and produced a result worth keeping. Our engine path runs the SAME
	7	erosion_filter.inc through erosion_wrapper_glsl.h's heightmap_gn. They differ
	8	only in the unit conversion handed to the filter:
	9
	10	erode_heightmap heightmap_gn
	11	vertical H 0.5(hi-lo), the TILE y_rangeEN_HAMP*0.5, a CONSTANT
	12	horizontal L gully/cell = 85.7 m M_TH = 60000 m
	13	slope factor L/H ~0.025 ~22.9
	14
	15	If ~900x more slope saturates the filter's slope-gated masks, the reference
	16	should FOLLOW the field's existing valleys and ours should not.
	17
	18	THE STATISTIC IS THE SAME FOR BOTH and is computed against the SOURCE field, so
	19	neither implementation can share a bug with its own judge:
	20
	21	corr(\|delta\|, concavity) concavity = -laplacian; positive = already a valley
	22
	23	THE SHUFFLE CONTROL must collapse to ~0 for both, or the number is measuring the
	24	window rather than the filter.
	25
	26	python compare_reference.py
	27	"""
	28	import argparse, os, subprocess
	29	import numpy as np
	30	from sweep_assumed_slope import render, MAP_PARAMS
	31	from stitch_block import read_hblk
	32
	33	HERE = os.path.dirname(os.path.abspath(__file__))
	34	DC = os.path.abspath(os.path.join(HERE, "..", "dc-terrain"))
	35	HBLK = os.path.abspath(os.path.join(HERE, "..", "hyper-tps", "assets",
	36	"terrain", "diffusion_gully_30km.hblk"))
	37	MAP = dict(MAP_PARAMS)
	38	MEAS_RMS, MEAS_STRENGTH = 0.227, MAP["en_strength"]
	39
	40
	41	def args_with(**over):
	42	d = dict(MAP_PARAMS); d.update(over)
	43	out = []
	44	for k, v in d.items():
	45	out += ["--set", k, str(v)]
	46	return out
	47
	48
	49	def stats(delta, concav, slope, rng):
	50	ae = np.abs(delta)
	51	return (float(np.sqrt((delta**2).mean())),
	52	float(np.corrcoef(ae.ravel(), concav.ravel())[0, 1]),
	53	float(np.corrcoef(ae.ravel(), slope.ravel())[0, 1]),
	54	float(np.corrcoef(rng.permutation(ae.ravel()), concav.ravel())[0, 1]))
	55
	56
	57	def main():
	58	ap = argparse.ArgumentParser()
	59	ap.add_argument("--ox", type=float, default=50000.0)
	60	ap.add_argument("--oz", type=float, default=34000.0)
	61	ap.add_argument("--n", type=int, default=512)
	62	ap.add_argument("--mpp", type=float, default=30.0)
	63	a = ap.parse_args()
	64	rng = np.random.default_rng(0)
	65
	66	# ---- the source window, at the field's own 30 m pitch ------------------
	67	arr, meta = read_hblk(HBLK)
	68	i0 = int(round((a.ox - meta["origin_x"]) / meta["mpp"]))
	69	j0 = int(round((a.oz - meta["origin_z"]) / meta["mpp"]))
	70	step = int(round(a.mpp / meta["mpp"])) or 1
	71	src = arr[j0:j0+a.nstep:step, i0:i0+a.nstep:step].astype(np.float32)
	72	assert src.shape == (a.n, a.n), src.shape
	73	src.tofile("/tmp/cr_in.f32")
	74
	75	h = src.astype(np.float64)
	76	sz, sx = np.gradient(h, a.mpp)
	77	slope = np.sqrt(sxsx + szsz)
	78	concav = -(np.gradient(np.gradient(h, a.mpp, axis=0), a.mpp, axis=0) +
	79	np.gradient(np.gradient(h, a.mpp, axis=1), a.mpp, axis=1))
	80	print(f"window {a.n*a.mpp/1000:.2f} km at ({a.ox:.0f},{a.oz:.0f}) @ {a.mpp} m/px")
	81	print(f" relief {h.min():.0f} .. {h.max():.0f} m\n")
	82
	83	rows = []
	84
	85	# ---- A. the REFERENCE, at its own defaults ----------------------------
	86	for gully, strength, oct_ in [(900.0, 0.22, 5), (60.0, 0.22, 7)]:
	87	r = subprocess.run(["/tmp/erode_heightmap", "--in", "/tmp/cr_in.f32",
	88	"--n", str(a.n), "--mpp", str(a.mpp),
	89	"--out", "/tmp/cr_ref.f32",
	90	"--gully-metres", str(gully),
	91	"--strength", str(strength),
	92	"--octaves", str(oct_)],
	93	cwd=DC, capture_output=True, text=True)
	94	if r.returncode != 0:
	95	print(r.stdout, r.stderr); raise SystemExit("erode_heightmap failed")
	96	out = np.fromfile("/tmp/cr_ref.f32", dtype="<f4",
	97	count=a.n*a.n).astype(np.float64).reshape(a.n, a.n)
	98	rows.append((f"reference gully {gully:.0f} m, str {strength}",) +
	99	stats(out - h, concav, slope, rng))
	100
	101	# ---- B. OUR path, at the map's settings and at 10 m depth --------------
	102	off = render("/tmp/erosion_render", a.ox, a.oz, a.n, a.mpp,
	103	"/tmp/cr_off.f32", args_with(en_strength=0))
	104	for label, st in [("ours map settings", MEAS_STRENGTH),
	105	("ours at 10 m depth", MEAS_STRENGTH*(10.0/MEAS_RMS))]:
	106	g = render("/tmp/erosion_render", a.ox, a.oz, a.n, a.mpp,
	107	f"/tmp/cr_our_{st:.4g}.f32", args_with(en_strength=st))
	108	rows.append((label,) + stats(g - off, concav, slope, rng))
	109
	110	print(f"{'':>34} {'rms m':>8} {'corr concav':>12} {'corr slope':>11} {'shuffle':>9}")
	111	for name, rms, cc, cs, sh in rows:
	112	print(f"{name:>34} {rms:>8.2f} {cc:>12.3f} {cs:>11.3f} {sh:>9.3f}")
	113
	114	ref_best = max(abs(r[2]) for r in rows[:2])
	115	our_best = max(abs(r[2]) for r in rows[2:])
	116	print(f"\nreference best \|corr concav\| {ref_best:.3f} ours {our_best:.3f}")
	117	if ref_best > 0.15 and our_best < 0.10:
	118	print("CONFIRMED: the reference follows the field's existing valleys and\n"
	119	" our path does not. The difference is the unit conversion.")
	120	elif ref_best < 0.15:
	121	print("NOT CONFIRMED: the reference does not score high either, so the\n"
	122	" unit conversion is NOT the explanation. Look elsewhere.")
	123
	124
	125	if __name__ == "__main__":
	126	main()

depth_strip.py +109 -0

…	…	@@ -0,0 +1,109 @@
	1	#!/usr/bin/env python3
	2	"""
	3	THE 60 m BAND AT THREE DEPTHS, side by side, plus the field with no erosion.
	4
	5	amplitude_law.py measured the law: incision is proportional to strength (log-log
	6	slope +1.000) and to erosion_scale (+1.003), and at the map's operating point it
	7	delivers 0.227 m rms / 0.392 m peak. So a target depth in metres is reachable by
	8	scaling strength alone, at the band the map already asks for.
	9
	10	THE TARGET IS RMS, because that is the number every measurement in this thread
	11	has quoted. Peak is printed beside it for every panel, so nothing here depends
	12	on which one you mean by "depth" -- both are on the label.
	13
	14	THE OFF PANEL IS THE CONTROL and it is not decoration: at the map's shipped
	15	0.227 m it should be nearly indistinguishable from the 3 m panel's neighbour,
	16	and if OFF and 3 m look identical then the difference being judged is not the
	17	erosion.
	18
	19	python depth_strip.py
	20	"""
	21	import argparse, os
	22	import numpy as np
	23	from sweep_assumed_slope import render, MAP_PARAMS, EROSION_H, GULLY_METRES
	24
	25	MAP = dict(MAP_PARAMS)
	26	MEAS_RMS, MEAS_STRENGTH = 0.227, MAP["en_strength"] # amplitude_law.py
	27	TARGETS = [3.0, 10.0, 30.0]
	28
	29
	30	def args_with(**over):
	31	d = dict(MAP_PARAMS); d.update(over)
	32	out = []
	33	for k, v in d.items():
	34	out += ["--set", k, str(v)]
	35	return out
	36
	37
	38	def tint(h, lo, hi, step):
	39	t = np.clip((h - lo) / max(hi - lo, 1e-6), 0, 1)
	40	stops = np.array([[0.05,0.12,0.28],[0.20,0.35,0.45],[0.30,0.42,0.22],
	41	[0.55,0.50,0.30],[0.60,0.55,0.52],[0.92,0.92,0.94]])
	42	idx = t * (len(stops) - 1)
	43	i0 = np.floor(idx).astype(int); i1 = np.minimum(i0 + 1, len(stops) - 1)
	44	f = (idx - i0)[..., None]
	45	rgb = stops[i0] * (1 - f) + stops[i1] * f
	46	gy, gx = np.gradient(h, step)
	47	n = 1.0 / np.sqrt(gxgx + gygy + 1.0)
	48	shade = np.clip(0.40 + 0.85 * ((-gx - gy) * n * 0.5 + 0.5), 0, 1.4)
	49	return np.clip(rgb * shade[..., None], 0, 1)
	50
	51
	52	def main():
	53	ap = argparse.ArgumentParser()
	54	ap.add_argument("--bin", default="/tmp/erosion_render")
	55	ap.add_argument("--ox", type=float, default=50000.0)
	56	ap.add_argument("--oz", type=float, default=34000.0)
	57	ap.add_argument("--n", type=int, default=768)
	58	ap.add_argument("--mpp", type=float, default=8.0)
	59	ap.add_argument("--outdir", default="/tmp/depthstrip")
	60	a = ap.parse_args()
	61	os.makedirs(a.outdir, exist_ok=True)
	62
	63	print(f"window {a.n*a.mpp/1000:.2f} km at ({a.ox:.0f}, {a.oz:.0f}) @ {a.mpp} m/px, "
	64	f"gully band {GULLY_METRES:g} m ({GULLY_METRES/a.mpp:.1f} px per gully)")
	65
	66	off = render(a.bin, a.ox, a.oz, a.n, a.mpp, "/tmp/ds_off.f32",
	67	args_with(en_strength=0))
	68
	69	panels = [("erosion OFF", off, 0.0, 0.0, 0.0)]
	70	shipped = render(a.bin, a.ox, a.oz, a.n, a.mpp, "/tmp/ds_ship.f32", args_with())
	71	e = shipped - off
	72	panels.append((f"SHIPPED (asks 30 m)", shipped, MEAS_STRENGTH,
	73	float(np.sqrt((e**2).mean())), float(np.abs(e).max())))
	74
	75	for tgt in TARGETS:
	76	s = MEAS_STRENGTH * (tgt / MEAS_RMS)
	77	h = render(a.bin, a.ox, a.oz, a.n, a.mpp, f"/tmp/ds_{tgt:g}.f32",
	78	args_with(en_strength=s))
	79	e = h - off
	80	panels.append((f"{tgt:g} m", h, s,
	81	float(np.sqrt((e**2).mean())), float(np.abs(e).max())))
	82
	83	print(f"\n{'panel':>20} {'strength':>12} {'rms m':>9} {'peak m':>9} {'d:w':>7}")
	84	for name, _, s, rms, pk in panels:
	85	dw = pk / GULLY_METRES if pk else 0.0
	86	print(f"{name:>20} {s:>12.6g} {rms:>9.3f} {pk:>9.3f} {dw:>7.3f}")
	87
	88	lo = min(p[1].min() for p in panels); hi = max(p[1].max() for p in panels)
	89	imgs = [tint(p[1], lo, hi, a.mpp) for p in panels]
	90	strip = np.concatenate(imgs, axis=1)
	91
	92	from PIL import Image, ImageDraw
	93	im = Image.fromarray((strip*255).astype(np.uint8))
	94	d = ImageDraw.Draw(im)
	95	for k, (name, _, s, rms, pk) in enumerate(panels):
	96	x = k * a.n + 8
	97	d.rectangle([x-6, 4, x + 250, 46], fill=(0, 0, 0))
	98	d.text((x, 8), f"{name}", fill=(255, 255, 255))
	99	d.text((x, 22), f"rms {rms:.2f} m peak {pk:.2f} m", fill=(210, 210, 210))
	100	d.text((x, 34), f"strength {s:.5g}", fill=(160, 160, 160))
	101	Image.fromarray((imgs[k]*255).astype(np.uint8)).save(
	102	os.path.join(a.outdir, f"panel{k}_{name.split()[0]}.png"))
	103	out = os.path.join(a.outdir, "strip.png")
	104	im.save(out)
	105	print(f"\nwrote {out} ({im.width}x{im.height}, shared tint {lo:.0f}..{hi:.0f} m)")
	106
	107
	108	if __name__ == "__main__":
	109	main()

follows_structure.py +114 -0

…	…	@@ -0,0 +1,114 @@
	1	#!/usr/bin/env python3
	2	"""
	3	DOES THE FILTER ENHANCE THE FIELD'S EXISTING LINES, OR ADD ITS OWN?
	4
	5	erosion_filter.inc:549 places gullies with phacelle_noise(p*freq, gully_slope,..):
	6	the PATTERN is noise, merely ORIENTED downhill. So a priori the filter adds a
	7	new gully system rather than deepening the drainage the diffusion field already
	8	has. This measures whether that is true, and whether assumed_slope_mix changes it.
	9
	10	WHY THIS RE-TESTS A KNOB ALREADY "REFUTED". sweep_assumed_slope.py found mix
	11	inert -- but it ran at the shipped operating point, where the erosion moves the
	12	surface 0.227 m. A knob that modulates an effect which is not happening reads as
	13	inert no matter what it does. That measurement proved nothing; this one runs at
	14	10 m, where there is something to modulate.
	15
	16	THE TWO CORRELATIONS, both against the SOURCE field (never against a render, so
	17	they cannot share a bug with what they judge):
	18
	19	corr(\|ero\|, slope) does the filter cut harder on steep ground
	20	corr(\|ero\|, concavity) does it cut harder where the field is ALREADY a
	21	valley. This is the "enhance existing lines" question.
	22	Concavity = -laplacian, so positive = valley floor.
	23
	24	A filter that enhances has concavity correlation well above zero. A filter that
	25	adds an independent pattern has it near zero however deep it cuts.
	26
	27	THE CONTROL is the correlation of the erosion-OFF field with itself-derived
	28	concavity, which is meaningless by construction, plus a SHUFFLED \|ero\| whose
	29	correlation must collapse to ~0. If the shuffle does not collapse, the statistic
	30	is picking up something structural in the window rather than in the filter.
	31
	32	python follows_structure.py [--depth 10]
	33	"""
	34	import argparse, os
	35	import numpy as np
	36	from sweep_assumed_slope import render, MAP_PARAMS
	37	from stitch_block import read_hblk
	38
	39	HERE = os.path.dirname(os.path.abspath(__file__))
	40	HBLK = os.path.abspath(os.path.join(HERE, "..", "hyper-tps", "assets",
	41	"terrain", "diffusion_gully_30km.hblk"))
	42	MAP = dict(MAP_PARAMS)
	43	MEAS_RMS, MEAS_STRENGTH = 0.227, MAP["en_strength"]
	44	MIXES = [1.0, 0.75, 0.5, 0.25, 0.0]
	45
	46
	47	def args_with(**over):
	48	d = dict(MAP_PARAMS); d.update(over)
	49	out = []
	50	for k, v in d.items():
	51	out += ["--set", k, str(v)]
	52	return out
	53
	54
	55	def main():
	56	ap = argparse.ArgumentParser()
	57	ap.add_argument("--bin", default="/tmp/erosion_render")
	58	ap.add_argument("--ox", type=float, default=50000.0)
	59	ap.add_argument("--oz", type=float, default=34000.0)
	60	ap.add_argument("--n", type=int, default=512)
	61	ap.add_argument("--mpp", type=float, default=8.0)
	62	ap.add_argument("--depth", type=float, default=10.0)
	63	a = ap.parse_args()
	64
	65	strength = MEAS_STRENGTH * (a.depth / MEAS_RMS)
	66	print(f"window {a.n*a.mpp/1000:.2f} km @ {a.mpp} m/px, depth {a.depth:g} m rms "
	67	f"(strength {strength:.5g})\n")
	68
	69	# SOURCE field, sampled independently of any render.
	70	arr, meta = read_hblk(HBLK)
	71	ii = ((np.arange(a.n)*a.mpp + a.ox) - meta["origin_x"]) / meta["mpp"]
	72	jj = ((np.arange(a.n)*a.mpp + a.oz) - meta["origin_z"]) / meta["mpp"]
	73	gi = np.clip(np.round(ii).astype(int), 0, meta["side"]-1)
	74	gj = np.clip(np.round(jj).astype(int), 0, meta["side"]-1)
	75	hsrc = arr[np.ix_(gj, gi)].astype(np.float64)
	76	sz, sx = np.gradient(hsrc, a.mpp)
	77	slope = np.sqrt(sxsx + szsz)
	78	lap = (np.gradient(np.gradient(hsrc, a.mpp, axis=0), a.mpp, axis=0) +
	79	np.gradient(np.gradient(hsrc, a.mpp, axis=1), a.mpp, axis=1))
	80	concav = -lap # positive = already a valley
	81
	82	off = render(a.bin, a.ox, a.oz, a.n, a.mpp, "/tmp/fs_off.f32",
	83	args_with(en_strength=0))
	84
	85	rng = np.random.default_rng(0)
	86	print(f"{'mix':>6} {'rms m':>8} {'corr(\|ero\|,slope)':>19} {'corr(\|ero\|,concav)':>20} "
	87	f"{'shuffled':>10}")
	88	out = []
	89	for m in MIXES:
	90	h = render(a.bin, a.ox, a.oz, a.n, a.mpp, f"/tmp/fs_{m}.f32",
	91	args_with(en_strength=strength, assumed_slope_mix=m))
	92	ae = np.abs(h - off)
	93	cs = np.corrcoef(ae.ravel(), slope.ravel())[0, 1]
	94	cc = np.corrcoef(ae.ravel(), concav.ravel())[0, 1]
	95	sh = np.corrcoef(rng.permutation(ae.ravel()), concav.ravel())[0, 1]
	96	print(f"{m:>6.2f} {np.sqrt(((h-off)**2).mean()):>8.2f} {cs:>19.3f} "
	97	f"{cc:>20.3f} {sh:>10.3f}")
	98	out.append((m, cs, cc, sh))
	99
	100	cc0, cc1 = out[0][2], out[-1][2]
	101	print(f"\nconcavity correlation: {cc0:+.3f} (mix 1.0) -> {cc1:+.3f} (mix 0.0)")
	102	print(f"shuffle control : {max(abs(o[3]) for o in out):.3f} "
	103	"(must be ~0 or the statistic is not measuring the filter)")
	104	best = max(out, key=lambda o: abs(o[2]))
	105	if abs(best[2]) < 0.15:
	106	print("\nREADS AS: the filter ADDS an independent pattern. It does not\n"
	107	" find the field's existing valleys at ANY mix setting.")
	108	else:
	109	print(f"\nREADS AS: at mix {best[0]:g} the filter DOES follow existing valleys "
	110	f"(corr {best[2]:+.3f}).")
	111
	112
	113	if __name__ == "__main__":
	114	main()

sweep_assumed_slope.py +159 -0

…	…	@@ -0,0 +1,159 @@
	1	#!/usr/bin/env python3
	2	"""
	3	SWEEP assumed_slope_mix over the diffusion heightfield.
	4
	5	WHY THIS KNOB. erosion_filter.inc:547 is
	6
	7	gully_slope = mix(slope2, slope2/slope_length * assumed_slope.x, assumed_slope.y)
	8
	9	At mix = 1.0 -- what we ship -- that collapses to slope2/slope_length * 0.7:
	10	the heightmap's slope vector keeps its DIRECTION and loses its MAGNITUDE to a
	11	constant. The filter's own comment (:228) says the real slope then "survives
	12	only through slope_length in the initial mask". So incision is the same depth
	13	on a valley wall as on a plain, which is what "the erosion doesn't follow the
	14	heightmap's large-scale structures" looks like.
	15
	16	At mix = 0 the real slope magnitude drives the gullies.
	17
	18	WHAT IS MEASURED, so the picture is not the only evidence:
	19	relief sd of the surface -- does the knob change the terrain at all
	20	\|erosion\| rms of (this render - the erosion-OFF render), in metres:
	21	how much the filter is actually moving the surface
	22	corr(\|ero\|,slope) correlation between per-pixel erosion magnitude and the
	23	UNDERLYING heightmap slope. THIS is the question: at mix 1.0
	24	it should be near zero (erosion ignores steepness), and it
	25	should climb as mix falls. If it does not move, the knob is
	26	not the cause and this sweep says so.
	27
	28	python sweep_assumed_slope.py [--mpp 15] [--n 512]
	29	"""
	30	import argparse, os, subprocess, sys
	31	import numpy as np
	32	from stitch_block import read_hblk
	33
	34	HERE = os.path.dirname(os.path.abspath(__file__))
	35	DC = os.path.abspath(os.path.join(HERE, "..", "dc-terrain"))
	36	HBLK = os.path.abspath(os.path.join(HERE, "..", "hyper-tps", "assets",
	37	"terrain", "diffusion_gully_30km.hblk"))
	38	MIXES = [1.0, 0.75, 0.5, 0.25, 0.0]
	39
	40	# THE MAP'S CONSTANTS, not the renderer's defaults.
	41	#
	42	# erosion_render's ErosionWrapperParams defaults are the filter's own, and they
	43	# put 633 m rms of incision on the field. map_dc_diffusion.gd asks for 30 m.
	44	# Sweeping a knob at the wrong operating point measures a terrain nobody plays,
	45	# so every value below is transcribed from scenes/maps/map_dc_diffusion.gd and
	46	# the two DERIVED ones are derived here the same way the map derives them.
	47	FLOOR_Y, CEILING_Y, TERRAIN_HEIGHT = -10900.0, 10100.0, 60000.0
	48	GULLY_METRES, DEPTH_METRES = 60.0, 30.0
	49	CELL_SCALE, OCTAVES, GAIN, HEIGHT_SCALE = 0.7, 7, 0.5, 5.0
	50	EROSION_H = (CEILING_Y - FLOOR_Y) * 0.25 * 0.5 # map :313
	51	MAP_PARAMS = [
	52	("floor_y", FLOOR_Y),
	53	("ceiling_y", CEILING_Y),
	54	("terrain_height", TERRAIN_HEIGHT),
	55	("en_height_scale", HEIGHT_SCALE),
	56	("en_strength", DEPTH_METRES / (2.0 * EROSION_H)), # map :602
	57	("en_erosion_scale", GULLY_METRES / (TERRAIN_HEIGHT * CELL_SCALE)),
	58	("en_cell_scale", CELL_SCALE),
	59	("en_octaves", OCTAVES),
	60	("en_gain", GAIN),
	61	]
	62
	63
	64	def map_args():
	65	out = []
	66	for k, v in MAP_PARAMS:
	67	out += ["--set", k, str(v)]
	68	return out
	69
	70
	71	def render(binpath, ox, oz, n, mpp, out, extra=()):
	72	cmd = [binpath, "--heightmap", HBLK, "--ox", str(ox), "--oz", str(oz),
	73	"--n", str(n), "--mpp", str(mpp), "--out", out] + list(extra)
	74	r = subprocess.run(cmd, cwd=DC, capture_output=True, text=True)
	75	if r.returncode != 0:
	76	print(r.stdout, r.stderr); raise SystemExit("renderer failed: " + " ".join(cmd))
	77	return np.fromfile(out, dtype="<f4", count=n*n).astype(np.float64).reshape(n, n)
	78
	79
	80	def tint(h, lo, hi, step):
	81	t = np.clip((h - lo) / max(hi - lo, 1e-6), 0, 1)
	82	stops = np.array([[0.05,0.12,0.28],[0.20,0.35,0.45],[0.30,0.42,0.22],
	83	[0.55,0.50,0.30],[0.60,0.55,0.52],[0.92,0.92,0.94]])
	84	idx = t * (len(stops) - 1)
	85	i0 = np.floor(idx).astype(int); i1 = np.minimum(i0 + 1, len(stops) - 1)
	86	f = (idx - i0)[..., None]
	87	rgb = stops[i0] * (1 - f) + stops[i1] * f
	88	gy, gx = np.gradient(h, step)
	89	nrm = 1.0 / np.sqrt(gxgx + gygy + 1.0)
	90	shade = np.clip(0.45 + 0.75 * ((-gx - gy) * nrm * 0.5 + 0.5), 0, 1.35)
	91	return np.clip(rgb * shade[..., None], 0, 1)
	92
	93
	94	def main():
	95	ap = argparse.ArgumentParser()
	96	ap.add_argument("--bin", default="/tmp/erosion_render")
	97	ap.add_argument("--ox", type=float, default=50000.0)
	98	ap.add_argument("--oz", type=float, default=34000.0)
	99	ap.add_argument("--n", type=int, default=512)
	100	ap.add_argument("--mpp", type=float, default=15.0)
	101	ap.add_argument("--out", default="/tmp/sweep_assumed_slope.png")
	102	ap.add_argument("--defaults", action="store_true",
	103	help="use erosion_render's own defaults instead of the map's")
	104	a = ap.parse_args()
	105
	106	base_args = [] if a.defaults else map_args()
	107	print("parameters: " + ("erosion_render DEFAULTS" if a.defaults
	108	else "map_dc_diffusion.gd (gully %g m, depth %g m, strength %.6g)"
	109	% (GULLY_METRES, DEPTH_METRES, DEPTH_METRES/(2.0*EROSION_H))))
	110	span = a.n * a.mpp
	111	print(f"window {span/1000:.2f} km at ({a.ox:.0f}, {a.oz:.0f}), "
	112	f"{a.n}x{a.n} @ {a.mpp} m/px\n")
	113
	114	# The UNDERLYING heightmap slope, from the source field itself -- not from
	115	# any render -- so it is independent of what the filter did.
	116	arr, meta = read_hblk(HBLK)
	117	ii = ((np.arange(a.n) * a.mpp + a.ox) - meta["origin_x"]) / meta["mpp"]
	118	jj = ((np.arange(a.n) * a.mpp + a.oz) - meta["origin_z"]) / meta["mpp"]
	119	gi = np.clip(np.round(ii).astype(int), 0, meta["side"]-1)
	120	gj = np.clip(np.round(jj).astype(int), 0, meta["side"]-1)
	121	hsrc = arr[np.ix_(gj, gi)].astype(np.float64)
	122	sz, sx = np.gradient(hsrc, a.mpp)
	123	src_slope = np.sqrt(sxsx + szsz)
	124
	125	# erosion OFF: the baseline the filter is displacing from.
	126	base = render(a.bin, a.ox, a.oz, a.n, a.mpp, "/tmp/sw_off.f32",
	127	base_args + ["--set", "en_strength", "0"])
	128
	129	print(f"{'mix':>6} {'relief sd':>11} {'\|erosion\| rms':>15} "
	130	f"{'max \|ero\|':>11} {'corr(\|ero\|,slope)':>19}")
	131	rows, imgs = [], []
	132	for m in MIXES:
	133	h = render(a.bin, a.ox, a.oz, a.n, a.mpp, f"/tmp/sw_{m}.f32",
	134	base_args + ["--set", "assumed_slope_mix", str(m)])
	135	ero = h - base
	136	aero = np.abs(ero)
	137	c = np.corrcoef(aero.ravel(), src_slope.ravel())[0, 1]
	138	print(f"{m:>6.2f} {h.std():>11.1f} {np.sqrt((ero**2).mean()):>15.2f} "
	139	f"{aero.max():>11.2f} {c:>19.3f}")
	140	rows.append((m, c)); imgs.append(h)
	141
	142	lo = min(x.min() for x in imgs); hi = max(x.max() for x in imgs)
	143	strip = np.concatenate([tint(x, lo, hi, a.mpp) for x in imgs], axis=1)
	144	try:
	145	from PIL import Image
	146	Image.fromarray((strip*255).astype(np.uint8)).save(a.out)
	147	print(f"\nwrote {a.out} (mix {', '.join(str(m) for m in MIXES)}, left to right)")
	148	except ImportError:
	149	print("\n(PIL missing — numbers only)")
	150
	151	d = rows[-1][1] - rows[0][1]
	152	print(f"\ncorr moves {rows[0][1]:+.3f} (mix 1.0) -> {rows[-1][1]:+.3f} (mix 0.0), "
	153	f"delta {d:+.3f}")
	154	print("READS AS: the knob " + ("DOES" if abs(d) > 0.05 else "DOES NOT") +
	155	" change how far the erosion follows the heightmap's slope.")
	156
	157
	158	if __name__ == "__main__":
	159	main()

terrain-lookdev commits / d72cb94

Changed files (6) +734 -0