Adding root polishing step to solve_quartic function.

Newton-Raphson step added to polish initial roots found by the solve_quartic
function. The core solve_quartic tolerance allows roots to drift from <0.0
to >0.0 values with the latter causing artifacts and additional root
filtering. This the very likely the reason for the long too large 1e-2
intersection depth value in blob.cpp now reduced to
MIN_ISECT_DEPTH_RETURNED.

As part of this change created a new FUDGE_FACTOR4(1e-8) constant DBL value
to replace previous use of SMALL_ENOUGH(1e-10) within solve_quartic. Looks
like the value had been smaller to make roots more accurate, but at the cost
of missing roots in some difficult equation cases. With the root polishing
can move back to a larger value so as to always get roots. Yes, this better
addresses most of what the already remove difficult_coeffs() function and
bump into polysolve was trying to do.

Author: wfpokorny

source/core/math/polynomialsolver.cpp

@@ -85,6 +85,15 @@ const DBL FUDGE_FACTOR2 = -1.0e-5;
 ///
 const DBL FUDGE_FACTOR3 = 1.0e-7;
 
+/// @var FUDGE_FACTOR4
+/// @brief const DBL value used in the active solve_quartic() function.
+///
+/// Roughly acts as @ref FUDGE_FACTOR2 and @ref FUDGE_FACTOR3 did for the
+/// original solve_quartic() versions but for the current solve_quartic() code.
+/// Value arrived at by running many scenes and settling on what worked best.
+///
+const DBL FUDGE_FACTOR4 = 1.0e-8;
+
 /// @var TWO_M_PI_3
 /// const DBL value used in solve_cubic() equal to 2.0 * pi / 3.0.
 ///
@@ -1391,7 +1400,7 @@ static int solve_quartic(const DBL *x, DBL *results)
 
     if (d1 < 0.0)
     {
-        if (d1 > -SMALL_ENOUGH)
+        if (d1 > -FUDGE_FACTOR4)
         {
             d1 = 0.0;
         }
@@ -1401,7 +1410,7 @@ static int solve_quartic(const DBL *x, DBL *results)
         }
     }
 
-    if (d1 < SMALL_ENOUGH)
+    if (d1 < FUDGE_FACTOR4)
     {
         d2 = z * z - r;
 
@@ -1461,6 +1470,38 @@ static int solve_quartic(const DBL *x, DBL *results)
         }
     }
 
+    // The quartic root finding method above somewhat inaccurate and more so at
+    // large scales. The follow code polishes the found roots so, for example,
+    // negative roots which otherwise drift positive and cause artifacts are
+    // corrected. A reason for the historically too large 1e-2 minimum intersection
+    // depth used with blobs that itself caused artifacts of other kinds. Further,
+    // the blob 1e-2 value was not large enough to filter all drift to 0+ roots in
+    // any case. Newton-Raphson method.
+    //
+    DBL pv, dpv, t, dt;
+    for (int c = 0; c < i; c++)
+    {
+        t = results[c];
+        for (int j = 0; j < 7; j++)
+        {
+            pv  = x[0] * t + x[1];
+            dpv = x[0];
+            for (int k=2; k<=4; k++)
+            {
+                dpv = dpv * t + pv;
+                pv  = pv * t + x[k];
+            }
+
+            dt = pv / dpv;
+            t -= dt;
+            if (fabs(dt) < SMALL_ENOUGH)
+            {
+                results[c] = t;
+                break;
+            }
+        }
+    }
+
     return(i);
 }
 #endif