Decompose splines with t flexibly updated

This commit focuses on decomposing a spline into fewer segments by
finding an optimal value of t that meets our target tolerance for
flatness while ensuring it is at least as large as the maximum distance
from the line segment connecting the starting and ending points of the
original curve.

The approach employs the bisection method to determine the value of t.
The idea behind this method is based on the properties of curves. By
choosing the point of maximum curvature as the junction between two line
segments during curve fitting, we ensure that the curvature of all
points on both segments is less than the maximum curvature of the
original curve. This strategy guarantees that the two line segments fit
the original curve as well as possible.

Additionally, the properties of these two segments are such that the
slope of the tangent line at all points in one segment is relatively
positive compared to the slope of the line connecting the starting and
ending points of the original curve, while the slope in the other
line segment is relatively negative. Consequently, the point where the
tangent line's slope is zero relative to the connecting line results in
the maximum distance.

This distance is recorded to ensure we maximize it as much as possible.
Although this approach requires more computation to decompose the
spline into fewer segments, the extra effort is justified for improved
fit.

See https://keithp.com/blogs/more-iterative-splines/ for details.
See also: sysprog21#2

src/spline.c

@@ -11,38 +11,38 @@ typedef struct _twin_spline {
 } twin_spline_t;
 
 /*
- * Linearly interpolate between points 'a' and 'b' with a shift factor.
- * The shift factor determines the position between 'a' and 'b'.
+ * Linearly interpolate between points 'a' and 'b' with a 't' factor.
+ * The 't' factor determines the position between 'a' and 'b'.
  * The result is stored in 'result'.
  */
 static void _lerp(twin_spoint_t *a,
                   twin_spoint_t *b,
-                  int shift,
+                  twin_dfixed_t t,
                   twin_spoint_t *result)
 {
-    result->x = a->x + ((b->x - a->x) >> shift);
-    result->y = a->y + ((b->y - a->y) >> shift);
+    result->x = a->x + (((twin_dfixed_t) (b->x - a->x) * t) >> 8);
+    result->y = a->y + (((twin_dfixed_t) (b->y - a->y) * t) >> 8);
 }
 
 /*
- * Perform the de Casteljau algorithm to split a spline at a given shift
+ * Perform the de Casteljau algorithm to split a spline at a given 't'
  * factor. The spline is split into two new splines 's1' and 's2'.
  */
 static void _de_casteljau(twin_spline_t *spline,
-                          int shift,
+                          twin_dfixed_t t,
                           twin_spline_t *s1,
                           twin_spline_t *s2)
 {
     twin_spoint_t ab, bc, cd;
     twin_spoint_t abbc, bccd;
     twin_spoint_t final;
 
-    _lerp(&spline->a, &spline->b, shift, &ab);
-    _lerp(&spline->b, &spline->c, shift, &bc);
-    _lerp(&spline->c, &spline->d, shift, &cd);
-    _lerp(&ab, &bc, shift, &abbc);
-    _lerp(&bc, &cd, shift, &bccd);
-    _lerp(&abbc, &bccd, shift, &final);
+    _lerp(&spline->a, &spline->b, t, &ab);
+    _lerp(&spline->b, &spline->c, t, &bc);
+    _lerp(&spline->c, &spline->d, t, &cd);
+    _lerp(&ab, &bc, t, &abbc);
+    _lerp(&bc, &cd, t, &bccd);
+    _lerp(&abbc, &bccd, t, &final);
 
     s1->a = spline->a;
     s1->b = ab;
@@ -93,14 +93,37 @@ static void _twin_spline_decompose(twin_path_t *path,
     _twin_path_sdraw(path, spline->a.x, spline->a.y);
 
     while (!is_flat(spline, tolerance_squared)) {
-        int shift = 1;
+        twin_dfixed_t hi = TWIN_SFIXED_ONE * TWIN_SFIXED_ONE, lo = 0;
+        twin_dfixed_t best = 0, t;
+        twin_dfixed_t best_error = 0, error;
         twin_spline_t left, right;
 
-        /* FIXME: Find the optimal shift value to decompose the spline */
-        do {
-            _de_casteljau(spline, shift, &left, &right);
-            shift++;
-        } while (!is_flat(&left, tolerance_squared));
+        while (true) {
+            t = (hi + lo) >> 1;
+            if (best == t)
+                break;
+
+            _de_casteljau(spline, t, &left, &right);
+
+            error = _twin_spline_error_squared(&left);
+            if (error < best_error)
+                break;
+
+            /* Left is close enough to the bezier curve */
+            if (error <= tolerance_squared) {
+                best_error = error;
+                best = t;
+                /* Try to find better t */
+                lo = t;
+            } else {
+                /* Find good t already and cannot find better t */
+                if (best)
+                    break;
+                hi = t;
+            }
+        }
+
+        _de_casteljau(spline, best, &left, &right);
 
         /* Draw the left segment */
         _twin_path_sdraw(path, left.d.x, left.d.y);