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, ParamCurve, ParamCurveArclen, 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 evaulated 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/// When a higher degree of optimization is desired (at considerably more runtime cost),
158/// consider [`fit_to_bezpath_opt`] instead.
159pub fn fit_to_bezpath(source: &impl ParamCurveFit, accuracy: f64) -> BezPath {
160 let mut path: BezPath = BezPath::new();
161 fit_to_bezpath_rec(source, range:0.0..1.0, accuracy, &mut path);
162 path
163}
164
165// Discussion question: possibly should take endpoint samples, to avoid
166// duplication of that work.
167fn fit_to_bezpath_rec(
168 source: &impl ParamCurveFit,
169 range: Range<f64>,
170 accuracy: f64,
171 path: &mut BezPath,
172) {
173 let start: f64 = range.start;
174 let end: f64 = range.end;
175 if let Some(t: f64) = source.break_cusp(range:start..end) {
176 fit_to_bezpath_rec(source, range:start..t, accuracy, path);
177 fit_to_bezpath_rec(source, range:t..end, accuracy, path);
178 } else if let Some((c: CubicBez, _)) = fit_to_cubic(source, range:start..end, accuracy) {
179 if path.is_empty() {
180 path.move_to(c.p0);
181 }
182 path.curve_to(c.p1, c.p2, c.p3);
183 } else {
184 // A smarter approach is possible than midpoint subdivision, but would be
185 // a significant increase in complexity.
186 let t: f64 = 0.5 * (start + end);
187 fit_to_bezpath_rec(source, range:start..t, accuracy, path);
188 fit_to_bezpath_rec(source, range:t..end, accuracy, path);
189 }
190}
191
192const N_SAMPLE: usize = 20;
193
194/// An acceleration structure for estimating curve distance
195struct CurveDist {
196 samples: ArrayVec<CurveFitSample, N_SAMPLE>,
197 arcparams: ArrayVec<f64, N_SAMPLE>,
198 range: Range<f64>,
199 /// A "spicy" curve is one with potentially extreme curvature variation,
200 /// so use arc length measurement for better accuracy.
201 spicy: bool,
202}
203
204impl CurveDist {
205 fn from_curve(source: &impl ParamCurveFit, range: Range<f64>) -> Self {
206 let step = (range.end - range.start) * (1.0 / (N_SAMPLE + 1) as f64);
207 let mut last_tan = None;
208 let mut spicy = false;
209 const SPICY_THRESH: f64 = 0.2;
210 let mut samples = ArrayVec::new();
211 for i in 0..N_SAMPLE + 2 {
212 let sample = source.sample_pt_tangent(range.start + i as f64 * step, 1.0);
213 if let Some(last_tan) = last_tan {
214 let cross = sample.tangent.cross(last_tan);
215 let dot = sample.tangent.dot(last_tan);
216 if cross.abs() > SPICY_THRESH * dot.abs() {
217 spicy = true;
218 }
219 }
220 last_tan = Some(sample.tangent);
221 if i > 0 && i < N_SAMPLE + 1 {
222 samples.push(sample);
223 }
224 }
225 CurveDist {
226 samples,
227 arcparams: Default::default(),
228 range,
229 spicy,
230 }
231 }
232
233 fn compute_arc_params(&mut self, source: &impl ParamCurveFit) {
234 const N_SUBSAMPLE: usize = 10;
235 let (start, end) = (self.range.start, self.range.end);
236 let dt = (end - start) * (1.0 / ((N_SAMPLE + 1) * N_SUBSAMPLE) as f64);
237 let mut arclen = 0.0;
238 for i in 0..N_SAMPLE + 1 {
239 for j in 0..N_SUBSAMPLE {
240 let t = start + dt * ((i * N_SUBSAMPLE + j) as f64 + 0.5);
241 let (_, deriv) = source.sample_pt_deriv(t);
242 arclen += deriv.hypot();
243 }
244 if i < N_SAMPLE {
245 self.arcparams.push(arclen);
246 }
247 }
248 let arclen_inv = arclen.recip();
249 for x in &mut self.arcparams {
250 *x *= arclen_inv;
251 }
252 }
253
254 /// Evaluate distance based on arc length parametrization
255 fn eval_arc(&self, c: CubicBez, acc2: f64) -> Option<f64> {
256 // TODO: this could perhaps be tuned.
257 const EPS: f64 = 1e-9;
258 let c_arclen = c.arclen(EPS);
259 let mut max_err2 = 0.0;
260 for (sample, s) in self.samples.iter().zip(&self.arcparams) {
261 let t = c.inv_arclen(c_arclen * s, EPS);
262 let err = sample.p.distance_squared(c.eval(t));
263 max_err2 = err.max(max_err2);
264 if max_err2 > acc2 {
265 return None;
266 }
267 }
268 Some(max_err2)
269 }
270
271 /// Evaluate distance to a cubic approximation.
272 ///
273 /// If distance exceeds stated accuracy, can return None. Note that
274 /// `acc2` is the square of the target.
275 ///
276 /// Returns the squre of the error, which is intended to be a good
277 /// approximation of the Fréchet distance.
278 fn eval_ray(&self, c: CubicBez, acc2: f64) -> Option<f64> {
279 let mut max_err2 = 0.0;
280 for sample in &self.samples {
281 let mut best = acc2 + 1.0;
282 for t in sample.intersect(c) {
283 let err = sample.p.distance_squared(c.eval(t));
284 best = best.min(err);
285 }
286 max_err2 = best.max(max_err2);
287 if max_err2 > acc2 {
288 return None;
289 }
290 }
291 Some(max_err2)
292 }
293
294 fn eval_dist(&mut self, source: &impl ParamCurveFit, c: CubicBez, acc2: f64) -> Option<f64> {
295 // Always compute cheaper distance, hoping for early-out.
296 let ray_dist = self.eval_ray(c, acc2)?;
297 if !self.spicy {
298 return Some(ray_dist);
299 }
300 if self.arcparams.is_empty() {
301 self.compute_arc_params(source);
302 }
303 self.eval_arc(c, acc2)
304 }
305}
306
307/// As described in [Simplifying Bézier paths], strictly optimizing for
308/// Fréchet distance can create bumps. The problem is curves with long
309/// control arms (distance from the control point to the corresponding
310/// endpoint). We mitigate that by applying a penalty as a multiplier to
311/// the measured error (approximate Fréchet distance). This is ReLU-like,
312/// with a value of 1.0 below the elbow, and a given slope above it. The
313/// values here have been determined empirically to give good results.
314///
315/// [Simplifying Bézier paths]:
316/// https://raphlinus.github.io/curves/2023/04/18/bezpath-simplify.html
317const D_PENALTY_ELBOW: f64 = 0.65;
318const D_PENALTY_SLOPE: f64 = 2.0;
319
320/// Fit a single cubic to a range of the source curve.
321///
322/// Returns the cubic segment and the square of the error.
323/// Discussion question: should this be a method on the trait instead?
324pub fn fit_to_cubic(
325 source: &impl ParamCurveFit,
326 range: Range<f64>,
327 accuracy: f64,
328) -> Option<(CubicBez, f64)> {
329 // TODO: handle case where chord is short; return `None` if subdivision
330 // will be useful (ie resulting chord is longer), otherwise simple short
331 // cubic (maybe raised line).
332
333 let start = source.sample_pt_tangent(range.start, 1.0);
334 let end = source.sample_pt_tangent(range.end, -1.0);
335 let d = end.p - start.p;
336 let th = d.atan2();
337 let chord2 = d.hypot2();
338 fn mod_2pi(th: f64) -> f64 {
339 let th_scaled = th * core::f64::consts::FRAC_1_PI * 0.5;
340 core::f64::consts::PI * 2.0 * (th_scaled - th_scaled.round())
341 }
342 let th0 = mod_2pi(start.tangent.atan2() - th);
343 let th1 = mod_2pi(th - end.tangent.atan2());
344
345 let (mut area, mut x, mut y) = source.moment_integrals(range.clone());
346 let (x0, y0) = (start.p.x, start.p.y);
347 let (dx, dy) = (d.x, d.y);
348 // Subtract off area of chord
349 area -= dx * (y0 + 0.5 * dy);
350 // `area` is signed area of closed curve segment.
351 // This quantity is invariant to translation and rotation.
352
353 // Subtract off moment of chord
354 let dy_3 = dy * (1. / 3.);
355 x -= dx * (x0 * y0 + 0.5 * (x0 * dy + y0 * dx) + dy_3 * dx);
356 y -= dx * (y0 * y0 + y0 * dy + dy_3 * dy);
357 // Translate start point to origin; convert raw integrals to moments.
358 x -= x0 * area;
359 y = 0.5 * y - y0 * area;
360 // Rotate into place (this also scales up by chordlength for efficiency).
361 let moment = d.x * x + d.y * y;
362 // `moment` is the chordlength times the x moment of the curve translated
363 // so its start point is on the origin, and rotated so its end point is on the
364 // x axis.
365
366 let chord2_inv = chord2.recip();
367 let unit_area = area * chord2_inv;
368 let mx = moment * chord2_inv.powi(2);
369 // `unit_area` is signed area scaled to unit chord; `mx` is scaled x moment
370
371 let chord = chord2.sqrt();
372 let aff = Affine::translate(start.p.to_vec2()) * Affine::rotate(th) * Affine::scale(chord);
373 let mut curve_dist = CurveDist::from_curve(source, range);
374 let acc2 = accuracy * accuracy;
375 let mut best_c = None;
376 let mut best_err2 = None;
377 for (cand, d0, d1) in cubic_fit(th0, th1, unit_area, mx) {
378 let c = aff * cand;
379 if let Some(err2) = curve_dist.eval_dist(source, c, acc2) {
380 fn scale_f(d: f64) -> f64 {
381 1.0 + (d - D_PENALTY_ELBOW).max(0.0) * D_PENALTY_SLOPE
382 }
383 let scale = scale_f(d0).max(scale_f(d1)).powi(2);
384 let err2 = err2 * scale;
385 if err2 < acc2 && best_err2.map(|best| err2 < best).unwrap_or(true) {
386 best_c = Some(c);
387 best_err2 = Some(err2);
388 }
389 }
390 }
391 match (best_c, best_err2) {
392 (Some(c), Some(err2)) => Some((c, err2)),
393 _ => None,
394 }
395}
396
397/// Returns curves matching area and moment, given unit chord.
398fn cubic_fit(th0: f64, th1: f64, area: f64, mx: f64) -> ArrayVec<(CubicBez, f64, f64), 4> {
399 // Note: maybe we want to take unit vectors instead of angle? Shouldn't
400 // matter much either way though.
401 let (s0, c0) = th0.sin_cos();
402 let (s1, c1) = th1.sin_cos();
403 let a4 = -9.
404 * c0
405 * (((2. * s1 * c1 * c0 + s0 * (2. * c1 * c1 - 1.)) * c0 - 2. * s1 * c1) * c0
406 - c1 * c1 * s0);
407 let a3 = 12.
408 * ((((c1 * (30. * area * c1 - s1) - 15. * area) * c0 + 2. * s0
409 - c1 * s0 * (c1 + 30. * area * s1))
410 * c0
411 + c1 * (s1 - 15. * area * c1))
412 * c0
413 - s0 * c1 * c1);
414 let a2 = 12.
415 * ((((70. * mx + 15. * area) * s1 * s1 + c1 * (9. * s1 - 70. * c1 * mx - 5. * c1 * area))
416 * c0
417 - 5. * s0 * s1 * (3. * s1 - 4. * c1 * (7. * mx + area)))
418 * c0
419 - c1 * (9. * s1 - 70. * c1 * mx - 5. * c1 * area));
420 let a1 = 16.
421 * (((12. * s0 - 5. * c0 * (42. * mx - 17. * area)) * s1
422 - 70. * c1 * (3. * mx - area) * s0
423 - 75. * c0 * c1 * area * area)
424 * s1
425 - 75. * c1 * c1 * area * area * s0);
426 let a0 = 80. * s1 * (42. * s1 * mx - 25. * area * (s1 - c1 * area));
427 // TODO: "roots" is not a good name for this variable, as it also contains
428 // the real part of complex conjugate pairs.
429 let mut roots = ArrayVec::<f64, 4>::new();
430 const EPS: f64 = 1e-12;
431 if a4.abs() > EPS {
432 let a = a3 / a4;
433 let b = a2 / a4;
434 let c = a1 / a4;
435 let d = a0 / a4;
436 if let Some(quads) = factor_quartic_inner(a, b, c, d, false) {
437 for (qc1, qc0) in quads {
438 let qroots = solve_quadratic(qc0, qc1, 1.0);
439 if qroots.is_empty() {
440 // Real part of pair of complex roots
441 roots.push(-0.5 * qc1);
442 } else {
443 roots.extend(qroots);
444 }
445 }
446 }
447 } else if a3.abs() > EPS {
448 roots.extend(solve_cubic(a0, a1, a2, a3));
449 } else {
450 roots.extend(solve_quadratic(a0, a1, a2));
451 }
452
453 let s01 = s0 * c1 + s1 * c0;
454 roots
455 .iter()
456 .filter_map(|&d0| {
457 let (d0, d1) = if d0 > 0.0 {
458 let d1 = (d0 * s0 - area * (10. / 3.)) / (0.5 * d0 * s01 - s1);
459 if d1 > 0.0 {
460 (d0, d1)
461 } else {
462 (s1 / s01, 0.0)
463 }
464 } else {
465 (0.0, s0 / s01)
466 };
467 // We could implement a maximum d value here.
468 if d0 >= 0.0 && d1 >= 0.0 {
469 Some((
470 CubicBez::new(
471 (0.0, 0.0),
472 (d0 * c0, d0 * s0),
473 (1.0 - d1 * c1, d1 * s1),
474 (1.0, 0.0),
475 ),
476 d0,
477 d1,
478 ))
479 } else {
480 None
481 }
482 })
483 .collect()
484}
485
486/// Generate a highly optimized Bézier path that fits the source curve.
487///
488/// This function is still experimental and the signature might change; it's possible
489/// it might become a method on the [`ParamCurveFit`] trait.
490///
491/// This function is considerably slower than [`fit_to_bezpath`], as it computes
492/// optimal subdivision points. Its result is expected to be very close to the optimum
493/// possible Bézier path for the source curve, in that it has a minimal number of curve
494/// segments, and a minimal error over all paths with that number of segments.
495///
496/// See [`fit_to_bezpath`] for an explanation of the `accuracy` parameter.
497pub fn fit_to_bezpath_opt(source: &impl ParamCurveFit, accuracy: f64) -> BezPath {
498 let mut cusps: Vec = Vec::new();
499 let mut path: BezPath = BezPath::new();
500 let mut t0: f64 = 0.0;
501 loop {
502 let t1: f64 = cusps.last().copied().unwrap_or(default:1.0);
503 match fit_to_bezpath_opt_inner(source, accuracy, range:t0..t1, &mut path) {
504 Some(t: f64) => cusps.push(t),
505 None => match cusps.pop() {
506 Some(t: f64) => t0 = t,
507 None => break,
508 },
509 }
510 }
511 path
512}
513
514/// Fit a range without cusps.
515///
516/// On Ok return, range has been added to the path. On Err return, report a cusp (and don't
517/// mutate path).
518fn fit_to_bezpath_opt_inner(
519 source: &impl ParamCurveFit,
520 accuracy: f64,
521 range: Range<f64>,
522 path: &mut BezPath,
523) -> Option<f64> {
524 if let Some(t) = source.break_cusp(range.clone()) {
525 return Some(t);
526 }
527 let err;
528 if let Some((c, err2)) = fit_to_cubic(source, range.clone(), accuracy) {
529 err = err2.sqrt();
530 if err < accuracy {
531 if range.start == 0.0 {
532 path.move_to(c.p0);
533 }
534 path.curve_to(c.p1, c.p2, c.p3);
535 return None;
536 }
537 } else {
538 err = 2.0 * accuracy;
539 }
540 let (mut t0, t1) = (range.start, range.end);
541 let mut n = 0;
542 let last_err;
543 loop {
544 n += 1;
545 match fit_opt_segment(source, accuracy, t0..t1) {
546 FitResult::ParamVal(t) => t0 = t,
547 FitResult::SegmentError(err) => {
548 last_err = err;
549 break;
550 }
551 FitResult::CuspFound(t) => return Some(t),
552 }
553 }
554 t0 = range.start;
555 const EPS: f64 = 1e-9;
556 let f = |x| fit_opt_err_delta(source, x, accuracy, t0..t1, n);
557 let k1 = 0.2 / accuracy;
558 let ya = -err;
559 let yb = accuracy - last_err;
560 let (_, x) = match solve_itp_fallible(f, 0.0, accuracy, EPS, 1, k1, ya, yb) {
561 Ok(x) => x,
562 Err(t) => return Some(t),
563 };
564 //println!("got fit with n={}, err={}", n, x);
565 let path_len = path.elements().len();
566 for i in 0..n {
567 let t1 = if i < n - 1 {
568 match fit_opt_segment(source, x, t0..range.end) {
569 FitResult::ParamVal(t1) => t1,
570 FitResult::SegmentError(_) => range.end,
571 FitResult::CuspFound(t) => {
572 path.truncate(path_len);
573 return Some(t);
574 }
575 }
576 } else {
577 range.end
578 };
579 let (c, _) = fit_to_cubic(source, t0..t1, accuracy).unwrap();
580 if i == 0 && range.start == 0.0 {
581 path.move_to(c.p0);
582 }
583 path.curve_to(c.p1, c.p2, c.p3);
584 t0 = t1;
585 if t0 == range.end {
586 // This is unlikely but could happen when not monotonic.
587 break;
588 }
589 }
590 None
591}
592
593fn measure_one_seg(source: &impl ParamCurveFit, range: Range<f64>, limit: f64) -> Option<f64> {
594 fit_to_cubic(source, range, accuracy:limit).map(|(_, err2: f64)| err2.sqrt())
595}
596
597enum FitResult {
598 /// The parameter (`t`) value that meets the desired accuracy.
599 ParamVal(f64),
600 /// Error of the measured segment.
601 SegmentError(f64),
602 /// The parameter value where a cusp was found.
603 CuspFound(f64),
604}
605
606fn fit_opt_segment(source: &impl ParamCurveFit, accuracy: f64, range: Range<f64>) -> FitResult {
607 if let Some(t: f64) = source.break_cusp(range.clone()) {
608 return FitResult::CuspFound(t);
609 }
610 let missing_err: f64 = accuracy * 2.0;
611 let err: f64 = measure_one_seg(source, range.clone(), accuracy).unwrap_or(default:missing_err);
612 if err <= accuracy {
613 return FitResult::SegmentError(err);
614 }
615 let (t0: f64, t1: f64) = (range.start, range.end);
616 let f: impl Fn(f64) -> Result = |x: f64| {
617 if let Some(t: f64) = source.break_cusp(range.clone()) {
618 return Err(t);
619 }
620 let err: f64 = measure_one_seg(source, t0..x, accuracy).unwrap_or(default:missing_err);
621 Ok(err - accuracy)
622 };
623 const EPS: f64 = 1e-9;
624 let k1: f64 = 2.0 / (t1 - t0);
625 match solve_itp_fallible(f, a:t0, b:t1, EPS, n0:1, k1, -accuracy, yb:err - accuracy) {
626 Ok((t1: f64, _)) => FitResult::ParamVal(t1),
627 Err(t: f64) => FitResult::CuspFound(t),
628 }
629}
630
631// Ok result is delta error (accuracy - error of last seg).
632// Err result is a cusp.
633fn fit_opt_err_delta(
634 source: &impl ParamCurveFit,
635 accuracy: f64,
636 limit: f64,
637 range: Range<f64>,
638 n: usize,
639) -> Result<f64, f64> {
640 let (mut t0: f64, t1: f64) = (range.start, range.end);
641 for _ in 0..n - 1 {
642 t0 = match fit_opt_segment(source, accuracy, range:t0..t1) {
643 FitResult::ParamVal(t0: f64) => t0,
644 // In this case, n - 1 will work, which of course means the error is highly
645 // non-monotonic. We should probably harvest that solution.
646 FitResult::SegmentError(_) => return Ok(accuracy),
647 FitResult::CuspFound(t: f64) => return Err(t),
648 }
649 }
650 let err: f64 = measure_one_seg(source, t0..t1, limit).unwrap_or(default:accuracy * 2.0);
651 Ok(accuracy - err)
652}
653