Find optimal shift value to decompose the spline

The previous method for decomposing a spline into fewer segments aimed
to identify the optimal value of 't', which corresponds to the location
of maximum curvature on the curve. However, using a bisection method to
locate this point increases the average number of '_de_casteljau' calls
per point, as demonstrated in Experiment 2. Even with the early stopping
method applied, the improvement is little.

Moreover, given the high cost of finding points with the bisection
method, we have decided to abandon the flexible update of 't' using
bisection. Therefore, we can reduce computation cost by only using
shift operation to determine the coefficients of lerp in the
'_de_casteljau' function.

Based on the previous experiments, we have identified that the primary
issue is the necessity to call the '_de_casteljau' function twice to
generate a point, even under the original settings. On average, the
number of _de_casteljau calls per point is 1.99. Ideally, we want the
rendering process to call the '_de_casteljau' function only once for
each point in every iteration, indicating that the first call is
redundant.

Also, we've observed that during the initial 60 percent of the rendering
process, the original method on average requires more than two shift
attempts to determine the optimal value in every iteration. Therefore,
we propose adjusting the initial 't' value from 0.5 to 0.25 by applying
an initial shift of 2.

Additionally, as the spline rendering process progresses to the later
stages, the amount of shift gradually decreases. Hence, we store the
amount of shift used in each iteration as a global variable to record
the value of the last shift change. This way allows us to use this value
directly in the next iteration, eliminating the need to start from an
initial shift of 2 again.

Furthermore, instead of merely decreasing 't' by adding the amount of
shift, we also reduce the amount of shift to a minimum of 1. Based on
Experiment 3, this limitation of the scope can decrease the average
number of function '_de_casteljau' calls per point and increase the
number of points when rendering. This means that these points do not
always have the maximum curvature, and some optimal points with maximum
curvature will be overlooked while limiting the scope of 't' to [0,
0.25].

In summary, we have chosen the Original (shift2) setting in this pull
request, as demonstrated in Experiment 1.

The modified implementation of font-edit, which utilizes fixed-point
arithmetic, serves as the evaluation testbed for the following
experiments:

Experiment 1:
Original - Average number of _de_casteljau calls per point: 1.99
Original - Average points per character: 18.89
Original (shift2) - Average number of _de_casteljau calls per point:
1.51
Original (shift2) - Average points per character: 18.98

Experiment 2:
Flexible - Average number of _de_casteljau calls per point: 4.53
Flexible - Average points per character: 16.30
Flexible (shift2) - Average number of _de_casteljau calls per point:
4.40
Flexible (shift2) - Average points per character: 21.16
Flexible (early stopping) - Average number of _de_casteljau calls per
point: 4.23
Flexible (early stopping) - Average points per character: 16.18
Flexible (early stopping) (shift2) - Average number of _de_casteljau
calls per point: 3.99
Flexible (early stopping) (shift2) - Average points per character: 21.09

Experiment 3:
Original (shift2) (limit scope) - Average number of _de_casteljau calls
per point: 1.18
Original (shift2) (limit scope) - Average points per character: 22.57

src/spline.c

@@ -11,8 +11,8 @@ 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 'shift' factor.
+ * The 'shift' factor determines the position between 'a' and 'b'.
  * The result is stored in 'result'.
  */
 static void _lerp(const twin_spoint_t *a,
@@ -25,17 +25,15 @@ static void _lerp(const twin_spoint_t *a,
 }
 
 /*
- * 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 'shift'
  * factor. The spline is split into two new splines 's1' and 's2'.
  */
 static void _de_casteljau(twin_spline_t *spline,
                           int shift,
                           twin_spline_t *s1,
                           twin_spline_t *s2)
 {
-    twin_spoint_t ab, bc, cd;
-    twin_spoint_t abbc, bccd;
-    twin_spoint_t final;
+    twin_spoint_t ab, bc, cd, abbc, bccd, final;
 
     _lerp(&spline->a, &spline->b, shift, &ab);
     _lerp(&spline->b, &spline->c, shift, &bc);
@@ -56,28 +54,31 @@ static void _de_casteljau(twin_spline_t *spline,
 }
 
 /*
- * Return an upper bound on the error (squared) that could result from
- * approximating a spline with a line segment connecting the two endpoints.
+ * Return an upper bound on the distance (squared) that could result from
+ * approximating a spline with a line segment connecting the two endpoints,
+ * which is based on the Convex Hull Property of Bézier Curves: The Bézier Curve
+ * lies completely in the convex hull of the given control points. Therefore, we
+ * can use control points B and C to approximate the actual spline.
  */
-static twin_dfixed_t _twin_spline_error_squared(twin_spline_t *spline)
+static twin_dfixed_t _twin_spline_distance_squared(twin_spline_t *spline)
 {
-    twin_dfixed_t berr, cerr;
+    twin_dfixed_t bdist, cdist;
 
-    berr = _twin_distance_to_line_squared(&spline->b, &spline->a, &spline->d);
-    cerr = _twin_distance_to_line_squared(&spline->c, &spline->a, &spline->d);
+    bdist = _twin_distance_to_line_squared(&spline->b, &spline->a, &spline->d);
+    cdist = _twin_distance_to_line_squared(&spline->c, &spline->a, &spline->d);
 
-    if (berr > cerr)
-        return berr;
-    return cerr;
+    if (bdist > cdist)
+        return bdist;
+    return cdist;
 }
 
 /*
- * Check if a spline is flat enough by comparing the error against the
+ * Check if a spline is flat enough by comparing the distance against the
  * tolerance.
  */
 static bool is_flat(twin_spline_t *spline, twin_dfixed_t tolerance_squared)
 {
-    return _twin_spline_error_squared(spline) <= tolerance_squared;
+    return _twin_spline_distance_squared(spline) <= tolerance_squared;
 }
 
 /*
@@ -91,16 +92,30 @@ static void _twin_spline_decompose(twin_path_t *path,
 {
     /* Draw starting point */
     _twin_path_sdraw(path, spline->a.x, spline->a.y);
-
+    /*
+     * It on average requires over two shift attempts per iteration to find the
+     * optimal value. To reduce redundancy in shift 1, adjust the initial 't'
+     * value from 0.5 to 0.25 by applying an initial shift of 2. As spline
+     * rendering progresses, the shift amount decreases. Store the last shift
+     * value as a global variable to use directly in the next iteration,
+     * avoiding a reset to an initial shift of 2.
+     */
+    int shift = 2;
     while (!is_flat(spline, tolerance_squared)) {
-        int shift = 1;
         twin_spline_t left, right;
 
-        /* FIXME: Find the optimal shift value to decompose the spline */
-        do {
+        while (true) {
             _de_casteljau(spline, shift, &left, &right);
+            if (is_flat(&left, tolerance_squared)) {
+                /* Limiting the scope of 't' may overlook optimal points with
+                 * maximum curvature. Therefore, dynamically reduce the shift
+                 * amount to a minimum of 1. */
+                if (shift > 1)
+                    shift--;
+                break;
+            }
             shift++;
-        } while (!is_flat(&left, tolerance_squared));
+        }
 
         /* Draw the left segment */
         _twin_path_sdraw(path, left.d.x, left.d.y);