kurbo/

offset.rs

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

//! Computation of offset curves of cubic Béziers.
//!
//! The main algorithm in this module is a new technique designed for robustness
//! and speed. The details are involved; hopefully there will be a paper.

use arrayvec::ArrayVec;

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

use crate::{
    common::{solve_itp, solve_quadratic},
    BezPath, CubicBez, ParamCurve, ParamCurveDeriv, PathSeg, Point, QuadBez, Vec2,
};

/// State used for computing an offset curve of a single cubic.
struct CubicOffset {
    /// The cubic being offset. This has been regularized.
    c: CubicBez,
    /// The derivative of `c`.
    q: QuadBez,
    /// The offset distance (same as the argument).
    d: f64,
    /// `c0 + c1 t + c2 t^2` is the cross product of second and first
    /// derivatives of the underlying cubic, multiplied by the offset.
    /// This is used for computing cusps on the offset curve.
    ///
    /// Note that given a curve `c(t)`, its signed curvature is
    /// `c''(t) x c'(t) / ||c'(t)||^3`. See also [`Self::cusp_sign`].
    c0: f64,
    c1: f64,
    c2: f64,
    /// The tolerance (same as the argument).
    tolerance: f64,
}

// We never let cusp values have an absolute value smaller than
// this. When a cusp is found, determine its sign and use this value.
const CUSP_EPSILON: f64 = 1e-12;

/// Number of points for least-squares fit and error evaluation.
///
/// This value is a tradeoff between accuracy and performance. The main risk in it
/// being to small is under-sampling the error and thus letting excessive error
/// slip through. That said, the "arc drawing" approach is designed to be robust
/// and not generate approximate results with narrow error peaks, even in near-cusp
/// "J" shape curves.
const N_LSE: usize = 8;

/// The proportion of transverse error that is blended in the least-squares logic.
const BLEND: f64 = 1e-3;

/// Maximum recursion depth.
///
/// Recursion is bounded to this depth, so the total number of subdivisions will
/// not exceed two to this power.
///
/// This is primarily a "belt and suspenders" robustness guard. In normal operation,
/// the recursion bound should never be reached, as accuracy improves very quickly
/// on subdivision. For unreasonably large coordinate values or small tolerances, it
/// is possible, and in those cases the result will be out of tolerance.
///
/// Perhaps should be configurable.
const MAX_DEPTH: usize = 8;

/// State local to a subdivision
struct OffsetRec {
    t0: f64,
    t1: f64,
    // unit tangent at t0
    utan0: Vec2,
    // unit tangent at t1
    utan1: Vec2,
    cusp0: f64,
    cusp1: f64,
    /// Recursion depth
    depth: usize,
}

/// Result of error evaluation
struct ErrEval {
    /// Maximum detected error
    err_squared: f64,
    /// Unit normals sampled uniformly across approximation
    unorms: [Vec2; N_LSE],
    /// Difference between point on source curve and normal from approximation.
    err_vecs: [Vec2; N_LSE],
}

/// Result of subdivision
struct SubdivisionPoint {
    /// Source curve t value at subdivision point
    t: f64,
    /// Unit tangent at subdivision point
    utan: Vec2,
}

/// Compute an approximate offset curve.
///
/// The parallel curve of `c` offset by `d` is written to the `result` path.
///
/// There is a fair amount of attention to robustness, but this method is not suitable
/// for degenerate cubics with entirely co-linear control points. Those cases should be
/// handled before calling this function, by replacing them with linear segments.
pub fn offset_cubic(c: CubicBez, d: f64, tolerance: f64, result: &mut BezPath) {
    result.truncate(0);
    // A tuning parameter for regularization. A value too large may distort the curve,
    // while a value too small may fail to generate smooth curves. This is a somewhat
    // arbitrary value, and should be revisited.
    const DIM_TUNE: f64 = 0.25;
    // We use regularization to perturb the curve to avoid *interior* zero-derivative
    // cusps. There is robustness logic in place to handle zero derivatives at the
    // endpoints.
    //
    // As a performance note, it might be a good idea to move regularization and
    // tangent determination to the caller, as those computations are the same for both
    // signs of `d`.
    let c_regularized = c.regularize_cusp(tolerance * DIM_TUNE);
    let co = CubicOffset::new(c_regularized, d, tolerance);
    let (tan0, tan1) = PathSeg::Cubic(c).tangents();
    let utan0 = tan0.normalize();
    let utan1 = tan1.normalize();
    let cusp0 = co.endpoint_cusp(co.q.p0, co.c0);
    let cusp1 = co.endpoint_cusp(co.q.p2, co.c0 + co.c1 + co.c2);
    result.move_to(c.p0 + d * utan0.turn_90());
    let rec = OffsetRec::new(0., 1., utan0, utan1, cusp0, cusp1, 0);
    co.offset_rec(&rec, result);
}

impl CubicOffset {
    /// Create a new curve from Bézier segment and offset.
    fn new(c: CubicBez, d: f64, tolerance: f64) -> Self {
        let q = c.deriv();
        let d2 = 2.0 * d;
        let p1xp0 = q.p1.to_vec2().cross(q.p0.to_vec2());
        let p2xp0 = q.p2.to_vec2().cross(q.p0.to_vec2());
        let p2xp1 = q.p2.to_vec2().cross(q.p1.to_vec2());
        CubicOffset {
            c,
            q,
            d,
            c0: d2 * p1xp0,
            c1: d2 * (p2xp0 - 2.0 * p1xp0),
            c2: d2 * (p2xp1 - p2xp0 + p1xp0),
            tolerance,
        }
    }

    /// Compute curvature of the source curve times offset plus 1.
    ///
    /// This quantity is called "cusp" because cusps appear in the offset curve
    /// where this value crosses zero. This is based on the geometric property
    /// that the offset curve has a cusp when the radius of curvature of the
    /// source curve is equal to the offset curve's distance.
    ///
    /// Note: there is a potential division by zero when the derivative vanishes.
    /// We avoid doing so for interior points by regularizing the cubic beforehand.
    /// We avoid doing so for endpoints by calling `endpoint_cusp` instead.
    fn cusp_sign(&self, t: f64) -> f64 {
        let ds2 = self.q.eval(t).to_vec2().hypot2();
        ((self.c2 * t + self.c1) * t + self.c0) / (ds2 * ds2.sqrt()) + 1.0
    }

    /// Compute cusp value of endpoint.
    ///
    /// This is a special case of [`Self::cusp_sign`]. For the start point, `tan` should be
    /// the start point tangent and `y` should be `c0`. For the end point, `tan` should be
    /// the end point tangent and `y` should be `c0 + c1 + c2`.
    ///
    /// This is just evaluating the polynomial at t=0 and t=1.
    ///
    /// See [`Self::cusp_sign`] for a description of what "cusp value" means.
    fn endpoint_cusp(&self, tan: Point, y: f64) -> f64 {
        // Robustness to avoid divide-by-zero when derivatives vanish
        const TAN_DIST_EPSILON: f64 = 1e-12;
        let tan_dist = tan.to_vec2().hypot().max(TAN_DIST_EPSILON);
        let rsqrt = 1.0 / tan_dist;
        y * (rsqrt * rsqrt * rsqrt) + 1.0
    }

    /// Primary entry point for recursive subdivision.
    ///
    /// At a high level, this method determines whether subdivision is necessary. If
    /// so, it determines a subdivision point and then recursively calls itself on
    /// both subdivisions. If not, it computes a single cubic Bézier to approximate
    /// the offset curve, and appends it to `result`.
    fn offset_rec(&self, rec: &OffsetRec, result: &mut BezPath) {
        // First, determine whether the offset curve contains a cusp. If the sign
        // of the cusp value (curvature times offset plus 1) is different at the
        // subdivision endpoints, then there is definitely a cusp inside. Find it and
        // subdivide there.
        //
        // Note that there's a possibility the curve has two (or, potentially, any
        // even number). We don't rigorously check for this case; if the measured
        // error comes in under the tolerance, we simply accept it. Otherwise, in
        // the common case we expect to detect a sign crossing from the new
        // subdivision point.
        if rec.cusp0 * rec.cusp1 < 0.0 {
            let a = rec.t0;
            let b = rec.t1;
            let s = rec.cusp1.signum();
            let f = |t| s * self.cusp_sign(t);
            let k1 = 0.2 / (b - a);
            const ITP_EPS: f64 = 1e-12;
            let t = solve_itp(f, a, b, ITP_EPS, 1, k1, s * rec.cusp0, s * rec.cusp1);
            // TODO(robustness): If we're unlucky, there will be 3 cusps between t0
            // and t1, and the solver will land on the middle one. In that case, the
            // derivative on cusp will be the opposite sign as expected.
            //
            // If we're even more unlucky, there is a second-order cusp, both zero
            // cusp value and zero derivative.
            let utan_t = self.q.eval(t).to_vec2().normalize();
            let cusp_t_minus = CUSP_EPSILON.copysign(rec.cusp0);
            let cusp_t_plus = CUSP_EPSILON.copysign(rec.cusp1);
            self.subdivide(rec, result, t, utan_t, cusp_t_minus, cusp_t_plus);
            return;
        }
        // We determine the first approximation to the offset curve.
        let (mut a, mut b) = self.draw_arc(rec);
        let dt = (rec.t1 - rec.t0) * (1.0 / (N_LSE + 1) as f64);
        // These represent t values on the source curve.
        let mut ts = core::array::from_fn(|i| rec.t0 + (i + 1) as f64 * dt);
        let mut c_approx = self.apply(rec, a, b);
        let err_init = self.eval_err(rec, c_approx, &mut ts);
        let mut err = err_init;
        // Number of least-squares refinement steps. More gives a smaller
        // error, but takes more time.
        const N_REFINE: usize = 2;
        for _ in 0..N_REFINE {
            if err.err_squared <= self.tolerance * self.tolerance {
                break;
            }
            let (a2, b2) = self.refine_least_squares(rec, a, b, &err);
            let c_approx2 = self.apply(rec, a2, b2);
            let err2 = self.eval_err(rec, c_approx2, &mut ts);
            if err2.err_squared >= err.err_squared {
                break;
            }
            err = err2;
            (a, b) = (a2, b2);
            c_approx = c_approx2;
        }
        if rec.depth < MAX_DEPTH && err.err_squared > self.tolerance * self.tolerance {
            let SubdivisionPoint { t, utan } = self.find_subdivision_point(rec);
            // TODO(robustness): if cusp is extremely near zero, then assign epsilon
            // with alternate signs based on derivative of cusp.
            let cusp = self.cusp_sign(t);
            self.subdivide(rec, result, t, utan, cusp, cusp);
        } else {
            result.curve_to(c_approx.p1, c_approx.p2, c_approx.p3);
        }
    }

    /// Recursively subdivide.
    ///
    /// In the case of subdividing at a cusp, the cusp value at the subdivision point
    /// is mathematically zero, but in those cases we treat it as a signed infinitesimal
    /// value representing the values at t minus epsilon and t plus epsilon.
    ///
    /// Note that unit tangents are passed down explicitly. In the general case, they
    /// are equal to the derivative (evaluated at that t value) normalized to unit
    /// length, but in cases where the derivative is near-zero, they are computed more
    /// robustly.
    fn subdivide(
        &self,
        rec: &OffsetRec,
        result: &mut BezPath,
        t: f64,
        utan_t: Vec2,
        cusp_t_minus: f64,
        cusp_t_plus: f64,
    ) {
        let rec0 = OffsetRec::new(
            rec.t0,
            t,
            rec.utan0,
            utan_t,
            rec.cusp0,
            cusp_t_minus,
            rec.depth + 1,
        );
        self.offset_rec(&rec0, result);
        let rec1 = OffsetRec::new(
            t,
            rec.t1,
            utan_t,
            rec.utan1,
            cusp_t_plus,
            rec.cusp1,
            rec.depth + 1,
        );
        self.offset_rec(&rec1, result);
    }

    /// Convert from (a, b) parameter space to the approximate cubic Bézier.
    ///
    /// The offset approximation can be considered `B(t) + d * D(t)`, where `D(t)`
    /// is roughly a unit vector in the direction of the unit normal of the source
    /// curve. (The word "roughly" is appropriate because transverse error may
    /// cancel out normal error, resulting in a lower error than either alone).
    /// The endpoints of `D(t)` must be the unit normals of the source curve, and
    /// the endpoint tangents of `D(t)` must tangent to the endpoint tangents of
    /// the source curve, to ensure G1 continuity.
    ///
    /// The (a, b) parameters refer to the magnitude of the vector from the endpoint
    /// to the corresponding control point in `D(t)`, the direction being determined
    /// by the unit tangent.
    ///
    /// When the candidate solution would lead to negative distance from the
    /// endpoint to the control point, that distance is clamped to zero. Otherwise
    /// such solutions should be considered invalid, and have the unpleasant
    /// property of sometimes passing error tolerance checks.
    fn apply(&self, rec: &OffsetRec, a: f64, b: f64) -> CubicBez {
        // wondering if p0 and p3 should be in rec
        // Scale factor from derivatives to displacements
        let s = (1. / 3.) * (rec.t1 - rec.t0);
        let p0 = self.c.eval(rec.t0) + self.d * rec.utan0.turn_90();
        let l0 = s * self.q.eval(rec.t0).to_vec2().length() + a * self.d;
        let mut p1 = p0;
        if l0 * rec.cusp0 > 0.0 {
            p1 += l0 * rec.utan0;
        }
        let p3 = self.c.eval(rec.t1) + self.d * rec.utan1.turn_90();
        let mut p2 = p3;
        let l1 = s * self.q.eval(rec.t1).to_vec2().length() - b * self.d;
        if l1 * rec.cusp1 > 0.0 {
            p2 -= l1 * rec.utan1;
        }
        CubicBez::new(p0, p1, p2, p3)
    }

    /// Compute arc approximation.
    ///
    /// This is called "arc drawing" because if we just look at the delta
    /// vector, it describes an arc from the initial unit normal to the final unit
    /// normal, with "as smooth as possible" parametrization. This approximation
    /// is not necessarily great, but is very robust, and in particular, accuracy
    /// does not degrade for J shaped near-cusp source curves or when the offset
    /// distance is large (with respect to the source curve arc length).
    ///
    /// It is a pretty good approximation overall and has very clean O(n^4) scaling.
    /// Its worst performance is on curves with a large cubic component, where it
    /// undershoots. The theory is that the least squares refinement improves those
    /// cases.
    fn draw_arc(&self, rec: &OffsetRec) -> (f64, f64) {
        // possible optimization: this can probably be done with vectors
        // rather than arctangent
        let th = rec.utan1.cross(rec.utan0).atan2(rec.utan1.dot(rec.utan0));
        let a = (2. / 3.) / (1.0 + (0.5 * th).cos()) * 2.0 * (0.5 * th).sin();
        let b = -a;
        (a, b)
    }

    /// Evaluate error and also refine t values.
    ///
    /// Returns evaluation of error including error vectors and (squared)
    /// maximum error.
    ///
    /// The vector of t values represents points on the source curve; the logic
    /// here is a Newton step to bring those points closer to the normal ray of
    /// the approximation.
    fn eval_err(&self, rec: &OffsetRec, c_approx: CubicBez, ts: &mut [f64; N_LSE]) -> ErrEval {
        let qa = c_approx.deriv();
        let mut err_squared = 0.0;
        let mut unorms = [Vec2::ZERO; N_LSE];
        let mut err_vecs = [Vec2::ZERO; N_LSE];
        for i in 0..N_LSE {
            let ta = (i + 1) as f64 * (1.0 / (N_LSE + 1) as f64);
            let mut t = ts[i];
            let p = self.c.eval(t);
            // Newton step to refine t value
            let pa = c_approx.eval(ta);
            let tana = qa.eval(ta).to_vec2();
            t += tana.dot(pa - p) / tana.dot(self.q.eval(t).to_vec2());
            t = t.max(rec.t0).min(rec.t1);
            ts[i] = t;
            let cusp = rec.cusp0.signum();
            let unorm = cusp * tana.normalize().turn_90();
            unorms[i] = unorm;
            let p_new = self.c.eval(t) + self.d * unorm;
            let err_vec = pa - p_new;
            err_vecs[i] = err_vec;
            let mut dist_err_squared = err_vec.length_squared();
            if !dist_err_squared.is_finite() {
                // A hack to make sure we reject bad refinements
                dist_err_squared = 1e12;
            }
            // Note: consider also incorporating angle error
            err_squared = dist_err_squared.max(err_squared);
        }
        ErrEval {
            err_squared,
            unorms,
            err_vecs,
        }
    }

    /// Refine an approximation, minimizing least squares error.
    ///
    /// Compute the approximation that minimizes least squares error, based on the given error
    /// evaluation.
    ///
    /// The effectiveness of this refinement varies. Basically, if the curve has a large cubic
    /// component, then the arc drawing will undershoot systematically and this refinement will
    /// reduce error considerably. In other cases, it will eventually converge to a local
    /// minimum, but convergence is slow.
    ///
    /// The `BLEND` parameter controls a tradeoff between robustness and speed of convergence.
    /// In the happy case, convergence is fast and not very sensitive to this parameter. If the
    /// parameter is too small, then in near-parabola cases the determinant will be small and
    /// the result not numerically stable.
    ///
    /// A value of 1.0 for `BLEND` corresponds to essentially the Hoschek method, minimizing
    /// Euclidean distance (which tends to over-anchor on the given t values). A value of 0 would
    /// minimize the dot product of error wrt the normal vector, ignoring the cross product
    /// component.
    ///
    /// A possible future direction would be to tune the parameter adaptively.
    fn refine_least_squares(&self, rec: &OffsetRec, a: f64, b: f64, err: &ErrEval) -> (f64, f64) {
        let mut aa = 0.0;
        let mut ab = 0.0;
        let mut ac = 0.0;
        let mut bb = 0.0;
        let mut bc = 0.0;
        for i in 0..N_LSE {
            let n = err.unorms[i];
            let err_vec = err.err_vecs[i];
            let c_n = err_vec.dot(n);
            let c_t = err_vec.cross(n);
            let a_n = A_WEIGHTS[i] * rec.utan0.dot(n);
            let a_t = A_WEIGHTS[i] * rec.utan0.cross(n);
            let b_n = B_WEIGHTS[i] * rec.utan1.dot(n);
            let b_t = B_WEIGHTS[i] * rec.utan1.cross(n);
            aa += a_n * a_n + BLEND * (a_t * a_t);
            ab += a_n * b_n + BLEND * a_t * b_t;
            ac += a_n * c_n + BLEND * a_t * c_t;
            bb += b_n * b_n + BLEND * (b_t * b_t);
            bc += b_n * c_n + BLEND * b_t * c_t;
        }
        let idet = 1.0 / (self.d * (aa * bb - ab * ab));
        let delta_a = idet * (ac * bb - ab * bc);
        let delta_b = idet * (aa * bc - ac * ab);
        (a - delta_a, b - delta_b)
    }

    /// Decide where to subdivide when error is exceeded.
    ///
    /// For curves not containing an inflection point, subdivide at the tangent
    /// bisecting the endpoint tangents. The logic is that for a near cusp in the
    /// source curve, you want the subdivided sections to be approximately
    /// circular arcs with progressively smaller angles.
    ///
    /// When there is an inflection point (or, more specifically, when the curve
    /// crosses its chord), bisecting the angle can lead to very lopsided arc
    /// lengths, so just subdivide by t in that case.
    fn find_subdivision_point(&self, rec: &OffsetRec) -> SubdivisionPoint {
        let t = 0.5 * (rec.t0 + rec.t1);
        let q_t = self.q.eval(t).to_vec2();
        let x0 = rec.utan0.cross(q_t).abs();
        let x1 = rec.utan1.cross(q_t).abs();
        const SUBDIVIDE_THRESH: f64 = 0.1;
        if x0 > SUBDIVIDE_THRESH * x1 && x1 > SUBDIVIDE_THRESH * x0 {
            let utan = q_t.normalize();
            return SubdivisionPoint { t, utan };
        }

        // Note: do we want to track p0 & p3 in rec, to avoid repeated eval?
        let chord = self.c.eval(rec.t1) - self.c.eval(rec.t0);
        if chord.cross(rec.utan0) * chord.cross(rec.utan1) < 0.0 {
            let tan = rec.utan0 + rec.utan1;
            if let Some(subdivision) =
                self.subdivide_for_tangent(rec.utan0, rec.t0, rec.t1, tan, false)
            {
                return subdivision;
            }
        }
        // Curve definitely has an inflection point
        // Try to subdivide based on integral of absolute curvature.

        // Tangents at recursion endpoints and inflection points.
        let mut tangents: ArrayVec<Vec2, 4> = ArrayVec::new();
        let mut ts: ArrayVec<f64, 4> = ArrayVec::new();
        tangents.push(rec.utan0);
        ts.push(rec.t0);
        for t in self.c.inflections() {
            if t > rec.t0 && t < rec.t1 {
                tangents.push(self.q.eval(t).to_vec2());
                ts.push(t);
            }
        }
        tangents.push(rec.utan1);
        ts.push(rec.t1);
        let mut arc_angles: ArrayVec<f64, 3> = ArrayVec::new();
        let mut sum = 0.0;
        for i in 0..tangents.len() - 1 {
            let tan0 = tangents[i];
            let tan1 = tangents[i + 1];
            let th = tan0.cross(tan1).atan2(tan0.dot(tan1));
            sum += th.abs();
            arc_angles.push(th);
        }
        let mut target = sum * 0.5;
        let mut i = 0;
        while arc_angles[i].abs() < target {
            target -= arc_angles[i].abs();
            i += 1;
        }
        let rotation = Vec2::from_angle(target.copysign(arc_angles[i]));
        let base = tangents[i];
        let tan = base.rotate_scale(rotation);
        let utan0 = if i == 0 { rec.utan0 } else { base.normalize() };
        self.subdivide_for_tangent(utan0, ts[i], ts[i + 1], tan, true)
            .unwrap()
    }

    /// Find a subdivision point, given a tangent vector.
    ///
    /// When subdividing by bisecting the angle (or, more generally, subdividing by
    /// the L1 norm of curvature when there are inflection points), we find the
    /// subdivision point by solving for the tangent matching, specifically the
    /// cross-product of the tangent and the curve's derivative being zero. For
    /// internal cusps, subdividing near the cusp is a good thing, but there is
    /// still a robustness concern for vanishing derivative at the endpoints.
    fn subdivide_for_tangent(
        &self,
        utan0: Vec2,
        t0: f64,
        t1: f64,
        tan: Vec2,
        force: bool,
    ) -> Option<SubdivisionPoint> {
        let mut t = 0.0;
        let mut n_soln = 0;
        // set up quadratic equation for matching tangents
        let z0 = tan.cross(self.q.p0.to_vec2());
        let z1 = tan.cross(self.q.p1.to_vec2());
        let z2 = tan.cross(self.q.p2.to_vec2());
        let c0 = z0;
        let c1 = 2.0 * (z1 - z0);
        let c2 = (z2 - z1) - (z1 - z0);
        for root in solve_quadratic(c0, c1, c2) {
            if root >= t0 && root <= t1 {
                t = root;
                n_soln += 1;
            }
        }
        if n_soln != 1 {
            if !force {
                return None;
            }
            // Numerical failure, try to subdivide at cusp; we pick the
            // smaller derivative.
            if self.q.eval(t0).to_vec2().length_squared()
                > self.q.eval(t1).to_vec2().length_squared()
            {
                t = t1;
            } else {
                t = t0;
            }
        }
        let q = self.q.eval(t).to_vec2();
        const UTAN_EPSILON: f64 = 1e-12;
        let utan = if n_soln == 1 && q.length_squared() >= UTAN_EPSILON {
            q.normalize()
        } else if tan.length_squared() >= UTAN_EPSILON {
            // Curve has a zero-derivative cusp but angles well defined
            tan.normalize()
        } else {
            // 180 degree U-turn, arbitrarily pick a direction.
            // If we get to this point, there will probably be a failure.
            utan0.turn_90()
        };
        Some(SubdivisionPoint { t, utan })
    }
}

impl OffsetRec {
    #[allow(clippy::too_many_arguments)]
    fn new(
        t0: f64,
        t1: f64,
        utan0: Vec2,
        utan1: Vec2,
        cusp0: f64,
        cusp1: f64,
        depth: usize,
    ) -> Self {
        OffsetRec {
            t0,
            t1,
            utan0,
            utan1,
            cusp0,
            cusp1,
            depth,
        }
    }
}

/// Compute Bézier weights for evenly subdivided t values.
const fn mk_a_weights(rev: bool) -> [f64; N_LSE] {
    let mut result = [0.0; N_LSE];
    let mut i = 0;
    while i < N_LSE {
        let t = (i + 1) as f64 / (N_LSE + 1) as f64;
        let mt = 1. - t;
        let ix = if rev { N_LSE - 1 - i } else { i };
        result[ix] = 3.0 * mt * t * mt;
        i += 1;
    }
    result
}

const A_WEIGHTS: [f64; N_LSE] = mk_a_weights(false);
const B_WEIGHTS: [f64; N_LSE] = mk_a_weights(true);

#[cfg(test)]
mod tests {
    use super::offset_cubic;
    use crate::{BezPath, CubicBez, PathEl};

    // This test tries combinations of parameters that have caused problems in the past.
    #[test]
    fn pathological_curves() {
        let curve = CubicBez {
            p0: (-1236.3746269978635, 152.17981429574826).into(),
            p1: (-1175.18662093517, 108.04721798590596).into(),
            p2: (-1152.142883879584, 105.76260301083356).into(),
            p3: (-1151.842639804639, 105.73040758939104).into(),
        };
        let offset = 3603.7267536453924;
        let accuracy = 0.1;

        let mut result = BezPath::new();
        offset_cubic(curve, offset, accuracy, &mut result);
        assert!(matches!(result.iter().next(), Some(PathEl::MoveTo(_))));

        let mut result = BezPath::new();
        offset_cubic(curve, offset, accuracy, &mut result);
        assert!(matches!(result.iter().next(), Some(PathEl::MoveTo(_))));
    }

    /// Cubic offset that used to trigger infinite recursion.
    #[test]
    fn infinite_recursion() {
        const TOLERANCE: f64 = 0.1;
        const OFFSET: f64 = -0.5;
        let c = CubicBez::new(
            (1096.2962962962963, 593.90243902439033),
            (1043.6213991769548, 593.90243902439033),
            (1030.4526748971193, 593.90243902439033),
            (1056.7901234567901, 593.90243902439033),
        );

        let mut result = BezPath::new();
        offset_cubic(c, OFFSET, TOLERANCE, &mut result);
    }

    #[test]
    fn test_cubic_offset_simple_line() {
        let cubic = CubicBez::new((0., 0.), (10., 0.), (20., 0.), (30., 0.));
        let mut result = BezPath::new();
        offset_cubic(cubic, 5., 1e-6, &mut result);
    }
}