Fix smooth_min

The old smooth_min(x, k) had the properties
  smooth_min(1, k)   = k
  smooth_min(inf, k) = 1

Now it's smooth_min(x, k, j) with the properties
  smooth_min(1, k, j)   = k
  smooth_min(inf, k, j) = j

1 files changed, 36 insertions, 24 deletions

diff --git a/tone.cginc b/tone.cginc
index 40965c5..8558cb8 100644
--- a/tone.cginc
+++ b/tone.cginc
@@ -21,33 +21,45 @@ float3 aces_filmic(float3 x) {
 // Nice properties:
 //  1. At x=0, the derivative is 1.
 //  2. No transcendental ops, and branchless.
-float3 smooth_min(float3 x, float k) {
-  // Derivation of `b` from `k`:
-  //  f(x, b) = b * x / (x + b)
-  //  We want f(1, b) = k.
-  //  In other words, we want the max value the function can take on [0, 1] to
-  //  be k.
-  //   k = f(1, b)
-  //     = b / (1 + b)
-  //   b = k * (1 + b)
-  //     = k + kb
-  //   1 = k/b + k
-  //   1 - k = k/b
-  //   1/(1-k) = b/k
-  //   b = k/(1-k)
-  float e = 1E-4;
-  k = min(1-e, k);
-  float b = k/(1-k);
-  return b * x / (x + b);
+
+// This original attempt at a smooth minimum function has a problem:
+//   f(x, k) = k * x / (x + k)
+// As k -> inf, f(x, k) -> 1. We want f(x, k) -> k.
+// Claude suggests this:
+//  f(x, a, j) = j * (1 + a) * x / (1 + ax)
+// At x=1:
+//  f(1, a, j) = j * (1 + a) / (1 + a)
+//             = j
+// At infinity, we know that:
+//  b * x / (x + b) -> 1
+// So:
+// f(x, a, j) = j * (1 + a) * x / (1 + ax)
+//            = (j + ja) * x / (1 + ax)
+//            = (jx + jax) / (1 + ax)
+//            = jx / (1 + ax) + jax / (1 + ax)
+//            = j (x / (1 + ax) + ax / (1 + ax))
+// At infinity, this becomes:
+//              j (1/a + 1)
+// So if we want the limit to be k:
+//  k       = j (1/a + 1)
+//  k/j     = 1/a + 1
+//  k/j - 1 = 1/a
+//  a       = 1 / (k/j - 1)
+
+// Smooth, analytic min function.
+// Guarantees that x <= k for all positive x. At x=1, returns j.
+// Caller must ensure that j < k.
+float smooth_min(float x, float j, float k) {
+  float a = 1 / (k / j - 1);
+  return j * (1 + a) * x / (1 + a * x);
 }
-float smooth_min(float x, float k) {
-  float e = 1E-4;
-  k = min(1-e, k);
-  float b = k/(1-k);
-  return b * x / (x + b);
+float3 smooth_min(float3 x, float j, float k) {
+  float a = 1 / (k / j - 1);
+  return j * (1 + a) * x / (1 + a * x);
 }
+
 float smooth_clamp(float x, float lo, float hi) {
-  return smooth_max(smooth_min(x, hi), lo);
+  return smooth_max(smooth_min(x, (lo + (hi - lo)/2), hi), lo);
 }
 
 #endif  // __TONE_INC