From 55d241a3891419cf0e23adc6cad9bd84d41098c6 Mon Sep 17 00:00:00 2001 From: Wei-Hsin Yeh Date: Thu, 8 Aug 2024 02:08:03 +0800 Subject: [PATCH] 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. Close #2 --- src/spline.c | 76 ++++++++++++++++++++++++++++++++++++---------------- 1 file changed, 53 insertions(+), 23 deletions(-) diff --git a/src/spline.c b/src/spline.c index 66203ed..cdb443f 100644 --- a/src/spline.c +++ b/src/spline.c @@ -11,25 +11,25 @@ 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) { @@ -37,12 +37,12 @@ static void _de_casteljau(twin_spline_t *spline, 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; @@ -56,10 +56,10 @@ static void _de_casteljau(twin_spline_t *spline, } /* - * Return an upper bound on the error (squared) that could result from + * Return an upper bound on the distance (squared) that could result from * approximating a spline with a line segment connecting the two endpoints. */ -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; @@ -72,12 +72,12 @@ static twin_dfixed_t _twin_spline_error_squared(twin_spline_t *spline) } /* - * 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; } /* @@ -93,14 +93,44 @@ 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 t_optimal = 0, t; + twin_dfixed_t max_distance = 0, distance; 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 (t_optimal == t) + break; + + _de_casteljau(spline, t, &left, &right); + + distance = _twin_spline_distance_squared(&left); + if (distance < max_distance) + break; + + /* The left segment is close enough to fit the original spline. */ + if (distance <= tolerance_squared) { + max_distance = distance; + t_optimal = t; + /* 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; + } + } + + _de_casteljau(spline, t_optimal, &left, &right); /* Draw the left segment */ _twin_path_sdraw(path, left.d.x, left.d.y);