Add faithful advanced erosion filter core

README.md

@@ -66,6 +66,11 @@ local `AetheriaMath` surface. The value-noise family follows public shader
 literature from Inigo Quilez's noise and mini-spline articles, expressed in
 CultMath names with tests instead of copied project-local helper drift.
 
+The advanced erosion filter is an MPL-2.0 adaptation of Rune Skovbo Johansen's
+published filter. Its faithful planar core is isolated in
+`AdvancedErosionFilter.cs` and `AdvancedErosionFilter.hlsl`; spherical mapping
+is a separate CultMath concern. See `THIRD-PARTY-NOTICES.md`.
+
 ## Build
 
 ```powershell

THIRD-PARTY-NOTICES.md

@@ -0,0 +1,21 @@
+# Third-Party Notices
+
+## Advanced Terrain Erosion Filter
+
+`src/CultMath/AdvancedErosionFilter.cs` and
+`shaders/AdvancedErosionFilter.hlsl` adapt the mathematical procedure from
+Rune Skovbo Johansen's **Advanced Terrain Erosion Filter**:
+
+- Copyright 2025 Rune Skovbo Johansen.
+- Upstream: <https://www.shadertoy.com/view/wXcfWn>
+- Explanation: <https://blog.runevision.com/2026/03/fast-and-gorgeous-erosion-filter.html>
+- License: Mozilla Public License 2.0,
+  <https://www.mozilla.org/MPL/2.0/>
+
+The CultMath files retain the MPL notice and separate the faithful planar
+filter from CultMath's spherical adaptation. Modifications include C#/HLSL
+type lowering, explicit bounded octave counts, finite zero-rounding behavior,
+and CultMath naming.
+
+The technique builds on earlier work credited by the upstream author to Fewes
+and Clay John. The hash form follows the public shader math of Inigo Quilez.

shaders/AdvancedErosionFilter.hlsl

@@ -0,0 +1,112 @@
+// 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/.
+// Mathematical procedure adapted from Rune Skovbo Johansen's Advanced Terrain
+// Erosion Filter (2025): https://www.shadertoy.com/view/wXcfWn
+
+#ifndef CULTMATH_ADVANCED_EROSION_FILTER_HLSL
+#define CULTMATH_ADVANCED_EROSION_FILTER_HLSL
+
+struct CultMathAdvancedErosionParameters
+{
+    float scale; float strength; float gully_weight; float detail;
+    float4 rounding; float4 onset; float2 assumed_slope;
+    float cell_scale; float normalization; int octaves; float lacunarity; float gain;
+};
+
+struct CultMathAdvancedErosionResult
+{
+    float3 delta; float magnitude; float ridge_map; float fade_target;
+};
+
+float2 cultmath_erosion_hash2(float2 input)
+{
+    float2 k = float2(0.3183099, 0.3678794);
+    float2 x = input * k + k.yx;
+    return -1.0 + 2.0 * frac(16.0 * k * frac(x.x * x.y * (x.x + x.y)));
+}
+
+float4 cultmath_phacelle_noise(float2 position, float2 normal_direction, float frequency, float offset_cycles, float normalization)
+{
+    float2 side_direction = float2(-normal_direction.y, normal_direction.x) * frequency * 6.28318530718;
+    float offset = offset_cycles * 6.28318530718;
+    float2 cell = floor(position);
+    float2 local = frac(position);
+    float2 phase = 0.0;
+    float weight_sum = 0.0;
+    [unroll] for (int x = -1; x <= 2; x++)
+    [unroll] for (int y = -1; y <= 2; y++)
+    {
+        float2 grid_offset = float2(x, y);
+        float2 grid_point = cell + grid_offset;
+        float2 random_offset = cultmath_erosion_hash2(grid_point) * 0.5;
+        float2 delta = local - grid_offset - random_offset;
+        float weight = max(0.0, exp(-dot(delta, delta) * 2.0) - 0.01111);
+        weight_sum += weight;
+        float wave = dot(delta, side_direction) + offset;
+        phase += float2(cos(wave), sin(wave)) * weight;
+    }
+    float2 interpolated = phase / max(weight_sum, 1.0e-20);
+    float magnitude = max(1.0 - saturate(normalization), length(interpolated));
+    return float4(interpolated / max(magnitude, 1.0e-20), side_direction);
+}
+
+float cultmath_erosion_pow_inverse(float value, float power) { return 1.0 - pow(1.0 - saturate(value), power); }
+float cultmath_erosion_ease_out(float value) { float inverse = 1.0 - saturate(value); return 1.0 - inverse * inverse; }
+float cultmath_erosion_smooth_start(float value, float smoothing)
+{
+    return smoothing <= 1.0e-10 || value >= smoothing ? value - 0.5 * smoothing : 0.5 * value * value / smoothing;
+}
+float2 cultmath_erosion_safe_normalize(float2 value)
+{
+    float value_length = length(value);
+    return abs(value_length) > 1.0e-10 ? value / value_length : value;
+}
+
+CultMathAdvancedErosionResult cultmath_advanced_erosion_filter(float2 position, float3 base_height_and_slope, float fade_target, CultMathAdvancedErosionParameters p)
+{
+    float strength = p.strength * p.scale;
+    fade_target = clamp(fade_target, -1.0, 1.0);
+    float3 input_value = base_height_and_slope;
+    float3 height_and_slope = input_value;
+    float frequency = 1.0 / max(p.scale * p.cell_scale, 1.0e-10);
+    float2 slope = height_and_slope.yz;
+    float slope_length = max(length(slope), 1.0e-10);
+    float magnitude = 0.0;
+    float rounding_multiplier = 1.0;
+    float input_rounding = lerp(p.rounding.y, p.rounding.x, saturate(fade_target + 0.5)) * p.rounding.z;
+    float combined_mask = cultmath_erosion_ease_out(cultmath_erosion_smooth_start(slope_length * p.onset.x, input_rounding * p.onset.x));
+    float ridge_mask = cultmath_erosion_ease_out(slope_length * p.onset.z);
+    float ridge_fade = fade_target;
+    float2 gully_slope = lerp(slope, slope / slope_length * p.assumed_slope.x, p.assumed_slope.y);
+
+    [loop] for (int octave = 0; octave < min(max(p.octaves, 0), 8); octave++)
+    {
+        float4 phacelle = cultmath_phacelle_noise(position * frequency, cultmath_erosion_safe_normalize(gully_slope), p.cell_scale, 0.25, p.normalization);
+        float2 derivative_direction = phacelle.zw * -frequency;
+        float sloping = abs(phacelle.y);
+        gully_slope += sign(phacelle.y) * derivative_direction * strength * p.gully_weight;
+        float3 gullies = float3(phacelle.x, phacelle.y * derivative_direction.x, phacelle.y * derivative_direction.y);
+        float3 faded = lerp(float3(fade_target, 0.0, 0.0), gullies * p.gully_weight, combined_mask);
+        height_and_slope += faded * strength;
+        magnitude += strength;
+        fade_target = faded.x;
+        float octave_rounding = lerp(p.rounding.y, p.rounding.x, saturate(phacelle.x + 0.5)) * rounding_multiplier;
+        float new_mask = cultmath_erosion_ease_out(cultmath_erosion_smooth_start(sloping * p.onset.y, octave_rounding * p.onset.y));
+        combined_mask = cultmath_erosion_pow_inverse(combined_mask, p.detail) * new_mask;
+        ridge_fade = lerp(ridge_fade, gullies.x, ridge_mask);
+        ridge_mask *= cultmath_erosion_ease_out(sloping * p.onset.w);
+        strength *= p.gain;
+        frequency *= p.lacunarity;
+        rounding_multiplier *= p.rounding.w;
+    }
+
+    CultMathAdvancedErosionResult result;
+    result.delta = height_and_slope - input_value;
+    result.magnitude = magnitude;
+    result.ridge_map = ridge_fade * (1.0 - ridge_mask);
+    result.fade_target = fade_target;
+    return result;
+}
+
+#endif

shaders/CultMath.hlsl

@@ -263,6 +263,8 @@ float cultmath_value_noise_texture_bicubic(Texture2D<float> noise_texture, Sampl
     return cultmath_catmullrom(row0, row1, row2, row3, local.y);
 }
 
+#include "AdvancedErosionFilter.hlsl"
+
 struct CultMathSphericalErosionParameters
 {
     float frequency; float amplitude; float lacunarity; float gain;

src/CultMath/AdvancedErosionFilter.cs

@@ -0,0 +1,129 @@
+// 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/.
+//
+// Mathematical procedure adapted from Rune Skovbo Johansen's Advanced Terrain
+// Erosion Filter (2025): https://www.shadertoy.com/view/wXcfWn
+
+namespace CultMath;
+
+public readonly record struct AdvancedErosionParameters(
+    float Scale,
+    float Strength,
+    float GullyWeight,
+    float Detail,
+    float4 Rounding,
+    float4 Onset,
+    float2 AssumedSlope,
+    float CellScale,
+    float Normalization,
+    int Octaves,
+    float Lacunarity,
+    float Gain)
+{
+    public static AdvancedErosionParameters Default => new(
+        0.15f, 0.22f, 0.5f, 1.5f,
+        new float4(0.1f, 0.0f, 0.1f, 2.0f),
+        new float4(1.25f, 1.25f, 2.8f, 1.5f),
+        new float2(0.7f, 1.0f),
+        0.7f, 0.5f, 5, 2.0f, 0.5f);
+}
+
+public readonly record struct AdvancedErosionResult(float3 Delta, float Magnitude, float RidgeMap, float FadeTarget);
+
+public static class AdvancedErosionFilter
+{
+    public static AdvancedErosionResult Sample(float2 position, float3 baseHeightAndSlope, float fadeTarget, AdvancedErosionParameters p)
+    {
+        var strength = p.Strength * p.Scale;
+        fadeTarget = math.clamp(fadeTarget, -1.0f, 1.0f);
+        var input = baseHeightAndSlope;
+        var heightAndSlope = input;
+        var frequency = 1.0f / math.max(p.Scale * p.CellScale, 1.0e-10f);
+        var slope = new float2(heightAndSlope.y, heightAndSlope.z);
+        var slopeLength = math.max(math.length(slope), 1.0e-10f);
+        var magnitude = 0.0f;
+        var roundingMultiplier = 1.0f;
+
+        var inputRounding = math.lerp(p.Rounding.y, p.Rounding.x, math.saturate(fadeTarget + 0.5f)) * p.Rounding.z;
+        var combinedMask = EaseOut(SmoothStart(slopeLength * p.Onset.x, inputRounding * p.Onset.x));
+        var ridgeMask = EaseOut(slopeLength * p.Onset.z);
+        var ridgeFade = fadeTarget;
+        var gullySlope = math.lerp(slope, slope / slopeLength * p.AssumedSlope.x, p.AssumedSlope.y);
+
+        var octaves = math.min(math.max(p.Octaves, 0), 8);
+        for (var octave = 0; octave < octaves; octave++)
+        {
+            var phacelle = Phacelle(position * frequency, SafeNormalize(gullySlope), p.CellScale, 0.25f, p.Normalization);
+            var derivativeDirection = new float2(phacelle.z, phacelle.w) * -frequency;
+            var sloping = math.abs(phacelle.y);
+            gullySlope += math.sign(phacelle.y) * derivativeDirection * strength * p.GullyWeight;
+
+            var gullies = new float3(phacelle.x, phacelle.y * derivativeDirection.x, phacelle.y * derivativeDirection.y);
+            var faded = math.lerp(new float3(fadeTarget, 0.0f, 0.0f), gullies * p.GullyWeight, combinedMask);
+            heightAndSlope += faded * strength;
+            magnitude += strength;
+            fadeTarget = faded.x;
+
+            var octaveRounding = math.lerp(p.Rounding.y, p.Rounding.x, math.saturate(phacelle.x + 0.5f)) * roundingMultiplier;
+            var newMask = EaseOut(SmoothStart(sloping * p.Onset.y, octaveRounding * p.Onset.y));
+            combinedMask = PowInverse(combinedMask, p.Detail) * newMask;
+
+            ridgeFade = math.lerp(ridgeFade, gullies.x, ridgeMask);
+            ridgeMask *= EaseOut(sloping * p.Onset.w);
+            strength *= p.Gain;
+            frequency *= p.Lacunarity;
+            roundingMultiplier *= p.Rounding.w;
+        }
+
+        return new AdvancedErosionResult(heightAndSlope - input, magnitude, ridgeFade * (1.0f - ridgeMask), fadeTarget);
+    }
+
+    public static float4 Phacelle(float2 position, float2 normalDirection, float frequency, float offsetCycles, float normalization)
+    {
+        var sideDirection = new float2(-normalDirection.y, normalDirection.x) * frequency * math.TAU;
+        var offset = offsetCycles * math.TAU;
+        var cell = math.floor(position);
+        var local = math.frac(position);
+        var phase = float2.zero;
+        var weightSum = 0.0f;
+
+        for (var x = -1; x <= 2; x++)
+        for (var y = -1; y <= 2; y++)
+        {
+            var gridOffset = new float2(x, y);
+            var gridPoint = cell + gridOffset;
+            var randomOffset = Hash2(gridPoint) * 0.5f;
+            var delta = local - gridOffset - randomOffset;
+            var weight = math.max(0.0f, math.exp(-math.dot(delta, delta) * 2.0f) - 0.01111f);
+            weightSum += weight;
+            var wave = math.dot(delta, sideDirection) + offset;
+            phase += new float2(math.cos(wave), math.sin(wave)) * weight;
+        }
+
+        var interpolated = phase / math.max(weightSum, 1.0e-20f);
+        var magnitude = math.max(1.0f - math.saturate(normalization), math.length(interpolated));
+        var normalized = interpolated / math.max(magnitude, 1.0e-20f);
+        return new float4(normalized, sideDirection.x, sideDirection.y);
+    }
+
+    private static float2 Hash2(float2 input)
+    {
+        var k = new float2(0.3183099f, 0.3678794f);
+        var x = input * k + new float2(k.y, k.x);
+        var scalar = math.frac(x.x * x.y * (x.x + x.y));
+        return -float2.one + 2.0f * math.frac(16.0f * k * scalar);
+    }
+
+    private static float2 SafeNormalize(float2 value)
+    {
+        var length = math.length(value);
+        return math.abs(length) > 1.0e-10f ? value / length : value;
+    }
+
+    private static float PowInverse(float value, float power) => 1.0f - math.pow(1.0f - math.saturate(value), power);
+    private static float EaseOut(float value) { var inverse = 1.0f - math.saturate(value); return 1.0f - inverse * inverse; }
+    private static float SmoothStart(float value, float smoothing) => smoothing <= 1.0e-10f || value >= smoothing
+        ? value - 0.5f * smoothing
+        : 0.5f * value * value / smoothing;
+}

tests/CultMath.Tests/AdvancedErosionFilterTests.cs

@@ -0,0 +1,44 @@
+using Xunit;
+
+namespace CultMath.Tests;
+
+public sealed class AdvancedErosionFilterTests
+{
+    [Fact]
+    public void DefaultReferenceFixtureIsStable()
+    {
+        var result = AdvancedErosionFilter.Sample(
+            new float2(0.42f, 0.31f),
+            new float3(0.5f, 0.1f, -0.2f),
+            0.0f,
+            AdvancedErosionParameters.Default);
+
+        Assert.Equal(0.0639375f, result.Magnitude, precision: 7);
+        Assert.True(float.IsFinite(result.Delta.x));
+        Assert.True(float.IsFinite(result.Delta.y));
+        Assert.True(float.IsFinite(result.Delta.z));
+        Assert.InRange(result.RidgeMap, -1.0f, 1.0f);
+    }
+
+    [Fact]
+    public void FlatInputAndZeroRoundingRemainFinite()
+    {
+        var parameters = AdvancedErosionParameters.Default with { Rounding = float4.zero };
+        var result = AdvancedErosionFilter.Sample(new float2(0.25f, 0.75f), float3.zero, 0.0f, parameters);
+
+        Assert.True(float.IsFinite(result.Delta.x));
+        Assert.True(float.IsFinite(result.RidgeMap));
+        Assert.True(float.IsFinite(result.FadeTarget));
+    }
+
+    [Fact]
+    public void PhacellePartialNormalizationIsFiniteAtEndpoints()
+    {
+        foreach (var normalization in new[] { 0.0f, 0.5f, 1.0f })
+        {
+            var sample = AdvancedErosionFilter.Phacelle(new float2(1.7f, -3.2f), new float2(0.6f, -0.8f), 0.7f, 0.25f, normalization);
+            Assert.True(float.IsFinite(sample.x));
+            Assert.True(float.IsFinite(sample.y));
+        }
+    }
+}

tests/CultMath.Tests/HlslMirrorTests.cs

@@ -7,11 +7,22 @@ public sealed class HlslMirrorTests
     [Fact]
     public void HlslMirrorPublishesSphericalErosionKernel()
     {
-        var source = File.ReadAllText(FindHlslMirror());
+        var source = File.ReadAllText(Path.Combine(GetRepositoryRoot(), "shaders", "CultMath.hlsl"));
         Assert.Contains("cultmath_spherical_erosion", source);
         Assert.Contains("CultMathSphericalErosionParameters", source);
     }
 
+    [Fact]
+    public void HlslMirrorPublishesFaithfulAdvancedErosionKernel()
+    {
+        var mirror = Path.Combine(GetRepositoryRoot(), "shaders", "CultMath.hlsl");
+        var source = File.ReadAllText(mirror);
+        Assert.Contains("AdvancedErosionFilter.hlsl", source);
+        var advanced = File.ReadAllText(Path.Combine(Path.GetDirectoryName(mirror)!, "AdvancedErosionFilter.hlsl"));
+        Assert.Contains("cultmath_phacelle_noise", advanced);
+        Assert.Contains("cultmath_advanced_erosion_filter", advanced);
+    }
+
     [Fact]
     public void HlslMirrorPublishesCultMathPrimitives()
     {

tests/CultMath.Tests/SphericalErosionTests.cs

@@ -1,3 +1,5 @@
+using Xunit;
+
 namespace CultMath.Tests;
 
 public sealed class SphericalErosionTests