vec3 mod289_v3(vec3 x) {
return x - floor(x * (1.0 / 289.0)) * 289.0;
}
vec4 mod289_v4(vec4 x) {
return x - floor(x * (1.0 / 289.0)) * 289.0;
}
vec4 permute(vec4 x) {
return mod289_v4(((x * 34.0) + 10.0) * x);
}
vec4 taylor_inv_sqrt(vec4 r) {
return 1.79284291400159 - 0.85373472095314 * r;
}
float simplex_noise_3d(vec3 v) {
const vec2 C = vec2(1.0 / 6.0, 1.0 / 3.0);
const vec4 D = vec4(0.0, 0.5, 1.0, 2.0);
// First corner
vec3 i = floor(v + dot(v, C.yyy));
vec3 x0 = v - i + dot(i, C.xxx);
// Other corners
vec3 g = step(x0.yzx, x0.xyz);
vec3 l = 1.0 - g;
vec3 i1 = min(g.xyz, l.zxy);
vec3 i2 = max(g.xyz, l.zxy);
vec3 x1 = x0 - i1 + C.xxx;
vec3 x2 = x0 - i2 + C.yyy;
vec3 x3 = x0 - D.yyy;
// Permutations
i = mod289_v3(i);
vec4 p = permute(permute(permute(
i.z + vec4(0.0, i1.z, i2.z, 1.0))
+ i.y + vec4(0.0, i1.y, i2.y, 1.0))
+ i.x + vec4(0.0, i1.x, i2.x, 1.0));
// Gradients: 7x7 points over a square, mapped onto an octahedron
float n_ = 0.142857142857; // 1.0 / 7.0
vec3 ns = n_ * D.wyz - D.xzx;
vec4 j = p - 49.0 * floor(p * ns.z * ns.z);
vec4 x_ = floor(j * ns.z);
vec4 y_ = floor(j - 7.0 * x_);
vec4 x = x_ * ns.x + ns.yyyy;
vec4 y = y_ * ns.x + ns.yyyy;
vec4 h = 1.0 - abs(x) - abs(y);
vec4 b0 = vec4(x.xy, y.xy);
vec4 b1 = vec4(x.zw, y.zw);
vec4 s0 = floor(b0) * 2.0 + 1.0;
vec4 s1 = floor(b1) * 2.0 + 1.0;
vec4 sh = -step(h, vec4(0.0));
vec4 a0 = b0.xzyw + s0.xzyw * sh.xxyy;
vec4 a1 = b1.xzyw + s1.xzyw * sh.zzww;
vec3 p0 = vec3(a0.xy, h.x);
vec3 p1 = vec3(a0.zw, h.y);
vec3 p2 = vec3(a1.xy, h.z);
vec3 p3 = vec3(a1.zw, h.w);
// Normalise gradients
vec4 norm = taylor_inv_sqrt(vec4(dot(p0, p0), dot(p1, p1), dot(p2, p2), dot(p3, p3)));
p0 *= norm.x;
p1 *= norm.y;
p2 *= norm.z;
p3 *= norm.w;
// Mix contributions from the four corners
vec4 m = max(0.5 - vec4(dot(x0, x0), dot(x1, x1), dot(x2, x2), dot(x3, x3)), 0.0);
m = m * m;
return 105.0 * dot(m * m, vec4(dot(p0, x0), dot(p1, x1), dot(p2, x2), dot(p3, x3)));
}
float noise_fbm(vec3 pos, int octaves, float roughness, float lacunarity) {
float value = 0.0;
float amplitude = 1.0;
float frequency = 1.0;
float max_amp = 0.0;
for (int i = 0; i < octaves; i++) {
value += amplitude * simplex_noise_3d(pos * frequency);
max_amp += amplitude;
amplitude *= roughness;
frequency *= lacunarity;
}
return value / max_amp;
}
// Voronoi noise â distance to nearest cell center
float noise_voronoi(vec3 pos) {
vec3 i_pos = floor(pos);
vec3 f_pos = fract(pos);
float min_dist = 1e10;
for (int x = -1; x <= 1; x++) {
for (int y = -1; y <= 1; y++) {
for (int z = -1; z <= 1; z++) {
vec3 offset = vec3(float(x), float(y), float(z));
vec3 cell = i_pos + offset;
// Hash-based random point in cell
vec3 r = fract(sin(vec3(
dot(cell, vec3(127.1, 311.7, 74.7)),
dot(cell, vec3(269.5, 183.3, 246.1)),
dot(cell, vec3(113.5, 271.9, 124.6))
)) * 43758.5453);
vec3 cell_pos = offset + r - f_pos;
float d = length(cell_pos);
min_dist = min(min_dist, d);
}
}
}
return min_dist;
}
// Musgrave Multifractal
float noise_musgrave_multi(vec3 pos, int octaves, float roughness, float lacunarity) {
float value = 1.0;
float frequency = 1.0;
for (int i = 0; i < octaves; i++) {
float n = simplex_noise_3d(pos * frequency);
value *= (n * roughness + 1.0);
frequency *= lacunarity;
}
return value;
}
// Musgrave Heterogeneous
float noise_musgrave_hetero(vec3 pos, int octaves, float roughness, float lacunarity) {
float value = simplex_noise_3d(pos) + 1.0;
float frequency = lacunarity;
float amplitude = roughness;
for (int i = 1; i < octaves; i++) {
float n = (simplex_noise_3d(pos * frequency) + 1.0) * amplitude * value;
value += n;
frequency *= lacunarity;
amplitude *= roughness;
}
return value;
}
// Musgrave Hybrid Multifractal
float noise_musgrave_hybrid(vec3 pos, int octaves, float roughness, float lacunarity) {
float value = simplex_noise_3d(pos) + 1.0;
float weight = value;
float frequency = lacunarity;
float amplitude = roughness;
for (int i = 1; i < octaves; i++) {
weight = min(weight, 1.0);
float n = (simplex_noise_3d(pos * frequency) + 1.0) * amplitude;
value += weight * n;
weight *= n;
frequency *= lacunarity;
amplitude *= roughness;
}
return value;
}
// Musgrave Ridged Multifractal
float noise_musgrave_ridged(vec3 pos, int octaves, float roughness, float lacunarity) {
float value = 0.0;
float frequency = 1.0;
float amplitude = 1.0;
float weight = 1.0;
for (int i = 0; i < octaves; i++) {
float n = abs(simplex_noise_3d(pos * frequency));
n = 1.0 - n;
n *= n;
n *= weight;
weight = clamp(n * 2.0, 0.0, 1.0);
value += n * amplitude;
frequency *= lacunarity;
amplitude *= roughness;
}
return value;
}
// =========================================================================
// Perlin noise â Blender's exact implementation (Jenkins Lookup3 hash)
// =========================================================================
#define PERLIN_ROT(x, k) (((x) << (k)) | ((x) >> (32u - (k))))
uint jenkins_hash3(uint kx, uint ky, uint kz) {
uint a, b, c;
a = b = c = 0xdeadbeefu + (3u << 2u) + 13u;
c += kz; b += ky; a += kx;
// hash_bit_final
c ^= b; c -= PERLIN_ROT(b, 14u);
a ^= c; a -= PERLIN_ROT(c, 11u);
b ^= a; b -= PERLIN_ROT(a, 25u);
c ^= b; c -= PERLIN_ROT(b, 16u);
a ^= c; a -= PERLIN_ROT(c, 4u);
b ^= a; b -= PERLIN_ROT(a, 14u);
c ^= b; c -= PERLIN_ROT(b, 24u);
return c;
}
float perlin_fade(float t) {
return t * t * t * (t * (t * 6.0 - 15.0) + 10.0);
}
float perlin_negate_if(float val, uint cond) {
return (cond != 0u) ? -val : val;
}
float perlin_noise_grad(uint hash, float x, float y, float z) {
uint h = hash & 15u;
float u = h < 8u ? x : y;
float vt = ((h == 12u) || (h == 14u)) ? x : z;
float v = h < 4u ? y : vt;
return perlin_negate_if(u, h & 1u) + perlin_negate_if(v, h & 2u);
}
float perlin_tri_mix(
float v0, float v1, float v2, float v3,
float v4, float v5, float v6, float v7,
float x, float y, float z) {
float x1 = 1.0 - x;
float y1 = 1.0 - y;
float z1 = 1.0 - z;
return z1 * (y1 * (v0 * x1 + v1 * x) + y * (v2 * x1 + v3 * x)) +
z * (y1 * (v4 * x1 + v5 * x) + y * (v6 * x1 + v7 * x));
}
float perlin_noise_3d(vec3 p) {
int X = int(floor(p.x)); float fx = p.x - floor(p.x);
int Y = int(floor(p.y)); float fy = p.y - floor(p.y);
int Z = int(floor(p.z)); float fz = p.z - floor(p.z);
float u = perlin_fade(fx);
float v = perlin_fade(fy);
float w = perlin_fade(fz);
return perlin_tri_mix(
perlin_noise_grad(jenkins_hash3(uint(X), uint(Y), uint(Z)), fx, fy, fz),
perlin_noise_grad(jenkins_hash3(uint(X+1), uint(Y), uint(Z)), fx-1, fy, fz),
perlin_noise_grad(jenkins_hash3(uint(X), uint(Y+1), uint(Z)), fx, fy-1, fz),
perlin_noise_grad(jenkins_hash3(uint(X+1), uint(Y+1), uint(Z)), fx-1, fy-1, fz),
perlin_noise_grad(jenkins_hash3(uint(X), uint(Y), uint(Z+1)), fx, fy, fz-1),
perlin_noise_grad(jenkins_hash3(uint(X+1), uint(Y), uint(Z+1)), fx-1, fy, fz-1),
perlin_noise_grad(jenkins_hash3(uint(X), uint(Y+1), uint(Z+1)), fx, fy-1, fz-1),
perlin_noise_grad(jenkins_hash3(uint(X+1), uint(Y+1), uint(Z+1)), fx-1, fy-1, fz-1),
u, v, w);
}
float perlin_signed_3d(vec3 p) {
// Precision correction for very large coordinates (matches Blender noise.cc).
// For |x| < 1000000, this is a no-op (c=0, fmod returns x unchanged).
// CRITICAL: Use sign-preserving fmod, NOT GLSL mod().
// GLSL mod(x, y) = x - y*floor(x/y) which wraps negatives to positive.
// C++ fmodf(x, y) preserves the sign of x. For x in [-100000, 100000],
// fmodf is identity. We replicate this with trunc-based modulo.
vec3 c = vec3(
abs(p.x) >= 1000000.0 ? 0.5 : 0.0,
abs(p.y) >= 1000000.0 ? 0.5 : 0.0,
abs(p.z) >= 1000000.0 ? 0.5 : 0.0);
// trunc-based fmod: x - trunc(x/y) * y (preserves sign like C++ fmodf)
p = p - trunc(p / 100000.0) * 100000.0 + c;
return perlin_noise_3d(p) * 0.9820;
}
// Perlin fBM (matches Blender exactly, including fractional detail interpolation)
float perlin_fbm(vec3 pos, int octaves, float roughness, float lacunarity, bool do_normalize) {
float fscale = 1.0;
float amp = 1.0;
float maxamp = 0.0;
float sum = 0.0;
float detail = float(octaves);
for (int i = 0; i <= octaves; i++) {
sum += perlin_signed_3d(pos * fscale) * amp;
maxamp += amp;
amp *= roughness;
fscale *= lacunarity;
}
// Fractional detail interpolation (matches CPU + Blender)
float rmd = detail - floor(detail);
if (rmd != 0.0) {
float t = perlin_signed_3d(pos * fscale);
float sum2 = sum + t * amp;
float maxamp2 = maxamp + amp;
return do_normalize ?
(1.0 - rmd) * (0.5 * sum / maxamp + 0.5) + rmd * (0.5 * sum2 / maxamp2 + 0.5) :
(1.0 - rmd) * sum + rmd * sum2;
}
return do_normalize ? 0.5 * sum / maxamp + 0.5 : sum;
}
// Perlin Multifractal
float perlin_multi_fractal(vec3 pos, int octaves, float roughness, float lacunarity) {
float value = 1.0;
float pwr = 1.0;
for (int i = 0; i <= octaves; i++) {
value *= (pwr * perlin_signed_3d(pos) + 1.0);
pwr *= roughness;
pos *= lacunarity;
}
return value;
}
// Perlin Hybrid Multifractal
float perlin_hybrid(vec3 pos, int octaves, float roughness, float lacunarity, float offset, float gain) {
float pwr = 1.0;
float value = 0.0;
float weight = 1.0;
for (int i = 0; (weight > 0.001) && (i <= octaves); i++) {
weight = min(weight, 1.0);
float signal = (perlin_signed_3d(pos) + offset) * pwr;
pwr *= roughness;
value += weight * signal;
weight *= gain * signal;
pos *= lacunarity;
}
return value;
}
// Perlin Ridged Multifractal
float perlin_ridged(vec3 pos, int octaves, float roughness, float lacunarity, float offset, float gain) {
float pwr = roughness;
float signal = offset - abs(perlin_signed_3d(pos));
signal *= signal;
float value = signal;
for (int i = 1; i <= octaves; i++) {
pos *= lacunarity;
float weight = clamp(signal * gain, 0.0, 1.0);
signal = offset - abs(perlin_signed_3d(pos));
signal *= signal;
signal *= weight;
value += signal * pwr;
pwr *= roughness;
}
return value;
}
// Perlin Hetero Terrain
float perlin_hetero(vec3 pos, int octaves, float roughness, float lacunarity, float offset) {
float pwr = roughness;
float value = offset + perlin_signed_3d(pos);
pos *= lacunarity;
for (int i = 1; i <= octaves; i++) {
float increment = (perlin_signed_3d(pos) + offset) * pwr * value;
value += increment;
pwr *= roughness;
pos *= lacunarity;
}
return value;
}
// =========================================================================
// Erosion terrain functions â iq/runevision/clayjohn
// =========================================================================
// Hash for erosion noise (same constants as iq/runevision)
vec2 erosion_hash(vec2 x) {
const float kx = 0.3183099, ky = 0.3678794;
x = vec2(x.x * kx + ky, x.y * ky + kx);
float inner = fract(x.x * x.y * (x.x + x.y));
return vec2(
-1.0 + 2.0 * fract(16.0 * kx * inner),
-1.0 + 2.0 * fract(16.0 * ky * inner));
}
// 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
vec3 fractal_noise_2d(vec2 p, float freq, int octaves, float lac, float gain) {
vec3 n = vec3(0.0);
float nf = freq;
float na = 1.0;
for (int i = 0; i < octaves; i++) {
vec3 d = grad_noised_2d(p * nf);
n += vec3(d.x * na, d.y * na * nf, d.z * na * nf);
na *= gain;
nf *= lac;
}
return n;
}
// Gullies: 16-cell weighted cosine waves with analytical derivatives
// Returns vec3(height, dh/dx, dh/dy)
vec3 gullies_2d(vec2 p, vec2 slope) {
const float PI = 3.14159265;
vec2 side_dir = vec2(-slope.y, slope.x) * 2.0 * PI;
vec2 p_int = floor(p);
vec2 p_frac = p - p_int;
vec3 has = vec3(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 vec_from_cell = p_frac - grid_offset - random_offset;
float sqr_dist = dot(vec_from_cell, vec_from_cell);
float weight = max(0.0, exp(-sqr_dist * 2.0) - 0.01111);
weight_sum += weight;
float wave_input = dot(vec_from_cell, side_dir);
float cos_w = cos(wave_input);
float sin_w = sin(wave_input);
has += vec3(cos_w, -sin_w * side_dir.x, -sin_w * side_dir.y) * weight;
}
}
return has / max(weight_sum, 1e-6);
}
// Multi-octave erosion with analytical derivative accumulation
// Returns the erosion DELTA (to add to input height+slope)
vec3 erosion_filter(vec2 p, vec3 has, float scale, float strength,
float slope_power, float cell_scale, int octaves, float gain, float lacunarity) {
vec3 input_has = has;
float freq = 1.0 / (scale * cell_scale);
strength *= scale;
for (int i = 0; i < octaves; i++) {
float sqr_len = has.y * has.y + has.z * has.z;
vec2 input_slope = vec2(0.0);
if (sqr_len > 1e-10) {
float factor = pow(sqr_len, 0.5 * (slope_power - 1.0));
input_slope = vec2(has.y, has.z) * factor;
}
vec3 g = gullies_2d(p * freq, input_slope * cell_scale);
has += vec3(g.x * strength, g.y * strength * freq, g.z * strength * freq);
strength *= gain;
freq *= lacunarity;
}
return has - input_has;
}
float magnitude_sum(int octaves, float gain) {
return (1.0 - pow(gain, float(octaves))) / (1.0 - gain);
}
float terrain_density(vec3 pos) {
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.1);
float gn_height = gn_result.x;
vec3 gn_slope = vec3(gn_result.y, 0.0, gn_result.z);
// Node 1: Erosion Noise
float HEIGHT_AMP = 0.25;
float EROSION_SCALE = 0.08333;
float WATER_HEIGHT = 0.465;
float EROSION_HEIGHT_OFFSET = -0.5;
float EROSION_LAC = 2.0;
vec2 erosion_uv = vec2(pos.x, pos.z);
// Prepare height + slope from upstream GradientNoise
float en_h = gn_height * HEIGHT_AMP;
vec3 en_n = vec3(en_h * 0.5 + 0.5, gn_slope.x * HEIGHT_AMP * 0.5, gn_slope.z * HEIGHT_AMP * 0.5);
// Water-level modulation of erosion strength
float en_sm = smoothstep(WATER_HEIGHT - 0.1, WATER_HEIGHT + 0.1, en_n.x);
float en_str = 0.16 * en_sm;
// Erosion with analytical derivatives
vec3 en_delta = erosion_filter(erosion_uv, en_n, EROSION_SCALE, en_str,
0.6, 1.0, int(min(5, pn_max_octaves)), 0.5, EROSION_LAC);
float en_erosion_mag = EROSION_SCALE * en_str * magnitude_sum(int(min(5, pn_max_octaves)), 0.5);
float en_offset = en_erosion_mag * EROSION_HEIGHT_OFFSET;
float en_final_height = en_n.x + en_delta.x + en_offset;
float en_erosion = (en_erosion_mag > 0.001) ? en_delta.x / en_erosion_mag : 0.0;
pos.y = en_final_height * 5.0;
// SDF: terrain_height - query_y. Positive below surface, negative above.
return pos.y - query_y;
}