// Copyright 2022 the Kurbo Authors
// SPDX-License-Identifier: Apache-2.0 OR MIT

//! An implementation of cubic Bézier curve fitting based on a quartic
//! solver making signed area and moment match the source curve.

use core::ops::Range;

use alloc::vec::Vec;

use arrayvec::ArrayVec;

use crate::{
    common::{
        factor_quartic_inner, solve_cubic, solve_itp_fallible, solve_quadratic,
        GAUSS_LEGENDRE_COEFFS_16,
    },
    Affine, BezPath, CubicBez, Line, ParamCurve, ParamCurveArclen, ParamCurveNearest, Point, Vec2,
};

#[cfg(not(feature = "std"))]
use crate::common::FloatFuncs;

/// The source curve for curve fitting.
///
/// This trait is a general representation of curves to be used as input to a curve
/// fitting process. It can represent source curves with cusps and corners, though
/// if the corners are known in advance, it may be better to run curve fitting on
/// subcurves bounded by the corners.
///
/// The trait primarily works by sampling the source curve and computing the position
/// and derivative at each sample. Those derivatives are then used for multiple
/// sub-tasks, including ensuring G1 continuity at subdivision points, computing the
/// area and moment of the curve for curve fitting, and casting rays for evaluation
/// of a distance metric to test accuracy.
///
/// A major motivation is computation of offset curves, which often have cusps, but
/// the presence and location of those cusps is not generally known. It is also
/// intended for conversion between curve types (for example, piecewise Euler spiral
/// or NURBS), and distortion effects such as perspective transform.
///
/// Note general similarities to [`ParamCurve`] but also important differences.
/// Instead of separate [`eval`] and evaluation of derivative, have a single
/// [`sample_pt_deriv`] method which can be more efficient and also handles cusps more
/// robustly. Also there is no method for subsegment, as that is not needed and
/// would be annoying to implement.
///
/// [`ParamCurve`]: crate::ParamCurve
/// [`eval`]: crate::ParamCurve::eval
/// [`sample_pt_deriv`]: ParamCurveFit::sample_pt_deriv
pub trait ParamCurveFit {
    /// Evaluate the curve and its tangent at parameter `t`.
    ///
    /// For a regular curve (one not containing a cusp or corner), the
    /// derivative is a good choice for the tangent vector and the `sign`
    /// parameter can be ignored. Otherwise, the `sign` parameter selects which
    /// side of the discontinuity the tangent will be sampled from.
    ///
    /// Generally `t` is in the range [0..1].
    fn sample_pt_tangent(&self, t: f64, sign: f64) -> CurveFitSample;

    /// Evaluate the point and derivative at parameter `t`.
    ///
    /// In curves with cusps, the derivative can go to zero.
    fn sample_pt_deriv(&self, t: f64) -> (Point, Vec2);

    /// Compute moment integrals.
    ///
    /// This method computes the integrals of y dx, x y dx, and y^2 dx over the
    /// length of this curve. From these integrals it is fairly straightforward
    /// to derive the moments needed for curve fitting.
    ///
    /// A default implementation is proved which does quadrature integration
    /// with Green's theorem, in terms of samples evaluated with
    /// [`sample_pt_deriv`].
    ///
    /// [`sample_pt_deriv`]: ParamCurveFit::sample_pt_deriv
    fn moment_integrals(&self, range: Range<f64>) -> (f64, f64, f64) {
        let t0 = 0.5 * (range.start + range.end);
        let dt = 0.5 * (range.end - range.start);
        let (a, x, y) = GAUSS_LEGENDRE_COEFFS_16
            .iter()
            .map(|(wi, xi)| {
                let t = t0 + xi * dt;
                let (p, d) = self.sample_pt_deriv(t);
                let a = wi * d.x * p.y;
                let x = p.x * a;
                let y = p.y * a;
                (a, x, y)
            })
            .fold((0.0, 0.0, 0.0), |(a0, x0, y0), (a1, x1, y1)| {
                (a0 + a1, x0 + x1, y0 + y1)
            });
        (a * dt, x * dt, y * dt)
    }

    /// Find a cusp or corner within the given range.
    ///
    /// If the range contains a corner or cusp, return it. If there is more
    /// than one such discontinuity, any can be reported, as the function will
    /// be called repeatedly after subdivision of the range.
    ///
    /// Do not report cusps at the endpoints of the range, as this may cause
    /// potentially infinite subdivision. In particular, when a cusp is reported
    /// and this method is called on a subdivided range bounded by the reported
    /// cusp, then the subsequent call should not report a cusp there.
    ///
    /// The definition of what exactly constitutes a cusp is somewhat loose.
    /// If a cusp is missed, then the curve fitting algorithm will attempt to
    /// fit the curve with a smooth curve, which is generally not a disaster
    /// will usually result in more subdivision. Conversely, it might be useful
    /// to report near-cusps, specifically points of curvature maxima where the
    /// curvature is large but still mathematically finite.
    fn break_cusp(&self, range: Range<f64>) -> Option<f64>;
}

/// A sample point of a curve for fitting.
pub struct CurveFitSample {
    /// A point on the curve at the sample location.
    pub p: Point,
    /// A vector tangent to the curve at the sample location.
    pub tangent: Vec2,
}

impl CurveFitSample {
    /// Intersect a ray orthogonal to the tangent with the given cubic.
    ///
    /// Returns a vector of `t` values on the cubic.
    fn intersect(&self, c: CubicBez) -> ArrayVec<f64, 3> {
        let p1 = 3.0 * (c.p1 - c.p0);
        let p2 = 3.0 * c.p2.to_vec2() - 6.0 * c.p1.to_vec2() + 3.0 * c.p0.to_vec2();
        let p3 = (c.p3 - c.p0) - 3.0 * (c.p2 - c.p1);
        let c0 = (c.p0 - self.p).dot(self.tangent);
        let c1 = p1.dot(self.tangent);
        let c2 = p2.dot(self.tangent);
        let c3 = p3.dot(self.tangent);
        solve_cubic(c0, c1, c2, c3)
            .into_iter()
            .filter(|t| (0.0..=1.0).contains(t))
            .collect()
    }
}

/// Generate a Bézier path that fits the source curve.
///
/// This function is still experimental and the signature might change; it's possible
/// it might become a method on the [`ParamCurveFit`] trait.
///
/// This function recursively subdivides the curve in half by the parameter when the
/// accuracy is not met. That gives a reasonably optimized result but not necessarily
/// the minimum number of segments.
///
/// In general, the resulting Bézier path should have a Fréchet distance less than
/// the provided `accuracy` parameter. However, this is not a rigorous guarantee, as
/// the error metric is computed approximately.
///
/// This function is intended for use when the source curve is piecewise continuous,
/// with the discontinuities reported by the cusp method. In applications (such as
/// stroke expansion) where this property may not hold, it is up to the client to
/// detect and handle such cases. Even so, best effort is made to avoid infinite
/// subdivision.
///
/// When a higher degree of optimization is desired (at considerably more runtime cost),
/// consider [`fit_to_bezpath_opt`] instead.
pub fn fit_to_bezpath(source: &impl ParamCurveFit, accuracy: f64) -> BezPath {
    let mut path = BezPath::new();
    fit_to_bezpath_rec(source, 0.0..1.0, accuracy, &mut path);
    path
}

// Discussion question: possibly should take endpoint samples, to avoid
// duplication of that work.
fn fit_to_bezpath_rec(
    source: &impl ParamCurveFit,
    range: Range<f64>,
    accuracy: f64,
    path: &mut BezPath,
) {
    let start = range.start;
    let end = range.end;
    let start_p = source.sample_pt_tangent(range.start, 1.0).p;
    let end_p = source.sample_pt_tangent(range.end, -1.0).p;
    if start_p.distance_squared(end_p) <= accuracy * accuracy {
        if let Some((c, _)) = try_fit_line(source, accuracy, range, start_p, end_p) {
            if path.is_empty() {
                path.move_to(c.p0);
            }
            path.curve_to(c.p1, c.p2, c.p3);
            return;
        }
    }
    let t = if let Some(t) = source.break_cusp(start..end) {
        t
    } else if let Some((c, _)) = fit_to_cubic(source, start..end, accuracy) {
        if path.is_empty() {
            path.move_to(c.p0);
        }
        path.curve_to(c.p1, c.p2, c.p3);
        return;
    } else {
        // A smarter approach is possible than midpoint subdivision, but would be
        // a significant increase in complexity.
        0.5 * (start + end)
    };
    if t == start || t == end {
        // infinite recursion, just draw a line
        let p1 = start_p.lerp(end_p, 1.0 / 3.0);
        let p2 = end_p.lerp(start_p, 1.0 / 3.0);
        if path.is_empty() {
            path.move_to(start_p);
        }
        path.curve_to(p1, p2, end_p);
        return;
    }
    fit_to_bezpath_rec(source, start..t, accuracy, path);
    fit_to_bezpath_rec(source, t..end, accuracy, path);
}

const N_SAMPLE: usize = 20;

/// An acceleration structure for estimating curve distance
struct CurveDist {
    samples: ArrayVec<CurveFitSample, N_SAMPLE>,
    arcparams: ArrayVec<f64, N_SAMPLE>,
    range: Range<f64>,
    /// A "spicy" curve is one with potentially extreme curvature variation,
    /// so use arc length measurement for better accuracy.
    spicy: bool,
}

impl CurveDist {
    fn from_curve(source: &impl ParamCurveFit, range: Range<f64>) -> Self {
        let step = (range.end - range.start) * (1.0 / (N_SAMPLE + 1) as f64);
        let mut last_tan = None;
        let mut spicy = false;
        const SPICY_THRESH: f64 = 0.2;
        let mut samples = ArrayVec::new();
        for i in 0..N_SAMPLE + 2 {
            let sample = source.sample_pt_tangent(range.start + i as f64 * step, 1.0);
            if let Some(last_tan) = last_tan {
                let cross = sample.tangent.cross(last_tan);
                let dot = sample.tangent.dot(last_tan);
                if cross.abs() > SPICY_THRESH * dot.abs() {
                    spicy = true;
                }
            }
            last_tan = Some(sample.tangent);
            if i > 0 && i < N_SAMPLE + 1 {
                samples.push(sample);
            }
        }
        CurveDist {
            samples,
            arcparams: Default::default(),
            range,
            spicy,
        }
    }

    fn compute_arc_params(&mut self, source: &impl ParamCurveFit) {
        const N_SUBSAMPLE: usize = 10;
        let (start, end) = (self.range.start, self.range.end);
        let dt = (end - start) * (1.0 / ((N_SAMPLE + 1) * N_SUBSAMPLE) as f64);
        let mut arclen = 0.0;
        for i in 0..N_SAMPLE + 1 {
            for j in 0..N_SUBSAMPLE {
                let t = start + dt * ((i * N_SUBSAMPLE + j) as f64 + 0.5);
                let (_, deriv) = source.sample_pt_deriv(t);
                arclen += deriv.hypot();
            }
            if i < N_SAMPLE {
                self.arcparams.push(arclen);
            }
        }
        let arclen_inv = arclen.recip();
        for x in &mut self.arcparams {
            *x *= arclen_inv;
        }
    }

    /// Evaluate distance based on arc length parametrization
    fn eval_arc(&self, c: CubicBez, acc2: f64) -> Option<f64> {
        // TODO: this could perhaps be tuned.
        const EPS: f64 = 1e-9;
        let c_arclen = c.arclen(EPS);
        let mut max_err2 = 0.0;
        for (sample, s) in self.samples.iter().zip(&self.arcparams) {
            let t = c.inv_arclen(c_arclen * s, EPS);
            let err = sample.p.distance_squared(c.eval(t));
            max_err2 = err.max(max_err2);
            if max_err2 > acc2 {
                return None;
            }
        }
        Some(max_err2)
    }

    /// Evaluate distance to a cubic approximation.
    ///
    /// If distance exceeds stated accuracy, can return `None`. Note that
    /// `acc2` is the square of the target.
    ///
    /// Returns the square of the error, which is intended to be a good
    /// approximation of the Fréchet distance.
    fn eval_ray(&self, c: CubicBez, acc2: f64) -> Option<f64> {
        let mut max_err2 = 0.0;
        for sample in &self.samples {
            let mut best = acc2 + 1.0;
            for t in sample.intersect(c) {
                let err = sample.p.distance_squared(c.eval(t));
                best = best.min(err);
            }
            max_err2 = best.max(max_err2);
            if max_err2 > acc2 {
                return None;
            }
        }
        Some(max_err2)
    }

    fn eval_dist(&mut self, source: &impl ParamCurveFit, c: CubicBez, acc2: f64) -> Option<f64> {
        // Always compute cheaper distance, hoping for early-out.
        let ray_dist = self.eval_ray(c, acc2)?;
        if !self.spicy {
            return Some(ray_dist);
        }
        if self.arcparams.is_empty() {
            self.compute_arc_params(source);
        }
        self.eval_arc(c, acc2)
    }
}

/// As described in [Simplifying Bézier paths], strictly optimizing for
/// Fréchet distance can create bumps. The problem is curves with long
/// control arms (distance from the control point to the corresponding
/// endpoint). We mitigate that by applying a penalty as a multiplier to
/// the measured error (approximate Fréchet distance). This is ReLU-like,
/// with a value of 1.0 below the elbow, and a given slope above it. The
/// values here have been determined empirically to give good results.
///
/// [Simplifying Bézier paths]:
/// https://raphlinus.github.io/curves/2023/04/18/bezpath-simplify.html
const D_PENALTY_ELBOW: f64 = 0.65;
const D_PENALTY_SLOPE: f64 = 2.0;

/// Try fitting a line.
///
/// This is especially useful for very short chords, in which the standard
/// cubic fit is not numerically stable. The tangents are not considered, so
/// it's useful in cusp and near-cusp situations where the tangents are not
/// reliable, as well.
///
/// Returns the line raised to a cubic and the error, if within tolerance.
fn try_fit_line(
    source: &impl ParamCurveFit,
    accuracy: f64,
    range: Range<f64>,
    start: Point,
    end: Point,
) -> Option<(CubicBez, f64)> {
    let acc2 = accuracy * accuracy;
    let chord_l = Line::new(start, end);
    const SHORT_N: usize = 7;
    let mut max_err2 = 0.0;
    let dt = (range.end - range.start) / (SHORT_N + 1) as f64;
    for i in 0..SHORT_N {
        let t = range.start + (i + 1) as f64 * dt;
        let p = source.sample_pt_deriv(t).0;
        let err2 = chord_l.nearest(p, accuracy).distance_sq;
        if err2 > acc2 {
            // Not in tolerance; likely subdivision will help.
            return None;
        }
        max_err2 = err2.max(max_err2);
    }
    let p1 = start.lerp(end, 1.0 / 3.0);
    let p2 = end.lerp(start, 1.0 / 3.0);
    let c = CubicBez::new(start, p1, p2, end);
    Some((c, max_err2))
}

/// Fit a single cubic to a range of the source curve.
///
/// Returns the cubic segment and the square of the error.
/// Discussion question: should this be a method on the trait instead?
pub fn fit_to_cubic(
    source: &impl ParamCurveFit,
    range: Range<f64>,
    accuracy: f64,
) -> Option<(CubicBez, f64)> {
    let start = source.sample_pt_tangent(range.start, 1.0);
    let end = source.sample_pt_tangent(range.end, -1.0);
    let d = end.p - start.p;
    let chord2 = d.hypot2();
    let acc2 = accuracy * accuracy;
    if chord2 <= acc2 {
        // Special case very short chords; try to fit a line.
        return try_fit_line(source, accuracy, range, start.p, end.p);
    }
    let th = d.atan2();
    fn mod_2pi(th: f64) -> f64 {
        let th_scaled = th * core::f64::consts::FRAC_1_PI * 0.5;
        core::f64::consts::PI * 2.0 * (th_scaled - th_scaled.round())
    }
    let th0 = mod_2pi(start.tangent.atan2() - th);
    let th1 = mod_2pi(th - end.tangent.atan2());

    let (mut area, mut x, mut y) = source.moment_integrals(range.clone());
    let (x0, y0) = (start.p.x, start.p.y);
    let (dx, dy) = (d.x, d.y);
    // Subtract off area of chord
    area -= dx * (y0 + 0.5 * dy);
    // `area` is signed area of closed curve segment.
    // This quantity is invariant to translation and rotation.

    // Subtract off moment of chord
    let dy_3 = dy * (1. / 3.);
    x -= dx * (x0 * y0 + 0.5 * (x0 * dy + y0 * dx) + dy_3 * dx);
    y -= dx * (y0 * y0 + y0 * dy + dy_3 * dy);
    // Translate start point to origin; convert raw integrals to moments.
    x -= x0 * area;
    y = 0.5 * y - y0 * area;
    // Rotate into place (this also scales up by chordlength for efficiency).
    let moment = d.x * x + d.y * y;
    // `moment` is the chordlength times the x moment of the curve translated
    // so its start point is on the origin, and rotated so its end point is on the
    // x axis.

    let chord2_inv = chord2.recip();
    let unit_area = area * chord2_inv;
    let mx = moment * chord2_inv.powi(2);
    // `unit_area` is signed area scaled to unit chord; `mx` is scaled x moment

    let chord = chord2.sqrt();
    let aff = Affine::translate(start.p.to_vec2()) * Affine::rotate(th) * Affine::scale(chord);
    let mut curve_dist = CurveDist::from_curve(source, range);
    let mut best_c = None;
    let mut best_err2 = None;
    for (cand, d0, d1) in cubic_fit(th0, th1, unit_area, mx) {
        let c = aff * cand;
        if let Some(err2) = curve_dist.eval_dist(source, c, acc2) {
            fn scale_f(d: f64) -> f64 {
                1.0 + (d - D_PENALTY_ELBOW).max(0.0) * D_PENALTY_SLOPE
            }
            let scale = scale_f(d0).max(scale_f(d1)).powi(2);
            let err2 = err2 * scale;
            if err2 < acc2 && best_err2.map(|best| err2 < best).unwrap_or(true) {
                best_c = Some(c);
                best_err2 = Some(err2);
            }
        }
    }
    match (best_c, best_err2) {
        (Some(c), Some(err2)) => Some((c, err2)),
        _ => None,
    }
}

/// Returns curves matching area and moment, given unit chord.
fn cubic_fit(th0: f64, th1: f64, area: f64, mx: f64) -> ArrayVec<(CubicBez, f64, f64), 4> {
    // Note: maybe we want to take unit vectors instead of angle? Shouldn't
    // matter much either way though.
    let (s0, c0) = th0.sin_cos();
    let (s1, c1) = th1.sin_cos();
    let a4 = -9.
        * c0
        * (((2. * s1 * c1 * c0 + s0 * (2. * c1 * c1 - 1.)) * c0 - 2. * s1 * c1) * c0
            - c1 * c1 * s0);
    let a3 = 12.
        * ((((c1 * (30. * area * c1 - s1) - 15. * area) * c0 + 2. * s0
            - c1 * s0 * (c1 + 30. * area * s1))
            * c0
            + c1 * (s1 - 15. * area * c1))
            * c0
            - s0 * c1 * c1);
    let a2 = 12.
        * ((((70. * mx + 15. * area) * s1 * s1 + c1 * (9. * s1 - 70. * c1 * mx - 5. * c1 * area))
            * c0
            - 5. * s0 * s1 * (3. * s1 - 4. * c1 * (7. * mx + area)))
            * c0
            - c1 * (9. * s1 - 70. * c1 * mx - 5. * c1 * area));
    let a1 = 16.
        * (((12. * s0 - 5. * c0 * (42. * mx - 17. * area)) * s1
            - 70. * c1 * (3. * mx - area) * s0
            - 75. * c0 * c1 * area * area)
            * s1
            - 75. * c1 * c1 * area * area * s0);
    let a0 = 80. * s1 * (42. * s1 * mx - 25. * area * (s1 - c1 * area));
    // TODO: "roots" is not a good name for this variable, as it also contains
    // the real part of complex conjugate pairs.
    let mut roots = ArrayVec::<f64, 4>::new();
    const EPS: f64 = 1e-12;
    if a4.abs() > EPS {
        let a = a3 / a4;
        let b = a2 / a4;
        let c = a1 / a4;
        let d = a0 / a4;
        if let Some(quads) = factor_quartic_inner(a, b, c, d, false) {
            for (qc1, qc0) in quads {
                let qroots = solve_quadratic(qc0, qc1, 1.0);
                if qroots.is_empty() {
                    // Real part of pair of complex roots
                    roots.push(-0.5 * qc1);
                } else {
                    roots.extend(qroots);
                }
            }
        }
    } else if a3.abs() > EPS {
        roots.extend(solve_cubic(a0, a1, a2, a3));
    } else if a2.abs() > EPS || a1.abs() > EPS || a0.abs() > EPS {
        roots.extend(solve_quadratic(a0, a1, a2));
    } else {
        return [(
            CubicBez::new((0.0, 0.0), (1. / 3., 0.0), (2. / 3., 0.0), (1., 0.0)),
            1f64 / 3.,
            1f64 / 3.,
        )]
        .into_iter()
        .collect();
    }

    let s01 = s0 * c1 + s1 * c0;
    roots
        .iter()
        .filter_map(|&d0| {
            let (d0, d1) = if d0 > 0.0 {
                let d1 = (d0 * s0 - area * (10. / 3.)) / (0.5 * d0 * s01 - s1);
                if d1 > 0.0 {
                    (d0, d1)
                } else {
                    (s1 / s01, 0.0)
                }
            } else {
                (0.0, s0 / s01)
            };
            // We could implement a maximum d value here.
            if d0 >= 0.0 && d1 >= 0.0 {
                Some((
                    CubicBez::new(
                        (0.0, 0.0),
                        (d0 * c0, d0 * s0),
                        (1.0 - d1 * c1, d1 * s1),
                        (1.0, 0.0),
                    ),
                    d0,
                    d1,
                ))
            } else {
                None
            }
        })
        .collect()
}

/// Generate a highly optimized Bézier path that fits the source curve.
///
/// This function is still experimental and the signature might change; it's possible
/// it might become a method on the [`ParamCurveFit`] trait.
///
/// This function is considerably slower than [`fit_to_bezpath`], as it computes
/// optimal subdivision points. Its result is expected to be very close to the optimum
/// possible Bézier path for the source curve, in that it has a minimal number of curve
/// segments, and a minimal error over all paths with that number of segments.
///
/// See [`fit_to_bezpath`] for an explanation of the `accuracy` parameter.
pub fn fit_to_bezpath_opt(source: &impl ParamCurveFit, accuracy: f64) -> BezPath {
    let mut cusps = Vec::new();
    let mut path = BezPath::new();
    let mut t0 = 0.0;
    loop {
        let t1 = cusps.last().copied().unwrap_or(1.0);
        match fit_to_bezpath_opt_inner(source, accuracy, t0..t1, &mut path) {
            Some(t) => cusps.push(t),
            None => match cusps.pop() {
                Some(t) => t0 = t,
                None => break,
            },
        }
    }
    path
}

/// Fit a range without cusps.
///
/// On Ok return, range has been added to the path. On Err return, report a cusp (and don't
/// mutate path).
fn fit_to_bezpath_opt_inner(
    source: &impl ParamCurveFit,
    accuracy: f64,
    range: Range<f64>,
    path: &mut BezPath,
) -> Option<f64> {
    if let Some(t) = source.break_cusp(range.clone()) {
        return Some(t);
    }
    let err;
    if let Some((c, err2)) = fit_to_cubic(source, range.clone(), accuracy) {
        err = err2.sqrt();
        if err < accuracy {
            if range.start == 0.0 {
                path.move_to(c.p0);
            }
            path.curve_to(c.p1, c.p2, c.p3);
            return None;
        }
    } else {
        err = 2.0 * accuracy;
    }
    let (mut t0, t1) = (range.start, range.end);
    let mut n = 0;
    let last_err;
    loop {
        n += 1;
        match fit_opt_segment(source, accuracy, t0..t1) {
            FitResult::ParamVal(t) => t0 = t,
            FitResult::SegmentError(err) => {
                last_err = err;
                break;
            }
            FitResult::CuspFound(t) => return Some(t),
        }
    }
    t0 = range.start;
    const EPS: f64 = 1e-9;
    let f = |x| fit_opt_err_delta(source, x, accuracy, t0..t1, n);
    let k1 = 0.2 / accuracy;
    let ya = -err;
    let yb = accuracy - last_err;
    let (_, x) = match solve_itp_fallible(f, 0.0, accuracy, EPS, 1, k1, ya, yb) {
        Ok(x) => x,
        Err(t) => return Some(t),
    };
    //println!("got fit with n={}, err={}", n, x);
    let path_len = path.elements().len();
    for i in 0..n {
        let t1 = if i < n - 1 {
            match fit_opt_segment(source, x, t0..range.end) {
                FitResult::ParamVal(t1) => t1,
                FitResult::SegmentError(_) => range.end,
                FitResult::CuspFound(t) => {
                    path.truncate(path_len);
                    return Some(t);
                }
            }
        } else {
            range.end
        };
        let (c, _) = fit_to_cubic(source, t0..t1, accuracy).unwrap();
        if i == 0 && range.start == 0.0 {
            path.move_to(c.p0);
        }
        path.curve_to(c.p1, c.p2, c.p3);
        t0 = t1;
        if t0 == range.end {
            // This is unlikely but could happen when not monotonic.
            break;
        }
    }
    None
}

fn measure_one_seg(source: &impl ParamCurveFit, range: Range<f64>, limit: f64) -> Option<f64> {
    fit_to_cubic(source, range, limit).map(|(_, err2)| err2.sqrt())
}

enum FitResult {
    /// The parameter (`t`) value that meets the desired accuracy.
    ParamVal(f64),
    /// Error of the measured segment.
    SegmentError(f64),
    /// The parameter value where a cusp was found.
    CuspFound(f64),
}

fn fit_opt_segment(source: &impl ParamCurveFit, accuracy: f64, range: Range<f64>) -> FitResult {
    if let Some(t) = source.break_cusp(range.clone()) {
        return FitResult::CuspFound(t);
    }
    let missing_err = accuracy * 2.0;
    let err = measure_one_seg(source, range.clone(), accuracy).unwrap_or(missing_err);
    if err <= accuracy {
        return FitResult::SegmentError(err);
    }
    let (t0, t1) = (range.start, range.end);
    let f = |x| {
        if let Some(t) = source.break_cusp(range.clone()) {
            return Err(t);
        }
        let err = measure_one_seg(source, t0..x, accuracy).unwrap_or(missing_err);
        Ok(err - accuracy)
    };
    const EPS: f64 = 1e-9;
    let k1 = 2.0 / (t1 - t0);
    match solve_itp_fallible(f, t0, t1, EPS, 1, k1, -accuracy, err - accuracy) {
        Ok((t1, _)) => FitResult::ParamVal(t1),
        Err(t) => FitResult::CuspFound(t),
    }
}

// Ok result is delta error (accuracy - error of last seg).
// Err result is a cusp.
fn fit_opt_err_delta(
    source: &impl ParamCurveFit,
    accuracy: f64,
    limit: f64,
    range: Range<f64>,
    n: usize,
) -> Result<f64, f64> {
    let (mut t0, t1) = (range.start, range.end);
    for _ in 0..n - 1 {
        t0 = match fit_opt_segment(source, accuracy, t0..t1) {
            FitResult::ParamVal(t0) => t0,
            // In this case, n - 1 will work, which of course means the error is highly
            // non-monotonic. We should probably harvest that solution.
            FitResult::SegmentError(_) => return Ok(accuracy),
            FitResult::CuspFound(t) => return Err(t),
        }
    }
    let err = measure_one_seg(source, t0..t1, limit).unwrap_or(accuracy * 2.0);
    Ok(accuracy - err)
}