// ============================================================================
// THE MAP WRAPPER, AS ONE STRING, EMITTED ONCE.
//
// erosion_filter.inc is the erosion algorithm and nothing else. Around it sits
// a WRAPPER -- the mountain layer, the relief modulation, the detail layer, the
// y-mapping and the floor/ceiling clip -- which is what turns the filter into
// a particular map. That wrapper existed twice inside test_erosion_parity.cpp
// (as GLSL text and as C++), and the moment a second GPU consumer appeared it
// was about to exist a third time.
//
// #45 spent a session reducing the FILTER from four copies to one. This is the
// same move one level out, for the WRAPPER. Every consumer that needs the
// wrapper as GLSL calls erosion_wrapper_glsl() and gets the same text.
//
// WHAT THIS IS NOT
//
// It is not the engine's wrapper. DCTerrainGeneratorProcNodes::_recompile_glsl
// assembles the shipping string and is unreachable without booting godot; the
// text below is test_erosion_parity.cpp's TRANSCRIPTION of it, annotated by
// file and line, and guarded only in the sense that the oracle checks the GLSL
// transcription against the C++ transcription. Hoisting the engine's own
// assembly out of godot is still the open item, and is what would make this
// text the shipping artifact rather than a faithful copy of it.
//
// So: this removes the SECOND-consumer duplication. It does not close the
// engine gap. Do not read it as more than that.
//
// ON THE DEFAULTS BELOW: they are a CONVENIENCE, not the map. They were typed
// by hand once and were wrong within hours -- MAP_FLOOR_Y moved from -1000 to
// -10900 in dc-terrain 156b2e4 and this file did not. A caller that wants the
// map must fill the struct from erosion_wrapper_params_from_map() in
// test_erosion_parity.cpp, which is the single place those constants live.
// Treat a default here as a starting point for a knob, never as ground truth.
// ============================================================================
#pragma once
#include <cstdio>
#include <string>
#include <vector>
struct ErosionWrapperParams {
// map_dc_erosion.gd -- the terrain field
float floor_y = -10900.0f; // FLOOR_Y (see note above)
float ceiling_y = 10100.0f; // CEILING_Y
float terrain_height = 60000.0f; // :48
float mountain_scale = 1.0f; // :51
float mountain_height = 0.55f; // :55
float mountain_gain = 0.68f;
float mountain_variance = 0.0f; // generator "mountain_variance"
// 0 = Musgrave height-weighted multifractal (mountain_variance is its knob)
// 1 = iq's gradient trick (mountain_damping is its knob)
int mountain_mode = 0;
float mountain_damping = 1.0f;
int mountain_octaves = 8; // :53
float relief_scale = 0.5f; // :60
float relief_floor = 0.08f; // :62
float detail_scale = 600.0f; // :79
float detail_height = 30.0f; // :80
float detail_gain = 0.5f; // :86
int detail_octaves = 4; // :85
float noise_lacunarity = 1.87f; // :92
// The LOD voxel size the Nyquist clamp is keyed to. 1.0 is LOD 0, where
// nothing is clamped. Raising it band-limits the generator, which is what
// a coarse LOD does and what any sampled measurement of it must do.
float voxel_size = 1.0f; // :69
// Graph node defaults, map_dc_erosion.gd:122-133
float gn_scale = 3.0f; // GradientNoise "Scale"
float gn_gain = 0.1f; // GradientNoise "Gain"
int gn_detail = 3; // GradientNoise "Detail"
float en_height_scale = 5.0f; // ErosionNoise "HeightScale"
float en_strength = 0.15f; // ErosionNoise "ErosionStr"
float en_slope_power = 0.6f; // accepted and IGNORED by the fade approach
float en_cell_scale = 0.7f; // ErosionNoise "CellScale"
float en_gain = 0.5f; // ErosionNoise "Gain"
int en_octaves = 7; // ErosionNoise "Octaves"
float en_erosion_scale = 0.15f; // generator "erosion_scale"
float en_fade_div = 4.8f; // altitude-fade divisor, wrapper literal
// The filter parameters that erosion_node_eval BAKES on the engine path,
// because GpuCompiler's push-constant block is seven floats and rounding
// plus onset alone need eight. A consumer that emits GLSL text has no such
// cap, so the lookdev tool can turn them. Defaults are runevision's
// reference values (lpmitchell/AdvancedTerrainErosion README, and the
// Shadertoy lineage), so leaving them alone reproduces the engine exactly.
//
// rounding.z is the one the blog post singles out: "the erosion function
// doesn't know about the sizes of the terrain features coming from [the
// input height function], so an input value is provided for tweaking the
// rounding of those." It is the input-feature-size control.
float en_height_amp = 0.25f;
float en_erosion_lac = 2.0f;
float en_gully_weight = 0.5f;
float en_detail = 1.5f;
float en_normalization = 0.5f;
float en_height_offset_a = -0.65f;
float en_height_offset_b = 0.0f;
float rounding_ridges = 0.1f; // rounding.x
float rounding_creases = 0.0f; // rounding.y
float rounding_input = 0.1f; // rounding.z <- input feature size
float rounding_per_oct = 2.0f; // rounding.w
float onset_input = 1.25f; // onset.x
float onset_per_oct = 1.25f; // onset.y
float onset_ridge_input = 2.8f; // onset.z
float onset_ridge_oct = 1.5f; // onset.w
float assumed_slope_mag = 0.7f; // assumed_slope.x
float assumed_slope_mix = 1.0f; // assumed_slope.y
// ---- THE FIELD'S OWN RELIEF, for the SLOPE only ----------------------
//
// The vertical unit for the slope is NOT the one for the height, and
// conflating them is why the filter under-reads the terrain.
//
// HEIGHT normalises by H = y_range*height_amp*0.5, and it must: that
// mapping is what puts the terrain at the right ALTITUDE, and the map's
// HeightScale is its exact inverse (map_dc_diffusion.gd says HeightScale
// MUST be 5.0 for this reason). Changing it shifts the whole world
// vertically -- the 2-6.5 m offset bug.
//
// SLOPE has no such constraint. It feeds slope_length, which gates every
// structure-sensitive mask in the filter, and what those masks want is
// steepness relative to the FIELD's own relief -- which is what
// erode_heightmap.cpp gives them by normalising per tile.
//
// y_range is the PLAYABLE vertical band (-10900..10100 = 21000 m). The
// field only ever spans -4829..5110 (measured over all 841 cached tiles),
// so normalising slope by the band under-reads it by about 2.1x.
//
// GLOBAL, NOT PER-TILE, and that is the whole point. erode_heightmap uses
// each tile's own min/max, so two adjacent tiles of different relief
// normalise differently and the join steps -- which is why this header
// removed per-tile normalisation in the first place. A constant taken from
// the WHOLE field is seam-free by construction, and stays seam-free when
// the single block becomes many.
//
// ⚠ gn.yz is therefore NOT the exact gradient of gn.x any more. That is
// safe HERE because the filter never integrates it: n.x feeds fade_target,
// n.yz feeds gully direction and the masks, and nothing reconstructs one
// from the other. Do not reuse this vec3 as a (value, gradient) pair.
float field_half_range = 0.0f; // 0 = fall back to y_range*0.5
// ---- EXTERNAL HEIGHTFIELD SOURCE (hyper-games #50 section 5.3) ---------
//
// When true, gn comes from a bound heightfield in METRES instead of from
// fractal_noise_2d + terrain_offset. The consumer must define
//
// float hblk_height_m(float x, float z);
//
// ABOVE the wrapper text. The split is deliberate: STORAGE differs per
// consumer (an SSBO in the standalone harness, another set/binding in the
// engine) but the UNIT CONVERSION must not, because getting that wrong is
// exactly the 60x wrapper bug of #45. The conversion lives here, once, and
// is derived rather than typed -- see erosion_heightmap_gn_glsl below.
//
// It also removes the PER-TILE normalisation erode_heightmap.cpp uses.
// That tool takes mid and H from the tile's own min/max, so two adjacent
// tiles of different relief normalise differently and the join between
// them steps. Here both come from the map's floor/ceiling, which are world
// constants -- so the eroded field is seam-free for the same reason the
// diffusion field underneath it is.
bool heightmap_source = false;
// Central-difference step in METRES. Should be the source field's own
// pixel pitch: smaller re-reads one texel and returns quantisation noise,
// larger smooths real slope away.
float heightmap_step = 30.0f;
};
// The constant block. Split out so a caller can see exactly which names the
// wrapper text below depends on.
inline std::string erosion_wrapper_consts_glsl(const ErosionWrapperParams& p) {
char cfg[4096];
snprintf(cfg, sizeof(cfg),
"const float M_FLOOR_Y=%.9g, M_CEIL_Y=%.9g, M_TH=%.9g;\n"
"const float M_MSCALE=%.9g, M_MHEIGHT=%.9g, M_MGAIN=%.9g;\n"
"const float M_MVAR=%.9g;\n"
"const int M_MMODE=%d;\n"
"const float M_MDAMP=%.9g;\n"
"const float M_VOXEL=%.9g;\n"
"const int M_MOCT=%d;\n"
"const float M_RSCALE=%.9g, M_RFLOOR=%.9g;\n"
"const float M_DSCALE=%.9g, M_DHEIGHT=%.9g, M_DGAIN=%.9g;\n"
"const int M_DOCT=%d;\n"
"const float M_LAC=%.9g;\n"
"const float GN_SCALE=%.9g, GN_GAIN=%.9g;\n"
"const int GN_DETAIL=%d;\n"
"const float EN_HS=%.9g, EN_STR=%.9g, EN_SP=%.9g, EN_CS=%.9g, EN_GAIN=%.9g;\n"
"const int EN_OCT=%d;\n"
"const float EN_ESCALE=%.9g;\n"
"const float EN_FADEDIV=%.9g;\n"
"const float EN_HAMP=%.9g, EN_ELAC=%.9g, EN_GW=%.9g, EN_DET=%.9g, EN_NORM=%.9g;\n"
"const vec2 EN_HOFF=vec2(%.9g, %.9g);\n"
"const vec4 EN_ROUND=vec4(%.9g, %.9g, %.9g, %.9g);\n"
"const vec4 EN_ONSET=vec4(%.9g, %.9g, %.9g, %.9g);\n"
"const vec2 EN_ASLOPE=vec2(%.9g, %.9g);\n"
"const float EN_FIELDH=%.9g;\n",
p.floor_y, p.ceiling_y, p.terrain_height,
p.mountain_scale, p.mountain_height, p.mountain_gain,
p.mountain_variance, p.mountain_mode, p.mountain_damping,
p.voxel_size, p.mountain_octaves,
p.relief_scale, p.relief_floor,
p.detail_scale, p.detail_height, p.detail_gain, p.detail_octaves,
p.noise_lacunarity,
p.gn_scale, p.gn_gain, p.gn_detail,
p.en_height_scale, p.en_strength, p.en_slope_power, p.en_cell_scale,
p.en_gain, p.en_octaves,
p.en_erosion_scale,
p.en_fade_div,
p.en_height_amp, p.en_erosion_lac, p.en_gully_weight, p.en_detail,
p.en_normalization,
p.en_height_offset_a, p.en_height_offset_b,
p.rounding_ridges, p.rounding_creases, p.rounding_input, p.rounding_per_oct,
p.onset_input, p.onset_per_oct, p.onset_ridge_input, p.onset_ridge_oct,
p.assumed_slope_mag, p.assumed_slope_mix,
p.field_half_range);
std::string out(cfg);
// Emitted ONLY in heightmap mode, so the default emission stays
// byte-identical to what every existing consumer already compiles.
if (p.heightmap_source) {
char hb[128];
snprintf(hb, sizeof(hb), "const float HBLK_STEP=%.9g;\n", p.heightmap_step);
out += hb;
}
return out;
}
// The wrapper itself. Requires erosion_filter.inc to have been emitted above
// it (it calls nyquist_octaves, fractal_*_2d_lod and erosion_node_eval).
//
// `gpu_wrapper` and `cpu_wrapper` are IDENTICAL text and that is deliberate:
// the parity oracle compares each against its own C++ twin, so emitting one
// and aliasing the other would silently halve what the oracle can catch.
// The heightfield -> filter-space conversion, DERIVED from
// erosion_node_eval_full rather than assumed. That function does
//
// n.x = gn.x * height_amp, then n.x = n.x*0.5 + 0.5
// eroded = n.x + delta + offset, return eroded * height_scale
//
// and the wrapper compares the result against
//
// y_mapped = (y - M_FLOOR_Y) * (EN_HS / y_range)
//
// Zeroing delta and offset (strength 0, EN_HOFF zero) and solving
// en.x == y_mapped for the y that reproduces the input height gives
//
// gn.x = (h - mid) / H mid = (M_FLOOR_Y + M_CEIL_Y) / 2
// H = y_range * EN_HAMP / 2
//
// which is erode_heightmap.cpp's form with WORLD constants in place of its
// per-tile min/max. The slope follows by the chain rule -- but NOT with M_TH
// as the horizontal unit, which is the bug this note used to describe.
//
// THE SLOPE MUST BE IN THE UNITS OF THE GULLIES IT STEERS. p is world/M_TH, but
// the filter rescales p by gully_band_scale internally, so the gullies live at
// M_TH * EN_ESCALE, not at M_TH. Handing the filter a slope measured in M_TH
// units over-scales it by 1/EN_ESCALE -- 700x at the map's settings. Every
// structure-sensitive term is gated on slope_length (combi_mask at the
// ease_out(smooth_start(slope_length * onset.x, ..)) line, and the ridge mask
// below it), so at 700x they saturate everywhere, the mask goes flat, and what
// survives is the phacelle noise: an added gully pattern with no reference to
// the field's own terrain. Measured corr(|erosion|, terrain slope) = 0.020,
// against 0.368 for erode_heightmap.cpp on the same field and band.
//
// M_TH * EN_ESCALE is not a correction factor, it is the reference's own
// horizontal unit: erode_heightmap uses L = gully_metres / cell_scale with
// erosion_scale 1, and 60000 * 0.00142857 = 85.7143 = 60 / 0.7 exactly.
//
// gn.yz = (dh/dworld) * M_TH * EN_ESCALE / H.
//
// gate_heightmap checks exactly this: at strength 0 with no height offset the
// surface must come back EQUAL to the input field. A unit error cannot survive
// that, which is the whole reason it is written as an identity.
// ---------------------------------------------------------------------------
// THE MIP PYRAMID.
//
// Every other layer in this generator band-limits itself per LOD through
// nyquist_octaves(..., voxel_scale, ..), dropping octaves as the voxel grows. A
// sampled heightfield has no such property: at LOD 6 the mesher point-samples a
// 30 m field at 64 m spacing and the surface jumps between neighbouring voxels.
// Dual contouring then faithfully meshes an aliased field, which is why the
// artefacts appear on COARSE chunks and the fine ones are clean.
//
// Level k texel j is the box average of level-0 texels [j*2^k, (j+1)*2^k),
// computed DIRECTLY from level 0 -- repeated halving walks the footprint off the
// power of two the sampler's half-texel shift assumes, and the error accumulates
// down the chain.
struct HeightBlockMips {
std::vector<float> data; // every level, concatenated
std::vector<int> off; // start index of each level
std::vector<int> side; // side length of each level
int levels = 0;
};
inline HeightBlockMips erosion_heightmap_build_mips(const float* level0, int n) {
HeightBlockMips m;
for (int k = 0;; k++) {
int s = n >> k;
if (s < 1) s = 1;
m.side.push_back(s);
m.off.push_back((int)m.data.size());
if (k == 0) {
m.data.insert(m.data.end(), level0, level0 + (size_t)n * (size_t)n);
} else {
const int f = 1 << k;
for (int j = 0; j < s; j++) {
for (int i = 0; i < s; i++) {
double acc = 0.0; int cnt = 0;
for (int dj = 0; dj < f; dj++) {
int sj = j * f + dj; if (sj >= n) break;
for (int di = 0; di < f; di++) {
int si = i * f + di; if (si >= n) break;
acc += (double)level0[(size_t)sj * (size_t)n + (size_t)si];
cnt++;
}
}
m.data.push_back(cnt ? (float)(acc / (double)cnt) : 0.0f);
}
}
}
if (s == 1) { m.levels = k + 1; break; }
}
return m;
}
inline std::string erosion_heightmap_levels_glsl(const HeightBlockMips& m) {
std::string s = "const int HBLK_LEVELS = " + std::to_string(m.levels) + ";\n";
s += "const int HBLK_OFF[" + std::to_string(m.levels) + "] = int[" +
std::to_string(m.levels) + "](";
for (int i = 0; i < m.levels; i++) { if (i) s += ", "; s += std::to_string(m.off[i]); }
s += ");\n";
s += "const int HBLK_SIDE[" + std::to_string(m.levels) + "] = int[" +
std::to_string(m.levels) + "](";
for (int i = 0; i < m.levels; i++) { if (i) s += ", "; s += std::to_string(m.side[i]); }
s += ");\n";
return s;
}
// ============================================================================
// THE TILE GRID — one buffer holding an INCOMPLETE, GROWING set of tiles
// ============================================================================
//
// The 900 km world is 119x119 = 14,161 tiles and 841 exist (5.9%), so the set
// is INCOMPLETE and it GROWS. A single stitched block cannot express that: it
// fixes membership at bake time, and a dense all-841 block CANNOT EVEN BE
// ALLOCATED (side 27,650 -> H = 248,868 at 4096 wide against a measured
// maxImageDimension2D of 32,768).
//
// So custom_data becomes [ INDEX TABLE | per-tile rects ]: nx*ny floats, one
// per world cell, holding the tile's start offset or -1 for absent. This is
// the SSBO layout, and it changes the CONTENTS of custom_data and the emitted
// GLSL but NOT the binding set -- a materially safer class of change than the
// texture-atlas swap, which was reverted.
//
// ⛔ AN ABSENT TILE IS A STUB, NOT A CLAMP. Clamping smears the nearest loaded
// row across the hole, which reads as terrain -- the silent-seam failure. A
// flat stub is visibly a stub.
//
// ⚠ GUTTERS DO NOT BITE THIS LAYOUT, and that is a property of the SSBO rather
// than a thing anyone chose. A Catmull-Rom tap that crosses a tile boundary
// resolves through the index table to the NEIGHBOUR'S OWN payload texel, so
// the stored 1 px apron is never read here and no margin is needed. The
// 2-texel-per-level gutter requirement is real for the paged ATLAS design,
// where a page is a rect and a tap cannot leave it.
//
// ⚠ AXIS SWAP, INHERITED. stitch_block.py: "tx indexes the FIRST array axis ->
// WORLD Z; ty indexes the SECOND -> WORLD X". Getting it backwards mirrors the
// world and still looks plausible, which is why the gate reads coordinates from
// filenames and has a control that swaps them deliberately.
struct HeightTileSrc {
int tx = 0, ty = 0; // grid coords; tx -> world Z, ty -> world X
int n = 0; // STORED side, including the apron
int apron = 1;
float mpp = 30.0f;
std::vector<float> h; // n*n, row = z, column = x
};
// THE WORLD, not the footprint. Fixed for a given octree + tile size, and the
// reason the emitted constants can stop changing: see erosion_heightmap_pack_tiles.
struct HeightWorld {
int tx0 = 0, ty0 = 0, nx = 0, ny = 0;
};
// Derive the world grid from the octree span and the tile span, both in metres.
// map_dc_diffusion.gd: root span = 28 * 2^(OCTREE_DEPTH-1) * VOXEL_SIZE =
// 917,504 m, and a tile is 256 px * 30 m = 7,680 m, so 119 tiles per axis and
// 119^2 = 14,161 cells -- the count the specs quote.
inline HeightWorld erosion_heightmap_world(double octree_span_m, double tile_span_m) {
HeightWorld w;
const int n = (int)(octree_span_m / tile_span_m);
w.nx = w.ny = n;
w.tx0 = w.ty0 = -(n / 2);
return w;
}
// Everything needed to ADDRESS the packed buffer, without the buffer. The GLSL
// emitter, the C++ twin and the engine all take this, so there is one
// description of the layout rather than three agreeing by hand.
struct HeightGridDesc {
int nx = 0, ny = 0; // cells along tx (Z) and ty (X)
int tx0 = 0, ty0 = 0; // grid origin in tile coords
int payload = 0; // interior side, apron excluded
int apron = 0;
int n = 0; // stored side = payload + 2*apron
int table = 0; // floats before the first tile rect = nx*ny
float mpp = 30.0f;
float stub = 0.0f;
bool active() const { return nx > 0 && ny > 0 && payload > 0; }
// World metres of GRID texel (0,0). HBLK_OX/OZ, so hblk_sample_level's
// world->texel mapping needs no change.
float ox() const { return (float)((double)ty0 * payload * mpp); }
float oz() const { return (float)((double)tx0 * payload * mpp); }
// The texel extent the sampler clamps to.
int side() const { return (nx > ny ? nx : ny) * payload; }
};
struct HeightGridPack {
std::vector<float> data; // [nx*ny index table | tile rects]
HeightGridDesc desc;
int tiles = 0;
int rejected = 0; // tiles refused for sitting outside the world window
bool ok = false;
};
// Pack a tile list into one buffer.
//
// ★ THE TABLE IS SIZED BY THE WORLD, NOT BY THE LOADED FOOTPRINT, and that is
// the whole difference between a tileset that LOADS and one that MUTATES.
//
// HBLK_G_NX/NY/TX0/TY0/TABLE are emitted as GLSL COMPILE-TIME CONSTANTS. Size
// them to the tiles you happen to hold and the first tile that arrives outside
// that bounding box changes them, which RECOMPILES THE DENSITY SHADER -- the
// defect tile-residency-and-generation.md section 1.2 records as "residency is
// a compile-time constant". Size them to the world and they are fixed forever:
// a tile arriving anywhere is then a pure DATA upload, and dc_gpu_thread.h
// already change-detects the shader and the data SEPARATELY (:262-282), so the
// no-recompile path exists and needs nothing new.
//
// The cost is one float per world cell: 119^2 = 14,161 = 55 KiB. Negligible,
// and exactly the page table tile-residency-and-generation.md section 2.4 asks
// for.
//
// ⛔ A TILE OUTSIDE THE WORLD WINDOW IS REFUSED, NOT WRAPPED. The generated
// cache reaches tx 77 (591 km) while a centred 119-tile window stops at 59, so
// out-of-window tiles REALLY EXIST. Aliasing one into an occupied cell would
// silently replace real terrain somewhere else in the world.
//
// Returns ok=false for a MIXED set rather than picking a winner: tiles of
// different side or pitch in one grid is a generator bug, and silently packing
// them would produce terrain that is wrong in a way no probe would name.
inline HeightGridPack erosion_heightmap_pack_tiles(
const std::vector<HeightTileSrc>& tiles, float stub,
const HeightWorld& world, bool swap_axes = false) {
HeightGridPack g;
HeightGridDesc &D = g.desc;
D.stub = stub;
if (tiles.empty() || world.nx <= 0 || world.ny <= 0) return g;
D.n = tiles[0].n;
D.apron = tiles[0].apron;
D.mpp = tiles[0].mpp;
D.payload = D.n - 2 * D.apron;
if (D.payload <= 0) return g;
for (size_t k = 0; k < tiles.size(); k++) {
const HeightTileSrc& t = tiles[k];
if (t.n != D.n || t.apron != D.apron || t.mpp != D.mpp) return g;
if ((int)t.h.size() != D.n * D.n) return g;
}
D.tx0 = world.tx0; D.ty0 = world.ty0;
D.nx = world.nx; D.ny = world.ny;
D.table = D.nx * D.ny;
// -1 is ABSENT. The table holds a tile INDEX, NOT a buffer offset, and that
// is not a style choice -- it is a precision one.
//
// custom_data is float32, which represents integers exactly only up to
// 2^24 = 16,777,216. A tile is 258^2 = 66,564 floats, so a table of
// OFFSETS goes inexact at 16,777,216 / 66,564 = the 252nd tile. The
// residency budget is 512 slots, so an offset table would start silently
// rounding to the wrong texel INSIDE the intended working set, and a
// rounded offset does not fail -- it reads a neighbouring row of a
// different tile.
//
// An INDEX stays small: 511 for the budget, 14,160 for the whole world,
// both exact by an enormous margin. The shader multiplies it back up in
// int, where 14,161 * 258^2 = 942,653,124 still fits int32 with room.
g.data.assign((size_t)D.nx * (size_t)D.ny, -1.0f);
for (size_t k = 0; k < tiles.size(); k++) {
const HeightTileSrc& t = tiles[k];
const int tx = swap_axes ? t.ty : t.tx;
const int ty = swap_axes ? t.tx : t.ty;
if (tx < D.tx0 || tx >= D.tx0 + D.nx || ty < D.ty0 || ty >= D.ty0 + D.ny) {
g.rejected++; // outside the world window -- refused, not wrapped
continue;
}
const size_t cell = (size_t)(tx - D.tx0) * (size_t)D.ny + (size_t)(ty - D.ty0);
g.data[cell] = (float)g.tiles;
g.data.insert(g.data.end(), t.h.begin(), t.h.end());
g.tiles++;
}
g.ok = true;
return g;
}
// The grid's constants. Emitted alongside HBLK_* rather than replacing them:
// HBLK_OX/OZ/MPP still say where texel (0,0) of the GRID sits, so
// hblk_sample_level's world -> texel mapping is untouched.
inline std::string erosion_heightmap_grid_glsl(const HeightGridDesc& g) {
std::string s;
s += "const int HBLK_G_NX = " + std::to_string(g.nx) + ";\n";
s += "const int HBLK_G_NY = " + std::to_string(g.ny) + ";\n";
s += "const int HBLK_G_TX0 = " + std::to_string(g.tx0) + ";\n";
s += "const int HBLK_G_TY0 = " + std::to_string(g.ty0) + ";\n";
s += "const int HBLK_G_PAYLOAD = " + std::to_string(g.payload) + ";\n";
s += "const int HBLK_G_APRON = " + std::to_string(g.apron) + ";\n";
s += "const int HBLK_G_TN = " + std::to_string(g.n) + ";\n";
s += "const int HBLK_G_TABLE = " + std::to_string(g.table) + ";\n";
s += "const float HBLK_G_STUB = " + std::to_string(g.stub) + ";\n";
return s;
}
// ---------------------------------------------------------------------------
// THE C++ TWIN OF THE PAGED FETCH.
//
// The engine meshes on the GPU but answers sample_density / get_distance /
// get_surface_normal / get_capsule_push from C++, so the tile grid needs BOTH.
// test_erosion_parity.cpp exists because those two drifted once and put 4 km of
// air between the surface you saw and the surface you stood on.
//
// So this sits beside the string that emits the GLSL, in one file, and the gate
// compares them on the same points. It is still a transcription -- GLSL text
// against C++ -- but it is a transcription of TWELVE LINES that a probe can
// check, not of a wrapper nobody can reach without booting the engine.
inline float erosion_heightmap_grid_texel(const float* d, const HeightGridDesc& g,
int i, int j) {
// Mirrors hblk_texel_l in erosion_heightmap_sampler_glsl(paged=true).
auto fdiv = [](int a, int b) { return (a >= 0) ? (a / b) : -(((-a) + b - 1) / b); };
const int cy = fdiv(i, g.payload);
const int cx = fdiv(j, g.payload);
if (cx < 0 || cx >= g.nx || cy < 0 || cy >= g.ny) return g.stub;
const int idx = (int)d[(size_t)cx * g.ny + cy];
if (idx < 0) return g.stub;
const int lx = i - cy * g.payload + g.apron;
const int lz = j - cx * g.payload + g.apron;
const size_t base = (size_t)g.table + (size_t)idx * g.n * g.n;
return d[base + (size_t)lz * g.n + lx];
}
// Bilinear, in world metres.
//
// ⚠ BILINEAR, NOT CATMULL-ROM, AND THAT IS THE EXISTING CPU/GPU DIVERGENCE, NOT
// A NEW ONE. Every CPU site in dc_terrain_generator_procnodes.cpp samples the
// heightfield bilinearly today while the GPU runs the 16-tap Catmull-Rom. This
// function reproduces the CPU behaviour exactly so that WIRING the grid changes
// addressing and nothing else -- folding a filter change into it would make any
// resulting difference unattributable. Closing that gap is its own change with
// its own measurement.
inline float erosion_heightmap_grid_sample(const float* d, const HeightGridDesc& g,
float x, float z) {
if (!d || !g.active()) return g.stub;
const float side = (float)(g.side() - 1);
float u = (x - g.ox()) / g.mpp;
float v = (z - g.oz()) / g.mpp;
u = u < 0.0f ? 0.0f : (u > side ? side : u);
v = v < 0.0f ? 0.0f : (v > side ? side : v);
const int i0 = (int)u, j0 = (int)v;
const float fu = u - (float)i0, fv = v - (float)j0;
const float h00 = erosion_heightmap_grid_texel(d, g, i0, j0);
const float h10 = erosion_heightmap_grid_texel(d, g, i0 + 1, j0);
const float h01 = erosion_heightmap_grid_texel(d, g, i0, j0 + 1);
const float h11 = erosion_heightmap_grid_texel(d, g, i0 + 1, j0 + 1);
const float a = h00 + (h10 - h00) * fu;
const float b = h01 + (h11 - h01) * fu;
return a + (b - a) * fv;
}
// The SAMPLER, shared. The consumer declares the storage -- hblk_data[] plus
// HBLK_N / HBLK_OX / HBLK_OZ / HBLK_MPP -- and this reads it.
//
// CATMULL-ROM, NOT BILINEAR, AND AN ANALYTIC GRADIENT. Both matter, and neither
// shows up in a lookdev.
//
// The lookdev samples the field at 30 m/px, the resolution the data is defined
// at, so its samples land exactly ON texel boundaries. dc-terrain samples the
// same field at 1 m voxels -- thirty times finer, through the middle of every
// cell. Bilinear is C0 there: it creases at every texel edge. Measured on the
// shipped block over 240 m at 1 m steps, the surface slope jumped by up to
// 0.1167 in a single step, and dual contouring turns each of those into a
// facet. Catmull-Rom is C1 and brings that to 0.0123, while still passing
// exactly through the samples (0.000000 m error at texel centres).
//
// The gradient is worse than the creasing. heightmap_gn used to take central
// differences at +/-30 m, which is not the gradient of the surface it also
// returns -- measured rms 0.052, max 0.11 apart, about 6 degrees. The erosion
// filter is SLOPE-AWARE and orients gullies to that vector, so it was carving
// channels aligned to a surface that was not the one underneath them. In the
// lookdev this is invisible: at 30 m spacing the central difference IS the
// sampled gradient. Differentiating the interpolant instead gives rms 0.00037.
//
// THE TEXEL FETCH IS THE ONLY THING THE TILE GRID CHANGES. hblk_sample_level's
// 16-tap Catmull-Rom and its analytic gradient are untouched, deliberately:
// they are the measured-correct part (0.0123 slope-jump against bilinear's
// 0.1167; gradient rms 0.00037 against a central difference's 0.052) and
// duplicating them to add addressing is how a "do not touch" invariant dies.
//
// paged = false emits exactly the text every current consumer already gets.
inline std::string erosion_heightmap_sampler_glsl(bool paged = false) {
// ONE TILE: i,j index a single square, out of range CLAMPS.
std::string fetch =
"float hblk_texel_l(int i, int j, int of, int sd) {\n"
" return custom_data[of + clamp(j,0,sd-1)*sd + clamp(i,0,sd-1)];\n"
"}\n"
"float hblk_texel(int i, int j) {\n"
" return custom_data[clamp(j,0,HBLK_N-1)*HBLK_N + clamp(i,0,HBLK_N-1)];\n"
"}\n";
if (paged) {
// THE TILE GRID: i,j are GRID-relative texels -- i along world X (ty),
// j along world Z (tx) -- resolved through the index table.
//
// A tap that leaves its tile lands in the NEIGHBOUR'S payload, not in
// an apron and not in a clamp, because every world texel resolves to
// the tile that OWNS it. That is why this layout needs no gutter.
//
// ⛔ A MISS RETURNS THE STUB. Clamping here would smear a loaded row
// across the hole and read as terrain -- the silent seam.
//
// of/sd are accepted and IGNORED: the grid supplies its own extent.
// The signature is kept so hblk_sample_level compiles unchanged.
fetch =
"int hblk_fdiv(int a, int b) {\n"
" return (a >= 0) ? (a / b) : -(((-a) + b - 1) / b);\n"
"}\n"
"float hblk_texel_l(int i, int j, int of, int sd) {\n"
" int cy = hblk_fdiv(i, HBLK_G_PAYLOAD);\n"
" int cx = hblk_fdiv(j, HBLK_G_PAYLOAD);\n"
" if (cx < 0 || cx >= HBLK_G_NX || cy < 0 || cy >= HBLK_G_NY) return HBLK_G_STUB;\n"
" int idx = int(custom_data[cx * HBLK_G_NY + cy]);\n"
" if (idx < 0) return HBLK_G_STUB;\n"
" int lx = i - cy * HBLK_G_PAYLOAD + HBLK_G_APRON;\n"
" int lz = j - cx * HBLK_G_PAYLOAD + HBLK_G_APRON;\n"
" int base = HBLK_G_TABLE + idx * HBLK_G_TN * HBLK_G_TN;\n"
" return custom_data[base + lz * HBLK_G_TN + lx];\n"
"}\n"
"float hblk_texel(int i, int j) { return hblk_texel_l(i, j, 0, 0); }\n";
}
// custom_data[] is the CONSUMER's buffer and the name is deliberate: the
// engine transports opaque floats and does not know they are elevations.
return fetch +
"float hblk_cr(float p0, float p1, float p2, float p3, float t) {\n"
" return 0.5*((2.0*p1) + (-p0+p2)*t + (2.0*p0-5.0*p1+4.0*p2-p3)*t*t\n"
" + (-p0+3.0*p1-3.0*p2+p3)*t*t*t);\n"
"}\n"
"float hblk_dcr(float p0, float p1, float p2, float p3, float t) {\n"
" return 0.5*((-p0+p2) + 2.0*(2.0*p0-5.0*p1+4.0*p2-p3)*t\n"
" + 3.0*(-p0+3.0*p1-3.0*p2+p3)*t*t);\n"
"}\n"
// returns (height in metres, dh/dx, dh/dz) -- the gradient of the very
// surface returned in .x, by construction rather than by finite difference.
// One level of the pyramid. Level l texel j covers level-0 texels
// [j*2^l, (j+1)*2^l), whose centre sits at j*2^l + (2^l-1)/2 -- the
// level-0 mapping, shifted then scaled. Level 0 reduces to the identity.
"vec3 hblk_sample_level(int l, float x, float z) {\n"
" float sc = float(1 << l);\n"
" int sd = HBLK_SIDE[l];\n"
" int of = HBLK_OFF[l];\n"
" float u = clamp((((x - HBLK_OX) / HBLK_MPP) - (sc-1.0)*0.5) / sc, 0.0, float(sd-1));\n"
" float v = clamp((((z - HBLK_OZ) / HBLK_MPP) - (sc-1.0)*0.5) / sc, 0.0, float(sd-1));\n"
" int i = int(floor(u)); int j = int(floor(v));\n"
" float fu = u - float(i); float fv = v - float(j);\n"
" float r[4]; float d[4];\n"
" for (int k = 0; k < 4; k++) {\n"
" int jj = j + k - 1;\n"
" float p0 = hblk_texel_l(i-1, jj, of, sd); float p1 = hblk_texel_l(i, jj, of, sd);\n"
" float p2 = hblk_texel_l(i+1, jj, of, sd); float p3 = hblk_texel_l(i+2, jj, of, sd);\n"
" r[k] = hblk_cr (p0, p1, p2, p3, fu);\n"
" d[k] = hblk_dcr(p0, p1, p2, p3, fu);\n"
" }\n"
" float pitch = HBLK_MPP * sc;\n"
" return vec3(hblk_cr (r[0],r[1],r[2],r[3], fv),\n"
" hblk_cr (d[0],d[1],d[2],d[3], fv) / pitch,\n"
" hblk_dcr(r[0],r[1],r[2],r[3], fv) / pitch);\n"
"}\n"
// FRACTIONAL, for the reason nyquist_octaves gives for its own taper: an
// octave fades in "rather than snapping it on and moving the surface
// underfoot". A hard mip switch would step the ground at every LOD line.
"vec3 hblk_sample(float x, float z) {\n"
" float lod = clamp(log2(max(voxel_scale, HBLK_MPP) / HBLK_MPP),\n"
" 0.0, float(HBLK_LEVELS - 1));\n"
" int l0 = int(floor(lod));\n"
" float f = lod - float(l0);\n"
" vec3 s0 = hblk_sample_level(l0, x, z);\n"
// ONE LEVEL WHEN THE BLEND IS ZERO, which is every voxel at or below
// the block's own pitch. lod = log2(max(voxel_scale, HBLK_MPP)/HBLK_MPP)
// is EXACTLY 0 for voxel_scale <= 30 m, so f is 0 and the second level
// was being sampled in full and then multiplied by nothing. Each level
// is a 16-tap Catmull-Rom, so that is 32 buffer reads per density
// sample where 16 will do -- and edge_intersect evaluates the density
// six times per edge for the normal, so 192 reads instead of 96.
//
// The mix is kept for the coarse LODs where f is genuinely non-zero;
// this only skips work that provably cannot affect the result.
" if (f <= 0.0) return s0;\n"
" int l1 = min(l0 + 1, HBLK_LEVELS - 1);\n"
" return mix(s0, hblk_sample_level(l1, x, z), f);\n"
"}\n"
"vec3 hblk_sample_unused(float x, float z) {\n"
" float u = clamp((x - HBLK_OX) / HBLK_MPP, 0.0, float(HBLK_N - 1));\n"
" float v = clamp((z - HBLK_OZ) / HBLK_MPP, 0.0, float(HBLK_N - 1));\n"
" int i = int(floor(u)); int j = int(floor(v));\n"
" float fu = u - float(i); float fv = v - float(j);\n"
" float r[4]; float d[4];\n"
" for (int k = 0; k < 4; k++) {\n"
" int jj = j + k - 1;\n"
" float p0 = hblk_texel(i-1, jj); float p1 = hblk_texel(i, jj);\n"
" float p2 = hblk_texel(i+1, jj); float p3 = hblk_texel(i+2, jj);\n"
" r[k] = hblk_cr (p0, p1, p2, p3, fu);\n"
" d[k] = hblk_dcr(p0, p1, p2, p3, fu);\n"
" }\n"
" float h = hblk_cr (r[0], r[1], r[2], r[3], fv);\n"
" float dhdx = hblk_cr (d[0], d[1], d[2], d[3], fv) / HBLK_MPP;\n"
" float dhdz = hblk_dcr(r[0], r[1], r[2], r[3], fv) / HBLK_MPP;\n"
" return vec3(h, dhdx, dhdz);\n"
"}\n";
}
inline std::string erosion_heightmap_gn_glsl() {
return
"vec3 heightmap_gn(float x, float z) {\n"
" float y_range = M_CEIL_Y - M_FLOOR_Y;\n"
" float H = y_range * EN_HAMP * 0.5;\n"
" float mid = M_FLOOR_Y + y_range * 0.5;\n"
// The slope's vertical unit: the FIELD's relief, not the playable band.
// EN_HAMP is folded in because the filter re-applies it to all three
// components, so this is what makes the slope arrive at field scale.
" float Hs = (EN_FIELDH > 0.0 ? EN_FIELDH : y_range * 0.5) * EN_HAMP;\n"
" vec3 s = hblk_sample(x, z);\n"
" return vec3((s.x - mid) / H, s.y * M_TH * EN_ESCALE / Hs, s.z * M_TH * EN_ESCALE / Hs);\n"
"}\n";
}
// The C++ TWIN of heightmap_gn, for the CPU evaluate path.
//
// It has to exist because collision samples the CPU field (DCTerrain::
// sample_density) while rendering samples the GPU one: if the two disagree the
// player falls through ground that is visibly there. It is a second expression
// of the same algebra, which is a duplication -- but the alternative, leaving
// the CPU on the fBm mountain while the GPU reads the model, is not a
// duplication, it is two different worlds.
//
// Kept adjacent to the GLSL it mirrors, and to the derivation above, so the
// "MUST stay identical" pairing is visible in one file rather than across two
// repositories. hblk_at is the caller's bilinear sampler in METRES.
template <typename SampleFn>
inline void erosion_heightmap_gn_cpu(float x, float z,
float floor_y, float ceil_y, float terrain_height,
float height_amp, float erosion_scale, float field_half_range,
float step_m, SampleFn hblk_at,
float& out_h, float& out_dx, float& out_dz) {
const float y_range = ceil_y - floor_y;
const float H = y_range * height_amp * 0.5f;
const float mid = floor_y + y_range * 0.5f;
const float e = step_m;
const float h = hblk_at(x, z);
const float dx = (hblk_at(x + e, z) - hblk_at(x - e, z)) / (2.0f * e);
const float dz = (hblk_at(x, z + e) - hblk_at(x, z - e)) / (2.0f * e);
out_h = (h - mid) / H;
// MUST match heightmap_gn above, including the erosion_scale factor: the
// player falls through ground that is visibly there if these diverge.
const float Hs = (field_half_range > 0.0f ? field_half_range
: y_range * 0.5f) * height_amp;
out_dx = dx * terrain_height * erosion_scale / Hs;
out_dz = dz * terrain_height * erosion_scale / Hs;
}
inline std::string erosion_wrapper_glsl(const ErosionWrapperParams& p) {
std::string s = erosion_wrapper_consts_glsl(p);
if (p.heightmap_source) s += erosion_heightmap_sampler_glsl() + erosion_heightmap_gn_glsl();
// The two lines that choose where the base field comes from. Spelled out
// as a variable rather than inline because the wrapper bodies below stay
// SEPARATELY written out -- the header's note that emitting one and
// aliasing the other would halve what the oracle catches applies to the
// wrapper functions, not to this shared fragment of both.
const std::string gn_block = p.heightmap_source
? std::string(
" vec3 gn = heightmap_gn(x, z);\n"
" vec2 en = erosion_node_eval_full(p, gn, EN_HS, EN_STR,\n")
: std::string(
" vec3 gn = fractal_noise_2d(p, GN_SCALE, GN_DETAIL, 2.0, GN_GAIN);\n"
" vec2 en = erosion_node_eval_full(p, gn + terrain_offset(x,z), EN_HS, EN_STR,\n");
// IN HEIGHTMAP MODE, EMIT NEITHER NOISE LAYER.
//
// gpu_wrapper does not call terrain_offset when the field comes from a
// heightmap, and detail_offset early-returns on M_DSCALE <= 0 -- so both were
// dead. Emitting them anyway put a multifractal mountain and a 4-octave fBm
// in the shader text of a map whose whole point is that it has neither, which
// is indefensible on a read-through even if the compiler strips it. The
// lookdev's terrain is the heightmap and the erosion filter; the shader
// should say so.
if (p.heightmap_source) {
s +=
"float detail_offset(float x, float z) { return 0.0; }\n";
} else {
s +=
"vec3 terrain_offset(float x, float z) {\n"
" float inv_scale = 1.0 / M_TH;\n"
" float mtn_lim = nyquist_octaves(M_TH / M_MSCALE, M_VOXEL, M_LAC);\n"
" vec3 mtn = (M_MMODE == 1)\n"
" ? fractal_gradient_2d_lod(vec2(x,z)*inv_scale*M_MSCALE, 1.0, M_MOCT, M_LAC, M_MGAIN, M_MDAMP, mtn_lim)\n"
" : fractal_multifractal_2d_lod(vec2(x,z)*inv_scale*M_MSCALE, 1.0, M_MOCT, M_LAC, M_MGAIN, M_MVAR, mtn_lim);\n"
" float relief_amp = 1.0;\n"
" if (M_RSCALE > 0.0) {\n"
" vec3 rel = fractal_noise_2d_lod(vec2(x,z)*inv_scale*M_RSCALE, 1.0, 2, 2.0, 0.5, nyquist_octaves(M_TH/M_RSCALE, M_VOXEL, 2.0));\n"
" float t = clamp(rel.x*0.5+0.5, 0.0, 1.0);\n"
" relief_amp = M_RFLOOR + (1.0-M_RFLOOR)*t;\n"
" }\n"
" float amp = M_MHEIGHT * relief_amp;\n"
" return vec3(mtn.x, mtn.y*M_MSCALE, mtn.z*M_MSCALE) * (amp * 8.0);\n"
"}\n"
"float detail_offset(float x, float z) {\n"
" if (M_DSCALE <= 0.0) return 0.0;\n"
" float inv_scale = 1.0 / M_TH;\n"
" float y_range = M_CEIL_Y - M_FLOOR_Y;\n"
" vec3 det = fractal_noise_2d_lod(vec2(x,z)*inv_scale*M_DSCALE, 1.0, M_DOCT, M_LAC, M_DGAIN, nyquist_octaves(M_TH/M_DSCALE, M_VOXEL, M_LAC));\n"
" return det.x * (M_DHEIGHT / y_range);\n"
"}\n"
"";
}
s +=
"float y_mapped_of(float y, float det_off) {\n"
" float y_range = M_CEIL_Y - M_FLOOR_Y;\n"
" return (y - M_FLOOR_Y) * (EN_HS / y_range) - det_off * EN_HS;\n"
"}\n"
"float clip_d(float d, float y) {\n"
" d = max(d, (M_FLOOR_Y - y) * 2.0);\n"
" d = min(d, (M_CEIL_Y - y) * 2.0);\n"
" return d;\n"
"}\n"
"float gpu_wrapper(float x, float y, float z) {\n"
" float inv_scale = 1.0 / M_TH;\n"
" float y_range = M_CEIL_Y - M_FLOOR_Y;\n"
" float y_mapped = y_mapped_of(y, detail_offset(x, z));\n"
" vec2 p = vec2(x*inv_scale, z*inv_scale);\n";
s += gn_block;
s +=
" EN_CS, EN_OCT, EN_GAIN, EN_FADEDIV, EN_ESCALE,\n"
" EN_HAMP, EN_ELAC, EN_GW, EN_DET, EN_NORM, EN_HOFF,\n"
// 999 = no fractional taper, which is bit-identical to the arithmetic
// before max_octaves existed. The lookdev renders one fixed grid and has
// no LOD to taper against; the engine is where that parameter earns its
// keep.
" EN_ROUND, EN_ONSET, EN_ASLOPE, 999.0);\n"
// NODE -> WRAPPER egress, matching ErosionNoiseNode's emitted snippet.
// This wrapper calls erosion_node_eval_full directly instead of going
// through the node, so without this line the lookdev and standalone
// paths would compute ridge and drop it while the engine path exports
// it -- the two-paths-disagree failure this file's header warns about.
// Assigns a global and nothing else; d is untouched.
" pn_out_ridge = en.y;\n"
" float d = en.x - y_mapped;\n"
" d *= (y_range / EN_HS);\n"
" return clip_d(d, y);\n"
"}\n"
"float cpu_wrapper(float x, float y, float z) {\n"
" float inv_scale = 1.0 / M_TH;\n"
" float y_range = M_CEIL_Y - M_FLOOR_Y;\n"
" float y_mapped = y_mapped_of(y, detail_offset(x, z));\n"
" vec2 p = vec2(x*inv_scale, z*inv_scale);\n";
s += gn_block;
s +=
" EN_CS, EN_OCT, EN_GAIN, EN_FADEDIV, EN_ESCALE,\n"
" EN_HAMP, EN_ELAC, EN_GW, EN_DET, EN_NORM, EN_HOFF,\n"
// 999 = no fractional taper, which is bit-identical to the arithmetic
// before max_octaves existed. The lookdev renders one fixed grid and has
// no LOD to taper against; the engine is where that parameter earns its
// keep.
" EN_ROUND, EN_ONSET, EN_ASLOPE, 999.0);\n"
// NODE -> WRAPPER egress, matching ErosionNoiseNode's emitted snippet.
// This wrapper calls erosion_node_eval_full directly instead of going
// through the node, so without this line the lookdev and standalone
// paths would compute ridge and drop it while the engine path exports
// it -- the two-paths-disagree failure this file's header warns about.
// Assigns a global and nothing else; d is untouched.
" pn_out_ridge = en.y;\n"
" float d = en.x - y_mapped;\n"
" d *= (y_range / EN_HS);\n"
" return clip_d(d, y);\n"
"}\n";
return s;
}