Ray-caster: avoid roundoff error in quadratic equation

app-proto/inversive-display/src/inversive.frag

@@ -110,23 +110,37 @@ taggedFrag sphere_shading(vecInv v, vec3 pt, vec3 base_color, int id) {
 
 // --- ray-casting ---
 
+// if `a/b` is less than this threshold, we approximate `a*u^2 + b*u + c` by
+// the linear function `b*u + c`
+const float DEG_THRESHOLD = 1e-9;
+
+// the depths, represented as multiples of `dir`, where the line generated by
+// `dir` hits the sphere represented by `v`. if both depths are positive, the
+// smaller one is returned in the first component. if only one depth is
+// positive, it could be returned in either component
 vec2 sphere_cast(vecInv v, vec3 dir) {
     float a = -v.lt.s * dot(dir, dir);
     float b = dot(v.sp, dir);
     float c = -v.lt.t;
     
-    float scale = -b/(2.*a);
     float adjust = 4.*a*c/(b*b);
-    
     if (adjust < 1.) {
-        float offset = sqrt(1. - adjust);
-        return vec2(
-            scale * (1. - offset),
-            scale * (1. + offset)
-        );
+        // as long as `b` is non-zero, the linear approximation of
+        //
+        //   a*u^2 + b*u + c
+        //
+        // at `u = 0` will reach zero at a finite depth `u_lin`. the root of the
+        // quadratic adjacent to `u_lin` is stored in `lin_root`. if both roots
+        // have the same sign, `lin_root` will be the one closer to `u = 0`
+        float square_rect_ratio = 1. + sqrt(1. - adjust);
+        float lin_root = -(2.*c)/b / square_rect_ratio;
+        if (abs(a) > DEG_THRESHOLD * abs(b)) {
+            return vec2(lin_root, -b/(2.*a) * square_rect_ratio);
         } else {
-        // these parameters describe points behind the camera, so the
-        // corresponding fragments won't be drawn
+            return vec2(lin_root, -1.);
+        }
+    } else {
+        // the line through `dir` misses the sphere completely
         return vec2(-1., -1.);
     }
 }