1// Copyright 2022 the Kurbo Authors
2// SPDX-License-Identifier: Apache-2.0 OR MIT
3
4//! An implementation of cubic Bézier curve fitting based on a quartic
5//! solver making signed area and moment match the source curve.
6
7use core::ops::Range;
8
9use alloc::vec::Vec;
10
11use arrayvec::ArrayVec;
12
13use crate::{
14 common::{
15 factor_quartic_inner, solve_cubic, solve_itp_fallible, solve_quadratic,
16 GAUSS_LEGENDRE_COEFFS_16,
17 },
18 Affine, BezPath, CubicBez, Line, ParamCurve, ParamCurveArclen, ParamCurveNearest, Point, Vec2,
19};
20
21#[cfg(not(feature = "std"))]
22use crate::common::FloatFuncs;
23
24/// The source curve for curve fitting.
25///
26/// This trait is a general representation of curves to be used as input to a curve
27/// fitting process. It can represent source curves with cusps and corners, though
28/// if the corners are known in advance, it may be better to run curve fitting on
29/// subcurves bounded by the corners.
30///
31/// The trait primarily works by sampling the source curve and computing the position
32/// and derivative at each sample. Those derivatives are then used for multiple
33/// sub-tasks, including ensuring G1 continuity at subdivision points, computing the
34/// area and moment of the curve for curve fitting, and casting rays for evaluation
35/// of a distance metric to test accuracy.
36///
37/// A major motivation is computation of offset curves, which often have cusps, but
38/// the presence and location of those cusps is not generally known. It is also
39/// intended for conversion between curve types (for example, piecewise Euler spiral
40/// or NURBS), and distortion effects such as perspective transform.
41///
42/// Note general similarities to [`ParamCurve`] but also important differences.
43/// Instead of separate [`eval`] and evaluation of derivative, have a single
44/// [`sample_pt_deriv`] method which can be more efficient and also handles cusps more
45/// robustly. Also there is no method for subsegment, as that is not needed and
46/// would be annoying to implement.
47///
48/// [`ParamCurve`]: crate::ParamCurve
49/// [`eval`]: crate::ParamCurve::eval
50/// [`sample_pt_deriv`]: ParamCurveFit::sample_pt_deriv
51pub trait ParamCurveFit {
52 /// Evaluate the curve and its tangent at parameter `t`.
53 ///
54 /// For a regular curve (one not containing a cusp or corner), the
55 /// derivative is a good choice for the tangent vector and the `sign`
56 /// parameter can be ignored. Otherwise, the `sign` parameter selects which
57 /// side of the discontinuity the tangent will be sampled from.
58 ///
59 /// Generally `t` is in the range [0..1].
60 fn sample_pt_tangent(&self, t: f64, sign: f64) -> CurveFitSample;
61
62 /// Evaluate the point and derivative at parameter `t`.
63 ///
64 /// In curves with cusps, the derivative can go to zero.
65 fn sample_pt_deriv(&self, t: f64) -> (Point, Vec2);
66
67 /// Compute moment integrals.
68 ///
69 /// This method computes the integrals of y dx, x y dx, and y^2 dx over the
70 /// length of this curve. From these integrals it is fairly straightforward
71 /// to derive the moments needed for curve fitting.
72 ///
73 /// A default implementation is proved which does quadrature integration
74 /// with Green's theorem, in terms of samples evaluated with
75 /// [`sample_pt_deriv`].
76 ///
77 /// [`sample_pt_deriv`]: ParamCurveFit::sample_pt_deriv
78 fn moment_integrals(&self, range: Range<f64>) -> (f64, f64, f64) {
79 let t0 = 0.5 * (range.start + range.end);
80 let dt = 0.5 * (range.end - range.start);
81 let (a, x, y) = GAUSS_LEGENDRE_COEFFS_16
82 .iter()
83 .map(|(wi, xi)| {
84 let t = t0 + xi * dt;
85 let (p, d) = self.sample_pt_deriv(t);
86 let a = wi * d.x * p.y;
87 let x = p.x * a;
88 let y = p.y * a;
89 (a, x, y)
90 })
91 .fold((0.0, 0.0, 0.0), |(a0, x0, y0), (a1, x1, y1)| {
92 (a0 + a1, x0 + x1, y0 + y1)
93 });
94 (a * dt, x * dt, y * dt)
95 }
96
97 /// Find a cusp or corner within the given range.
98 ///
99 /// If the range contains a corner or cusp, return it. If there is more
100 /// than one such discontinuity, any can be reported, as the function will
101 /// be called repeatedly after subdivision of the range.
102 ///
103 /// Do not report cusps at the endpoints of the range, as this may cause
104 /// potentially infinite subdivision. In particular, when a cusp is reported
105 /// and this method is called on a subdivided range bounded by the reported
106 /// cusp, then the subsequent call should not report a cusp there.
107 ///
108 /// The definition of what exactly constitutes a cusp is somewhat loose.
109 /// If a cusp is missed, then the curve fitting algorithm will attempt to
110 /// fit the curve with a smooth curve, which is generally not a disaster
111 /// will usually result in more subdivision. Conversely, it might be useful
112 /// to report near-cusps, specifically points of curvature maxima where the
113 /// curvature is large but still mathematically finite.
114 fn break_cusp(&self, range: Range<f64>) -> Option<f64>;
115}
116
117/// A sample point of a curve for fitting.
118pub struct CurveFitSample {
119 /// A point on the curve at the sample location.
120 pub p: Point,
121 /// A vector tangent to the curve at the sample location.
122 pub tangent: Vec2,
123}
124
125impl CurveFitSample {
126 /// Intersect a ray orthogonal to the tangent with the given cubic.
127 ///
128 /// Returns a vector of `t` values on the cubic.
129 fn intersect(&self, c: CubicBez) -> ArrayVec<f64, 3> {
130 let p1: Vec2 = 3.0 * (c.p1 - c.p0);
131 let p2: Vec2 = 3.0 * c.p2.to_vec2() - 6.0 * c.p1.to_vec2() + 3.0 * c.p0.to_vec2();
132 let p3: Vec2 = (c.p3 - c.p0) - 3.0 * (c.p2 - c.p1);
133 let c0: f64 = (c.p0 - self.p).dot(self.tangent);
134 let c1: f64 = p1.dot(self.tangent);
135 let c2: f64 = p2.dot(self.tangent);
136 let c3: f64 = p3.dot(self.tangent);
137 solve_cubicimpl Iterator(c0, c1, c2, c3)
138 .into_iter()
139 .filter(|t: &f64| (0.0..=1.0).contains(item:t))
140 .collect()
141 }
142}
143
144/// Generate a Bézier path that fits the source curve.
145///
146/// This function is still experimental and the signature might change; it's possible
147/// it might become a method on the [`ParamCurveFit`] trait.
148///
149/// This function recursively subdivides the curve in half by the parameter when the
150/// accuracy is not met. That gives a reasonably optimized result but not necessarily
151/// the minimum number of segments.
152///
153/// In general, the resulting Bézier path should have a Fréchet distance less than
154/// the provided `accuracy` parameter. However, this is not a rigorous guarantee, as
155/// the error metric is computed approximately.
156///
157/// This function is intended for use when the source curve is piecewise continuous,
158/// with the discontinuities reported by the cusp method. In applications (such as
159/// stroke expansion) where this property may not hold, it is up to the client to
160/// detect and handle such cases. Even so, best effort is made to avoid infinite
161/// subdivision.
162///
163/// When a higher degree of optimization is desired (at considerably more runtime cost),
164/// consider [`fit_to_bezpath_opt`] instead.
165pub fn fit_to_bezpath(source: &impl ParamCurveFit, accuracy: f64) -> BezPath {
166 let mut path: BezPath = BezPath::new();
167 fit_to_bezpath_rec(source, range:0.0..1.0, accuracy, &mut path);
168 path
169}
170
171// Discussion question: possibly should take endpoint samples, to avoid
172// duplication of that work.
173fn fit_to_bezpath_rec(
174 source: &impl ParamCurveFit,
175 range: Range<f64>,
176 accuracy: f64,
177 path: &mut BezPath,
178) {
179 let start = range.start;
180 let end = range.end;
181 let start_p = source.sample_pt_tangent(range.start, 1.0).p;
182 let end_p = source.sample_pt_tangent(range.end, -1.0).p;
183 if start_p.distance_squared(end_p) <= accuracy * accuracy {
184 if let Some((c, _)) = try_fit_line(source, accuracy, range, start_p, end_p) {
185 if path.is_empty() {
186 path.move_to(c.p0);
187 }
188 path.curve_to(c.p1, c.p2, c.p3);
189 return;
190 }
191 }
192 let t = if let Some(t) = source.break_cusp(start..end) {
193 t
194 } else if let Some((c, _)) = fit_to_cubic(source, start..end, accuracy) {
195 if path.is_empty() {
196 path.move_to(c.p0);
197 }
198 path.curve_to(c.p1, c.p2, c.p3);
199 return;
200 } else {
201 // A smarter approach is possible than midpoint subdivision, but would be
202 // a significant increase in complexity.
203 0.5 * (start + end)
204 };
205 if t == start || t == end {
206 // infinite recursion, just draw a line
207 let p1 = start_p.lerp(end_p, 1.0 / 3.0);
208 let p2 = end_p.lerp(start_p, 1.0 / 3.0);
209 if path.is_empty() {
210 path.move_to(start_p);
211 }
212 path.curve_to(p1, p2, end_p);
213 return;
214 }
215 fit_to_bezpath_rec(source, start..t, accuracy, path);
216 fit_to_bezpath_rec(source, t..end, accuracy, path);
217}
218
219const N_SAMPLE: usize = 20;
220
221/// An acceleration structure for estimating curve distance
222struct CurveDist {
223 samples: ArrayVec<CurveFitSample, N_SAMPLE>,
224 arcparams: ArrayVec<f64, N_SAMPLE>,
225 range: Range<f64>,
226 /// A "spicy" curve is one with potentially extreme curvature variation,
227 /// so use arc length measurement for better accuracy.
228 spicy: bool,
229}
230
231impl CurveDist {
232 fn from_curve(source: &impl ParamCurveFit, range: Range<f64>) -> Self {
233 let step = (range.end - range.start) * (1.0 / (N_SAMPLE + 1) as f64);
234 let mut last_tan = None;
235 let mut spicy = false;
236 const SPICY_THRESH: f64 = 0.2;
237 let mut samples = ArrayVec::new();
238 for i in 0..N_SAMPLE + 2 {
239 let sample = source.sample_pt_tangent(range.start + i as f64 * step, 1.0);
240 if let Some(last_tan) = last_tan {
241 let cross = sample.tangent.cross(last_tan);
242 let dot = sample.tangent.dot(last_tan);
243 if cross.abs() > SPICY_THRESH * dot.abs() {
244 spicy = true;
245 }
246 }
247 last_tan = Some(sample.tangent);
248 if i > 0 && i < N_SAMPLE + 1 {
249 samples.push(sample);
250 }
251 }
252 CurveDist {
253 samples,
254 arcparams: Default::default(),
255 range,
256 spicy,
257 }
258 }
259
260 fn compute_arc_params(&mut self, source: &impl ParamCurveFit) {
261 const N_SUBSAMPLE: usize = 10;
262 let (start, end) = (self.range.start, self.range.end);
263 let dt = (end - start) * (1.0 / ((N_SAMPLE + 1) * N_SUBSAMPLE) as f64);
264 let mut arclen = 0.0;
265 for i in 0..N_SAMPLE + 1 {
266 for j in 0..N_SUBSAMPLE {
267 let t = start + dt * ((i * N_SUBSAMPLE + j) as f64 + 0.5);
268 let (_, deriv) = source.sample_pt_deriv(t);
269 arclen += deriv.hypot();
270 }
271 if i < N_SAMPLE {
272 self.arcparams.push(arclen);
273 }
274 }
275 let arclen_inv = arclen.recip();
276 for x in &mut self.arcparams {
277 *x *= arclen_inv;
278 }
279 }
280
281 /// Evaluate distance based on arc length parametrization
282 fn eval_arc(&self, c: CubicBez, acc2: f64) -> Option<f64> {
283 // TODO: this could perhaps be tuned.
284 const EPS: f64 = 1e-9;
285 let c_arclen = c.arclen(EPS);
286 let mut max_err2 = 0.0;
287 for (sample, s) in self.samples.iter().zip(&self.arcparams) {
288 let t = c.inv_arclen(c_arclen * s, EPS);
289 let err = sample.p.distance_squared(c.eval(t));
290 max_err2 = err.max(max_err2);
291 if max_err2 > acc2 {
292 return None;
293 }
294 }
295 Some(max_err2)
296 }
297
298 /// Evaluate distance to a cubic approximation.
299 ///
300 /// If distance exceeds stated accuracy, can return `None`. Note that
301 /// `acc2` is the square of the target.
302 ///
303 /// Returns the square of the error, which is intended to be a good
304 /// approximation of the Fréchet distance.
305 fn eval_ray(&self, c: CubicBez, acc2: f64) -> Option<f64> {
306 let mut max_err2 = 0.0;
307 for sample in &self.samples {
308 let mut best = acc2 + 1.0;
309 for t in sample.intersect(c) {
310 let err = sample.p.distance_squared(c.eval(t));
311 best = best.min(err);
312 }
313 max_err2 = best.max(max_err2);
314 if max_err2 > acc2 {
315 return None;
316 }
317 }
318 Some(max_err2)
319 }
320
321 fn eval_dist(&mut self, source: &impl ParamCurveFit, c: CubicBez, acc2: f64) -> Option<f64> {
322 // Always compute cheaper distance, hoping for early-out.
323 let ray_dist = self.eval_ray(c, acc2)?;
324 if !self.spicy {
325 return Some(ray_dist);
326 }
327 if self.arcparams.is_empty() {
328 self.compute_arc_params(source);
329 }
330 self.eval_arc(c, acc2)
331 }
332}
333
334/// As described in [Simplifying Bézier paths], strictly optimizing for
335/// Fréchet distance can create bumps. The problem is curves with long
336/// control arms (distance from the control point to the corresponding
337/// endpoint). We mitigate that by applying a penalty as a multiplier to
338/// the measured error (approximate Fréchet distance). This is ReLU-like,
339/// with a value of 1.0 below the elbow, and a given slope above it. The
340/// values here have been determined empirically to give good results.
341///
342/// [Simplifying Bézier paths]:
343/// https://raphlinus.github.io/curves/2023/04/18/bezpath-simplify.html
344const D_PENALTY_ELBOW: f64 = 0.65;
345const D_PENALTY_SLOPE: f64 = 2.0;
346
347/// Try fitting a line.
348///
349/// This is especially useful for very short chords, in which the standard
350/// cubic fit is not numerically stable. The tangents are not considered, so
351/// it's useful in cusp and near-cusp situations where the tangents are not
352/// reliable, as well.
353///
354/// Returns the line raised to a cubic and the error, if within tolerance.
355fn try_fit_line(
356 source: &impl ParamCurveFit,
357 accuracy: f64,
358 range: Range<f64>,
359 start: Point,
360 end: Point,
361) -> Option<(CubicBez, f64)> {
362 let acc2: f64 = accuracy * accuracy;
363 let chord_l: Line = Line::new(p0:start, p1:end);
364 const SHORT_N: usize = 7;
365 let mut max_err2: f64 = 0.0;
366 let dt: f64 = (range.end - range.start) / (SHORT_N + 1) as f64;
367 for i: usize in 0..SHORT_N {
368 let t: f64 = range.start + (i + 1) as f64 * dt;
369 let p: Point = source.sample_pt_deriv(t).0;
370 let err2: f64 = chord_l.nearest(p, accuracy).distance_sq;
371 if err2 > acc2 {
372 // Not in tolerance; likely subdivision will help.
373 return None;
374 }
375 max_err2 = err2.max(max_err2);
376 }
377 let p1: Point = start.lerp(other:end, t:1.0 / 3.0);
378 let p2: Point = end.lerp(other:start, t:1.0 / 3.0);
379 let c: CubicBez = CubicBez::new(p0:start, p1, p2, p3:end);
380 Some((c, max_err2))
381}
382
383/// Fit a single cubic to a range of the source curve.
384///
385/// Returns the cubic segment and the square of the error.
386/// Discussion question: should this be a method on the trait instead?
387pub fn fit_to_cubic(
388 source: &impl ParamCurveFit,
389 range: Range<f64>,
390 accuracy: f64,
391) -> Option<(CubicBez, f64)> {
392 let start = source.sample_pt_tangent(range.start, 1.0);
393 let end = source.sample_pt_tangent(range.end, -1.0);
394 let d = end.p - start.p;
395 let chord2 = d.hypot2();
396 let acc2 = accuracy * accuracy;
397 if chord2 <= acc2 {
398 // Special case very short chords; try to fit a line.
399 return try_fit_line(source, accuracy, range, start.p, end.p);
400 }
401 let th = d.atan2();
402 fn mod_2pi(th: f64) -> f64 {
403 let th_scaled = th * core::f64::consts::FRAC_1_PI * 0.5;
404 core::f64::consts::PI * 2.0 * (th_scaled - th_scaled.round())
405 }
406 let th0 = mod_2pi(start.tangent.atan2() - th);
407 let th1 = mod_2pi(th - end.tangent.atan2());
408
409 let (mut area, mut x, mut y) = source.moment_integrals(range.clone());
410 let (x0, y0) = (start.p.x, start.p.y);
411 let (dx, dy) = (d.x, d.y);
412 // Subtract off area of chord
413 area -= dx * (y0 + 0.5 * dy);
414 // `area` is signed area of closed curve segment.
415 // This quantity is invariant to translation and rotation.
416
417 // Subtract off moment of chord
418 let dy_3 = dy * (1. / 3.);
419 x -= dx * (x0 * y0 + 0.5 * (x0 * dy + y0 * dx) + dy_3 * dx);
420 y -= dx * (y0 * y0 + y0 * dy + dy_3 * dy);
421 // Translate start point to origin; convert raw integrals to moments.
422 x -= x0 * area;
423 y = 0.5 * y - y0 * area;
424 // Rotate into place (this also scales up by chordlength for efficiency).
425 let moment = d.x * x + d.y * y;
426 // `moment` is the chordlength times the x moment of the curve translated
427 // so its start point is on the origin, and rotated so its end point is on the
428 // x axis.
429
430 let chord2_inv = chord2.recip();
431 let unit_area = area * chord2_inv;
432 let mx = moment * chord2_inv.powi(2);
433 // `unit_area` is signed area scaled to unit chord; `mx` is scaled x moment
434
435 let chord = chord2.sqrt();
436 let aff = Affine::translate(start.p.to_vec2()) * Affine::rotate(th) * Affine::scale(chord);
437 let mut curve_dist = CurveDist::from_curve(source, range);
438 let mut best_c = None;
439 let mut best_err2 = None;
440 for (cand, d0, d1) in cubic_fit(th0, th1, unit_area, mx) {
441 let c = aff * cand;
442 if let Some(err2) = curve_dist.eval_dist(source, c, acc2) {
443 fn scale_f(d: f64) -> f64 {
444 1.0 + (d - D_PENALTY_ELBOW).max(0.0) * D_PENALTY_SLOPE
445 }
446 let scale = scale_f(d0).max(scale_f(d1)).powi(2);
447 let err2 = err2 * scale;
448 if err2 < acc2 && best_err2.map(|best| err2 < best).unwrap_or(true) {
449 best_c = Some(c);
450 best_err2 = Some(err2);
451 }
452 }
453 }
454 match (best_c, best_err2) {
455 (Some(c), Some(err2)) => Some((c, err2)),
456 _ => None,
457 }
458}
459
460/// Returns curves matching area and moment, given unit chord.
461fn cubic_fit(th0: f64, th1: f64, area: f64, mx: f64) -> ArrayVec<(CubicBez, f64, f64), 4> {
462 // Note: maybe we want to take unit vectors instead of angle? Shouldn't
463 // matter much either way though.
464 let (s0, c0) = th0.sin_cos();
465 let (s1, c1) = th1.sin_cos();
466 let a4 = -9.
467 * c0
468 * (((2. * s1 * c1 * c0 + s0 * (2. * c1 * c1 - 1.)) * c0 - 2. * s1 * c1) * c0
469 - c1 * c1 * s0);
470 let a3 = 12.
471 * ((((c1 * (30. * area * c1 - s1) - 15. * area) * c0 + 2. * s0
472 - c1 * s0 * (c1 + 30. * area * s1))
473 * c0
474 + c1 * (s1 - 15. * area * c1))
475 * c0
476 - s0 * c1 * c1);
477 let a2 = 12.
478 * ((((70. * mx + 15. * area) * s1 * s1 + c1 * (9. * s1 - 70. * c1 * mx - 5. * c1 * area))
479 * c0
480 - 5. * s0 * s1 * (3. * s1 - 4. * c1 * (7. * mx + area)))
481 * c0
482 - c1 * (9. * s1 - 70. * c1 * mx - 5. * c1 * area));
483 let a1 = 16.
484 * (((12. * s0 - 5. * c0 * (42. * mx - 17. * area)) * s1
485 - 70. * c1 * (3. * mx - area) * s0
486 - 75. * c0 * c1 * area * area)
487 * s1
488 - 75. * c1 * c1 * area * area * s0);
489 let a0 = 80. * s1 * (42. * s1 * mx - 25. * area * (s1 - c1 * area));
490 // TODO: "roots" is not a good name for this variable, as it also contains
491 // the real part of complex conjugate pairs.
492 let mut roots = ArrayVec::<f64, 4>::new();
493 const EPS: f64 = 1e-12;
494 if a4.abs() > EPS {
495 let a = a3 / a4;
496 let b = a2 / a4;
497 let c = a1 / a4;
498 let d = a0 / a4;
499 if let Some(quads) = factor_quartic_inner(a, b, c, d, false) {
500 for (qc1, qc0) in quads {
501 let qroots = solve_quadratic(qc0, qc1, 1.0);
502 if qroots.is_empty() {
503 // Real part of pair of complex roots
504 roots.push(-0.5 * qc1);
505 } else {
506 roots.extend(qroots);
507 }
508 }
509 }
510 } else if a3.abs() > EPS {
511 roots.extend(solve_cubic(a0, a1, a2, a3));
512 } else if a2.abs() > EPS || a1.abs() > EPS || a0.abs() > EPS {
513 roots.extend(solve_quadratic(a0, a1, a2));
514 } else {
515 return [(
516 CubicBez::new((0.0, 0.0), (1. / 3., 0.0), (2. / 3., 0.0), (1., 0.0)),
517 1f64 / 3.,
518 1f64 / 3.,
519 )]
520 .into_iter()
521 .collect();
522 }
523
524 let s01 = s0 * c1 + s1 * c0;
525 roots
526 .iter()
527 .filter_map(|&d0| {
528 let (d0, d1) = if d0 > 0.0 {
529 let d1 = (d0 * s0 - area * (10. / 3.)) / (0.5 * d0 * s01 - s1);
530 if d1 > 0.0 {
531 (d0, d1)
532 } else {
533 (s1 / s01, 0.0)
534 }
535 } else {
536 (0.0, s0 / s01)
537 };
538 // We could implement a maximum d value here.
539 if d0 >= 0.0 && d1 >= 0.0 {
540 Some((
541 CubicBez::new(
542 (0.0, 0.0),
543 (d0 * c0, d0 * s0),
544 (1.0 - d1 * c1, d1 * s1),
545 (1.0, 0.0),
546 ),
547 d0,
548 d1,
549 ))
550 } else {
551 None
552 }
553 })
554 .collect()
555}
556
557/// Generate a highly optimized Bézier path that fits the source curve.
558///
559/// This function is still experimental and the signature might change; it's possible
560/// it might become a method on the [`ParamCurveFit`] trait.
561///
562/// This function is considerably slower than [`fit_to_bezpath`], as it computes
563/// optimal subdivision points. Its result is expected to be very close to the optimum
564/// possible Bézier path for the source curve, in that it has a minimal number of curve
565/// segments, and a minimal error over all paths with that number of segments.
566///
567/// See [`fit_to_bezpath`] for an explanation of the `accuracy` parameter.
568pub fn fit_to_bezpath_opt(source: &impl ParamCurveFit, accuracy: f64) -> BezPath {
569 let mut cusps: Vec = Vec::new();
570 let mut path: BezPath = BezPath::new();
571 let mut t0: f64 = 0.0;
572 loop {
573 let t1: f64 = cusps.last().copied().unwrap_or(default:1.0);
574 match fit_to_bezpath_opt_inner(source, accuracy, range:t0..t1, &mut path) {
575 Some(t: f64) => cusps.push(t),
576 None => match cusps.pop() {
577 Some(t: f64) => t0 = t,
578 None => break,
579 },
580 }
581 }
582 path
583}
584
585/// Fit a range without cusps.
586///
587/// On Ok return, range has been added to the path. On Err return, report a cusp (and don't
588/// mutate path).
589fn fit_to_bezpath_opt_inner(
590 source: &impl ParamCurveFit,
591 accuracy: f64,
592 range: Range<f64>,
593 path: &mut BezPath,
594) -> Option<f64> {
595 if let Some(t) = source.break_cusp(range.clone()) {
596 return Some(t);
597 }
598 let err;
599 if let Some((c, err2)) = fit_to_cubic(source, range.clone(), accuracy) {
600 err = err2.sqrt();
601 if err < accuracy {
602 if range.start == 0.0 {
603 path.move_to(c.p0);
604 }
605 path.curve_to(c.p1, c.p2, c.p3);
606 return None;
607 }
608 } else {
609 err = 2.0 * accuracy;
610 }
611 let (mut t0, t1) = (range.start, range.end);
612 let mut n = 0;
613 let last_err;
614 loop {
615 n += 1;
616 match fit_opt_segment(source, accuracy, t0..t1) {
617 FitResult::ParamVal(t) => t0 = t,
618 FitResult::SegmentError(err) => {
619 last_err = err;
620 break;
621 }
622 FitResult::CuspFound(t) => return Some(t),
623 }
624 }
625 t0 = range.start;
626 const EPS: f64 = 1e-9;
627 let f = |x| fit_opt_err_delta(source, x, accuracy, t0..t1, n);
628 let k1 = 0.2 / accuracy;
629 let ya = -err;
630 let yb = accuracy - last_err;
631 let (_, x) = match solve_itp_fallible(f, 0.0, accuracy, EPS, 1, k1, ya, yb) {
632 Ok(x) => x,
633 Err(t) => return Some(t),
634 };
635 //println!("got fit with n={}, err={}", n, x);
636 let path_len = path.elements().len();
637 for i in 0..n {
638 let t1 = if i < n - 1 {
639 match fit_opt_segment(source, x, t0..range.end) {
640 FitResult::ParamVal(t1) => t1,
641 FitResult::SegmentError(_) => range.end,
642 FitResult::CuspFound(t) => {
643 path.truncate(path_len);
644 return Some(t);
645 }
646 }
647 } else {
648 range.end
649 };
650 let (c, _) = fit_to_cubic(source, t0..t1, accuracy).unwrap();
651 if i == 0 && range.start == 0.0 {
652 path.move_to(c.p0);
653 }
654 path.curve_to(c.p1, c.p2, c.p3);
655 t0 = t1;
656 if t0 == range.end {
657 // This is unlikely but could happen when not monotonic.
658 break;
659 }
660 }
661 None
662}
663
664fn measure_one_seg(source: &impl ParamCurveFit, range: Range<f64>, limit: f64) -> Option<f64> {
665 fit_to_cubic(source, range, accuracy:limit).map(|(_, err2: f64)| err2.sqrt())
666}
667
668enum FitResult {
669 /// The parameter (`t`) value that meets the desired accuracy.
670 ParamVal(f64),
671 /// Error of the measured segment.
672 SegmentError(f64),
673 /// The parameter value where a cusp was found.
674 CuspFound(f64),
675}
676
677fn fit_opt_segment(source: &impl ParamCurveFit, accuracy: f64, range: Range<f64>) -> FitResult {
678 if let Some(t: f64) = source.break_cusp(range.clone()) {
679 return FitResult::CuspFound(t);
680 }
681 let missing_err: f64 = accuracy * 2.0;
682 let err: f64 = measure_one_seg(source, range.clone(), accuracy).unwrap_or(default:missing_err);
683 if err <= accuracy {
684 return FitResult::SegmentError(err);
685 }
686 let (t0: f64, t1: f64) = (range.start, range.end);
687 let f: impl Fn(f64) -> Result = |x: f64| {
688 if let Some(t: f64) = source.break_cusp(range.clone()) {
689 return Err(t);
690 }
691 let err: f64 = measure_one_seg(source, t0..x, accuracy).unwrap_or(default:missing_err);
692 Ok(err - accuracy)
693 };
694 const EPS: f64 = 1e-9;
695 let k1: f64 = 2.0 / (t1 - t0);
696 match solve_itp_fallible(f, a:t0, b:t1, EPS, n0:1, k1, -accuracy, yb:err - accuracy) {
697 Ok((t1: f64, _)) => FitResult::ParamVal(t1),
698 Err(t: f64) => FitResult::CuspFound(t),
699 }
700}
701
702// Ok result is delta error (accuracy - error of last seg).
703// Err result is a cusp.
704fn fit_opt_err_delta(
705 source: &impl ParamCurveFit,
706 accuracy: f64,
707 limit: f64,
708 range: Range<f64>,
709 n: usize,
710) -> Result<f64, f64> {
711 let (mut t0: f64, t1: f64) = (range.start, range.end);
712 for _ in 0..n - 1 {
713 t0 = match fit_opt_segment(source, accuracy, range:t0..t1) {
714 FitResult::ParamVal(t0: f64) => t0,
715 // In this case, n - 1 will work, which of course means the error is highly
716 // non-monotonic. We should probably harvest that solution.
717 FitResult::SegmentError(_) => return Ok(accuracy),
718 FitResult::CuspFound(t: f64) => return Err(t),
719 }
720 }
721 let err: f64 = measure_one_seg(source, t0..t1, limit).unwrap_or(default:accuracy * 2.0);
722 Ok(accuracy - err)
723}
724