#!/usr/bin/env python3
"""
THE GULLY BAND, swept, with DEPTH HELD CONSTANT.
The depth strip settled depth at 10 m. What it could not show is that all seven
octaves sit BELOW the band: at 60 m they run 60, 30, 15, 7.5, 3.8, 1.9, 0.94 m,
so nothing the filter draws is larger than a 60 m gully and the result reads as
one frequency of texture laid over the diffusion field rather than as a valley
system. That is the "not multifractal enough" complaint, and it is a property of
where the band sits, not of how deep the cuts are.
DEPTH IS HELD, WHICH IS THE WHOLE POINT. amplitude_law.py measured incision
proportional to strength AND to erosion_scale, so widening the band at fixed
strength deepens everything and the panels would differ in two ways at once.
Strength is therefore scaled by the inverse of the band, and the achieved rms
is printed per panel so you can see it held.
THE OCTAVE LADDER is printed for each band -- that is the thing being chosen.
python band_strip.py [--depth 10]
"""
import argparse, os
import numpy as np
from sweep_assumed_slope import (render, MAP_PARAMS, GULLY_METRES,
TERRAIN_HEIGHT, CELL_SCALE)
MAP = dict(MAP_PARAMS)
MEAS_RMS, MEAS_STRENGTH = 0.227, MAP["en_strength"] # amplitude_law.py, at 60 m
BANDS = [60.0, 120.0, 180.0, 360.0]
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 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)
n = 1.0 / np.sqrt(gx*gx + gy*gy + 1.0)
shade = np.clip(0.40 + 0.85 * ((-gx - gy) * n * 0.5 + 0.5), 0, 1.4)
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=768)
ap.add_argument("--mpp", type=float, default=8.0)
ap.add_argument("--depth", type=float, default=10.0, help="target rms, metres")
ap.add_argument("--octaves", type=int, default=int(MAP["en_octaves"]))
ap.add_argument("--outdir", default="/tmp/bandstrip")
a = ap.parse_args()
os.makedirs(a.outdir, exist_ok=True)
print(f"window {a.n*a.mpp/1000:.2f} km at ({a.ox:.0f}, {a.oz:.0f}) @ {a.mpp} m/px")
print(f"depth held at {a.depth:g} m rms, {a.octaves} octaves\n")
off = render(a.bin, a.ox, a.oz, a.n, a.mpp, "/tmp/bs_off.f32",
args_with(en_strength=0))
panels = [("erosion OFF", off, 0.0, 0.0, 0.0, "")]
for band in BANDS:
esc = band / (TERRAIN_HEIGHT * CELL_SCALE)
# depth ~ strength * erosion_scale, so hold the product.
s = MEAS_STRENGTH * (a.depth / MEAS_RMS) * (GULLY_METRES / band)
h = render(a.bin, a.ox, a.oz, a.n, a.mpp, f"/tmp/bs_{band:g}.f32",
args_with(en_strength=s, en_erosion_scale=esc,
en_octaves=a.octaves))
e = h - off
ladder = " ".join(f"{band/(2**k):.0f}" for k in range(a.octaves))
panels.append((f"band {band:g} m", h, s,
float(np.sqrt((e**2).mean())), float(np.abs(e).max()),
ladder))
print(f"{'panel':>14} {'strength':>11} {'rms m':>8} {'peak m':>8} {'d:w':>7} octaves (m)")
for name, _, s, rms, pk, lad in panels:
band = float(name.split()[1]) if "band" in name else 1.0
print(f"{name:>14} {s:>11.5g} {rms:>8.2f} {pk:>8.2f} "
f"{(pk/band if 'band' in name else 0):>7.3f} {lad}")
lo = min(p[1].min() for p in panels); hi = max(p[1].max() for p in panels)
imgs = [tint(p[1], lo, hi, a.mpp) for p in panels]
strip = np.concatenate(imgs, axis=1)
from PIL import Image, ImageDraw
im = Image.fromarray((strip*255).astype(np.uint8))
d = ImageDraw.Draw(im)
for k, (name, _, s, rms, pk, lad) in enumerate(panels):
x = k * a.n + 8
d.rectangle([x-6, 4, x + 300, 50], fill=(0, 0, 0))
d.text((x, 8), name, fill=(255, 255, 255))
d.text((x, 22), f"rms {rms:.2f} m peak {pk:.2f} m", fill=(210, 210, 210))
d.text((x, 36), f"oct: {lad}", fill=(150, 150, 150))
Image.fromarray((imgs[k]*255).astype(np.uint8)).save(
os.path.join(a.outdir, f"band{k}.png"))
out = os.path.join(a.outdir, "strip.png")
im.save(out)
print(f"\nwrote {out} ({im.width}x{im.height}, shared tint {lo:.0f}..{hi:.0f} m)")
if __name__ == "__main__":
main()