// =========================================================================
// HYPSO CURVE GATE — is the land/sea transfer curve monotonic, C1, and
// actually bimodal?
//
// The curve is the thing that turns a unimodal fBm into continents and
// ocean basins. Three ways it can be wrong and look fine in a screenshot:
//
// NOT MONOTONIC the terrain folds -- two input heights map to the same
// output, which dual contouring renders as a shelf that
// is not in the field.
// NOT C1 a slope discontinuity at the waterline, i.e. a crease
// ring around every coast. Invisible in a height plot.
// NOT BIMODAL the whole exercise fails silently: you get a water
// plane through unimodal noise, which is the thing this
// replaces.
//
// The derivative is checked against a central difference because the
// wrapper uses it for the chain rule on the gradient the erosion filter
// steers by -- an analytic derivative that is merely plausible would carve
// gullies aligned to a surface that is not there.
//
// Build (no engine, no godot-cpp, no scons, nothing into bin/):
// g++ -O2 -std=c++17 -I ../proc-nodes/src/shared \
// -o /tmp/hypso_gate tests/test_hypso_curve.cpp
// =========================================================================
#include <cstdio>
#include <cmath>
#include <vector>
#include <string>
#include <algorithm>
#include <cstdlib>
// The .inc is written in the GLSL/C++ common subset and names GLSL builtins,
// which exist only inside dcerosion -- erosion_cpp_shim.h supplies them there
// and pulls erosion_filter.inc in the same way (shim :128). Including this at
// global scope instead is what "'clamp' was not declared in this scope" means.
#include "erosion_cpp_shim.h"
namespace dcerosion {
#include "hypso_curve.inc"
}
using dcerosion::hypso_curve;
using dcerosion::vec2;
using dcerosion::vec3;
using dcerosion::fractal_multifractal_2d_lod;
// map_dc_erosion's intended defaults, in NORMALISED units.
static const float SEA = 0.00f;
static const float SHELF_W = 0.12f; // gentle shelf below sea
static const float SLOPE_W = 0.05f; // steep drop half-width
static const float ABYSS = -0.62f;
static const float PLAT_SLOPE = 0.13f; // continental platform: gentle
static const float MTN_ONSET = 0.35f;
static const float MTN_WIDTH = 0.10f;
static const float MTN_SLOPE = 1.30f; // mountains: steep
static const float PEAK_ONSET = 0.60f;
static const float PEAK_GAIN = 0.30f;
static vec2 C(float h, float pg = PEAK_GAIN) {
return hypso_curve(h, SEA, SHELF_W, SLOPE_W, ABYSS, PLAT_SLOPE,
MTN_ONSET, MTN_WIDTH, MTN_SLOPE, PEAK_ONSET, pg);
}
int main() {
// ---- 0. DOES IT COMPILE AS GLSL? ------------------------------------
//
// THE CHECK THIS GATE WAS MISSING. Everything below compiles the .inc as
// C++ and says nothing about the shader -- and the first version named a
// local `out`, which is an ordinary identifier in C++ and a RESERVED
// STORAGE QUALIFIER in GLSL. The C++ gate passed, the shader died with
// "unexpected OUT", the GPU silently fell back to the CPU mesher, and the
// only symptom a human saw was that terrain generation got slow.
//
// A dual-compiled file needs BOTH compilers or it is only half gated.
{
const char* src = "/tmp/hypso_gate_glsl.comp";
FILE* f = fopen(src, "w");
if (!f) { printf(" GLSL compile : SKIPPED (cannot write %s)\n", src); }
else {
fputs("#version 450\nlayout(local_size_x=1) in;\n", f);
FILE* g = fopen("../proc-nodes/src/shared/hypso_curve.inc", "r");
if (!g) { fclose(f); printf(" GLSL compile : SKIPPED (.inc not found)\n"); }
else {
int c; while ((c = fgetc(g)) != EOF) fputc(c, f);
fclose(g);
fputs("\nlayout(set=0,binding=0) buffer B { float v[]; };\n"
"void main() { v[0] = hypso_curve(v[1],v[2],v[3],v[4],v[5],"
"v[6],v[7],v[8],v[9],v[10],v[11]).x; }\n", f);
fclose(f);
int rc = system("glslc -fshader-stage=compute /tmp/hypso_gate_glsl.comp "
"-o /dev/null 2>/tmp/hypso_gate_glsl.err");
printf(" GLSL compile : %s\n", rc == 0 ? "ok" : "FAIL");
if (rc != 0) system("head -3 /tmp/hypso_gate_glsl.err");
}
}
}
printf("HYPSO CURVE GATE\n");
printf(" sea %.2f shelf_w %.2f slope_w %.2f abyss %.2f plat %.2f |"
" mtn %.2f/%.2f/%.2f | peak %.2f/%.2f\n\n",
SEA, SHELF_W, SLOPE_W, ABYSS, PLAT_SLOPE, MTN_ONSET, MTN_WIDTH,
MTN_SLOPE, PEAK_ONSET, PEAK_GAIN);
// ---- 1. MONOTONIC ---------------------------------------------------
// The real condition is dout/dh >= 0 everywhere -- that is monotonicity by
// construction. Comparing consecutive OUTPUTS in float32 also reports
// decreases of ~1e-5 that are accumulated rounding in the blend, not folds;
// a genuine fold is orders larger. So report both, and judge on the
// derivative while showing the float noise for what it is.
const int N = 20001;
float prev = C(-1.0f).x, worst_back = 0.0f; int backs = 0;
double min_d = 1e30; float min_at = 0.0f;
for (int i = 1; i < N; i++) {
float h = -1.0f + 2.0f * (float)i / (N - 1);
vec2 v = C(h);
if (v.x < prev) { backs++; worst_back = std::max(worst_back, prev - v.x); }
if (v.y < min_d) { min_d = v.y; min_at = h; }
prev = v.x;
}
printf(" min dout/dh : %+.5f at h=%+.3f %s\n", min_d, min_at,
min_d >= 0.0 ? "ok -- monotonic by construction" : "FAIL -- folds");
printf(" float32 backsteps: %d, worst %.2e (%s)\n", backs, worst_back,
worst_back < 1e-4 ? "rounding, not a fold" : "REAL FOLD");
// ---- 2. DERIVATIVE vs CENTRAL DIFFERENCE ----------------------------
double worst_d = 0.0; float worst_at = 0.0f;
const float eps = 1e-4f;
for (int i = 0; i <= 4000; i++) {
float h = -1.0f + 2.0f * (float)i / 4000.0f;
double fd = (C(h + eps).x - C(h - eps).x) / (2.0 * eps);
double an = C(h).y;
double e = fabs(fd - an);
if (e > worst_d) { worst_d = e; worst_at = h; }
}
printf(" d/dh vs finite : worst %.5f at h=%+.3f %s\n",
worst_d, worst_at, worst_d < 2e-2 ? "ok" : "FAIL");
// ---- 3. C1, BY SCALING -----------------------------------------------
// The previous version compared adjacent derivative samples against a
// fixed threshold, which measures |f''|*dh, not a discontinuity -- a
// smoothstep of width w has |f''| ~ 6/w^2, so it "failed" for being
// curved. A CONTINUOUS derivative's max adjacent step HALVES when the
// sampling halves; a jump does not move at all.
auto max_step = [&](int n) {
double worst = 0.0; float p = C(-1.0f).y;
for (int i = 1; i <= n; i++) {
float h = -1.0f + 2.0f * (float)i / (float)n;
float d = C(h).y;
worst = std::max(worst, (double)fabs(d - p));
p = d;
}
return worst;
};
double s1 = max_step(20000), s2 = max_step(40000);
double ratio = s2 > 0 ? s1 / s2 : 0.0;
printf(" C1 (scaling) : step %.5f -> %.5f, ratio %.2f %s\n",
s1, s2, ratio, (ratio > 1.7 && ratio < 2.3) ? "ok -- continuous"
: "FAIL -- jump");
// ---- 4. BIMODAL, against the null ------------------------------------
// Input: a triangular-ish distribution standing in for summed octaves.
// What matters is that the OUTPUT histogram has two modes where the
// INPUT has one, so the same input is measured through both.
const int S = 400000, B = 24;
std::vector<int> hin(B, 0), hout(B, 0);
unsigned s = 12345u;
auto rnd = [&]() { s = s * 1664525u + 1013904223u;
return (float)((s >> 8) & 0xFFFFFF) / (float)0xFFFFFF; };
float omin = C(-1.0f).x, omax = C(1.0f).x;
for (int i = 0; i < S; i++) {
float h = (rnd() + rnd() + rnd()) / 1.5f - 1.0f; // ~triangular on [-1,1]
float o = C(h).x;
int bi = (int)((h + 1.0f) * 0.5f * (B - 1) + 0.5f);
int bo = (int)((o - omin) / (omax - omin) * (B - 1) + 0.5f);
hin[std::min(std::max(bi,0),B-1)]++;
hout[std::min(std::max(bo,0),B-1)]++;
}
// Two modes with a real dip between them, not just "two bumps".
auto two_modes = [&](const std::vector<int>& hh, double& dip) {
int n = (int)hh.size(), i1 = 0;
for (int i = 0; i < n; i++) if (hh[i] > hh[i1]) i1 = i;
int i2 = -1;
for (int i = 0; i < n; i++) {
if (abs(i - i1) < 3) continue;
bool loc = (i == 0 || hh[i] >= hh[i-1]) && (i == n-1 || hh[i] >= hh[i+1]);
if (loc && (i2 < 0 || hh[i] > hh[i2])) i2 = i;
}
if (i2 < 0) { dip = 0; return 0; }
int a = std::min(i1, i2), b = std::max(i1, i2), tr = hh[a];
for (int i = a; i <= b; i++) tr = std::min(tr, hh[i]);
dip = 1.0 - (double)tr / (double)hh[i2]; // 1 = clean separation
return 2;
};
double dip_in = 0, dip_out = 0;
int m_in = two_modes(hin, dip_in), m_out = two_modes(hout, dip_out);
printf("\n output range : %.3f .. %.3f (waterline at %.4f)\n",
omin, omax, C(SEA).x);
printf(" input : %d mode(s), dip %.2f\n", m_in, dip_in);
printf(" output : %d mode(s), dip %.2f %s\n\n", m_out, dip_out,
(m_out == 2 && dip_out > 0.35) ? "ok -- BIMODAL"
: "FAIL -- not separated");
int pk = *std::max_element(hout.begin(), hout.end());
for (int i = 0; i < B; i++) {
float lo = omin + (omax - omin) * i / (float)(B - 1);
printf(" %+6.2f |%s\n", lo, std::string(hout[i] * 46 / pk, '#').c_str());
}
// ---- 5. PEAK LIFT TOUCHES ONLY THE TAIL ------------------------------
double worst_low = 0.0;
for (int i = 0; i <= 2000; i++) {
float h = -1.0f + (PEAK_ONSET + 1.0f) * (float)i / 2000.0f;
worst_low = std::max(worst_low, (double)fabs(C(h, PEAK_GAIN).x - C(h, 0.0f).x));
}
printf("\n peak_gain effect below onset: %.3e %s\n", worst_low,
worst_low < 1e-6 ? "ok -- lowlands untouched" : "FAIL");
printf(" peak_gain effect at h=1.0 : %.3f (asked for %.3f)\n",
C(1.0f, PEAK_GAIN).x - C(1.0f, 0.0f).x, PEAK_GAIN);
// ---- 6. LAND FRACTION, against the REAL mountain field ---------------
//
// The histogram above uses a synthetic input, which is fine for shape but
// useless for "how much of the world is land" -- that depends on the
// actual distribution the curve is fed. erosion_filter.inc carries
// fractal_multifractal_2d_lod, the very function map_dc_erosion's wrapper
// calls, so this samples THAT rather than a stand-in.
//
// map_dc_erosion.gd's constants (MOUNTAIN_SCALE is post-0.1).
const float TH = 60000.0f, MSCALE = 0.1f, LAC = 1.87f, MGAIN = 0.68f;
const int MOCT = 8;
const float VARIANCE = 0.0f; // generator default
const float HALF = 458800.0f; // depth-16 world, +/- 458.8 km
const int G = 420; // 176400 columns
printf("\n --- LAND FRACTION over the real mountain field ---\n");
printf(" terrain_height %.0f mountain_scale %.2f %d oct lac %.2f"
" gain %.2f\n", TH, MSCALE, MOCT, LAC, MGAIN);
std::vector<float> raw;
raw.reserve((size_t)G * G);
for (int j = 0; j < G; j++) {
for (int i = 0; i < G; i++) {
float x = -HALF + 2.0f * HALF * i / (G - 1);
float z = -HALF + 2.0f * HALF * j / (G - 1);
vec3 m = fractal_multifractal_2d_lod(
vec2(x / TH * MSCALE, z / TH * MSCALE), 1.0f,
MOCT, LAC, MGAIN, VARIANCE, 999.0f);
raw.push_back(m.x);
}
}
float rmin = raw[0], rmax = raw[0];
for (float v : raw) { rmin = std::min(rmin, v); rmax = std::max(rmax, v); }
printf(" raw multifractal range: %+.3f .. %+.3f (curve assumes ~[-1,1])\n",
rmin, rmax);
auto land_frac = [&](float sea_v) {
float wl = hypso_curve(sea_v, sea_v, SHELF_W, SLOPE_W, ABYSS, PLAT_SLOPE,
MTN_ONSET, MTN_WIDTH, MTN_SLOPE,
PEAK_ONSET, PEAK_GAIN).x;
size_t n = 0;
for (float v : raw) {
float o = hypso_curve(v, sea_v, SHELF_W, SLOPE_W, ABYSS, PLAT_SLOPE,
MTN_ONSET, MTN_WIDTH, MTN_SLOPE,
PEAK_ONSET, PEAK_GAIN).x;
if (o > wl) n++;
}
return 100.0 * (double)n / (double)raw.size();
};
printf("\n sea land %% (Earth is ~29%%)\n");
for (int k = -6; k <= 6; k++) {
float sv = 0.05f * k;
printf(" %+5.2f %5.1f%s\n", sv, land_frac(sv),
(k == 0) ? " <- current" : "");
}
printf("\n at sea=%.2f: land %.1f%%, mountains (h>%.2f) %.1f%%\n",
SEA, land_frac(SEA), MTN_ONSET,
100.0 * (double)std::count_if(raw.begin(), raw.end(),
[&](float v){ return v > MTN_ONSET; }) / raw.size());
// ---- 7. THE SHELF, which is the whole point of this revision --------
printf("\n --- COAST PROFILE ---\n");
float wl = C(SEA).x;
printf(" waterline (curve units) : %+.4f\n", wl);
float brk = C(SEA - SHELF_W + SLOPE_W).x;
printf(" depth at the shelf break : %+.4f (%.0f m at ~7790 m/unit)\n",
brk - wl, (brk - wl) * 7790.0f);
// gradient just above sea vs at the steepest point
float g_land = C(SEA + 0.02f).y;
float g_max = 0.0f; float g_at = 0.0f;
for (int i = 0; i <= 4000; i++) {
float h = -1.0f + 2.0f * (float)i / 4000.0f;
if (C(h).y > g_max) { g_max = C(h).y; g_at = h; }
}
printf(" gradient at the shore : %.3f\n", g_land);
printf(" gradient at the steepest pt : %.3f at h=%+.3f (%s sea)\n",
g_max, g_at, g_at < SEA ? "BELOW -- submerged, good" : "ABOVE -- a cliff");
printf(" steepest / shore ratio : %.1fx %s\n", g_max / g_land,
g_at < SEA ? "(the drop is under water)" : "(THE DROP IS THE COAST)");
return 0;
}