OverviewFilesCommitsIssuesDocsReleasesPeers

dc-terrain commits / bd097c5

dc-terrain: wire the hypso curve, mask detail by slope, export the slope
GENERATOR. The wrapper feeds the mountain layer through hypso_curve and carries its derivative through the chain rule, so the erosion filter is handed a coastline-shaped surface rather than one it has to invent. The curve is driven by mountain_height, NOT relief_amp: relief varies 12.5x spatially, so keying sea level off it would move the waterline per region. DETAIL IS MASKED BY SLOPE. It was applied flat, so the same 100 m mottling sat on flats, cliff faces, ridges and gully floors alike -- and because the material reads the NORMAL, that mottling propagated into the colour, which is why the terrain read as crumpled foil rather than as landforms. Physically it is soil and scree: it accumulates on low slopes and is stripped from steep faces. The gradient it masks against is free. pn_out_grad is converted to d(world y)/d(world x) HERE, where inv_scale and the y-scale are already literals, rather than by a consumer copying 1/60000 and 4200 into a shader by hand -- the multi-stage rescale this file has already got wrong once. In sqrt height space the inverse is not derived and it is zeroed rather than scaled by a constant that would be quietly wrong: a consumer sees FLAT, which is visible. GATES. test_hypso_curve checks the curve compiles as BOTH languages -- it is dual-compiled, and naming a local `out` compiled fine as C++ while killing the shader, which showed up only as "slower to generate" when the GPU fell back to the CPU mesher. test_analytic_gradient isolates the analytic derivative per stage and names the culprit: grad_noised_2d 0.019 deg, fractal_noise_2d_lod 0.019, fractal_multifractal_2d_lod 2.15, erosion_node_eval_grad 29.2. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
bd097c5686fdac06b007b33ebc11b51aa5207c98
by plankatron <fleeky@gmail.com> · 2026-09-07T17:28:11.000Z

Changed files (5) +812 -6

src/dc_terrain_generator_procnodes.cpp +265 -6
……@@ -15,6 +15,13 @@
1515// dependency on proc-nodes' classes -- the ClassDB reflection in the header
1616// comment still stands. It is a shared math header and nothing more.
1717#include "erosion_cpp_shim.h"
18// The curve as C++, inside the shim's namespace so it sees the GLSL builtins
19// -- the same arrangement the shim uses for erosion_filter.inc. The shader
20// compiles this identical text via HYPSO_CURVE_GLSL, so get_waterline_norm()
21// below cannot drift from what the GPU computes.
22namespace dcerosion {
23#include "hypso_curve.inc"
24}
1825#include "erosion_wrapper_glsl.h"
1926#include <godot_cpp/classes/file_access.hpp>
2027
……@@ -91,6 +98,46 @@ void DCTerrainGeneratorProcNodes::_bind_methods() {
9198 ClassDB::bind_method(D_METHOD("set_crease_rounding", "value"), &DCTerrainGeneratorProcNodes::set_crease_rounding);
9299 ClassDB::bind_method(D_METHOD("get_crease_rounding"), &DCTerrainGeneratorProcNodes::get_crease_rounding);
93100 ADD_PROPERTY(PropertyInfo(Variant::FLOAT, "crease_rounding"), "set_crease_rounding", "get_crease_rounding");
101 ClassDB::bind_method(D_METHOD("set_detail_slope_lo", "value"), &DCTerrainGeneratorProcNodes::set_detail_slope_lo);
102 ClassDB::bind_method(D_METHOD("get_detail_slope_lo"), &DCTerrainGeneratorProcNodes::get_detail_slope_lo);
103 ADD_PROPERTY(PropertyInfo(Variant::FLOAT, "detail_slope_lo"), "set_detail_slope_lo", "get_detail_slope_lo");
104 ClassDB::bind_method(D_METHOD("set_detail_slope_hi", "value"), &DCTerrainGeneratorProcNodes::set_detail_slope_hi);
105 ClassDB::bind_method(D_METHOD("get_detail_slope_hi"), &DCTerrainGeneratorProcNodes::get_detail_slope_hi);
106 ADD_PROPERTY(PropertyInfo(Variant::FLOAT, "detail_slope_hi"), "set_detail_slope_hi", "get_detail_slope_hi");
107 ClassDB::bind_method(D_METHOD("set_hypso_enabled", "value"), &DCTerrainGeneratorProcNodes::set_hypso_enabled);
108 ClassDB::bind_method(D_METHOD("get_hypso_enabled"), &DCTerrainGeneratorProcNodes::get_hypso_enabled);
109 ADD_PROPERTY(PropertyInfo(Variant::FLOAT, "hypso_enabled"), "set_hypso_enabled", "get_hypso_enabled");
110 ClassDB::bind_method(D_METHOD("set_hypso_sea", "value"), &DCTerrainGeneratorProcNodes::set_hypso_sea);
111 ClassDB::bind_method(D_METHOD("get_hypso_sea"), &DCTerrainGeneratorProcNodes::get_hypso_sea);
112 ADD_PROPERTY(PropertyInfo(Variant::FLOAT, "hypso_sea"), "set_hypso_sea", "get_hypso_sea");
113 ClassDB::bind_method(D_METHOD("set_hypso_shelf_width", "value"), &DCTerrainGeneratorProcNodes::set_hypso_shelf_width);
114 ClassDB::bind_method(D_METHOD("get_hypso_shelf_width"), &DCTerrainGeneratorProcNodes::get_hypso_shelf_width);
115 ADD_PROPERTY(PropertyInfo(Variant::FLOAT, "hypso_shelf_width"), "set_hypso_shelf_width", "get_hypso_shelf_width");
116 ClassDB::bind_method(D_METHOD("set_hypso_slope_width", "value"), &DCTerrainGeneratorProcNodes::set_hypso_slope_width);
117 ClassDB::bind_method(D_METHOD("get_hypso_slope_width"), &DCTerrainGeneratorProcNodes::get_hypso_slope_width);
118 ADD_PROPERTY(PropertyInfo(Variant::FLOAT, "hypso_slope_width"), "set_hypso_slope_width", "get_hypso_slope_width");
119 ClassDB::bind_method(D_METHOD("set_hypso_abyss", "value"), &DCTerrainGeneratorProcNodes::set_hypso_abyss);
120 ClassDB::bind_method(D_METHOD("get_hypso_abyss"), &DCTerrainGeneratorProcNodes::get_hypso_abyss);
121 ADD_PROPERTY(PropertyInfo(Variant::FLOAT, "hypso_abyss"), "set_hypso_abyss", "get_hypso_abyss");
122 ClassDB::bind_method(D_METHOD("set_hypso_plat_slope", "value"), &DCTerrainGeneratorProcNodes::set_hypso_plat_slope);
123 ClassDB::bind_method(D_METHOD("get_hypso_plat_slope"), &DCTerrainGeneratorProcNodes::get_hypso_plat_slope);
124 ADD_PROPERTY(PropertyInfo(Variant::FLOAT, "hypso_plat_slope"), "set_hypso_plat_slope", "get_hypso_plat_slope");
125 ClassDB::bind_method(D_METHOD("set_hypso_mtn_onset", "value"), &DCTerrainGeneratorProcNodes::set_hypso_mtn_onset);
126 ClassDB::bind_method(D_METHOD("get_hypso_mtn_onset"), &DCTerrainGeneratorProcNodes::get_hypso_mtn_onset);
127 ADD_PROPERTY(PropertyInfo(Variant::FLOAT, "hypso_mtn_onset"), "set_hypso_mtn_onset", "get_hypso_mtn_onset");
128 ClassDB::bind_method(D_METHOD("set_hypso_mtn_width", "value"), &DCTerrainGeneratorProcNodes::set_hypso_mtn_width);
129 ClassDB::bind_method(D_METHOD("get_hypso_mtn_width"), &DCTerrainGeneratorProcNodes::get_hypso_mtn_width);
130 ADD_PROPERTY(PropertyInfo(Variant::FLOAT, "hypso_mtn_width"), "set_hypso_mtn_width", "get_hypso_mtn_width");
131 ClassDB::bind_method(D_METHOD("set_hypso_mtn_slope", "value"), &DCTerrainGeneratorProcNodes::set_hypso_mtn_slope);
132 ClassDB::bind_method(D_METHOD("get_hypso_mtn_slope"), &DCTerrainGeneratorProcNodes::get_hypso_mtn_slope);
133 ADD_PROPERTY(PropertyInfo(Variant::FLOAT, "hypso_mtn_slope"), "set_hypso_mtn_slope", "get_hypso_mtn_slope");
134 ClassDB::bind_method(D_METHOD("set_hypso_peak_onset", "value"), &DCTerrainGeneratorProcNodes::set_hypso_peak_onset);
135 ClassDB::bind_method(D_METHOD("get_hypso_peak_onset"), &DCTerrainGeneratorProcNodes::get_hypso_peak_onset);
136 ADD_PROPERTY(PropertyInfo(Variant::FLOAT, "hypso_peak_onset"), "set_hypso_peak_onset", "get_hypso_peak_onset");
137 ClassDB::bind_method(D_METHOD("set_hypso_peak_gain", "value"), &DCTerrainGeneratorProcNodes::set_hypso_peak_gain);
138 ClassDB::bind_method(D_METHOD("get_hypso_peak_gain"), &DCTerrainGeneratorProcNodes::get_hypso_peak_gain);
139 ADD_PROPERTY(PropertyInfo(Variant::FLOAT, "hypso_peak_gain"), "set_hypso_peak_gain", "get_hypso_peak_gain");
140 ClassDB::bind_method(D_METHOD("get_waterline_norm"), &DCTerrainGeneratorProcNodes::get_waterline_norm);
94141 ClassDB::bind_method(D_METHOD("set_ridge_rounding", "value"), &DCTerrainGeneratorProcNodes::set_ridge_rounding);
95142 ClassDB::bind_method(D_METHOD("get_ridge_rounding"), &DCTerrainGeneratorProcNodes::get_ridge_rounding);
96143 ADD_PROPERTY(PropertyInfo(Variant::FLOAT, "ridge_rounding"), "set_ridge_rounding", "get_ridge_rounding");
……@@ -316,6 +363,107 @@ void DCTerrainGeneratorProcNodes::set_crease_rounding(float p_val) {
316363 emit_changed();
317364}
318365
366float DCTerrainGeneratorProcNodes::get_waterline_norm() const {
367 // Sea level in the curve's OUTPUT units. The curve is monotonic, so the
368 // waterline is simply its value at `sea` -- no bisection, no second
369 // constant to fall out of step with the first.
370 return dcerosion::hypso_curve(hypso_sea_, hypso_sea_, hypso_shelf_width_,
371 hypso_slope_width_, hypso_abyss_, hypso_plat_slope_, hypso_mtn_onset_,
372 hypso_mtn_width_, hypso_mtn_slope_, hypso_peak_onset_,
373 hypso_peak_gain_).x;
374}
375
376void DCTerrainGeneratorProcNodes::set_detail_slope_lo(float p_val) {
377 if (detail_slope_lo_ == p_val) return;
378 detail_slope_lo_ = p_val;
379 _recompile_glsl();
380 emit_changed();
381}
382
383void DCTerrainGeneratorProcNodes::set_detail_slope_hi(float p_val) {
384 if (detail_slope_hi_ == p_val) return;
385 detail_slope_hi_ = p_val;
386 _recompile_glsl();
387 emit_changed();
388}
389
390void DCTerrainGeneratorProcNodes::set_hypso_enabled(float p_val) {
391 if (hypso_enabled_ == p_val) return;
392 hypso_enabled_ = p_val;
393 _recompile_glsl();
394 emit_changed();
395}
396
397void DCTerrainGeneratorProcNodes::set_hypso_sea(float p_val) {
398 if (hypso_sea_ == p_val) return;
399 hypso_sea_ = p_val;
400 _recompile_glsl();
401 emit_changed();
402}
403
404void DCTerrainGeneratorProcNodes::set_hypso_shelf_width(float p_val) {
405 if (hypso_shelf_width_ == p_val) return;
406 hypso_shelf_width_ = p_val;
407 _recompile_glsl();
408 emit_changed();
409}
410
411void DCTerrainGeneratorProcNodes::set_hypso_slope_width(float p_val) {
412 if (hypso_slope_width_ == p_val) return;
413 hypso_slope_width_ = p_val;
414 _recompile_glsl();
415 emit_changed();
416}
417
418void DCTerrainGeneratorProcNodes::set_hypso_abyss(float p_val) {
419 if (hypso_abyss_ == p_val) return;
420 hypso_abyss_ = p_val;
421 _recompile_glsl();
422 emit_changed();
423}
424
425void DCTerrainGeneratorProcNodes::set_hypso_plat_slope(float p_val) {
426 if (hypso_plat_slope_ == p_val) return;
427 hypso_plat_slope_ = p_val;
428 _recompile_glsl();
429 emit_changed();
430}
431
432void DCTerrainGeneratorProcNodes::set_hypso_mtn_onset(float p_val) {
433 if (hypso_mtn_onset_ == p_val) return;
434 hypso_mtn_onset_ = p_val;
435 _recompile_glsl();
436 emit_changed();
437}
438
439void DCTerrainGeneratorProcNodes::set_hypso_mtn_width(float p_val) {
440 if (hypso_mtn_width_ == p_val) return;
441 hypso_mtn_width_ = p_val;
442 _recompile_glsl();
443 emit_changed();
444}
445
446void DCTerrainGeneratorProcNodes::set_hypso_mtn_slope(float p_val) {
447 if (hypso_mtn_slope_ == p_val) return;
448 hypso_mtn_slope_ = p_val;
449 _recompile_glsl();
450 emit_changed();
451}
452
453void DCTerrainGeneratorProcNodes::set_hypso_peak_onset(float p_val) {
454 if (hypso_peak_onset_ == p_val) return;
455 hypso_peak_onset_ = p_val;
456 _recompile_glsl();
457 emit_changed();
458}
459
460void DCTerrainGeneratorProcNodes::set_hypso_peak_gain(float p_val) {
461 if (hypso_peak_gain_ == p_val) return;
462 hypso_peak_gain_ = p_val;
463 _recompile_glsl();
464 emit_changed();
465}
466
319467void DCTerrainGeneratorProcNodes::set_ridge_rounding(float p_val) {
320468 if (ridge_rounding_ == p_val) return;
321469 ridge_rounding_ = p_val;
……@@ -1120,6 +1268,49 @@ void DCTerrainGeneratorProcNodes::_recompile_glsl() {
11201268 const bool want_mountain = has_noise_lib && mountain_height_ != 0.0f;
11211269 const bool want_detail = has_noise_lib && detail_scale_ > 0.0f;
11221270
1271 // THE MOUNTAIN ASSIGNMENT, with the land/sea curve folded in.
1272 //
1273 // The curve goes HERE -- on the mountain layer, BEFORE the erosion filter
1274 // -- because that is the order the landscape happens in: continents and
1275 // ocean basins are the shape, erosion is what runs over it. Applying it
1276 // after erosion would stretch gullies wherever the curve is steep.
1277 //
1278 // ⚠ THE GRADIENT MUST BE REMAPPED TOO. The filter is slope-aware and
1279 // steers gullies by the gradient the wrapper hands it, so a height that
1280 // has been remapped and a gradient that has not describes a surface that
1281 // does not exist. hypso_curve returns dout/dh for exactly this.
1282 //
1283 // Disabled emits the ORIGINAL line, byte for byte, so `hypso_enabled 0`
1284 // is not "the curve with neutral parameters" -- it is the previous shader.
1285 String mtn_assign = (hypso_enabled_ != 0.0f)
1286 ? String(
1287 " vec2 hc = hypso_curve(mtn.x, ")
1288 + String::num(hypso_sea_, 6) + ", " + String::num(hypso_shelf_width_, 6) + ", "
1289 + String::num(hypso_slope_width_, 6) + ", "
1290 + String::num(hypso_abyss_, 6) + ", " + String::num(hypso_plat_slope_, 6) + ",\n"
1291 " "
1292 + String::num(hypso_mtn_onset_, 6) + ", " + String::num(hypso_mtn_width_, 6) + ", "
1293 + String::num(hypso_mtn_slope_, 6) + ", " + String::num(hypso_peak_onset_, 6) + ", "
1294 + String::num(hypso_peak_gain_, 6) + ");\n"
1295 // ⚠ mtn_amp IS NOT USED HERE, and that is deliberate.
1296 //
1297 // mtn_amp = mountain_height * relief_amp, and relief_amp varies
1298 // SPATIALLY from relief_floor (0.08) to 1.0 -- a 12.5x swing.
1299 // Scaling the curve's output by it would put sea level at a
1300 // different world Y in every region, and a single water plane
1301 // then cannot match the shoreline anywhere. The curve defines
1302 // where the sea is; that cannot be a function of position.
1303 //
1304 // So the hypso path uses mountain_height ALONE. The cost is that
1305 // relief_scale/relief_floor stop modulating how mountainous a
1306 // region is while the curve is on. Re-coupling them belongs on
1307 // the curve's INPUT (or on mtn_slope/peak_gain), not on its
1308 // output, and is deliberately not done here.
1309 " pn_terrain_offset = vec3(hc.x, mtn.y * " + mscale_str + " * hc.y, mtn.z * "
1310 + mscale_str + " * hc.y) * " + mheight_str + " * 8.0;\n"
1311 : String(" pn_terrain_offset = vec3(mtn.x, mtn.y * ") + mscale_str
1312 + ", mtn.z * " + mscale_str + ") * mtn_amp * 8.0;\n";
1313
11231314 String adapted =
11241315 "// Octave limit — fixed high value for LOD-consistent surface position.\n"
11251316 "// Nyquist clamping shifts the erosion surface between LOD levels, causing\n"
……@@ -1179,7 +1370,7 @@ void DCTerrainGeneratorProcNodes::_recompile_glsl() {
11791370 " // zero erosion strength.\n"
11801371 + (hblk_n_ > 0
11811372 ? String(" pn_terrain_offset = heightmap_gn(x, z);\n")
1182 : String(" pn_terrain_offset = vec3(mtn.x, mtn.y * ") + mscale_str + ", mtn.z * " + mscale_str + ") * mtn_amp * 8.0;\n")
1373 : mtn_assign)
11831374 +
11841375 " // DETAIL rides ON TOP of the eroded surface and is deliberately NOT folded\n"
11851376 " // into pn_terrain_offset: its gradient carries a factor of detail_scale, so\n"
……@@ -1221,7 +1412,25 @@ void DCTerrainGeneratorProcNodes::_recompile_glsl() {
12211412 + (want_detail ? String(
12221413 " if (" + dscale_str + " > 0.0) {\n"
12231414 " vec3 det = fractal_noise_2d_lod(vec2(x, z) * inv_scale * " + dscale_str + ", 1.0, " + doct_str + ", " + lac_str + ", " + dgain_str + ", det_lim);\n"
1224 " det_offset = det.x * " + dfrac_str + ";\n"
1415 // DETAIL BELONGS ON GENTLE GROUND, NOT EVERYWHERE.
1416 //
1417 // It was applied flat, so the same 100 m mottling sat on flats,
1418 // cliff faces, ridges and gully floors alike -- and because the
1419 // material reads the NORMAL, that mottling propagated straight
1420 // into the colour, which is why the terrain read as crumpled foil
1421 // rather than as landforms.
1422 //
1423 // Physically it is soil, scree and regolith: it ACCUMULATES on low
1424 // slopes and is stripped from steep faces. Masking by slope is the
1425 // same statement.
1426 //
1427 // The gradient is free -- pn_terrain_offset.yz is the mountain
1428 // slope computed a few lines above, in gn units per unit p. No
1429 // extra evaluation, no new field.
1430 " float det_slope = length(vec2(pn_terrain_offset.y, pn_terrain_offset.z));\n"
1431 " float det_mask = 1.0 - smoothstep(" + String::num(detail_slope_lo_, 6)
1432 + ", " + String::num(detail_slope_hi_, 6) + ", det_slope);\n"
1433 " det_offset = det.x * " + dfrac_str + " * det_mask;\n"
12251434 " }\n"
12261435 ) : String())
12271436 + (height_space_sqrt_
……@@ -1234,8 +1443,30 @@ void DCTerrainGeneratorProcNodes::_recompile_glsl() {
12341443 + ""
12351444 " float d = proc_density(vec3(x * inv_scale, y_mapped, z * inv_scale));\n"
12361445 " d *= " + inv_yscale_str + ";\n"
1237 " // Clip: solid below floor_y, air above ceiling_y.\n"
1238 " d = max(d, (" + floor_str + " - y) * 2.0);\n"
1446 // pn_out_grad LEAVES proc_density IN NODE SPACE: d(mapped height)/d(p),
1447 // with p = world / terrain_height. A consumer that wants a surface
1448 // normal needs d(world y)/d(world x), and the two constants that close
1449 // that gap are emitted right here as literals. Converting in the
1450 // consumer would mean copying 1/60000 and 4200 into a shader by hand --
1451 // exactly the multi-stage rescale this file has already got wrong once
1452 // (the waterline, off by 3.4x). Do it where the constants live.
1453 //
1454 // ⚠ det_offset's slope is NOT in here. Detail is subtracted from
1455 // y_mapped and its gradient is deliberately not folded into
1456 // pn_terrain_offset (see the DETAIL comment above: the direction flips
1457 // several times inside one gully cell). The normal this yields is the
1458 // ERODED MOUNTAIN's, which is the structure that survives to distance.
1459 + (height_space_sqrt_
1460 ? String(
1461 " // sqrt height space: the mapped->world conversion is not linear and\n"
1462 " // the inverse is not derived here. Zeroed rather than scaled by a\n"
1463 " // constant that would be quietly wrong -- a consumer sees FLAT, which\n"
1464 " // is visible, instead of a normal tilted by an unaudited factor.\n"
1465 " pn_out_grad = vec2(0.0);\n")
1466 : String(
1467 " pn_out_grad *= inv_scale * ") + inv_yscale_str + ";\n")
1468 + String(" // Clip: solid below floor_y, air above ceiling_y.\n"
1469 " d = max(d, (") + floor_str + " - y) * 2.0);\n"
12391470 " d = min(d, (" + ceil_str + " - y) * 2.0);\n"
12401471 " return d;\n"
12411472 "}\n";
……@@ -1459,7 +1690,25 @@ float DCTerrainGeneratorProcNodes::_eval_at(const EvalParams &p,
14591690 float mtn_amp = mtn_height * relief_amp;
14601691 // The mountain's HEIGHT AND GRADIENT, in gn units, to hand to the erosion —
14611692 // the pn_terrain_offset the GLSL sets.
1462 vec3 terrain_offset = vec3(mtn.x, mtn.y * mtn_scale, mtn.z * mtn_scale) * (mtn_amp * 8.0f);
1693 // THE LAND/SEA CURVE, mirroring the GLSL emission exactly. Without this
1694 // the CPU answers from the un-remapped fBm while the GPU meshes the
1695 // remapped one -- the meshed surface and the collidable surface differ,
1696 // which is the failure this function's own note below was written about.
1697 //
1698 // mountain_height alone, NOT mtn_amp: relief_amp varies spatially and
1699 // scaling the curve by it would move sea level per region. Same reason,
1700 // same line, as the emission site.
1701 vec3 terrain_offset;
1702 if (hypso_enabled_ != 0.0f) {
1703 dcerosion::vec2 hc = dcerosion::hypso_curve(mtn.x,
1704 hypso_sea_, hypso_shelf_width_, hypso_slope_width_, hypso_abyss_, hypso_plat_slope_,
1705 hypso_mtn_onset_, hypso_mtn_width_, hypso_mtn_slope_,
1706 hypso_peak_onset_, hypso_peak_gain_);
1707 terrain_offset = vec3(hc.x, mtn.y * mtn_scale * hc.y,
1708 mtn.z * mtn_scale * hc.y) * (mtn_height * 8.0f);
1709 } else {
1710 terrain_offset = vec3(mtn.x, mtn.y * mtn_scale, mtn.z * mtn_scale) * (mtn_amp * 8.0f);
1711 }
14631712
14641713 // THE SAME SUBSTITUTION THE GLSL MAKES.
14651714 //
……@@ -1491,7 +1740,17 @@ float DCTerrainGeneratorProcNodes::_eval_at(const EvalParams &p,
14911740 float detail_frac = (y_range > 0.01f) ? (detail_height_ / y_range) : 0.0f;
14921741 vec3 det = dcerosion::fractal_noise_2d_lod(vec2(x, z) * inv_scale * detail_scale_,
14931742 1.0f, detail_octaves_, lac, detail_gain_, det_lim);
1494 det_offset = det.x * detail_frac;
1743 // MIRRORS THE GLSL. See the note at the emission site: detail is
1744 // masked by slope so it accumulates on gentle ground and is stripped
1745 // from steep faces. terrain_offset is pn_terrain_offset's twin.
1746 float det_slope = std::sqrt(terrain_offset.y * terrain_offset.y +
1747 terrain_offset.z * terrain_offset.z);
1748 float dm_t = (detail_slope_hi_ > detail_slope_lo_)
1749 ? (det_slope - detail_slope_lo_) / (detail_slope_hi_ - detail_slope_lo_)
1750 : 0.0f;
1751 dm_t = dm_t < 0.0f ? 0.0f : (dm_t > 1.0f ? 1.0f : dm_t);
1752 float det_mask = 1.0f - (dm_t * dm_t * (3.0f - 2.0f * dm_t));
1753 det_offset = det.x * detail_frac * det_mask;
14951754 }
14961755
14971756 // Map world Y to proc_density space
src/dc_terrain_generator_procnodes.h +49 -0
……@@ -231,6 +231,25 @@ class DCTerrainGeneratorProcNodes : public DCTerrainGenerator {
231231 float height_offset_follows_altitude_ = 0.0f;
232232 float crease_rounding_ = 0.1f;
233233 float ridge_rounding_ = 0.0f;
234
235 // Detail-layer slope mask: full detail below lo, none above hi.
236 // Units are gn-per-p, the same as pn_terrain_offset.yz.
237 float detail_slope_lo_ = 0.30f;
238 float detail_slope_hi_ = 1.20f;
239
240 // --- hypsometric land/sea curve (hypso_curve.inc) ---
241 float hypso_enabled_ = 0.0f; // 0 = off, and off is BIT-IDENTICAL to before this landed
242 float hypso_sea_ = 0.20f; // input at the middle of the continental slope; sets land fraction (0.20 -> 33.7%)
243 float hypso_shelf_width_ = 0.12f; // gentle shelf below sea
244 float hypso_slope_width_ = 0.05f; // steep drop half-width // half-width of the slope; smaller = sharper coast
245 float hypso_abyss_ = -0.62f; // OUTPUT value of the abyssal plain
246 float hypso_plat_slope_ = 0.13f; // continental platform gradient (gentle)
247 float hypso_mtn_onset_ = 0.35f; // where the mountain front begins
248 float hypso_mtn_width_ = 0.10f; // half-width of that front
249 float hypso_mtn_slope_ = 1.30f; // mountain gradient (steep)
250 float hypso_peak_onset_ = 0.60f; // where the summit lift begins
251 float hypso_peak_gain_ = 0.30f; // output added at h = 1
252
234253 float input_feature_rounding_ = 0.1f;
235254 float rounding_per_octave_ = 2.0f;
236255 float gully_onset_slope_ = 1.25f;
……@@ -408,6 +427,36 @@ public:
408427 float get_height_offset_follows_altitude() const { return height_offset_follows_altitude_; }
409428 void set_crease_rounding(float p_val);
410429 float get_crease_rounding() const { return crease_rounding_; }
430
431 void set_detail_slope_lo(float p_val);
432 float get_detail_slope_lo() const { return detail_slope_lo_; }
433 void set_detail_slope_hi(float p_val);
434 float get_detail_slope_hi() const { return detail_slope_hi_; }
435 void set_hypso_enabled(float p_val);
436 float get_hypso_enabled() const { return hypso_enabled_; }
437 void set_hypso_sea(float p_val);
438 float get_hypso_sea() const { return hypso_sea_; }
439 void set_hypso_shelf_width(float p_val);
440 float get_hypso_shelf_width() const { return hypso_shelf_width_; }
441 void set_hypso_slope_width(float p_val);
442 float get_hypso_slope_width() const { return hypso_slope_width_; }
443 void set_hypso_abyss(float p_val);
444 float get_hypso_abyss() const { return hypso_abyss_; }
445 void set_hypso_plat_slope(float p_val);
446 float get_hypso_plat_slope() const { return hypso_plat_slope_; }
447 void set_hypso_mtn_onset(float p_val);
448 float get_hypso_mtn_onset() const { return hypso_mtn_onset_; }
449 void set_hypso_mtn_width(float p_val);
450 float get_hypso_mtn_width() const { return hypso_mtn_width_; }
451 void set_hypso_mtn_slope(float p_val);
452 float get_hypso_mtn_slope() const { return hypso_mtn_slope_; }
453 void set_hypso_peak_onset(float p_val);
454 float get_hypso_peak_onset() const { return hypso_peak_onset_; }
455 void set_hypso_peak_gain(float p_val);
456 float get_hypso_peak_gain() const { return hypso_peak_gain_; }
457 // Sea level in OUTPUT units -- hypso_curve(sea). The curve is
458 // monotonic so this is exact; do NOT also hand-set a water height.
459 float get_waterline_norm() const;
411460 void set_ridge_rounding(float p_val);
412461 float get_ridge_rounding() const { return ridge_rounding_; }
413462 void set_input_feature_rounding(float p_val);
tests/test_analytic_gradient.cpp +193 -0
……@@ -0,0 +1,193 @@
1// =========================================================================
2// ANALYTIC GRADIENT GATE — which stage loses the derivative?
3//
4// pn_out_grad carries erosion_node_eval_grad's analytic slope to the material
5// so a detail normal costs ONE density evaluation instead of six (measured on
6// an RTX 2070 at 1080p: 1.57 ms vs 16.4 ms fullscreen). Measured against
7// central differences it agreed within 5 deg for only 22.7% of samples --
8// with no transposition, no sign slip, and a magnitude ratio scattered about
9// 1.0 rather than offset from it. So a stage is dropping a term.
10//
11// Each stage is checked against a central difference of its OWN output, in
12// isolation, so the answer is a stage and not a vibe:
13//
14// grad_noised_2d analytic by construction -- if THIS is wrong
15// nothing downstream can be right
16// fractal_noise_2d_lod a sum of exact gradients; should be exact
17// fractal_multifractal_2d_lod weights each octave by accumulated elevation
18// and DROPS d(weight)/dp -- its own comment says
19// so. Prime suspect.
20// erosion_node_eval_grad fed a gn whose gradient is EXACT (a closed-form
21// sinusoid), so any error here is the filter's.
22//
23// Build (no engine, no godot-cpp, no scons):
24// g++ -O2 -std=c++17 -I ../proc-nodes/src/shared \
25// -o /tmp/grad_gate tests/test_analytic_gradient.cpp
26// =========================================================================
27#include <cstdio>
28#include <cmath>
29#include <vector>
30#include <algorithm>
31
32#include "erosion_cpp_shim.h"
33
34using dcerosion::vec2;
35using dcerosion::vec3;
36using dcerosion::vec4;
37
38// Angular error between the analytic surface normal and a central-difference
39// one, for a height field h(p) with gradient g. Normals rather than raw
40// gradients because a normal is what the material consumes, and a fixed
41// gradient error matters far less on a cliff than on a flat.
42static double normal_angle_deg(double gx_a, double gz_a, double gx_c, double gz_c) {
43 double ax = -gx_a, ay = 1.0, az = -gz_a;
44 double cx = -gx_c, cy = 1.0, cz = -gz_c;
45 double la = std::sqrt(ax*ax + ay*ay + az*az);
46 double lc = std::sqrt(cx*cx + cy*cy + cz*cz);
47 double d = (ax*cx + ay*cy + az*cz) / (la * lc);
48 d = std::max(-1.0, std::min(1.0, d));
49 return std::acos(d) * 180.0 / M_PI;
50}
51
52struct Stats { double mean, p50, p90, worst; int n; };
53
54static Stats summarize(std::vector<double>& e) {
55 std::sort(e.begin(), e.end());
56 Stats s{};
57 s.n = (int)e.size();
58 double sum = 0; for (double v : e) sum += v;
59 s.mean = sum / e.size();
60 s.p50 = e[e.size() / 2];
61 s.p90 = e[(size_t)(e.size() * 0.9)];
62 s.worst = e.back();
63 return s;
64}
65
66static void report(const char* name, Stats s, double tol_p90) {
67 const char* verdict = (s.p90 <= tol_p90) ? "ok" : "LOSES IT";
68 printf(" %-30s mean %6.3f p50 %6.3f p90 %6.3f max %7.3f deg %s\n",
69 name, s.mean, s.p50, s.p90, s.worst, verdict);
70}
71
72// Central difference of a scalar field sampled through a lambda.
73template <typename F>
74static void central(F f, double x, double z, double e, double& gx, double& gz) {
75 gx = (f(x + e, z) - f(x - e, z)) / (2.0 * e);
76 gz = (f(x, z + e) - f(x, z - e)) / (2.0 * e);
77}
78
79int main() {
80 printf("ANALYTIC GRADIENT GATE — analytic vs central difference, per stage\n");
81 // eps small enough to be a local derivative, large enough that float32
82 // noise evaluated in double-precision differencing does not go to ULP soup.
83 const double E = 1e-4;
84 const int N = 64;
85 const double SPAN = 4.0; // in p units (world / terrain_height)
86
87 std::vector<double> e_base, e_fbm, e_mf, e_ero;
88
89 for (int j = 0; j < N; j++) {
90 for (int i = 0; i < N; i++) {
91 double x = -SPAN * 0.5 + SPAN * (i + 0.5) / N;
92 double z = -SPAN * 0.5 + SPAN * (j + 0.5) / N;
93 double gx, gz;
94
95 // --- 1. base gradient noise ---------------------------------
96 {
97 vec3 r = dcerosion::grad_noised_2d(vec2((float)x, (float)z));
98 central([](double a, double b) {
99 return (double)dcerosion::grad_noised_2d(vec2((float)a, (float)b)).x;
100 }, x, z, E, gx, gz);
101 e_base.push_back(normal_angle_deg(r.y, r.z, gx, gz));
102 }
103
104 // --- 2. plain fBm -------------------------------------------
105 {
106 vec3 r = dcerosion::fractal_noise_2d_lod(vec2((float)x, (float)z), 1.0f, 6, 1.87f, 0.5f, 999.0f);
107 central([](double a, double b) {
108 return (double)dcerosion::fractal_noise_2d_lod(vec2((float)a, (float)b), 1.0f, 6, 1.87f, 0.5f, 999.0f).x;
109 }, x, z, E, gx, gz);
110 e_fbm.push_back(normal_angle_deg(r.y, r.z, gx, gz));
111 }
112
113 // --- 3. elevation-weighted multifractal ---------------------
114 {
115 vec3 r = dcerosion::fractal_multifractal_2d_lod(vec2((float)x, (float)z), 1.0f, 6, 1.87f, 0.5f, 0.5f, 999.0f);
116 central([](double a, double b) {
117 return (double)dcerosion::fractal_multifractal_2d_lod(vec2((float)a, (float)b), 1.0f, 6, 1.87f, 0.5f, 0.5f, 999.0f).x;
118 }, x, z, E, gx, gz);
119 e_mf.push_back(normal_angle_deg(r.y, r.z, gx, gz));
120 }
121
122 // --- 4. erosion filter, fed an EXACT-gradient input ----------
123 // gn is a closed-form sinusoid so its derivative carries no error
124 // of its own; whatever shows up here belongs to erosion_filter.
125 {
126 auto gn_at = [](double a, double b) {
127 double h = 0.30 * std::sin(3.0 * a) * std::cos(3.0 * b);
128 double dx = 0.90 * std::cos(3.0 * a) * std::cos(3.0 * b);
129 double dz = -0.90 * std::sin(3.0 * a) * std::sin(3.0 * b);
130 return vec3((float)h, (float)dx, (float)dz);
131 };
132 auto ero = [&](double a, double b) {
133 vec3 gn = gn_at(a, b);
134 return dcerosion::erosion_node_eval_grad(
135 vec2((float)a, (float)b), gn,
136 5.0f, 0.150000006f, 0.699999988f, 7, 0.5f,
137 4.8f, 0.15f,
138 0.25f, 2.0f, 0.5f, 1.5f, 0.5f,
139 vec2(-0.65f, 0.0f),
140 vec4(0.1f, 0.0f, 0.1f, 2.0f),
141 vec4(1.25f, 1.25f, 2.8f, 1.5f),
142 vec2(0.7f, 1.0f),
143 999.0f);
144 };
145 vec4 r = ero(x, z);
146 central([&](double a, double b) { return (double)ero(a, b).x; }, x, z, E, gx, gz);
147 e_ero.push_back(normal_angle_deg(r.z, r.w, gx, gz));
148 }
149 }
150 }
151
152 printf(" %d samples per stage, eps %.0e, span %.1f p-units\n\n", N * N, E, SPAN);
153 // grad_noised_2d and fBm are exact by construction, so anything above a
154 // rounding-level p90 there is a real defect, not a tolerance to be tuned.
155 report("grad_noised_2d", summarize(e_base), 0.5);
156 report("fractal_noise_2d_lod", summarize(e_fbm), 0.5);
157 report("fractal_multifractal_2d_lod", summarize(e_mf), 0.5);
158 report("erosion_node_eval_grad", summarize(e_ero), 5.0);
159
160 printf("\n erosion error vs OCTAVE COUNT -- combi_mask and fade_target are\n");
161 printf(" REPLACED each octave from that octave's own noise, and neither\n");
162 printf(" derivative is carried. If that is the mechanism, error grows here.\n");
163 for (int oct : {1, 2, 3, 5, 7}) {
164 std::vector<double> e;
165 for (int j = 0; j < N; j++) {
166 for (int i = 0; i < N; i++) {
167 double x = -SPAN * 0.5 + SPAN * (i + 0.5) / N;
168 double z = -SPAN * 0.5 + SPAN * (j + 0.5) / N;
169 auto gn_at = [](double a, double b) {
170 return vec3((float)(0.30 * std::sin(3.0 * a) * std::cos(3.0 * b)),
171 (float)(0.90 * std::cos(3.0 * a) * std::cos(3.0 * b)),
172 (float)(-0.90 * std::sin(3.0 * a) * std::sin(3.0 * b)));
173 };
174 auto ero = [&](double a, double b) {
175 return dcerosion::erosion_node_eval_grad(
176 vec2((float)a, (float)b), gn_at(a, b),
177 5.0f, 0.150000006f, 0.699999988f, oct, 0.5f,
178 4.8f, 0.15f, 0.25f, 2.0f, 0.5f, 1.5f, 0.5f,
179 vec2(-0.65f, 0.0f), vec4(0.1f, 0.0f, 0.1f, 2.0f),
180 vec4(1.25f, 1.25f, 2.8f, 1.5f), vec2(0.7f, 1.0f), 999.0f);
181 };
182 vec4 r = ero(x, z);
183 double gx, gz;
184 central([&](double a, double b) { return (double)ero(a, b).x; }, x, z, E, gx, gz);
185 e.push_back(normal_angle_deg(r.z, r.w, gx, gz));
186 }
187 }
188 Stats s = summarize(e);
189 printf(" %d octave(s) : mean %6.3f p50 %6.3f p90 %6.3f deg\n",
190 oct, s.mean, s.p50, s.p90);
191 }
192 return 0;
193}
tests/test_erosion_parity.cpp +7 -0
……@@ -516,6 +516,13 @@ static std::string build_shader(const std::string& inc) {
516516 // Verified byte-for-byte at the hoist: 33840 bytes, checksum 1254537005,
517517 // before and after.
518518 ErosionWrapperParams wp = erosion_wrapper_params_from_map();
519 // The noise library declares pn_out_ridge (gpu_compiler.cpp) and the
520 // wrapper assigns it, but a standalone shader has no library -- the same
521 // contract that makes the engine prepend "int pn_max_octaves = 99;".
522 // Without this, glslc fails and this harness still prints a healthy
523 // CPU-vs-CPU table, which is how the omission hides.
524 s += "float pn_out_ridge = 0.0;\n";
525 s += "vec2 pn_out_grad = vec2(0.0);\n";
519526 s += erosion_wrapper_glsl(wp);
520527
521528 s +=
tests/test_hypso_curve.cpp +298 -0
……@@ -0,0 +1,298 @@
1// =========================================================================
2// HYPSO CURVE GATE — is the land/sea transfer curve monotonic, C1, and
3// actually bimodal?
4//
5// The curve is the thing that turns a unimodal fBm into continents and
6// ocean basins. Three ways it can be wrong and look fine in a screenshot:
7//
8// NOT MONOTONIC the terrain folds -- two input heights map to the same
9// output, which dual contouring renders as a shelf that
10// is not in the field.
11// NOT C1 a slope discontinuity at the waterline, i.e. a crease
12// ring around every coast. Invisible in a height plot.
13// NOT BIMODAL the whole exercise fails silently: you get a water
14// plane through unimodal noise, which is the thing this
15// replaces.
16//
17// The derivative is checked against a central difference because the
18// wrapper uses it for the chain rule on the gradient the erosion filter
19// steers by -- an analytic derivative that is merely plausible would carve
20// gullies aligned to a surface that is not there.
21//
22// Build (no engine, no godot-cpp, no scons, nothing into bin/):
23// g++ -O2 -std=c++17 -I ../proc-nodes/src/shared \
24// -o /tmp/hypso_gate tests/test_hypso_curve.cpp
25// =========================================================================
26#include <cstdio>
27#include <cmath>
28#include <vector>
29#include <string>
30#include <algorithm>
31#include <cstdlib>
32
33// The .inc is written in the GLSL/C++ common subset and names GLSL builtins,
34// which exist only inside dcerosion -- erosion_cpp_shim.h supplies them there
35// and pulls erosion_filter.inc in the same way (shim :128). Including this at
36// global scope instead is what "'clamp' was not declared in this scope" means.
37#include "erosion_cpp_shim.h"
38namespace dcerosion {
39#include "hypso_curve.inc"
40}
41using dcerosion::hypso_curve;
42using dcerosion::vec2;
43using dcerosion::vec3;
44using dcerosion::fractal_multifractal_2d_lod;
45
46// map_dc_erosion's intended defaults, in NORMALISED units.
47static const float SEA = 0.00f;
48static const float SHELF_W = 0.12f; // gentle shelf below sea
49static const float SLOPE_W = 0.05f; // steep drop half-width
50static const float ABYSS = -0.62f;
51static const float PLAT_SLOPE = 0.13f; // continental platform: gentle
52static const float MTN_ONSET = 0.35f;
53static const float MTN_WIDTH = 0.10f;
54static const float MTN_SLOPE = 1.30f; // mountains: steep
55static const float PEAK_ONSET = 0.60f;
56static const float PEAK_GAIN = 0.30f;
57
58static vec2 C(float h, float pg = PEAK_GAIN) {
59 return hypso_curve(h, SEA, SHELF_W, SLOPE_W, ABYSS, PLAT_SLOPE,
60 MTN_ONSET, MTN_WIDTH, MTN_SLOPE, PEAK_ONSET, pg);
61}
62
63int main() {
64 // ---- 0. DOES IT COMPILE AS GLSL? ------------------------------------
65 //
66 // THE CHECK THIS GATE WAS MISSING. Everything below compiles the .inc as
67 // C++ and says nothing about the shader -- and the first version named a
68 // local `out`, which is an ordinary identifier in C++ and a RESERVED
69 // STORAGE QUALIFIER in GLSL. The C++ gate passed, the shader died with
70 // "unexpected OUT", the GPU silently fell back to the CPU mesher, and the
71 // only symptom a human saw was that terrain generation got slow.
72 //
73 // A dual-compiled file needs BOTH compilers or it is only half gated.
74 {
75 const char* src = "/tmp/hypso_gate_glsl.comp";
76 FILE* f = fopen(src, "w");
77 if (!f) { printf(" GLSL compile : SKIPPED (cannot write %s)\n", src); }
78 else {
79 fputs("#version 450\nlayout(local_size_x=1) in;\n", f);
80 FILE* g = fopen("../proc-nodes/src/shared/hypso_curve.inc", "r");
81 if (!g) { fclose(f); printf(" GLSL compile : SKIPPED (.inc not found)\n"); }
82 else {
83 int c; while ((c = fgetc(g)) != EOF) fputc(c, f);
84 fclose(g);
85 fputs("\nlayout(set=0,binding=0) buffer B { float v[]; };\n"
86 "void main() { v[0] = hypso_curve(v[1],v[2],v[3],v[4],v[5],"
87 "v[6],v[7],v[8],v[9],v[10],v[11]).x; }\n", f);
88 fclose(f);
89 int rc = system("glslc -fshader-stage=compute /tmp/hypso_gate_glsl.comp "
90 "-o /dev/null 2>/tmp/hypso_gate_glsl.err");
91 printf(" GLSL compile : %s\n", rc == 0 ? "ok" : "FAIL");
92 if (rc != 0) system("head -3 /tmp/hypso_gate_glsl.err");
93 }
94 }
95 }
96
97 printf("HYPSO CURVE GATE\n");
98 printf(" sea %.2f shelf_w %.2f slope_w %.2f abyss %.2f plat %.2f |"
99 " mtn %.2f/%.2f/%.2f | peak %.2f/%.2f\n\n",
100 SEA, SHELF_W, SLOPE_W, ABYSS, PLAT_SLOPE, MTN_ONSET, MTN_WIDTH,
101 MTN_SLOPE, PEAK_ONSET, PEAK_GAIN);
102
103 // ---- 1. MONOTONIC ---------------------------------------------------
104 // The real condition is dout/dh >= 0 everywhere -- that is monotonicity by
105 // construction. Comparing consecutive OUTPUTS in float32 also reports
106 // decreases of ~1e-5 that are accumulated rounding in the blend, not folds;
107 // a genuine fold is orders larger. So report both, and judge on the
108 // derivative while showing the float noise for what it is.
109 const int N = 20001;
110 float prev = C(-1.0f).x, worst_back = 0.0f; int backs = 0;
111 double min_d = 1e30; float min_at = 0.0f;
112 for (int i = 1; i < N; i++) {
113 float h = -1.0f + 2.0f * (float)i / (N - 1);
114 vec2 v = C(h);
115 if (v.x < prev) { backs++; worst_back = std::max(worst_back, prev - v.x); }
116 if (v.y < min_d) { min_d = v.y; min_at = h; }
117 prev = v.x;
118 }
119 printf(" min dout/dh : %+.5f at h=%+.3f %s\n", min_d, min_at,
120 min_d >= 0.0 ? "ok -- monotonic by construction" : "FAIL -- folds");
121 printf(" float32 backsteps: %d, worst %.2e (%s)\n", backs, worst_back,
122 worst_back < 1e-4 ? "rounding, not a fold" : "REAL FOLD");
123
124 // ---- 2. DERIVATIVE vs CENTRAL DIFFERENCE ----------------------------
125 double worst_d = 0.0; float worst_at = 0.0f;
126 const float eps = 1e-4f;
127 for (int i = 0; i <= 4000; i++) {
128 float h = -1.0f + 2.0f * (float)i / 4000.0f;
129 double fd = (C(h + eps).x - C(h - eps).x) / (2.0 * eps);
130 double an = C(h).y;
131 double e = fabs(fd - an);
132 if (e > worst_d) { worst_d = e; worst_at = h; }
133 }
134 printf(" d/dh vs finite : worst %.5f at h=%+.3f %s\n",
135 worst_d, worst_at, worst_d < 2e-2 ? "ok" : "FAIL");
136
137 // ---- 3. C1, BY SCALING -----------------------------------------------
138 // The previous version compared adjacent derivative samples against a
139 // fixed threshold, which measures |f''|*dh, not a discontinuity -- a
140 // smoothstep of width w has |f''| ~ 6/w^2, so it "failed" for being
141 // curved. A CONTINUOUS derivative's max adjacent step HALVES when the
142 // sampling halves; a jump does not move at all.
143 auto max_step = [&](int n) {
144 double worst = 0.0; float p = C(-1.0f).y;
145 for (int i = 1; i <= n; i++) {
146 float h = -1.0f + 2.0f * (float)i / (float)n;
147 float d = C(h).y;
148 worst = std::max(worst, (double)fabs(d - p));
149 p = d;
150 }
151 return worst;
152 };
153 double s1 = max_step(20000), s2 = max_step(40000);
154 double ratio = s2 > 0 ? s1 / s2 : 0.0;
155 printf(" C1 (scaling) : step %.5f -> %.5f, ratio %.2f %s\n",
156 s1, s2, ratio, (ratio > 1.7 && ratio < 2.3) ? "ok -- continuous"
157 : "FAIL -- jump");
158
159 // ---- 4. BIMODAL, against the null ------------------------------------
160 // Input: a triangular-ish distribution standing in for summed octaves.
161 // What matters is that the OUTPUT histogram has two modes where the
162 // INPUT has one, so the same input is measured through both.
163 const int S = 400000, B = 24;
164 std::vector<int> hin(B, 0), hout(B, 0);
165 unsigned s = 12345u;
166 auto rnd = [&]() { s = s * 1664525u + 1013904223u;
167 return (float)((s >> 8) & 0xFFFFFF) / (float)0xFFFFFF; };
168 float omin = C(-1.0f).x, omax = C(1.0f).x;
169 for (int i = 0; i < S; i++) {
170 float h = (rnd() + rnd() + rnd()) / 1.5f - 1.0f; // ~triangular on [-1,1]
171 float o = C(h).x;
172 int bi = (int)((h + 1.0f) * 0.5f * (B - 1) + 0.5f);
173 int bo = (int)((o - omin) / (omax - omin) * (B - 1) + 0.5f);
174 hin[std::min(std::max(bi,0),B-1)]++;
175 hout[std::min(std::max(bo,0),B-1)]++;
176 }
177 // Two modes with a real dip between them, not just "two bumps".
178 auto two_modes = [&](const std::vector<int>& hh, double& dip) {
179 int n = (int)hh.size(), i1 = 0;
180 for (int i = 0; i < n; i++) if (hh[i] > hh[i1]) i1 = i;
181 int i2 = -1;
182 for (int i = 0; i < n; i++) {
183 if (abs(i - i1) < 3) continue;
184 bool loc = (i == 0 || hh[i] >= hh[i-1]) && (i == n-1 || hh[i] >= hh[i+1]);
185 if (loc && (i2 < 0 || hh[i] > hh[i2])) i2 = i;
186 }
187 if (i2 < 0) { dip = 0; return 0; }
188 int a = std::min(i1, i2), b = std::max(i1, i2), tr = hh[a];
189 for (int i = a; i <= b; i++) tr = std::min(tr, hh[i]);
190 dip = 1.0 - (double)tr / (double)hh[i2]; // 1 = clean separation
191 return 2;
192 };
193 double dip_in = 0, dip_out = 0;
194 int m_in = two_modes(hin, dip_in), m_out = two_modes(hout, dip_out);
195
196 printf("\n output range : %.3f .. %.3f (waterline at %.4f)\n",
197 omin, omax, C(SEA).x);
198 printf(" input : %d mode(s), dip %.2f\n", m_in, dip_in);
199 printf(" output : %d mode(s), dip %.2f %s\n\n", m_out, dip_out,
200 (m_out == 2 && dip_out > 0.35) ? "ok -- BIMODAL"
201 : "FAIL -- not separated");
202 int pk = *std::max_element(hout.begin(), hout.end());
203 for (int i = 0; i < B; i++) {
204 float lo = omin + (omax - omin) * i / (float)(B - 1);
205 printf(" %+6.2f |%s\n", lo, std::string(hout[i] * 46 / pk, '#').c_str());
206 }
207
208 // ---- 5. PEAK LIFT TOUCHES ONLY THE TAIL ------------------------------
209 double worst_low = 0.0;
210 for (int i = 0; i <= 2000; i++) {
211 float h = -1.0f + (PEAK_ONSET + 1.0f) * (float)i / 2000.0f;
212 worst_low = std::max(worst_low, (double)fabs(C(h, PEAK_GAIN).x - C(h, 0.0f).x));
213 }
214 printf("\n peak_gain effect below onset: %.3e %s\n", worst_low,
215 worst_low < 1e-6 ? "ok -- lowlands untouched" : "FAIL");
216 printf(" peak_gain effect at h=1.0 : %.3f (asked for %.3f)\n",
217 C(1.0f, PEAK_GAIN).x - C(1.0f, 0.0f).x, PEAK_GAIN);
218 // ---- 6. LAND FRACTION, against the REAL mountain field ---------------
219 //
220 // The histogram above uses a synthetic input, which is fine for shape but
221 // useless for "how much of the world is land" -- that depends on the
222 // actual distribution the curve is fed. erosion_filter.inc carries
223 // fractal_multifractal_2d_lod, the very function map_dc_erosion's wrapper
224 // calls, so this samples THAT rather than a stand-in.
225 //
226 // map_dc_erosion.gd's constants (MOUNTAIN_SCALE is post-0.1).
227 const float TH = 60000.0f, MSCALE = 0.1f, LAC = 1.87f, MGAIN = 0.68f;
228 const int MOCT = 8;
229 const float VARIANCE = 0.0f; // generator default
230 const float HALF = 458800.0f; // depth-16 world, +/- 458.8 km
231 const int G = 420; // 176400 columns
232
233 printf("\n --- LAND FRACTION over the real mountain field ---\n");
234 printf(" terrain_height %.0f mountain_scale %.2f %d oct lac %.2f"
235 " gain %.2f\n", TH, MSCALE, MOCT, LAC, MGAIN);
236
237 std::vector<float> raw;
238 raw.reserve((size_t)G * G);
239 for (int j = 0; j < G; j++) {
240 for (int i = 0; i < G; i++) {
241 float x = -HALF + 2.0f * HALF * i / (G - 1);
242 float z = -HALF + 2.0f * HALF * j / (G - 1);
243 vec3 m = fractal_multifractal_2d_lod(
244 vec2(x / TH * MSCALE, z / TH * MSCALE), 1.0f,
245 MOCT, LAC, MGAIN, VARIANCE, 999.0f);
246 raw.push_back(m.x);
247 }
248 }
249 float rmin = raw[0], rmax = raw[0];
250 for (float v : raw) { rmin = std::min(rmin, v); rmax = std::max(rmax, v); }
251 printf(" raw multifractal range: %+.3f .. %+.3f (curve assumes ~[-1,1])\n",
252 rmin, rmax);
253
254 auto land_frac = [&](float sea_v) {
255 float wl = hypso_curve(sea_v, sea_v, SHELF_W, SLOPE_W, ABYSS, PLAT_SLOPE,
256 MTN_ONSET, MTN_WIDTH, MTN_SLOPE,
257 PEAK_ONSET, PEAK_GAIN).x;
258 size_t n = 0;
259 for (float v : raw) {
260 float o = hypso_curve(v, sea_v, SHELF_W, SLOPE_W, ABYSS, PLAT_SLOPE,
261 MTN_ONSET, MTN_WIDTH, MTN_SLOPE,
262 PEAK_ONSET, PEAK_GAIN).x;
263 if (o > wl) n++;
264 }
265 return 100.0 * (double)n / (double)raw.size();
266 };
267
268 printf("\n sea land %% (Earth is ~29%%)\n");
269 for (int k = -6; k <= 6; k++) {
270 float sv = 0.05f * k;
271 printf(" %+5.2f %5.1f%s\n", sv, land_frac(sv),
272 (k == 0) ? " <- current" : "");
273 }
274 printf("\n at sea=%.2f: land %.1f%%, mountains (h>%.2f) %.1f%%\n",
275 SEA, land_frac(SEA), MTN_ONSET,
276 100.0 * (double)std::count_if(raw.begin(), raw.end(),
277 [&](float v){ return v > MTN_ONSET; }) / raw.size());
278 // ---- 7. THE SHELF, which is the whole point of this revision --------
279 printf("\n --- COAST PROFILE ---\n");
280 float wl = C(SEA).x;
281 printf(" waterline (curve units) : %+.4f\n", wl);
282 float brk = C(SEA - SHELF_W + SLOPE_W).x;
283 printf(" depth at the shelf break : %+.4f (%.0f m at ~7790 m/unit)\n",
284 brk - wl, (brk - wl) * 7790.0f);
285 // gradient just above sea vs at the steepest point
286 float g_land = C(SEA + 0.02f).y;
287 float g_max = 0.0f; float g_at = 0.0f;
288 for (int i = 0; i <= 4000; i++) {
289 float h = -1.0f + 2.0f * (float)i / 4000.0f;
290 if (C(h).y > g_max) { g_max = C(h).y; g_at = h; }
291 }
292 printf(" gradient at the shore : %.3f\n", g_land);
293 printf(" gradient at the steepest pt : %.3f at h=%+.3f (%s sea)\n",
294 g_max, g_at, g_at < SEA ? "BELOW -- submerged, good" : "ABOVE -- a cliff");
295 printf(" steepest / shore ratio : %.1fx %s\n", g_max / g_land,
296 g_at < SEA ? "(the drop is under water)" : "(THE DROP IS THE COAST)");
297 return 0;
298}