// ============================================================================
// ERODE A HEIGHTMAP — run runevision's filter on ANY input height field, in
// particular a terrain-diffusion tile.
//
// WHY
//
// #46 §4 measured that the diffusion model produces landscape (coherent
// catchments, beta steepening 4.18 -> 5.10) where our fBm mountain layer
// produces texture (flat beta). The blog is explicit that the filter "can be
// applied on top of any height function". So: let the model supply the large
// forms it is demonstrably good at, and let the filter supply detail below its
// 30 m floor. This is #46 §4b option 1, made runnable.
//
// It needs no GPU descriptor work. A heightmap transform is ONE erosion
// evaluation per pixel with no surface search, so the CPU shim is enough --
// unlike dump_erosion_heightfield, which pays 572 density evaluations per
// column bisecting for a surface.
//
// UNITS, which is the whole difficulty
//
// The filter is dimensionless. Its p is world position over a horizontal unit
// L, so the coarsest gully wavelength is L * erosion_scale * cell_scale metres.
// Its height is normalised by a vertical unit H, and its SLOPE must be in the
// same space: dh/dp = (dh/dx_world) * L / H. Get that factor wrong and the
// gullies come out oriented correctly but scaled absurdly, which is exactly the
// 60x wrapper bug #45 found.
//
// THE IDENTITY GATE
//
// With --strength 0 the filter contributes nothing, so the output must equal
// the input to float precision. --check runs that and refuses to proceed if it
// does not hold. A unit error that survives this cannot exist.
//
// Build:
// cd dc-terrain
// g++ -O2 -std=c++17 -I ../proc-nodes/src/shared -o /tmp/erode_heightmap \
// ../terrain-lookdev/erode_heightmap.cpp
//
// Run:
// /tmp/erode_heightmap --check
// /tmp/erode_heightmap --in tile.f32 --n 512 --mpp 30 --out eroded.f32 \
// --gully-metres 900 --strength 0.22 --octaves 5
// ============================================================================
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <vector>
#include <string>
#include "erosion_cpp_shim.h"
using dcerosion::vec2;
using dcerosion::vec3;
using dcerosion::vec4;
struct Opts {
int n = 512;
float mpp = 30.0f;
// The coarsest gully, in METRES. Everything else about the horizontal
// layout is derived from it, because a raw "scale" number is unreadable --
// the same lesson the lookdev panel relabelling just learned.
float gully_metres = 900.0f;
float strength = 0.22f; // runevision's reference
int octaves = 5; // runevision's reference
float gain = 0.5f;
float cell_scale = 0.7f;
float lacunarity = 2.0f;
float gully_weight = 0.5f;
float detail = 1.5f;
float normalization = 0.5f;
float rounding_input = 0.1f;
float height_offset_a = -0.65f;
// Vertical normalisation. The filter's strength is relative to a height of
// order 1, so H should be the input's own relief, not an arbitrary number.
// 0 means "use the tile's (max-min)/2".
float vscale = 0.0f;
// DEPTH : WIDTH, the scale-invariant way to ask for erosion.
//
// `strength` is a fraction of H, and H is the input's own relief -- so the
// same strength means different METRES on every tile, and the geometry it
// produces depends on the gully band it is paired with. Measured on a
// 2.3 km-relief tile with a 240 m coarsest gully:
//
// strength 0.22 (runevision's) -> depth:width 2.19
// strength 0.08 -> 0.80
// strength 0.02 -> 0.20
//
// Real valleys run about 0.05-0.25. runevision's 0.22 is right for HIS
// input amplitude of 0.125 and absurd for ours, which is why it had to be
// wound down by eye every time the tile changed.
//
// depth = 2 * strength * H, so strength = ratio * gully_metres / (2*H).
// Setting the ratio pins the GEOMETRY and lets strength follow.
float depth_ratio = 0.0f; // 0 = use --strength directly
};
// One evaluation. Kept separate from the loop so --check can call it directly.
static float erode_point(const Opts& o, float h_norm, float dhdx, float dhdz,
float px, float pz, float L, float H, float* ridge) {
// world slope -> filter slope: multiply by L/H, matching p = world/L and
// height = metres/H.
vec3 gn;
gn.x = h_norm;
gn.y = dhdx * (L / H);
gn.z = dhdz * (L / H);
// height_amp 1 and height_scale 1 leave the filter's own [0,1] remap as the
// only vertical transform, so the inverse below is exact.
vec2 en = dcerosion::erosion_node_eval_full(
vec2(px / L, pz / L), gn,
/* height_scale */ 1.0f,
o.strength,
o.cell_scale, o.octaves, o.gain,
/* fade_div */ 1.0f,
/* erosion_scale */ o.gully_metres / (L * o.cell_scale),
/* height_amp */ 1.0f,
o.lacunarity, o.gully_weight, o.detail, o.normalization,
vec2(o.height_offset_a, 0.0f),
vec4(0.1f, 0.0f, o.rounding_input, 2.0f),
vec4(1.25f, 1.25f, 2.8f, 1.5f),
vec2(0.7f, 1.0f),
// max_octaves: the FRACTIONAL octave limit erosion_node_eval_full
// gained after this tool was last built. The engine sets it per
// chunk from the voxel size; this tool has no voxel, so 999 = no
// taper, which is what the call did before the parameter existed.
999.0f);
if (ridge) *ridge = en.y;
// invert the remap: the filter returned n.x*0.5+0.5 plus its delta.
return (en.x - 0.5f) * 2.0f;
}
int main(int argc, char** argv) {
Opts o;
const char* in_path = nullptr;
const char* out_path = "/tmp/eroded.f32";
const char* ridge_path = nullptr;
bool check = false;
for (int i = 1; i < argc; i++) {
const char* a = argv[i];
auto next = [&]() { return argv[++i]; };
if (!strcmp(a, "--in") && i+1 < argc) in_path = next();
else if (!strcmp(a, "--out") && i+1 < argc) out_path = next();
else if (!strcmp(a, "--ridge")&& i+1 < argc) ridge_path = next();
else if (!strcmp(a, "--n") && i+1 < argc) o.n = atoi(next());
else if (!strcmp(a, "--mpp") && i+1 < argc) o.mpp = (float)atof(next());
else if (!strcmp(a, "--gully-metres") && i+1<argc) o.gully_metres = (float)atof(next());
else if (!strcmp(a, "--strength") && i+1<argc) o.strength = (float)atof(next());
else if (!strcmp(a, "--octaves") && i+1<argc) o.octaves = atoi(next());
else if (!strcmp(a, "--gain") && i+1<argc) o.gain = (float)atof(next());
else if (!strcmp(a, "--cell-scale")&& i+1<argc) o.cell_scale = (float)atof(next());
else if (!strcmp(a, "--rounding-input") && i+1<argc) o.rounding_input = (float)atof(next());
else if (!strcmp(a, "--vscale") && i+1<argc) o.vscale = (float)atof(next());
else if (!strcmp(a, "--offset") && i+1<argc) o.height_offset_a = (float)atof(next());
else if (!strcmp(a, "--depth-ratio") && i+1<argc) o.depth_ratio = (float)atof(next());
else if (!strcmp(a, "--check")) check = true;
else { fprintf(stderr, "unknown arg %s\n", a); return 2; }
}
// ---- the identity gate -------------------------------------------------
if (check) {
Opts z = o; z.strength = 0.0f; z.height_offset_a = 0.0f;
double worst = 0.0;
for (int k = 0; k < 200; k++) {
float hn = -1.0f + 2.0f * (float)k / 199.0f;
float g = erode_point(z, hn, 0.013f, -0.007f,
100.0f * k, 70.0f * k, 6000.0f, 500.0f, nullptr);
double d = fabs((double)g - (double)hn);
if (d > worst) worst = d;
}
printf("IDENTITY GATE — strength 0 and no height offset must be a no-op\n");
printf(" worst |out - in| over 200 heights : %.3e\n", worst);
bool ok = worst < 1e-5;
printf(" VERDICT: %s\n", ok ? "unit conversion is exact."
: "WRONG. The height mapping does not invert; "
"do not trust any eroded tile.");
if (!ok || !in_path) return ok ? 0 : 1;
}
if (!in_path) { fprintf(stderr, "need --in\n"); return 2; }
const int n = o.n;
std::vector<float> h((size_t)n * n);
FILE* f = fopen(in_path, "rb");
if (!f) { perror("open --in"); return 1; }
if (fread(h.data(), sizeof(float), h.size(), f) != h.size()) {
fprintf(stderr, "short read: expected %zu floats\n", h.size()); return 1;
}
fclose(f);
float lo = h[0], hi = h[0];
for (float v : h) { if (v < lo) lo = v; if (v > hi) hi = v; }
const float mid = 0.5f * (lo + hi);
const float H = (o.vscale > 0.0f) ? o.vscale
: fmaxf(0.5f * (hi - lo), 1.0f);
// L is fixed by the gully band we asked for, so the number the user gives
// is in metres and everything else follows.
const float L = o.gully_metres / o.cell_scale;
// Derive strength from the requested depth:width, so the same number gives
// the same geometry on any tile and at any gully band.
if (o.depth_ratio > 0.0f) {
o.strength = o.depth_ratio * o.gully_metres / (2.0f * H);
printf(" depth:width %.3f requested -> strength %.4f (H = %.0f m)\n",
o.depth_ratio, o.strength, H);
}
printf("erode heightmap: %dx%d @ %g m/px = %.2f km across\n",
n, n, o.mpp, n * o.mpp / 1000.0f);
printf(" input relief %.1f .. %.1f m (mid %.1f, H %.1f)\n", lo, hi, mid, H);
printf(" gullies %.0f m -> %.0f m (%d oct, lac %.2f)\n",
o.gully_metres, o.gully_metres / powf(o.lacunarity, (float)(o.octaves - 1)),
o.octaves, o.lacunarity);
printf(" strength %.3g rounding_input %.3g offset %.3g\n",
o.strength, o.rounding_input, o.height_offset_a);
std::vector<float> out(h.size()), ridge(h.size());
for (int j = 0; j < n; j++) {
for (int i = 0; i < n; i++) {
// central differences, clamped at the border
int im = (i > 0) ? i - 1 : i, ip = (i < n - 1) ? i + 1 : i;
int jm = (j > 0) ? j - 1 : j, jp = (j < n - 1) ? j + 1 : j;
float dhdx = (h[(size_t)j * n + ip] - h[(size_t)j * n + im]) /
(o.mpp * (float)(ip - im));
float dhdz = (h[(size_t)jp * n + i] - h[(size_t)jm * n + i]) /
(o.mpp * (float)(jp - jm));
float hn = (h[(size_t)j * n + i] - mid) / H;
float r = 0.0f;
float g = erode_point(o, hn, dhdx, dhdz,
(float)i * o.mpp, (float)j * o.mpp, L, H, &r);
out[(size_t)j * n + i] = mid + g * H;
ridge[(size_t)j * n + i] = r;
}
}
float olo = out[0], ohi = out[0];
for (float v : out) { if (v < olo) olo = v; if (v > ohi) ohi = v; }
printf(" output relief %.1f .. %.1f m\n", olo, ohi);
f = fopen(out_path, "wb");
if (!f) { perror("open --out"); return 1; }
fwrite(out.data(), sizeof(float), out.size(), f);
fclose(f);
printf(" wrote %s\n", out_path);
if (ridge_path) {
f = fopen(ridge_path, "wb");
if (f) { fwrite(ridge.data(), sizeof(float), ridge.size(), f); fclose(f);
printf(" wrote %s (ridge: -1 creases, +1 ridges)\n", ridge_path); }
}
return 0;
}