From d9b6b808f795722383bc7751210b7addff703121 Mon Sep 17 00:00:00 2001
From: Bruno Herbelin <bruno.herbelin@gmail.com>
Date: Thu, 6 Jul 2023 22:20:48 +0200
Subject: [PATCH] BugFix Rendering Mask Ellipse

---
 rsc/shaders/mask_elipse.fs | 158 ++++++++++++++++++++++++++++++-------
 1 file changed, 130 insertions(+), 28 deletions(-)

diff --git a/rsc/shaders/mask_elipse.fs b/rsc/shaders/mask_elipse.fs
index d6ea066..eb8174f 100644
--- a/rsc/shaders/mask_elipse.fs
+++ b/rsc/shaders/mask_elipse.fs
@@ -15,38 +15,140 @@ uniform vec4 uv;
 uniform vec2  size;            // size of the mask area
 uniform float blur;            // percent of blur
 
-// See: http://www.iquilezles.org/www/articles/ellipsoids/ellipsoids.htm
-float sdEllipsoid( in vec2 p, in vec2 r )
+float msign(in float x) { return (x<0.0)?-1.0:1.0; }
+
+// MIT License
+// Copyright 2021 Henrik Dick
+// https://www.shadertoy.com/view/slS3Rw
+float sdEllipse(vec2 p, vec2 e)
 {
-    float k0 = length(p/r);
-    float k1 = length(p/(r*r));
-    return k0*(k0-1.0)/k1;
-}
+    float x = p.x;
+    float y = p.y;
+    float ax = abs(p.x);
+    float ay = abs(p.y);
+    float a = e.x;
+    float b = e.y;
+    float aa = e.x*e.x;
+    float bb = e.y*e.y;
 
-float sdEllipse( in vec2 p, in vec2 e )
-{
-    p = abs( p );
+    vec2 closest = vec2(0.0);
+    int iterations = 0;
 
-    if( e.x<e.y ) { p = p.yx; e = e.yx; }
-
-    vec2 r = e*e;
-    vec2 z = p/e;
-    vec2 n = r*z;
-
-    float g = dot(z,z) - 1.0;
-    float s0 = z.y - 1.0;
-    float s1 = (g<0.0) ? 0.0 : length( n )/r.y - 1.0;
-    float s = 0.0;
-    for( int i=0; i<64; i++ )
-    {
-        s = 0.5*(s0+s1);
-        vec2 ratio = n/(r.y*s+r);
-        g = dot(ratio,ratio) - 1.0;
-        if( g>0.0 ) s0=s; else s1=s;
+    // edge special case, handle as AABB
+    if (a * b <= 1e-15) {
+        closest = clamp(p, -e, e);
+        return length(closest-p);
     }
 
-    vec2 q = p*r/(r.y*s+r);
-    return length( p-q ) * (((p.y-q.y)<0.0)?-1.0:1.0);
+    // this epsilon will guarantee float precision result
+    // (error<1e-6) for degenerate cases
+    float epsilon = 1e-3;
+    float diff = bb - aa;
+    if (a < b) {
+        if (ax <= epsilon * a) {
+            if (ay * b < diff) {
+                float yc = bb * y / diff;
+                float xc = a * sqrt(1.0 - yc * yc / bb);
+                closest = vec2(xc, yc);
+                return -length(closest-p);
+            }
+            closest = vec2(x, b * msign(y));
+            return ay - b;
+        }
+        else if (ay <= epsilon * b) {
+            closest = vec2(a * msign(x), y);
+            return ax - a;
+        }
+    }
+    else {
+        if (ay <= epsilon * b) {
+            if (ax * a < -diff) {
+                float xc = aa * x / -diff;
+                float yc = b * sqrt(1.0 - xc * xc / aa);
+                closest = vec2(xc, yc);
+                return -length(closest-p);
+            }
+            closest = vec2(a * msign(x), y);
+            return ax - a;
+        }
+        else if (ax <= epsilon * a) {
+            closest = vec2(x, b * msign(y));
+            return ay - b;
+        }
+    }
+
+    float rx = x / a;
+    float ry = y / b;
+    float inside = rx*rx + ry*ry - 1.0;
+
+    // get lower/upper bound for parameter t
+    float s2 = sqrt(2.0);
+    float tmin = max(a * ax - aa, b * ay - bb);
+    float tmax = max(s2 * a * ax - aa, s2 * b * ay - bb);
+
+    float xx = x * x * aa;
+    float yy = y * y * bb;
+    float rxx = rx * rx;
+    float ryy = ry * ry;
+    float t;
+    if (inside < 0.0) {
+        tmax = min(tmax, 0.0);
+        if (ryy < 1.0)
+            tmin = max(tmin, sqrt(xx / (1.0 - ryy)) - aa);
+        if (rxx < 1.0)
+            tmin = max(tmin, sqrt(yy / (1.0 - rxx)) - bb);
+        t = tmin * 0.95;
+    }
+    else {
+        tmin = max(tmin, 0.0);
+        if (ryy < 1.0)
+            tmax = min(tmax, sqrt(xx / (1.0 - ryy)) - aa);
+        if (rxx < 1.0)
+            tmax = min(tmax, sqrt(yy / (1.0 - rxx)) - bb);
+        t = tmin;//2.0 * tmin * tmax / (tmin + tmax);
+    }
+    t = clamp(t, tmin, tmax);
+
+    int newton_steps = 12;
+    if (tmin >= tmax) {
+        t = tmin;
+        newton_steps = 0;
+    }
+
+    // iterate, most of the time 3 iterations are sufficient.
+    // bisect/newton hybrid
+    int i;
+    for (i = 0; i < newton_steps; i++) {
+        float at = aa + t;
+        float bt = bb + t;
+        float abt = at * bt;
+        float xxbt = xx * bt;
+        float yyat = yy * at;
+
+        float f0 = xxbt * bt + yyat * at - abt * abt;
+        float f1 = 2.0 * (xxbt + yyat - abt * (bt + at));
+        // bisect
+        if (f0 < 0.0)
+            tmax = t;
+        else if (f0 > 0.0)
+            tmin = t;
+        // newton iteration
+        float newton = f0 / abs(f1);
+        newton = clamp(newton, tmin-t, tmax-t);
+        newton = min(newton, a*b*2.0);
+        t += newton;
+
+        float absnewton = abs(newton);
+        if (absnewton < 1e-6 * (abs(t) + 0.1) || tmin >= tmax)
+            break;
+    }
+    iterations = i;
+
+    closest = vec2(x * a / (aa + t), y * b / (bb + t));
+    // this normalization is a tradeoff in precision types
+    closest = normalize(closest);
+    closest *= e;
+    return length(closest-p) * msign(inside);
 }
 
 void main()
@@ -56,7 +158,7 @@ void main()
 
     float d = sdEllipse( uv, vec2(size.x * iResolution.x/iResolution.y, size.y)  );
 
-    vec3 col = vec3(1.0- sign(d));
+    vec3 col = vec3(1.0- msign(d));
     col *= 1.0 - exp( -600.0/ (blur * 1000.f + 1.0) * abs(d));
 
     FragColor = vec4( col, 1.0 );