Share planetary radial refinement semantics

shaders/CultMath.hlsl

@@ -337,4 +337,23 @@ CultMathSphericalErosionSample cultmath_spherical_erosion(float3 unit_position,
     return result;
 }
 
+float3 cultmath_planetary_radial_refinement_step(
+    float3 position,
+    float3 ray_direction,
+    float3 center,
+    float target_radius,
+    float minimum_derivative,
+    float maximum_correction)
+{
+    float3 local=position-center;
+    float radius=length(local);
+    float3 radial=local/max(radius,0.0001);
+    float error=radius-target_radius;
+    float derivative=dot(ray_direction,radial);
+    float safe_minimum=max(minimum_derivative,0.0001);
+    float safe_derivative=abs(derivative)>safe_minimum?derivative:(derivative<0?-safe_minimum:safe_minimum);
+    float correction=clamp(-error/safe_derivative,-abs(maximum_correction),abs(maximum_correction));
+    return position+ray_direction*correction;
+}
+
 #endif

src/CultMath/PlanetaryRadialRefinement.cs

@@ -0,0 +1,28 @@
+namespace CultMath;
+
+public static class PlanetaryRadialRefinement
+{
+    public static float3 Step(
+        float3 position,
+        float3 rayDirection,
+        float3 center,
+        float targetRadius,
+        float minimumDerivative = 0.2f,
+        float maximumCorrection = 0.08f)
+    {
+        var local = position - center;
+        var radius = math.length(local);
+        var radial = local / math.max(radius, 0.0001f);
+        var error = radius - targetRadius;
+        var derivative = math.dot(rayDirection, radial);
+        var safeMinimum = math.max(minimumDerivative, 0.0001f);
+        var safeDerivative = math.abs(derivative) > safeMinimum
+            ? derivative
+            : derivative < 0.0f ? -safeMinimum : safeMinimum;
+        var correction = math.clamp(
+            -error / safeDerivative,
+            -math.abs(maximumCorrection),
+            math.abs(maximumCorrection));
+        return position + rayDirection * correction;
+    }
+}

tests/CultMath.Tests/PlanetaryRadialRefinementTests.cs

@@ -0,0 +1,29 @@
+using Xunit;
+
+namespace CultMath.Tests;
+
+public sealed class PlanetaryRadialRefinementTests
+{
+    [Fact]
+    public void FourBoundedStepsConvergeFromTheRasterCaptureShell()
+    {
+        var center=new float3(0.3f,-0.2f,0.7f);
+        var radial=math.normalize(new float3(0.2f,-0.4f,0.9f));
+        var targetRadius=3.85f;
+        var target=center+radial*targetRadius;
+        var origin=target+math.normalize(radial+new float3(0.12f,-0.05f,0.02f))*6.0f;
+        var position=center+radial*(targetRadius+0.30f);
+        var ray=math.normalize(position-origin);
+        for(var step=0;step<4;step++)
+            position=PlanetaryRadialRefinement.Step(position,ray,center,targetRadius);
+        Assert.InRange(math.abs(math.length(position-center)-targetRadius),0.0f,1.0e-4f);
+    }
+
+    [Fact]
+    public void StepClampsWeakDerivativeCorrections()
+    {
+        var position=new float3(0,0,4.1f);
+        var result=PlanetaryRadialRefinement.Step(position,new float3(1,0,0),float3.zero,3.8f);
+        Assert.InRange(math.length(result-position),0.07999f,0.08001f);
+    }
+}