// GENERATED by tools/gen_density_shader.py from /tmp/erosion_grad.glsl
// DO NOT EDIT -- regenerate. That script documents what is transformed and why.
//
// This is the SAME TEXT the GPU mesher compiles for the density field, so the
// material cannot drift from the geometry: there is one source, and it is
// DCTerrainGeneratorProcNodes::get_compiled_glsl().
shader_type spatial;
render_mode unshaded;
#include "res://shaders/terrain_material.gdshaderinc"
uniform float force_voxel_scale = 0.0;
uniform int probe_taps = 4;
// PREVIEW. >0 maps UV across a world patch of this many metres, so the probe
// shows actual terrain instead of the two-metre cube it is painted on. Without
// it the probe samples a 4 m patch of a 60 km world -- one flat spot, one flat
// colour, which is correct output and a useless picture.
uniform float probe_world_span = 0.0;
uniform vec2 probe_world_origin = vec2(0.0);
uniform float probe_world_y = 0.0;
// 0 = normal as RGB (what a gate wants). 1 = lambert against a fixed sun,
// which is what actually reads as relief to a human eye.
uniform int probe_view = 1;
uniform vec3 probe_sun = vec3(0.45, 0.55, 0.7);
// Surface elevation band, for the material's altitude term. Defaults match
// map_dc_erosion.gd's SURFACE_RANGE.
uniform vec2 probe_field_range = vec2(-7800.0, 10600.0);
uniform float probe_curv_gain = 0.5;
// Octave limit â fixed high value for LOD-consistent surface position.
// Nyquist clamping shifts the erosion surface between LOD levels, causing
// chunks to disappear during refinement. Consistent octaves = consistent surface.
// =========================================================================
// erosion_filter.inc â THE erosion filter. One copy. Compiled BOTH as GLSL
// and as C++ from this same text.
//
// WHY THIS FILE EXISTS.
//
// This algorithm existed four times: the GLSL library string in
// proc-nodes/src/gpu/gpu_compiler.cpp, its C++ twin in
// proc-nodes/src/nodes/noise_nodes.cpp, a THIRD C++ copy in
// dc-terrain/src/dc_terrain_generator_procnodes.cpp, and a fourth in
// hyper-tps/shaders/erosion_terrain.gdshader. They were not kept in step.
// The GPU wrapper called it with scale=0.08333 and the dc-terrain CPU
// wrapper with scale=5.0 â a 60x difference in both gully frequency and
// gully amplitude â so the surface the mesher drew and the surface
// sample_density/collision reported were different functions. Nothing
// caught it because nothing compared them.
//
// This is dc_shader_source.h's invariant one level down:
//
// The GPU and the CPU evaluate the SAME TEXT.
//
// Not two implementations kept in step by review. One file, included by
// both. A divergence is now a compile error or an oracle failure, not a
// thing you find eighteen months later in a collision bug.
//
// HOW THE DUAL COMPILATION WORKS.
//
// The body below is written in the subset GLSL and C++ agree on: vec2/vec3
// constructors, member access, arithmetic operators, and the builtins
// floor/fract/dot/max/min/exp/cos/sin/pow. In a shader this text is
// injected verbatim and every one of those is a GLSL builtin. In C++ it is
// #included after erosion_cpp_shim.h, which supplies the same names with
// the same componentwise semantics. There is no #ifdef in here on purpose âÂÂ
// a branch is a place for the two languages to disagree.
//
// LICENSING -- THIS FILE IS MPL-2.0. It was MIT until the advanced filter
// landed; erosion_filter and phacelle_noise are runevision's rewrite and
// carry MPL-2.0, so this file does too.
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this file,
// You can obtain one at https://mozilla.org/MPL/2.0/.
//
// Advanced Terrain Erosion Filter + Phacelle Noise:
// (c) 2025 Rune Skovbo Johansen, MPL-2.0
// https://www.shadertoy.com/view/wXcfWn
//
// MPL-2.0 is FILE-level copyleft, which is part of why the algorithm lives
// in its own file rather than pasted into gpu_compiler.cpp: the obligation
// travels with these functions. The generated erosion_filter_glsl.gen.h is
// this text verbatim and carries the same terms.
//
// Still MIT, and still present:
//
// Gradient noise + derivatives (grad_noised_2d): (c) 2017 Inigo Quilez,
// MIT, https://www.shadertoy.com/view/XdXBRH
// Eroded terrain noise: (c) 2020 Clay John, MIT,
// https://www.shadertoy.com/view/MtGcWh
// Terrain erosion noise: (c) 2023 Fewes, MIT,
// https://www.shadertoy.com/view/7ljcRW
// Clean terrain erosion filter:(c) 2025 Rune Skovbo Johansen, MIT,
// https://www.shadertoy.com/view/33cXW8
//
// PROVENANCE. The body started as a byte-identical extraction of
// gpu_compiler.cpp lines 378-490 (verified by diff, not by eye) and has
// since changed in two places, each with its reasoning written at the
// function it changed:
//
// erosion_hash float trick -> integer lattice hash, because the float
// one was unstable at the coordinates this project uses,
// which made GPU/CPU parity impossible in principle.
// gullies_2d and the frequency approach -> phacelle_noise and the fade
// erosion_filter approach: runevision's 2026 rewrite.
//
// grad_noised_2d and fractal_noise_2d are still the extracted text.
// =========================================================================
// =========================================================================
// Erosion terrain functions â iq/runevision/clayjohn
// =========================================================================
// Hash for erosion noise (same constants as iq/runevision)
// Integer lattice hash. Both call sites -- grad_noised_2d's floor(p) corners
// and gullies_2d's p_int + grid_offset -- pass EXACT integers, so this hashes
// the integer directly and is bit-identical on every implementation that has
// 32-bit unsigned arithmetic. GLSL and C++ both define unsigned wraparound.
//
// WHAT IT REPLACED, AND WHY. The previous hash was iq's float trick:
//
// inner = fract(x.x * x.y * (x.x + x.y)); return -1 + 2*fract(16*k*inner);
//
// That product passes 2^24 quickly -- at lattice coord 100 it is already
// ~3.5e5 with float32 ULP coarse enough to matter, and by 1000 it is ~1e8 --
// after which fract() has no resolution left. The failure is not a gentle
// loss of quality: fract is a step function, so a 1-ULP difference anywhere
// upstream lands on the far side of an integer boundary and the hash returns
// an UNRELATED gradient. Measured against a float64 reference, the fraction
// of lattice points whose hash flips this way:
//
// lattice coord 1 10 100 1000 1e6
// flipped 0.54% 0.83% 83.86% 99.49% noise dies entirely
//
// (At 1e6 the product exceeds float32 entirely, fract returns 0 everywhere,
// and the hash degenerates to a constant -- the terrain goes flat.)
//
// map_dc_erosion's detail layer runs at scale 200, which puts it at lattice
// coord 100 at x=10km and ~654 by its fourth octave: squarely in the >84%
// regime. That is why the GPU and the CPU disagreed by 23 m on IDENTICAL
// source, and why single-sourcing alone could never have made them agree.
// Measured: turning the detail layer off dropped worst |glsl - c++| from
// 23.2251 m to 0.0171 m. See dc-terrain/tests/test_erosion_parity.cpp.
//
// The new ceiling is 2^24 (16.7M) -- the point where a float stops
// representing consecutive integers at all, so floor() itself would fail
// first. At 1 m voxels that is 16,777 km, against a 57.3 km world.
//
// The mix is pcg2d (Jarzynski & Olano, "Hash Functions for GPU Rendering",
// JCGT 2020), written with scalar uints rather than uvec2 because GLSL
// `inout` and C++ references are not the same thing and this file must
// compile as both. Output is masked to 24 bits before the float conversion:
// 2^24-1 is exactly representable in float32 and the scale is a power of
// two, so the conversion is exact rather than merely round-to-nearest on
// both sides.
//
// THIS CHANGES EVERY SEED. The world it generates is a different world from
// the one the float hash generated. That was the accepted cost of making the
// meshed surface and the collidable surface the same surface.
vec2 erosion_hash(vec2 x) {
uint vx = uint(int(x.x));
uint vy = uint(int(x.y));
vx = vx * 1664525u + 1013904223u;
vy = vy * 1664525u + 1013904223u;
vx += vy * 1664525u;
vy += vx * 1664525u;
vx = vx ^ (vx >> 16u);
vy = vy ^ (vy >> 16u);
vx += vy * 1664525u;
vy += vx * 1664525u;
vx = vx ^ (vx >> 16u);
vy = vy ^ (vy >> 16u);
return vec2(
-1.0 + 2.0 * (float(vx >> 8u) * (1.0 / 16777216.0)),
-1.0 + 2.0 * (float(vy >> 8u) * (1.0 / 16777216.0)));
}
// iq's gradient noise with analytical derivatives (2D)
// Returns vec3(height, dh/dx, dh/dy)
vec3 grad_noised_2d(vec2 p) {
vec2 i = floor(p);
vec2 f = p - i;
vec2 u = f*f*f*(f*(f*6.0-15.0)+10.0);
vec2 du = 30.0*f*f*(f*(f-2.0)+1.0);
vec2 ga = erosion_hash(vec2(i.x, i.y));
vec2 gb = erosion_hash(vec2(i.x+1.0, i.y));
vec2 gc = erosion_hash(vec2(i.x, i.y+1.0));
vec2 gd = erosion_hash(vec2(i.x+1.0, i.y+1.0));
float va = dot(ga, f);
float vb = dot(gb, f - vec2(1.0, 0.0));
float vc = dot(gc, f - vec2(0.0, 1.0));
float vd = dot(gd, f - vec2(1.0, 1.0));
float height = va + u.x*(vb-va) + u.y*(vc-va) + u.x*u.y*(va-vb-vc+vd);
vec2 deriv = ga + (gb-ga)*u.x + (gc-ga)*u.y + (ga-gb-gc+gd)*u.x*u.y +
vec2(du.x * (u.y*(va-vb-vc+vd) + vb-va),
du.y * (u.x*(va-vb-vc+vd) + vc-va));
return vec3(height, deriv.x, deriv.y);
}
// fBM of grad_noised_2d with derivative accumulation
// NYQUIST OCTAVE LIMIT, as a FRACTIONAL count.
//
// An octave whose wavelength falls below about twice the voxel size cannot be
// represented by the mesher: it aliases, and in a dual-contouring mesher that
// shows up as a field of thin vertical spikes. This generator had no protection
// at all -- the wrapper's mountain/relief/detail layers passed literal octave
// counts, and the graph nodes' min(N, pn_max_octaves) was disarmed by setting
// pn_max_octaves = 99. The aliasing is therefore always present and merely
// below notice at low mountain_scale; raising mountain_scale shortens every
// mountain octave together and drives the fine ones under the voxel size.
//
// Returns how many octaves survive, FRACTIONALLY. That matters: a hard integer
// cutoff makes an octave appear all at once when a chunk refines, which moves
// the surface and pops geometry between LOD levels -- the documented reason
// octave clamping was turned off here in the first place. A fractional limit
// lets the caller fade the last octave in and out continuously instead.
//
// voxel <= 0 means "caller does not know / does not care" and disables the
// limit. DCTerrain::sample_density passes 0 deliberately (dc_terrain.cpp:2841),
// so collision always evaluates the full field.
float nyquist_octaves(float base_wavelength, float voxel, float lac) {
if (voxel <= 0.0 || lac <= 1.0) {
return 999.0;
}
float ratio = base_wavelength / max(2.0 * voxel, 1e-6);
if (ratio <= 1.0) {
return 0.0;
}
return log(ratio) / log(lac);
}
// ELEVATION-WEIGHTED fBm -- a multifractal, and the answer to "all the mountain
// ranges look the same".
//
// Plain fBm gives every location the same amount of every octave, so the terrain
// is statistically HOMOGENEOUS: equally rugged in the basins and on the peaks.
// That sameness is what reads as "no variance between ranges" no matter how the
// octave count or persistence are tuned -- those change how rough it is
// everywhere, not where the roughness goes.
//
// A multifractal makes each octave's amplitude depend on the terrain already
// accumulated, so high ground becomes rugged and low ground stays smooth. Real
// ranges do this: jagged peaks above smooth valley floors.
//
// WHY NOT MUSGRAVE'S HETERO TERRAIN, which is the textbook choice and is already
// in this library as perlin_hetero: it MULTIPLIES the running value each octave,
// so its output range depends on offset, roughness and octave count together and
// grows without a fixed bound. Blending that against fBm would need a
// normalisation measured per parameter set, and an unbounded field one slider
// away from a user is a NaN waiting to happen. Weighting the AMPLITUDE instead
// keeps fBm's envelope exactly, so `variance` is a safe 0..1 dial and 0 is
// bit-identical to fractal_noise_2d.
//
// THE DERIVATIVE IS APPROXIMATE. d(weight)/dp is dropped -- only the noise term
// is differentiated. That is deliberate and it is safe HERE because this
// gradient is consumed for gully DIRECTION, which erosion_filter normalises
// (normalize_safe) and then largely discards anyway: assumed_slope mixes at 1.0,
// so the real slope survives only through slope_length in the initial mask. Do
// not reuse this derivative anywhere that needs it to be exact.
// iq's GRADIENT TRICK, from "fbm with derivatives"
// (https://iquilezles.org/articles/morenoise/). Each octave is damped by the
// gradient ACCUMULATED so far, so steep ground loses its fine octaves while
// flats keep them -- a cheap stand-in for water running faster on slopes.
//
// This is NOT what fractal_multifractal_2d_lod below does. That one weights by
// accumulated ELEVATION (Musgrave: high ground gets more detail). Same shape of
// idea, different axis, and the two are worth having side by side.
//
// Following iq exactly, the accumulator is the RAW per-octave derivative, with
// neither the amplitude nor the frequency folded in, and the weight is
// 1/(1 + damping * dot(d, d)). damping = 0 reduces this to plain fBm.
//
// CAVEAT, and it is the same one the multifractal already carries: the returned
// .yz is no longer the exact gradient of the returned .x, because the weight is
// itself position-dependent and its own derivative is not included. The erosion
// filter uses the input slope only for gully DIRECTION and replaces its
// magnitude (assumed_slope), so this matters less here than it would elsewhere.
vec3 fractal_gradient_2d_lod(vec2 p, float freq, int octaves, float lac,
float gain, float damping, float max_octaves) {
vec3 n = vec3(0.0);
vec2 dsum = vec2(0.0);
float nf = freq;
float na = 1.0;
for (int i = 0; i < octaves; i++) {
float lw = clamp(max_octaves - float(i), 0.0, 1.0);
if (lw <= 0.0) {
na *= gain;
nf *= lac;
continue;
}
vec3 d = grad_noised_2d(p * nf);
dsum += vec2(d.y, d.z);
float w = 1.0 / (1.0 + damping * dot(dsum, dsum));
n += vec3(d.x * na * w * lw, d.y * na * nf * w * lw, d.z * na * nf * w * lw);
na *= gain;
nf *= lac;
}
return n;
}
vec3 fractal_multifractal_2d_lod(vec2 p, float freq, int octaves, float lac,
float gain, float variance, float max_octaves) {
vec3 n = vec3(0.0);
float nf = freq;
float na = 1.0;
for (int i = 0; i < octaves; i++) {
float lw = clamp(max_octaves - float(i), 0.0, 1.0);
if (lw <= 0.0) {
na *= gain;
nf *= lac;
continue;
}
// Weight from the elevation SO FAR. The 3.0 maps a typical fBm value
// range of about +/-0.33 onto 0..1; below that the weight floors and the
// octave is suppressed, which is the point.
float w = mix(1.0, clamp(0.5 + n.x * 3.0, 0.0, 1.0), variance);
vec3 d = grad_noised_2d(p * nf);
n += vec3(d.x * na * w * lw, d.y * na * nf * w * lw, d.z * na * nf * w * lw);
na *= gain;
nf *= lac;
}
return n;
}
vec3 fractal_multifractal_2d(vec2 p, float freq, int octaves, float lac,
float gain, float variance) {
return fractal_multifractal_2d_lod(p, freq, octaves, lac, gain, variance, 999.0);
}
vec3 fractal_noise_2d_lod(vec2 p, float freq, int octaves, float lac, float gain,
float max_octaves) {
vec3 n = vec3(0.0);
float nf = freq;
float na = 1.0;
for (int i = 0; i < octaves; i++) {
// Fractional taper: full weight while the octave is comfortably above
// Nyquist, fading to zero as it crosses. Continuous in max_octaves, so a
// chunk refining to a finer LOD grows the octave in rather than snapping
// it on and moving the surface under the player.
float w = clamp(max_octaves - float(i), 0.0, 1.0);
if (w > 0.0) {
vec3 d = grad_noised_2d(p * nf);
n += vec3(d.x * na * w, d.y * na * nf * w, d.z * na * nf * w);
}
na *= gain;
nf *= lac;
}
return n;
}
// Unlimited fBm. Kept so existing callers -- GradientNoiseNode's emitted GLSL,
// ErosionNoiseNode's mesh path -- are untouched by the LOD work.
vec3 fractal_noise_2d(vec2 p, float freq, int octaves, float lac, float gain) {
return fractal_noise_2d_lod(p, freq, octaves, lac, gain, 999.0);
}
// -------------------------------------------------------------------------
// THE ADVANCED FILTER â runevision's 2026 rewrite, MPL-2.0.
//
// What was here before was the FREQUENCY APPROACH: gullies_2d took an
// unnormalised slope and built its stripe direction as
//
// side_dir = vec2(-slope.y, slope.x) * 2.0 * PI
//
// so stripe FREQUENCY was proportional to slope magnitude, and at zero slope
// the stripes became infinitely thick. That is the 2018 clayjohn / 2023 Fewes
// lineage. It has two consequences runevision documents and we had both:
// valleys get bulges at the bottom and cannot be made V-shaped, and the only
// way to hide it is to fade erosion out at low altitude -- which is exactly
// what the WATER_HEIGHT smoothstep in erosion_node_eval was doing. We
// inherited a workaround for a defect and shipped it as a feature.
//
// The fade approach keeps stripe width CONSTANT (phacelle_noise takes a
// NORMALISED direction plus a separate freq) and instead fades the stripes
// towards a caller-supplied fade target as steepness approaches zero. That
// single change is what makes the rest possible:
//
// - normalized gullies: because side_dir is applied AFTER interpolation
// rather than inside the 16-cell loop, the interpolated (cos, sin) pair
// can be treated as a point on a circle and renormalised, giving gullies
// a consistent magnitude and therefore crisp ridges. Partial (k=2, i.e.
// normalization 0.5) because full normalisation produces spiky loop
// artefacts where the waves cancel.
// - straight gullies: gully_slope accumulates sign(ph.y) rather than ph.y,
// as if the gullies were extruded triangle waves, so smaller gullies
// branch off larger ones cleanly instead of curling along them. The
// OUTPUT derivative keeps the true sine; only the internal one is faked,
// which is why there are two slope accumulators.
// - stacked fading: each octave's mask multiplies into a combi-mask
// (combi_mask = pow_inv(combi_mask, detail) * new_mask) and its faded
// gullies become the next octave's fade target, so large ridges survive
// the small octaves instead of dissolving into them.
// - the ridge map: falls out for free, tracked in parallel without the
// rounding and gully-weight terms that would otherwise smear it.
//
// LICENSING. erosion_filter and phacelle_noise below are
// Advanced Terrain Erosion Filter / Phacelle Noise
// (c) 2025 Rune Skovbo Johansen, MPL-2.0
// https://www.shadertoy.com/view/wXcfWn
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this file,
// You can obtain one at https://mozilla.org/MPL/2.0/.
//
// Translated with reference to Luke Mitchell's burstable C# port
// (github.com/lpmitchell/AdvancedTerrainErosion, MPL-2.0 for these two
// functions). grad_noised_2d remains iq's, MIT; see the header.
//
// NO `out` PARAMETERS AND NO CONSTRUCTORS BELOW. GLSL `out` is not a C++
// reference and GLSL's positional struct constructor is not C++ aggregate
// initialisation, so the result comes back as a struct whose fields are
// assigned one at a time. That is the subset both languages agree on.
// -------------------------------------------------------------------------
float clamp01(float t) {
return clamp(t, 0.0, 1.0);
}
// Flip, raise to the power, flip back. Controls how aggressively higher
// frequency gullies are restricted to steeper ground: lower = more detail.
float pow_inv(float t, float power) {
return 1.0 - pow(1.0 - clamp01(t), power);
}
// Mirror of sqrt around the diagonal. The obvious shaping function here is
// pow(slope, 0.5), and it produces visible folds at peaks and valleys because
// the sqrt curve starts VERTICAL -- the tiniest slant already gets near-full
// erosion. This starts gently and mostly removes them.
float ease_out(float t) {
float v = 1.0 - clamp01(t);
return 1.0 - v * v;
}
// Variable-size ease-in chained onto the shaping function; this is what
// "rounding" of ridges and creases actually is. Quadratic below `smoothing`,
// linear above, meeting with matching value and slope.
//
// The division is safe for the values this file passes: t is always a
// magnitude (>= 0) and the t >= smoothing branch catches smoothing == 0.
float smooth_start(float t, float smoothing) {
if (t >= smoothing) {
return t - 0.5 * smoothing;
}
return 0.5 * t * t / smoothing;
}
// normalize() is undefined at zero length in GLSL. gully_slope genuinely can
// reach zero (flat ground, first octave), so it needs the guard.
vec2 normalize_safe(vec2 v) {
float l = length(v);
if (l > 1e-20) {
return v / l;
}
return vec2(0.0);
}
// PHACELLE NOISE â phase + cell. A stripe pattern aligned to norm_dir, built
// by interpolating cosine and sine waves from a 4x4 neighbourhood of randomly
// offset cell points.
//
// norm_dir MUST be normalised. This is the whole difference from the
// old gullies_2d, which took a raw slope and let its magnitude
// set the frequency.
// freq stripes per cell. Keep near 1; large values distort.
// offset phase offset, 1.0 = a full cycle.
// normalization degree of renormalisation, 0..1. At 0.5, raw magnitudes
// below 0.5 are left alone and everything above is pulled to
// 1 -- the k=2 partial normalisation that gives consistent
// gully magnitude without the loop artefacts full
// normalisation produces where the waves cancel exactly.
//
// Returns (cos, sin) normalised, plus side_dir, which the caller multiplies
// onto the sine to get the cosine's derivatives. Postponing that multiply out
// of the loop is what makes the normalisation possible at all: side_dir is
// constant across all 16 samples, so it can be factored out.
vec4 phacelle_noise(vec2 p, vec2 norm_dir, float freq, float offset, float normalization) {
// PN_TAU, not TAU: Godot shading language RESERVES TAU as a builtin, and
// this .inc is compiled as a fog shader as well as GLSL 450 and C++.
// Measured: `SHADER ERROR: Redefinition of 'TAU'` from a FogVolume.
const float PN_TAU = 6.28318530717958647;
vec2 side_dir = vec2(-norm_dir.y, norm_dir.x) * freq * PN_TAU;
offset *= PN_TAU;
vec2 p_int = floor(p);
vec2 p_frac = p - p_int;
vec2 phase_dir = vec2(0.0);
float weight_sum = 0.0;
for (int i = -1; i <= 2; i++) {
for (int j = -1; j <= 2; j++) {
vec2 grid_offset = vec2(float(i), float(j));
vec2 grid_point = p_int + grid_offset;
vec2 random_offset = erosion_hash(grid_point) * 0.5;
vec2 v = p_frac - grid_offset - random_offset;
// Bell weight, 1 at distance 0 and ~0 at 1.5. The 0.01111 makes it
// exactly 0 there, which removes faint grid-line artefacts: with
// random offsets up to 0.5, 1.5 is the closest a cell point OUTSIDE
// the 4x4 neighbourhood can be.
float sqr_dist = dot(v, v);
float weight = max(0.0, exp(-sqr_dist * 2.0) - 0.01111);
weight_sum += weight;
float wave_input = dot(v, side_dir) + offset;
phase_dir += vec2(cos(wave_input), sin(wave_input)) * weight;
}
}
vec2 interpolated = phase_dir / max(weight_sum, 1e-6);
float magnitude = sqrt(dot(interpolated, interpolated));
magnitude = max(1.0 - normalization, magnitude);
return vec4(interpolated.x / magnitude, interpolated.y / magnitude,
side_dir.x, side_dir.y);
}
struct ErosionResult {
vec3 delta; // height and slope delta, to add to the input
float magnitude; // summed per-octave strength; scales the height offset
float ridge; // -1 in creases, +1 on ridges. Free; drives materials.
};
// rounding = (ridges, creases, input scale, per-octave multiplier)
// onset = (input, per-octave, ridge-map input, ridge-map per-octave)
// assumed_slope = (pretend magnitude, how much to mix it in)
//
// ASSUMED SLOPE deserves a note. The input height function's slope goes to
// zero at its own peaks and valleys, where the eroded result should be just as
// steep as anywhere else, so using the real slope makes the first octave
// dominate exactly where it should not. runevision tried several fixes and
// settled on pretending the input has a fixed slope everywhere. With the
// reference mix of 1.0 the real slope is not used for gully direction at all,
// which has a useful side effect here: this filter is insensitive to the
// slope-magnitude mismatches that the CPU/GPU divergence was made of.
// max_octaves is FRACTIONAL, and it is the only way to vary octave count with
// LOD without stepping the ground.
//
// Truncating the loop at an integer moves the surface bodily: measured on
// map_dc_diffusion's settings, 7 -> 5 octaves shifts the mean height 20.4 m and
// 7 -> 4 shifts it 48.1 m. Two things cause that and they partly cancel, which
// is why the obvious fix makes it worse -- pinning `magnitude` to its limit
// takes 7 -> 4 from +48 m to -280 m. The octaves are not a truncatable series:
// each one sets `fade_target` and multiplies into `combi_mask` for the next, so
// dropping one changes where every finer octave is allowed to act.
//
// The taper is the same one fractal_noise_2d_lod already uses for the wrapper's
// fBm layers: weight the last partial octave by its fraction rather than
// including or excluding it whole. Everything that accumulates -- the height,
// the slope, the magnitude that drives the offset, the gully steering -- is
// weighted together, so an octave fades in continuously.
//
// 999 (the default every existing caller passes) makes every weight exactly 1.0
// and the arithmetic bit-identical to before the parameter existed.
ErosionResult erosion_filter(vec2 p, vec3 height_and_slope, float fade_target,
float strength, float gully_weight, float detail,
vec4 rounding, vec4 onset, vec2 assumed_slope,
float scale, int octaves, float lacunarity,
float gain, float cell_scale, float normalization,
float max_octaves) {
strength *= scale;
fade_target = clamp(fade_target, -1.0, 1.0);
vec3 input_has = height_and_slope;
float freq = 1.0 / (scale * cell_scale);
vec2 slope2 = vec2(height_and_slope.y, height_and_slope.z);
float slope_length = max(length(slope2), 1e-10);
float magnitude = 0.0;
float rounding_mult = 1.0;
float rounding_for_input = mix(rounding.y, rounding.x, clamp01(fade_target + 0.5)) * rounding.z;
// The combined accumulating mask: initially the input slope, then each
// octave's own slope multiplied in.
float combi_mask = ease_out(smooth_start(slope_length * onset.x, rounding_for_input * onset.x));
// The ridge map runs the same machinery WITHOUT rounding or gully weight,
// because both of those smear exactly the lines it is trying to draw.
float ridge_combi_mask = ease_out(slope_length * onset.z);
float ridge_fade_target = fade_target;
vec2 gully_slope = mix(slope2, slope2 / slope_length * assumed_slope.x, assumed_slope.y);
for (int i = 0; i < octaves; i++) {
vec4 ph = phacelle_noise(p * freq, normalize_safe(gully_slope), cell_scale, 0.25, normalization);
// p was multiplied by freq, so the derivative is too. Negated because
// these are slope directions that point DOWNhill.
vec2 ph_zw = vec2(ph.z, ph.w) * (-freq);
// The octave's own weight: 1 for a full octave, its fraction for the
// last partial one, 0 beyond. Applied to everything this iteration
// contributes, so the octave arrives continuously rather than at once.
float ow = clamp(max_octaves - float(i), 0.0, 1.0);
float sw = strength * ow;
float sloping = abs(ph.y);
// STRAIGHT GULLIES: sign(), not the sine itself. Pretending the gully
// slope is constant top to bottom -- an extruded triangle wave rather
// than a sine -- is what stops small gullies running parallel along
// large ones before branching. Note this feeds the INTERNAL slope only.
gully_slope += sign(ph.y) * ph_zw * sw * gully_weight;
vec3 gullies = vec3(ph.x, ph.y * ph_zw.x, ph.y * ph_zw.y);
vec3 faded = mix(vec3(fade_target, 0.0, 0.0), gullies * gully_weight, combi_mask);
height_and_slope += faded * sw;
magnitude += sw;
// STACKED FADING: this octave's faded gullies become the next octave's
// fade target, and its slope mask multiplies into the combi-mask. Each
// octave therefore restricts where the next one may act at all.
// STRUCTURAL TERMS BLEND BY THE SAME WEIGHT. fade_target and combi_mask
// are not accumulated, they are REPLACED each octave -- so weighting
// only the additive terms would still let a partial octave re-shape
// every finer one the moment it appears. Mixing from the previous value
// is what makes the whole chain continuous.
fade_target = mix(fade_target, faded.x, ow);
float rounding_for_octave = mix(rounding.y, rounding.x, clamp01(ph.x + 0.5)) * rounding_mult;
float new_mask = ease_out(smooth_start(sloping * onset.y, rounding_for_octave * onset.y));
combi_mask = mix(combi_mask, pow_inv(combi_mask, detail) * new_mask, ow);
ridge_fade_target = mix(ridge_fade_target, gullies.x, ridge_combi_mask * ow);
ridge_combi_mask *= mix(1.0, ease_out(sloping * onset.w), ow);
strength *= gain;
freq *= lacunarity;
rounding_mult *= rounding.w;
}
ErosionResult r;
r.delta = height_and_slope - input_has;
r.magnitude = magnitude;
// The last octave would otherwise appear with much thicker lines than the
// rest, never having been filtered through a mask; fading it towards
// neutral is what turns a fade target into a ridge map.
r.ridge = ridge_fade_target * (1.0 - ridge_combi_mask);
return r;
}
// =========================================================================
// ErosionNoiseNode's wrapper, hoisted.
//
// The erosion_filter above was duplicated three times; THIS was duplicated
// twice, and the second copy is where the 60x argument bug actually lived.
// The remap into the Shadertoy's [0,1] height space, the water-line strength
// gate and the height offset existed once as a GLSL string
// (ErosionNoiseNode::compile_glsl) and once as C++ (ProcNodesBridge::evaluate),
// and the C++ one had none of them. Fixing the arguments would have patched
// that instance; hoisting the wrapper removes the place the next one can
// appear.
//
// gn is the upstream GradientNoiseNode's vec3(height, dh/dx, dh/dz).
// GradientNoiseNode emits gn_slope = vec3(gn_result.y, 0.0, gn_result.z), so
// its .x and .z are gn.y and gn.z respectively -- that indirection is why the
// two 0.5 factors below look asymmetric. They are not.
//
// Returns vec2(height * height_scale, erosion), where erosion is the
// delta normalised by the octave-sum magnitude: -1 in creases, +1 on ridges,
// which is what the colour shader consumes.
//
// The five constants are ErosionNoiseNode's, not the caller's
// (noise_nodes.cpp:1229-1233). They are deliberately NOT parameters: they
// were literals in both former copies, and promoting them to arguments now
// would just move the divergence to the call sites.
// The node wrapper with EVERY constant exposed. erosion_node_eval below is
// this function with runevision's reference values supplied, and is what the
// engine calls; this one exists because the lookdev tool needs to turn those
// values into knobs and the alternative was a second copy of the wrapper.
//
// The five that used to be baked are baked for a REASON on the engine path --
// GpuCompiler's push-constant block is param0..param6, seven floats, a hard cap
// (gpu_compiler.cpp, get_shader_template), and rounding and onset alone would
// need eight. That cap does not apply to a consumer that emits GLSL text, which
// is why the lookdev renderer can reach them and the node still cannot.
//
// rounding.z is the one runevision singles out in the blog post: "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, it has been pinned
// at 0.1 throughout, and terrain_height has since moved 20 km -> 60 km.
// EVERY KNOB, NAMED FOR WHAT IT DOES.
//
// The old signature carried `detail`, `normalization`, `rounding` and `onset` --
// four names that say nothing about their effect, three of them packed into
// vec4s whose components had no names at all. Tuning meant reading the filter
// body to find out which component you were changing. runevision's article
// describes each of these in words; these are those words.
//
// erosion_strength how much the filter displaces the input
// gully_depth scales the gully contribution down
// fine_gully_slope_restriction LOWER values keep high-frequency gullies to
// steeper slopes -- the "detail only on the
// sloped parts" control
// gully_normalize_threshold only gullies above this get normalised, which
// is what keeps peaks pointy
// crease_rounding rounding applied in creases
// ridge_rounding rounding applied on ridges
// input_feature_rounding rounding for the INPUT height function's own
// feature size -- the one runevision singles
// out, and the one that must move when
// terrain_height does
// rounding_per_octave rounding multiplier compounded per octave
// gully_onset_slope slope at which gullies begin at all
// octave_onset_slope per-octave slope gate
// ridge_onset_slope the same for the ridge map
// ridge_octave_onset per-octave ridge gate
// assumed_slope the slope magnitude to pretend the input has
// assumed_slope_mix 0 = use the real slope, 1 = fully assumed
// height_offset flat offset, times the octave magnitude
// height_offset_follows_altitude mixes that offset toward -fade_target
// input_height_amplitude scales the incoming field into filter space
// octave_lacunarity frequency step per octave
// peak_height_divisor what counts as "a peak" for the fade target
// gully_band_scale coarsest gully wavelength, as a fraction of
// terrain_height
// cell_size phacelle cell size; too large distorts, too
// small goes grainy
vec2 erosion_node_eval_full(vec2 p, vec3 gn, float height_scale, float erosion_strength,
float cell_size, int octaves, float gain,
float peak_height_divisor, float gully_band_scale,
float input_height_amplitude, float octave_lacunarity, float gully_depth,
float fine_gully_slope_restriction,
float gully_normalize_threshold, vec2 height_offset_pair,
vec4 rounding_set, vec4 onset_set, vec2 assumed_slope_pair,
float max_octaves) {
vec3 n = vec3(gn.x * input_height_amplitude, gn.y * input_height_amplitude,
gn.z * input_height_amplitude);
float fade_target = clamp(n.x / (input_height_amplitude * peak_height_divisor), -1.0, 1.0);
n = vec3(n.x * 0.5 + 0.5, n.y * 0.5, n.z * 0.5);
ErosionResult h = erosion_filter(p, n, fade_target,
erosion_strength, gully_depth, fine_gully_slope_restriction,
rounding_set, onset_set, assumed_slope_pair,
gully_band_scale, octaves, octave_lacunarity, gain, cell_size,
gully_normalize_threshold, max_octaves);
float offset = mix(height_offset_pair.x, -fade_target, height_offset_pair.y) * h.magnitude;
float eroded = n.x + h.delta.x + offset;
return vec2(eroded * height_scale, h.ridge);
}
// THE SAME EVALUATION, KEEPING THE SLOPE IT ALREADY COMPUTED.
//
// erosion_filter accumulates height AND both slope components into one vec3 and
// returns all three as `delta`; every caller then used `.x` and discarded `.y`
// and `.z`. Meanwhile edge_intersect.glsl reconstructs the very same gradient by
// central differences -- SIX more full evaluations of this filter per edge, each
// one seven octaves of phacelle noise. Measured on the real device: the octave
// loop takes the pipeline from 491 us to 3795 us per 32^3 grid, and
// edge_intersect is 91% of the total.
//
// So this returns what was already there. The slope is in the same space as the
// height: d(filter output)/dp where p = world / terrain_height, and n's own
// components carried the 0.5 factors from the [-1,1] -> [0,1] remap, which is
// why .y and .z are halved on the way in and are NOT re-scaled here -- the
// caller converts to world units exactly as it converts the height.
//
// vec4 rather than a second function with an out-parameter: `out` is a GLSL
// keyword the C++ shim would have to imitate, and this file is compiled as
// both. Returning a wider vector is the one thing both languages spell the same.
//
// .x = eroded height * height_scale (identical to erosion_node_eval_full)
// .y = ridge (identical)
// .z = d(eroded)/dp.x (NEW -- was computed and dropped)
// .w = d(eroded)/dp.z (NEW)
vec4 erosion_node_eval_grad(vec2 p, vec3 gn, float height_scale, float erosion_strength,
float cell_size, int octaves, float gain,
float peak_height_divisor, float gully_band_scale,
float input_height_amplitude, float octave_lacunarity, float gully_depth,
float fine_gully_slope_restriction,
float gully_normalize_threshold, vec2 height_offset_pair,
vec4 rounding_set, vec4 onset_set, vec2 assumed_slope_pair,
float max_octaves) {
vec3 n = vec3(gn.x * input_height_amplitude, gn.y * input_height_amplitude,
gn.z * input_height_amplitude);
float fade_target = clamp(n.x / (input_height_amplitude * peak_height_divisor), -1.0, 1.0);
n = vec3(n.x * 0.5 + 0.5, n.y * 0.5, n.z * 0.5);
ErosionResult h = erosion_filter(p, n, fade_target,
erosion_strength, gully_depth, fine_gully_slope_restriction,
rounding_set, onset_set, assumed_slope_pair,
gully_band_scale, octaves, octave_lacunarity, gain, cell_size,
gully_normalize_threshold, max_octaves);
float offset = mix(height_offset_pair.x, -fade_target, height_offset_pair.y) * h.magnitude;
float eroded = n.x + h.delta.x + offset;
// The offset is a function of fade_target, which is a function of the INPUT
// height, so strictly it carries a gradient of its own. It is dropped here
// for the same reason the wrapper drops d(relief_amp)/dp: with
// height_offset_follows_altitude at its default 0 the term is a constant
// times magnitude and its gradient is exactly zero. Raise that knob and this
// becomes an approximation -- stated here rather than discovered later.
return vec4(eroded * height_scale, h.ridge,
(n.y + h.delta.y) * height_scale,
(n.z + h.delta.z) * height_scale);
}
// runevision's reference values, for callers that do not tune.
//
// slope_power IS GONE from this signature. It was accepted and ignored -- the
// frequency approach's exponent on slope magnitude, which the fade approach
// does not have -- and a parameter that silently does nothing is worse than no
// parameter at all. Callers drop the argument.
vec2 erosion_node_eval(vec2 p, vec3 gn, float height_scale, float erosion_strength,
float cell_size, int octaves, float gain,
float peak_height_divisor, float gully_band_scale, float max_octaves) {
// Base-noise shaping, unchanged from the Shadertoy lineage.
// The advanced filter's parameters. These are runevision's reference
// defaults, and they are CONSTANTS here rather than node sockets because
// GpuCompiler's push-constant block is param0..param6 -- seven floats,
// hard cap (gpu_compiler.cpp, get_shader_template). rounding and onset
// alone would need eight. Promoting them would mean a second uniform
// buffer through the whole node pipeline; baking them is what the five
// constants above them have always done.
// EROSION_SCALE is a PARAMETER now, not a constant. It sets the gully band:
// the coarsest gully wavelength is terrain_height * erosion_scale *
// cell_scale, in world metres, halving per octave. Baked, it made the band a
// fixed FRACTION of the terrain scale -- so tripling terrain_height from
// 20 km to 60 km to buy view distance tripled the gully wavelength with it,
// 2100 m to 6300 m, and the coarsest octave became a smear across a whole
// mountain rather than a channel down its side. The two bands have no reason
// to be locked together.
// POINTY PEAKS; normalising the
// gullies otherwise rounds them off
// Height offset, mix(a, -fade_target, b) * magnitude. b = 0 makes it a
// flat -0.65 of the accumulated magnitude; raise b to tie the offset to
// altitude instead.
// slope_power is ACCEPTED AND IGNORED. It was the frequency approach's
// exponent on slope magnitude, and the fade approach has no such term --
// stripe frequency no longer depends on slope at all. The socket stays so
// existing graphs keep loading; wiring it to something would be worse than
// leaving it dead, because a parameter that silently does nothing else is
// harder to find than one that does nothing.
// THE FADE TARGET: -1 in valleys, +1 on peaks. Overshoot is fine; the
// filter clamps. Derived from altitude here, which is what runevision's
// reference does and is not robust in general -- a high valley reads as a
// peak. Curvature would be principled but needs second derivatives from
// every upstream node.
//
// fade_div sets what counts as "a peak". runevision uses 0.6 for a height
// field of order 1. It is a PARAMETER rather than that constant because
// the caller may hand this function a combined terrain field several times
// larger -- dc-terrain folds its mountain layer into gn, which is ~8x the
// bare gradient noise -- and 0.6 against that saturates the fade target to
// +/-1 everywhere, which is "every point is a peak" and erases the fade
// approach's whole reason to exist. Scale it with the field.
// [-1,1] -> [0,1], as the Shadertoy does.
// NO WATER-LINE GATE. The old wrapper multiplied strength by
// smoothstep(WATER_HEIGHT +/- 0.1, n.x), inherited from the Shadertoys,
// where it exists to hide the frequency approach's bulges at valley
// bottoms. The fade approach does not produce them, so the gate has
// nothing to hide and was actively wrong here: WATER_HEIGHT 0.465 lives in
// the Shadertoy's [0,1] height space and has no relationship to this
// world's floor_y/ceiling_y.
return erosion_node_eval_full(p, gn, height_scale, erosion_strength,
cell_size, octaves, gain, peak_height_divisor, gully_band_scale,
/* input_height_amplitude */ 0.25,
/* octave_lacunarity */ 2.0,
/* gully_depth */ 0.5,
/* fine_gully_slope_restriction */ 1.5,
/* gully_normalize_threshold */ 0.5,
/* height_offset, follows_alt */ vec2(-0.65, 0.0),
/* crease/ridge/input/per_oct */ vec4(0.1, 0.0, 0.1, 2.0),
/* gully/octave/ridge/ridge_oct */ vec4(1.25, 1.25, 2.8, 1.5),
/* assumed_slope, mix */ vec2(0.7, 1.0),
/* fractional octave taper */ max_octaves);
}
// =========================================================================
// hypso_curve.inc â the analytic land/sea transfer curve.
//
// ONE TEXT, compiled as GLSL (injected into the density shader) and as C++
// (through erosion_cpp_shim.h, for the oracle). Same rule as
// erosion_filter.inc: only the subset both languages agree on.
//
// WHAT IT IS FOR. A plain fBm has a UNIMODAL height distribution, so
// intersecting it with a water plane at any level gives a coastline of
// uniform fractal character everywhere -- no shelves, no basins, and no
// reason for land to clump. Earth's hypsography is BIMODAL: ~85% of the
// surface sits in two narrow bands, because continental crust and oceanic
// crust have different densities. This curve puts that shape into the
// height field directly, with no real-world data.
//
// WHY A BLEND AND NOT A PIECEWISE JOIN. The obvious construction --
// `h < sea ? ocean(h) : land(h)` -- is C0 but not C1: the slope jumps at
// the join, which is EXACTLY the waterline, and a slope discontinuity along
// every coastline renders as a crease ring. Blending two regimes with a
// smoothstep is C1 by construction because smoothstep is.
//
// WHERE THE BIMODALITY COMES FROM. Not from a threshold. Output values
// cluster where the curve is FLAT (low |dhv/dh|), so a flat abyssal
// plateau and a land ramp separated by a narrow STEEP section produce two
// modes with few samples between them -- the continental slope. `shelf` is
// the half-width of that steep section: smaller = sharper coast = more
// strongly bimodal.
//
// THE WATERLINE IS DERIVED. The curve is monotonic, so sea level in OUTPUT
// space is simply hypso_curve(sea).x. Do not also hand-set a water height;
// read it from here, or the two disagree and the shoreline renders where
// the terrain is not.
//
// PEAK LIFT is additive on the upper tail only, so it raises summits
// without steepening mid-slopes or deepening valleys. That matters beyond
// looks: total relief is bounded by the physics band (character movement
// degrades with |y|), and lifting only the tail is how peak height stops
// competing with valley depth for the same budget.
//
// h input height, nominally [-1, 1]
// sea input value at the middle of the continental slope
// shelf half-width of that slope, in input units (>0)
// abyss OUTPUT value of the abyssal plateau (negative)
// land_slope output units per input unit above sea
// peak_onset input value where the summit lift begins
// peak_gain output units added at h = 1
//
// Returns (hv, dhv/dh). The derivative is not decoration: the erosion
// filter is slope-aware and the wrapper hands it a gradient, so remapping
// the height without remapping the gradient by the chain rule would carve
// gullies aligned to a surface that is no longer there.
// =========================================================================
// WHY THREE REGIMES AND NOT TWO. A first version had an abyssal plateau and
// a linear land ramp, and the gate measured its output histogram as ONE mode
// plus a smear -- because a ramp has no mode. Earth's bimodality is two
// PLATEAUS: the continental platform (~0..+1 km) and the abyssal plain
// (~-4..-5 km), with the continental slope between them. So the land side
// needs a low-gradient shelf of its own, and the mountains rise ABOVE that.
//
// abyss --[continental slope]-- platform --[range front]-- mountains
// flat steep gentle steeper + peak lift
//
// Two smoothsteps composing three lines. Both are C1, so the whole is.
// smoothstep's derivative. GLSL has smoothstep but not its derivative, and
// writing 6t(1-t)/(hi-lo) at three call sites is how one of them ends up
// missing the 1/(hi-lo).
// â DO NOT NAME A LOCAL `out`. It is a GLSL storage qualifier and a perfectly
// ordinary identifier in C++, so a C++-only test compiles it and the SHADER
// dies with "syntax error, unexpected OUT" -- which is what happened: the GPU
// fell back to the CPU mesher and the only symptom was that terrain generation
// got slow. Same class as erosion_filter.inc's PN_TAU note. `in`, `inout`,
// `filter`, `sample` and `active` are reserved the same way.
float hypso_dsmooth(float lo, float hi, float x) {
float d = hi - lo;
if (d <= 0.0) return 0.0;
float t = clamp((x - lo) / d, 0.0, 1.0);
return 6.0 * t * (1.0 - t) / d;
}
// h input height, nominally [-1, 1]
// sea input value that IS sea level; the curve returns 0 there
// shelf_width how far below sea the gentle shelf runs before dropping
// slope_width half-width of the steep drop (smaller = sharper break)
// abyss OUTPUT value of the abyssal plain (negative)
// plat_slope output per input on the continental platform (small)
// mtn_onset input value where the mountain front begins
// mtn_width half-width of that front
// mtn_slope output per input in the mountains (larger)
// peak_onset input value where the summit lift begins
// peak_gain output added at h = 1
//
// Returns (hv, dhv/dh). The derivative is not decoration: the erosion
// filter is slope-aware and the wrapper hands it a gradient, so remapping
// the height without remapping the gradient by the chain rule would carve
// gullies aligned to a surface that is no longer there.
vec2 hypso_curve(float h, float sea, float shelf_width, float slope_width,
float abyss,
float plat_slope, float mtn_onset, float mtn_width,
float mtn_slope, float peak_onset, float peak_gain) {
// --- regime 1: the continental platform, a gentle line through (sea, 0)
float plat = plat_slope * (h - sea);
// --- the continental slope, CENTRED BELOW SEA LEVEL.
//
// â THIS IS WHY THE FIRST VERSION WAS A CLIFF. The blend was centred ON
// `sea`, so the steep section straddled the waterline and land met deep
// water with nothing in between -- a wall, roughly 50x steeper than the
// platform behind it. That steepness IS the bimodality, so it cannot be
// removed; it has to be MOVED.
//
// On Earth the steep part is submerged: land -> beach -> a gentle shelf
// ~100 m deep and tens of km wide -> shelf break -> continental slope ->
// abyssal plain. Putting the blend `shelf_width` BELOW sea leaves the
// gentle platform line running down through the waterline, and that
// stretch is the continental shelf.
//
// Requires shelf_width > slope_width, or the drop reaches back up through
// sea level and the shelf does not exist.
float slope_mid = sea - shelf_width;
float c0 = slope_mid - slope_width;
float c1 = slope_mid + slope_width;
float w = smoothstep(c0, c1, h);
float dw = hypso_dsmooth(c0, c1, h);
float base = abyss + (plat - abyss) * w;
float dbase = (plat - abyss) * dw + w * plat_slope;
// --- the mountain front, ADDED not blended to. See the note above:
// blending toward a steeper LINE folds the terrain below its onset,
// measured at dhv/dh = -0.112. An additive ramp cannot fold.
float a = mtn_onset - mtn_width;
float b = mtn_onset + mtn_width;
float sm = smoothstep(a, b, h);
float dsm = hypso_dsmooth(a, b, h);
// (h-a)*smoothstep: exactly 0 below a (smoothstep clamps), asymptotically
// linear with unit slope above b, C1 throughout, never decreasing.
float ramp = (h - a) * sm;
float dramp = sm + (h - a) * dsm;
float extra = mtn_slope - plat_slope; // the slope the range ADDS
float hv = base + extra * ramp;
float dhv = dbase + extra * dramp;
// --- summit lift, upper tail only. Same additive shape, same reason.
if (peak_gain != 0.0) {
hv += peak_gain * smoothstep(peak_onset, 1.0, h);
dhv += peak_gain * hypso_dsmooth(peak_onset, 1.0, h);
}
return vec2(hv, dhv);
}
float proc_density(vec3 pos, inout int pn_max_octaves, inout vec3 pn_terrain_offset, inout float pn_fade_div, inout float pn_erosion_scale, inout float pn_out_ridge, inout vec2 pn_out_grad, inout float pn_ero_input_height_amplitude, inout float pn_ero_octave_lacunarity, inout float pn_ero_gully_depth, inout float pn_ero_fine_gully_slope_restriction, inout float pn_ero_gully_normalize_threshold, inout float pn_ero_height_offset, inout float pn_ero_height_offset_follows_altitude, inout float pn_ero_crease_rounding, inout float pn_ero_ridge_rounding, inout float pn_ero_input_feature_rounding, inout float pn_ero_rounding_per_octave, inout float pn_ero_gully_onset_slope, inout float pn_ero_octave_onset_slope, inout float pn_ero_ridge_onset_slope, inout float pn_ero_ridge_octave_onset, inout float pn_ero_assumed_slope, inout float pn_ero_assumed_slope_mix, inout float pn_ero_max_octaves) {
float query_y = pos.y; // Save original query height before nodes modify pos
vec3 normal = vec3(0.0, 1.0, 0.0);
// Node 0: Gradient Noise
vec3 gn_result = fractal_noise_2d(vec2(pos.x, pos.z), 3.0, int(min(3, pn_max_octaves)), 2.0, 0.100000001);
float gn_height = gn_result.x;
vec3 gn_slope = vec3(gn_result.y, 0.0, gn_result.z);
// Node 1: Erosion Noise
vec4 en_result = erosion_node_eval_grad(
vec2(pos.x, pos.z),
vec3(gn_height, gn_slope.x, gn_slope.z) + pn_terrain_offset,
5.0, 0.150000006, 0.699999988, int(min(7, pn_max_octaves)), 0.5,
pn_fade_div, pn_erosion_scale,
pn_ero_input_height_amplitude, pn_ero_octave_lacunarity,
pn_ero_gully_depth, pn_ero_fine_gully_slope_restriction,
pn_ero_gully_normalize_threshold,
vec2(pn_ero_height_offset, pn_ero_height_offset_follows_altitude),
vec4(pn_ero_crease_rounding, pn_ero_ridge_rounding,
pn_ero_input_feature_rounding, pn_ero_rounding_per_octave),
vec4(pn_ero_gully_onset_slope, pn_ero_octave_onset_slope,
pn_ero_ridge_onset_slope, pn_ero_ridge_octave_onset),
vec2(pn_ero_assumed_slope, pn_ero_assumed_slope_mix),
pn_ero_max_octaves);
float en_erosion = en_result.y;
pn_out_ridge = en_erosion;
// SECOND NODE -> WRAPPER egress, same pattern as pn_out_ridge above.
// erosion_filter already accumulated both slope components; _full returned
// vec2 and threw them away, so every consumer that wanted a normal rebuilt
// it with SIX more full evaluations. _grad is the same evaluation with a
// wider return -- .x and .y are bit-identical to _full's, which is why this
// swap cannot move the terrain. Measured: 6 taps 16.4 ms vs 1 tap 1.57 ms
// fullscreen at 1080p, which is the whole reason this egress exists.
pn_out_grad = vec2(en_result.z, en_result.w);
pos.y = en_result.x;
// SDF: terrain_height - query_y. Positive below surface, negative above.
return pos.y - query_y;
}
float evaluate_density(float x, float y, float z, float voxel_scale, inout float pn_out_ridge, inout vec2 pn_out_grad) {
int pn_max_octaves = 0;
vec3 pn_terrain_offset = vec3(0.0);
float pn_fade_div = 0.6;
float pn_erosion_scale = 0.15;
float pn_ero_input_height_amplitude = 0.25;
float pn_ero_octave_lacunarity = 2.0;
float pn_ero_gully_depth = 0.5;
float pn_ero_fine_gully_slope_restriction = 1.5;
float pn_ero_gully_normalize_threshold = 0.5;
float pn_ero_height_offset = -0.65;
float pn_ero_height_offset_follows_altitude = 0.0;
float pn_ero_crease_rounding = 0.1;
float pn_ero_ridge_rounding = 0.0;
float pn_ero_input_feature_rounding = 0.1;
float pn_ero_rounding_per_octave = 2.0;
float pn_ero_gully_onset_slope = 1.25;
float pn_ero_octave_onset_slope = 1.25;
float pn_ero_ridge_onset_slope = 2.8;
float pn_ero_ridge_octave_onset = 1.5;
float pn_ero_assumed_slope = 0.7;
float pn_ero_assumed_slope_mix = 1.0;
float pn_ero_max_octaves = 999.0;
float inv_scale = 1.0 / 60000.0;
pn_max_octaves = 99;
// MOUNTAIN LAYER, height AND gradient. The gradient used to be computed and
// thrown away: mtn_offset was added to y_mapped AFTER the erosion, so the
// filter never saw the range it was carving and aligned its gullies to gn
// alone -- three octaves at gain 0.1, i.e. one smooth 6.7 km blob. It read
// as noise laid over a mountain rather than erosion of one.
// MULTIFRACTAL. variance 0 is bit-identical to fractal_noise_2d, so the
// default costs nothing; above 0 each octave is weighted by the elevation
// accumulated so far, putting the roughness on high ground. Plain fBm is
// statistically homogeneous -- equally rugged in the basins and on the
// peaks -- which is what reads as every range looking the same.
// NYQUIST LIMITS. voxel_scale is the LOD's voxel size in metres and is
// declared by BOTH shaders this body is injected into (density_fill.glsl
// and edge_intersect.glsl), above the injection point. One cycle of a
// layer spans terrain_height/scale world metres, so that is its base
// wavelength. Fractional, so an octave fades in as a chunk refines
// rather than snapping on and moving the surface underfoot.
float mtn_lim = nyquist_octaves(60000.0 / max(0.1, 1e-6), voxel_scale, 1.87);
float rel_lim = nyquist_octaves(60000.0 / max(0.5, 1e-6), voxel_scale, 2.0);
float det_lim = nyquist_octaves(60000.0 / max(600.0, 1e-6), voxel_scale, 1.87);
vec3 mtn = fractal_multifractal_2d_lod(vec2(x, z) * inv_scale * 0.1, 1.0, 8, 1.87, 0.68, 0.0, mtn_lim);
float relief_amp = 1.0;
if (0.5 > 0.0) {
vec3 rel = fractal_noise_2d_lod(vec2(x, z) * inv_scale * 0.5, 1.0, 2, 2.0, 0.5, rel_lim);
float t = clamp(rel.x * 0.5 + 0.5, 0.0, 1.0);
relief_amp = 0.08 + (1.0 - 0.08) * t;
}
float mtn_amp = 0.55 * relief_amp;
// mtn is evaluated at p * mountain_scale, so d/dp = mtn.yz * mountain_scale.
// d(relief_amp)/dp is dropped: relief is a 40 km envelope, orders coarser
// than anything gully direction responds to. The 8.0 is exact, not fitted --
// gn reaches the surface as gn.x * HEIGHT_AMP * 0.5 = gn.x * 0.125, so
// dividing the offset by 0.125 reproduces the identical base surface at
// zero erosion strength.
vec2 hc = hypso_curve(mtn.x, 0.2, 0.12, 0.05, -0.62, 0.13,
0.35, 0.1, 1.3, 0.6, 0.3);
pn_terrain_offset = vec3(hc.x, mtn.y * 0.1 * hc.y, mtn.z * 0.1 * hc.y) * 0.55 * 8.0;
// DETAIL rides ON TOP of the eroded surface and is deliberately NOT folded
// into pn_terrain_offset: its gradient carries a factor of detail_scale, so
// inside one gully cell the direction flips several times and the stripes go
// chaotic -- measured 77.6 m of relief across a 2 m step, i.e. vertical
// noise. Detail is what water runs over, not what steers it.
pn_fade_div = 4.8;
pn_erosion_scale = 0.15;
pn_ero_input_height_amplitude = 0.25;
pn_ero_octave_lacunarity = 2.0;
pn_ero_gully_depth = 0.5;
pn_ero_fine_gully_slope_restriction = 1.5;
pn_ero_gully_normalize_threshold = 0.5;
pn_ero_height_offset = -0.649999976;
pn_ero_height_offset_follows_altitude = 0.0;
pn_ero_crease_rounding = 0.100000001;
pn_ero_ridge_rounding = 0.0;
pn_ero_input_feature_rounding = 0.100000001;
pn_ero_rounding_per_octave = 2.0;
pn_ero_gully_onset_slope = 1.25;
pn_ero_octave_onset_slope = 1.25;
pn_ero_ridge_onset_slope = 2.799999952;
pn_ero_ridge_octave_onset = 1.5;
pn_ero_assumed_slope = 0.699999988;
pn_ero_assumed_slope_mix = 1.0;
pn_ero_max_octaves = nyquist_octaves(6300.0, voxel_scale, 2.0);
float det_offset = 0.0;
if (600.0 > 0.0) {
vec3 det = fractal_noise_2d_lod(vec2(x, z) * inv_scale * 600.0, 1.0, 4, 1.87, 0.5, det_lim);
float det_slope = length(vec2(pn_terrain_offset.y, pn_terrain_offset.z));
float det_mask = 1.0 - smoothstep(0.3, 1.2, det_slope);
det_offset = det.x * 0.00142857 * det_mask;
}
float y_mapped = (y - -9100.0) * 0.000238095 - det_offset * 5.0;
float d = proc_density(vec3(x * inv_scale, y_mapped, z * inv_scale), pn_max_octaves, pn_terrain_offset, pn_fade_div, pn_erosion_scale, pn_out_ridge, pn_out_grad, pn_ero_input_height_amplitude, pn_ero_octave_lacunarity, pn_ero_gully_depth, pn_ero_fine_gully_slope_restriction, pn_ero_gully_normalize_threshold, pn_ero_height_offset, pn_ero_height_offset_follows_altitude, pn_ero_crease_rounding, pn_ero_ridge_rounding, pn_ero_input_feature_rounding, pn_ero_rounding_per_octave, pn_ero_gully_onset_slope, pn_ero_octave_onset_slope, pn_ero_ridge_onset_slope, pn_ero_ridge_octave_onset, pn_ero_assumed_slope, pn_ero_assumed_slope_mix, pn_ero_max_octaves);
d *= 4200.0;
pn_out_grad *= inv_scale * 4200.0;
// Clip: solid below floor_y, air above ceiling_y.
d = max(d, (-9100.0 - y) * 2.0);
d = min(d, (11900.0 - y) * 2.0);
return d;
}
void fragment() {
vec3 wp = (INV_VIEW_MATRIX * vec4(VERTEX, 1.0)).xyz;
if (probe_world_span > 0.0) {
vec2 w = probe_world_origin + (UV - 0.5) * probe_world_span;
wp = vec3(w.x, probe_world_y, w.y);
}
float vs = force_voxel_scale > 0.0
? force_voxel_scale
: (probe_world_span > 0.0 ? probe_world_span / 1024.0
: max(length(fwidth(wp)), 0.01));
float ridge = 0.0;
vec2 grad = vec2(0.0);
float curv_k = 0.5;
vec3 n;
if (probe_taps == 1) {
evaluate_density(wp.x, wp.y, wp.z, vs, ridge, grad);
n = normalize(vec3(-grad.x, 1.0, -grad.y));
} else {
// FOUR taps, not six. d = H(x,z) - y, so dd/dy is analytically exact
// and constant -- sampling it was two evaluations spent rediscovering
// a number the wrapper already emits. Only x and z need differencing.
float e = vs;
float dc0 = evaluate_density(wp.x, wp.y, wp.z, vs, ridge, grad);
float dxp = evaluate_density(wp.x + e, wp.y, wp.z, vs, ridge, grad);
float dxm = evaluate_density(wp.x - e, wp.y, wp.z, vs, ridge, grad);
float dzp = evaluate_density(wp.x, wp.y, wp.z + e, vs, ridge, grad);
float dzm = evaluate_density(wp.x, wp.y, wp.z - e, vs, ridge, grad);
n = normalize(vec3(-(dxp - dxm), 2.0 * e, -(dzp - dzm)));
// erosion_filter's OWN crease/ridge map, not a numerical second
// derivative. A Laplacian at pixel scale is swamped by the 100 m detail
// layer -- second differences amplify exactly what detail has most of --
// and came out as salt-and-pepper everywhere det_mask was not zero.
curv_k = clamp(0.5 + probe_curv_gain * ridge, 0.0, 1.0);
}
float lam = clamp(dot(n, normalize(probe_sun)), 0.0, 1.0);
if (probe_view == 4) {
// Curvature alone: black creases, white ridges. Saturated to pure black
// and white means the gain is too high to carry information.
ALBEDO = vec3(curv_k);
} else if (probe_view == 0) {
ALBEDO = n * 0.5 + 0.5;
} else if (probe_view == 1) {
ALBEDO = vec3(0.05 + 0.95 * lam);
} else {
// d = H(x,z) - y in world metres, so the SURFACE elevation under this
// sample is just y + d. One extra evaluation, and it is what both the
// material's altitude term and the hypsometric ramp key off.
float dc = evaluate_density(wp.x, wp.y, wp.z, vs, ridge, grad);
float h = wp.y + dc;
float ny = clamp(abs(n.y), 0.001, 1.0);
float slope = sqrt(max(1.0 - ny * ny, 0.0)) / ny;
float t = clamp((h - probe_field_range.x)
/ max(probe_field_range.y - probe_field_range.x, 1.0), 0.0, 1.0);
// Curvature now comes from the field's Laplacian, the same way the
// map's shader derives it -- so this preview shows the real material.
vec3 base = probe_view == 2 ? lookdev_material(t, slope, curv_k)
: hypso_ramp(h);
ALBEDO = base * (0.15 + 0.85 * lam);
}
}