Find optimal shift value to decompose the spline

De Casteljau's algorithm is a common approach to rendering spline. There
are many studies that aim to improve the spline. Here, we only consider
three metrics to evaluate which method is more suitable. The first
metric is the number of generated points when rendering a curve. The
second metric is the number of '_de_casteljau' function calls. The third
metric is the computation cost of the '_de_casteljau' function.

In order to generate fewer points when rendering a curve, we first apply
a flexible update 't' that employs the bisection method to identify the
positions that are best suited for curve fitting. Most of these points
have maximum curvature on the segment. When we try to find the point at
which curvature is as maximum as possible, we will reduce the scope of
the unhandled segment of spline. Therefore, we can generate fewer points
while rendering an unhandled segment of spline. However, the drawback of
this method is the high cost of finding points and the high cost of the
bisection method.

Next, to solve the drawback mentioned above, we apply a flexible update
't' with early stopping, which aims to stop the call of the
'_de_casteljau' function. However, the improvement is little.

Based on the previous experiments, we can see that the root of the
problem is that there is a must to call '_de_casteljau' function twice
to generate a point even in the original setting, which has the most
fewer numbers of call '_de_casteljau'.

The ideal method is that the rendering process only needs to call
'_de_casteljau' once to generate a point. To achieve this goal, we
reduce the initial 't' from 0.5 to 0.25 because the average number of
'_de_casteljau' calls per point is 1.99. This means that the first call
of the function '_de_casteljau' is redundant. In the experiment, we can
see that the average points per character of the original method that
shift 2 is almost the same with the original method.

Experiment :
Use font-edit as the evaluation testbed to do experiments as follows:

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
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

After considering the computation cost of the '_de_casteljau' function
and the increase in average number of _de_casteljau calls per point, we
abandon the method flexible update 't' with the bisection method. Then,
Choose to apply the original method with shift 2 finally.
Take t = 0.25 for example (3/4) * A + (1/4) * B
Original: A + (B - A) >> 4
We can directly use (B - A) to get the coefficients 1/4 and 3/4.
Bisection: we need to multiply (B - A) with 0.25.
Also, the bisection has many possible values of t rather than the
original method. When using the original method, the '_de_casteljau'
function can all use shift operations because the t is always the power
of 2.

Author: weihsinyeh

src/spline.c

@@ -11,38 +11,36 @@ typedef struct _twin_spline {
 } twin_spline_t;
 
 /*
- * Linearly interpolate between points 'a' and 'b' with a 't' factor.
- * The 't' 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,
                   const twin_spoint_t *b,
-                  twin_dfixed_t t,
+                  int shift,
                   twin_spoint_t *result)
 {
-    result->x = a->x + twin_dfixed_to_sfixed(twin_dfixed_mul(
-                           twin_sfixed_to_dfixed(b->x - a->x), t));
-    result->y = a->y + twin_dfixed_to_sfixed(twin_dfixed_mul(
-                           twin_sfixed_to_dfixed(b->y - a->y), t));
+    result->x = a->x + ((b->x - a->x) >> shift);
+    result->y = a->y + ((b->y - a->y) >> shift);
 }
 
 /*
- * Perform the de Casteljau algorithm to split a spline at a given 't'
+ * 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,
-                          twin_dfixed_t t,
+                          int shift,
                           twin_spline_t *s1,
                           twin_spline_t *s2)
 {
     twin_spoint_t ab, bc, cd, abbc, bccd, 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);
+    _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);
 
     s1->a = spline->a;
     s1->b = ab;
@@ -64,14 +62,14 @@ static void _de_casteljau(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;
 }
 
 /*
@@ -94,49 +92,20 @@ static void _twin_spline_decompose(twin_path_t *path,
 {
     /* Draw starting point */
     _twin_path_sdraw(path, spline->a.x, spline->a.y);
-
+    int shift = 2;
     while (!is_flat(spline, tolerance_squared)) {
-        twin_dfixed_t hi = TWIN_SFIXED_ONE * TWIN_SFIXED_ONE, lo = 0,
-                      t_optimal = 0, t, max_distance = 0, left_dist, right_dist;
         twin_spline_t left, right;
 
         while (true) {
-            t = (hi + lo) >> 1;
-            if (t_optimal == t)
+            _de_casteljau(spline, shift, &left, &right);
+            if (is_flat(&left, tolerance_squared)) {
+                if (shift > 1)
+                    shift--;
                 break;
-
-            _de_casteljau(spline, t, &left, &right);
-
-            left_dist = _twin_spline_distance_squared(&left);
-            right_dist = _twin_spline_distance_squared(&right);
-
-            /* The left segment is close enough to fit the original spline. */
-            if (left_dist <= tolerance_squared) {
-                max_distance = left_dist;
-                t_optimal = t;
-                if (right_dist <= tolerance_squared)
-                    break;
-                /*
-                 * Try to find a better point such that the line segment
-                 * connecting it to the previous point can fit the spline, while
-                 * maximizing the distance from the spline as much as possible.
-                 */
-                lo = t;
-            } else {
-                /*
-                 * The line segment connecting the previous point to the
-                 * currently tested point with value 't' cannot fit the original
-                 * spline. The point with value 't_optimal' is sufficient to fit
-                 * the original spline.
-                 */
-                if (t_optimal)
-                    break;
-                hi = t;
-            }
+            } else
+                shift++;
         }
 
-        _de_casteljau(spline, t_optimal, &left, &right);
-
         /* Draw the left segment */
         _twin_path_sdraw(path, left.d.x, left.d.y);