ogeom_geom/

fit.rs

//! Fitting a B-spline to points, to a stated error target.
//!
//! Interpolation and fixed-count approximation live in `ogeom-algo`, and both
//! take a control-point count from the caller. This module chooses it: fit,
//! measure where the fit is worst, refine the knots exactly there, and repeat
//! until the error target is met.
//!
//! A fit that silently picks its own resolution and reports success is the shape of answer
//! that gets trusted. So the result here carries the error actually reached and
//! whether the target was met, and a fit that ran out of room says so rather
//! than rounding "close" up to "done".
//!
//! # Where the knots go
//!
//! Refinement is *where the error is*, not everywhere. Splitting every span
//! doubles the control points per round and most of them buy nothing. A curve
//! that is straight for most of its length and tight in one corner needs its
//! knots in the corner. Each round measures the error per span and splits only
//! the spans that exceed the target, so the knot density ends up tracking the
//! curvature, which is where it belongs.
//!
//! # What the error measures
//!
//! A fit through bare points ([`fit_points`], [`fit_points_at`],
//! [`fit_surface_grid`], [`fit_surface_scattered`] and their siblings) has
//! nothing to measure against but the points, and its error is the distance
//! at the points only. Between them the fit runs as close to the sampled
//! shape as a smooth spline through those samples does, which for smooth
//! data shrinks with the sample spacing, but no number here bounds it. A
//! caller holding the true geometry measures there, or fits with
//! [`fit_curve_sampled`] or [`fit_surface_sampled`], whose error covers the
//! points between the samples too and is no less than the widest miss.
//!
//! # Who this is for
//!
//! The marching intersector, first: a traced branch is a polyline with a stated
//! chord tolerance, and downstream code wants a curve, so the polyline is fitted
//! to the same tolerance and the result is as good as the trace. But nothing
//! here knows about intersections. It fits points, in three dimensions or two,
//! which is also what a digitized profile or an imported polyline needs.

use ogeom_core::{OgeomResult, Tolerances, ogeom_bail};
use ogeom_math::{KnotVector, Point, Point2};

use crate::curve::BSplineCurve;
use crate::curve2d::BSpline2d;

mod seeded;

pub use seeded::{Ended, ScatteredFit, ScatteredOptions, Target, fit_surface_scattered_over};

#[allow(
    clippy::cast_precision_loss,
    reason = "sample counts are far below 2^52"
)]
fn precise(n: usize) -> f64 {
    n as f64
}

/// What a fit produced.
#[derive(Debug, Clone, PartialEq)]
pub struct Fitted<C> {
    /// The curve.
    pub curve: C,
    /// The largest distance from any input point to the curve at its
    /// parameter. A fit through bare points measures at those points only;
    /// a sampled fit measures between its samples as well (see the module
    /// documentation).
    pub error: f64,
    /// Whether the error target was met.
    ///
    /// A fit can run out of room (more control points than points to fit
    /// solves nothing, since at that ratio least squares *is* interpolation),
    /// and then this is `false` and `error` says how close it got. Reported
    /// rather than rounded up to success.
    pub met: bool,
}

/// Fit a spline through 3D points, refining until `tolerance` is met.
///
/// The first and last points are honoured exactly: they are where the curve
/// joins whatever comes next, and a fit that drifts at its ends produces gaps
/// at every junction built on it. A closed polyline (first point repeated at
/// the end) therefore comes back closed.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if there are
/// fewer than two distinct points or the tolerance is not a positive distance.
pub fn fit_points(
    points: &[Point],
    degree: usize,
    tolerance: f64,
    tol: Tolerances,
) -> OgeomResult<Fitted<BSplineCurve>> {
    let (knots, control, error, met) = fit::<3>(
        &points.iter().map(|p| [p.x, p.y, p.z]).collect::<Vec<_>>(),
        degree,
        tolerance,
        false,
        tol,
    )?;
    let curve = BSplineCurve::new(
        knots,
        control
            .into_iter()
            .map(|c| Point::new(c[0], c[1], c[2]))
            .collect(),
        tol,
    )?;
    Ok(Fitted { curve, error, met })
}

/// Fit a *fair* curve: least squares over the points, pulled toward
/// minimum bending energy by a smoothing weight.
///
/// The energy is the squared second difference of the control polygon (the
/// discrete bending of the curve's own skeleton), added to the normal
/// equations as `λ·DᵀD`. At `λ = 0` this is plain least squares. As `λ`
/// grows the curve trades closeness for straightness, which is the batten
/// a drafter flexes through points. Both ends interpolate their points
/// exactly, whatever the weight. The reported error is the honest maximum
/// distance from the inputs, which *rises* with `λ`: fairness is spent
/// closeness, and the number says how much was spent.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if
/// fewer than two distinct points arrive, the control budget cannot carry
/// the degree, or the weight is not finite and non-negative.
pub fn fit_points_faired(
    points: &[Point],
    degree: usize,
    controls: usize,
    smoothing: f64,
    tol: Tolerances,
) -> OgeomResult<Fitted<BSplineCurve>> {
    if !smoothing.is_finite() || smoothing < 0.0 {
        ogeom_bail!(
            Construction,
            "a smoothing weight of {smoothing} is not a weight"
        );
    }
    if controls < degree + 1 {
        ogeom_bail!(
            Construction,
            "{controls} control points cannot carry degree {degree}"
        );
    }
    let m = points.len();
    if m < 2 {
        ogeom_bail!(Construction, "a fair curve needs at least two points");
    }

    // Chord-length parameters over the clamped knot domain.
    let knots = KnotVector::clamped_uniform(degree, controls)?;
    let (lo, hi) = (knots.domain_start(), knots.domain_end());
    let mut cumulative = vec![0.0f64; m];
    for i in 1..m {
        cumulative[i] = cumulative[i - 1] + points[i].distance(points[i - 1]);
    }
    let total = cumulative[m - 1];
    if total <= tol.confusion() {
        ogeom_bail!(
            Construction,
            "the points coincide; there is no curve to fair"
        );
    }
    let parameters: Vec<f64> = cumulative
        .iter()
        .map(|c| lo + (hi - lo) * (c / total))
        .collect();

    // The full stiffness K = AᵀA + λ·DᵀD over every control point, then the
    // end points pinned by moving their columns to the right-hand side.
    let n = controls;
    let mut a = nalgebra::DMatrix::<f64>::zeros(m, n);
    for (k, &u) in parameters.iter().enumerate() {
        let span = knots.span_unchecked(u);
        let basis = knots.basis(span, u);
        let first = span - degree;
        for (j, b) in basis.iter().enumerate() {
            a[(k, first + j)] = *b;
        }
    }
    let mut d = nalgebra::DMatrix::<f64>::zeros(n.saturating_sub(2), n);
    for r in 0..n.saturating_sub(2) {
        d[(r, r)] = 1.0;
        d[(r, r + 1)] = -2.0;
        d[(r, r + 2)] = 1.0;
    }
    let k_full = a.transpose() * &a + d.transpose() * &d * smoothing;

    let free: Vec<usize> = (1..n - 1).collect();
    let mut reduced = nalgebra::DMatrix::<f64>::zeros(free.len(), free.len());
    for (ri, &i) in free.iter().enumerate() {
        for (rj, &j) in free.iter().enumerate() {
            reduced[(ri, rj)] = k_full[(i, j)];
        }
    }
    let decomposition = reduced.lu();

    let mut control = vec![Point::ORIGIN; n];
    control[0] = points[0];
    control[n - 1] = points[m - 1];
    for axis in 0..3 {
        let b_data =
            nalgebra::DVector::from_iterator(m, points.iter().map(|p| [p.x, p.y, p.z][axis]));
        let full_rhs = a.transpose() * &b_data;
        let mut rhs = nalgebra::DVector::<f64>::zeros(free.len());
        for (ri, &i) in free.iter().enumerate() {
            rhs[ri] = full_rhs[i]
                - k_full[(i, 0)] * [points[0].x, points[0].y, points[0].z][axis]
                - k_full[(i, n - 1)] * [points[m - 1].x, points[m - 1].y, points[m - 1].z][axis];
        }
        let Some(solved) = decomposition.solve(&rhs) else {
            ogeom_bail!(
                Numeric,
                "the faired system is singular; fewer controls or more points"
            );
        };
        for (ri, &i) in free.iter().enumerate() {
            match axis {
                0 => control[i].x = solved[ri],
                1 => control[i].y = solved[ri],
                _ => control[i].z = solved[ri],
            }
        }
    }

    let curve = BSplineCurve::new(knots, control, tol)?;
    let mut error = 0.0f64;
    {
        use crate::traits::Curve3d as _;
        for (point, &u) in points.iter().zip(&parameters) {
            error = error.max(curve.point_at(u, tol)?.distance(*point));
        }
    }
    Ok(Fitted {
        curve,
        error,
        met: true,
    })
}

/// Fit a smoothly closed loop: as [`fit_points`], with the join C1.
///
/// The input must be a loop (the first point repeated at the end), and the
/// ends are honoured exactly as always. Beyond that, the tangent leaving the
/// join is constrained to equal the tangent arriving at it, eliminated
/// exactly inside the least-squares solve rather than patched on after, so a
/// surface built over the loop shows no crease at the seam.
///
/// The constraint spends shape freedom at the join: the last movable control
/// point follows the first, so a loop fits with slightly more knots than the
/// same data fitted open, and the tightest reachable error is a little
/// higher. That is the price of the join being smooth, and it is why this is
/// a separate entry rather than a change to [`fit_points`].
///
/// # Errors
///
/// As [`fit_points`], and additionally if the input is not a closed loop.
pub fn fit_points_closed(
    points: &[Point],
    degree: usize,
    tolerance: f64,
    tol: Tolerances,
) -> OgeomResult<Fitted<BSplineCurve>> {
    let Some((first, last)) = points.first().zip(points.last()) else {
        ogeom_bail!(Construction, "a closed fit needs points");
    };
    if !first.is_equal(*last, tol) {
        ogeom_bail!(
            Construction,
            "a closed fit needs a loop: the first point repeated at the end"
        );
    }
    let (knots, control, error, met) = fit::<3>(
        &points.iter().map(|p| [p.x, p.y, p.z]).collect::<Vec<_>>(),
        degree,
        tolerance,
        true,
        tol,
    )?;
    let curve = BSplineCurve::new(
        knots,
        control
            .into_iter()
            .map(|c| Point::new(c[0], c[1], c[2]))
            .collect(),
        tol,
    )?;
    Ok(Fitted { curve, error, met })
}

/// Fit a spline through 2D points, refining until `tolerance` is met.
///
/// The planar twin of [`fit_points`], for pcurves: an intersection curve lives
/// on both surfaces, and each face needs it in its own parameter space.
///
/// # Errors
///
/// As [`fit_points`].
pub fn fit_points_2d(
    points: &[Point2],
    degree: usize,
    tolerance: f64,
    tol: Tolerances,
) -> OgeomResult<Fitted<BSpline2d>> {
    let (knots, control, error, met) = fit::<2>(
        &points.iter().map(|p| [p.x, p.y]).collect::<Vec<_>>(),
        degree,
        tolerance,
        false,
        tol,
    )?;
    let curve = BSpline2d::new(
        knots,
        control
            .into_iter()
            .map(|c| Point2::new(c[0], c[1]))
            .collect(),
        tol,
    )?;
    Ok(Fitted { curve, error, met })
}

/// Fit one curve living in three spaces at once: a 3D curve and its two
/// parameter-space images, as a single seven-dimensional fit.
///
/// One parameterization, one knot vector, one correction: the three results
/// are same-parameter *by construction*, which separate fits cannot promise:
/// each fit's parameter correction drifts its parameterization independently,
/// and the drift is invisible to every per-fit error measure, so a pcurve
/// fitted apart can report a tiny error and still evaluate far from its own
/// curve. The reported error bounds the worst deviation across all seven
/// coordinates, so it bounds each space's deviation too.
///
/// # Errors
///
/// As [`fit_points`], and the three inputs must be equally long.
#[allow(clippy::type_complexity)]
pub fn fit_points_joint(
    points: &[Point],
    on_a: &[Point2],
    on_b: &[Point2],
    degree: usize,
    tolerance: f64,
    tol: Tolerances,
) -> OgeomResult<(Fitted<BSplineCurve>, BSpline2d, BSpline2d)> {
    fit_points_joint_inner(points, on_a, on_b, degree, tolerance, None, tol)
}

/// As [`fit_points_joint`], with a closed loop's join made C1.
///
/// The seven-dimensional twin of [`fit_points_closed`]: when the trace is a
/// loop (the first sample repeated at the end in every space), the join's
/// tangent constraint is eliminated inside the shared solve, so the curve
/// *and both pcurves* cross their seam smoothly, still same-parameter by
/// construction. An input that is not closed in all seven coordinates fits
/// as [`fit_points_joint`] would. The constraint simply never engages.
///
/// # Errors
///
/// As [`fit_points_joint`].
#[allow(clippy::type_complexity)]
pub fn fit_points_joint_closed(
    points: &[Point],
    on_a: &[Point2],
    on_b: &[Point2],
    degree: usize,
    tolerance: f64,
    tol: Tolerances,
) -> OgeomResult<(Fitted<BSplineCurve>, BSpline2d, BSpline2d)> {
    fit_points_joint_inner(points, on_a, on_b, degree, tolerance, Some(7), tol)
}

/// As [`fit_points_joint_closed`], for a loop closed in space whose pcurves
/// may end a whole period from where they began: a section winding once
/// round a periodic surface. The join is C1 in every coordinate, which is
/// right where the chart runs at one speed across its seam, as a periodic
/// surface's does. A patch that merely meets itself need not.
///
/// # Errors
///
/// As [`fit_points_joint`].
#[allow(clippy::type_complexity)]
pub fn fit_points_joint_winding(
    points: &[Point],
    on_a: &[Point2],
    on_b: &[Point2],
    degree: usize,
    tolerance: f64,
    tol: Tolerances,
) -> OgeomResult<(Fitted<BSplineCurve>, BSpline2d, BSpline2d)> {
    fit_points_joint_inner(points, on_a, on_b, degree, tolerance, Some(3), tol)
}

/// As [`fit_points_joint`], at chord-length parameters held fixed, its
/// knots starting at the parameters of the samples `stations` names.
///
/// A trace sampled far more finely than a cubic needs (a marched section,
/// stepped by its chord's sag) fits in one round where the caller knows
/// where the spans go: the knots start there and are refined only where a
/// sample misses. Every sample is fitted and measured. The error is the
/// widest seven-coordinate distance of a sample from the fit at its own
/// parameter, which bounds its distance from the fit anywhere. Where
/// `closed`, the ends are pinned and the join is C1 in every coordinate,
/// as [`fit_points_joint_winding`] makes it.
///
/// # Errors
///
/// As [`fit_points_joint`]; also where two consecutive samples coincide in
/// all seven coordinates, which leaves no parameter between them.
#[allow(
    clippy::type_complexity,
    clippy::too_many_arguments,
    reason = "one trace in three spaces, its stations and the fit's terms"
)]
pub fn fit_points_joint_from(
    points: &[Point],
    on_a: &[Point2],
    on_b: &[Point2],
    stations: &[usize],
    closed: bool,
    degree: usize,
    tolerance: f64,
    tol: Tolerances,
) -> OgeomResult<(Fitted<BSplineCurve>, BSpline2d, BSpline2d)> {
    joint_from(
        points, on_a, on_b, stations, closed, degree, tolerance, false, tol,
    )
}

/// As [`fit_points_joint_from`], splitting a span only where samples lie
/// strictly inside it on both sides of the new knot.
///
/// The refinement splits a span between its two middle samples. Where a
/// span holds two samples, one on its own start knot, that leaves a span
/// whose only sample sits on its knot: nothing holds the cubic inside it,
/// and the fit can pass every sample and swing between two of them by far
/// more than their spacing, a turn round a drum in its chart image between
/// samples a fraction of a millimetre apart. Here such a span is left
/// whole and the fit stops there, its residual stated as it is.
///
/// # Errors
///
/// As [`fit_points_joint_from`].
#[allow(
    clippy::type_complexity,
    clippy::too_many_arguments,
    reason = "one trace in three spaces, its stations and the fit's terms"
)]
pub fn fit_points_joint_from_held(
    points: &[Point],
    on_a: &[Point2],
    on_b: &[Point2],
    stations: &[usize],
    closed: bool,
    degree: usize,
    tolerance: f64,
    tol: Tolerances,
) -> OgeomResult<(Fitted<BSplineCurve>, BSpline2d, BSpline2d)> {
    joint_from(
        points, on_a, on_b, stations, closed, degree, tolerance, true, tol,
    )
}

#[allow(
    clippy::type_complexity,
    clippy::too_many_arguments,
    reason = "one trace in three spaces, its stations and the fit's terms"
)]
fn joint_from(
    points: &[Point],
    on_a: &[Point2],
    on_b: &[Point2],
    stations: &[usize],
    closed: bool,
    degree: usize,
    tolerance: f64,
    within: bool,
    tol: Tolerances,
) -> OgeomResult<(Fitted<BSplineCurve>, BSpline2d, BSpline2d)> {
    if points.len() != on_a.len() || points.len() != on_b.len() {
        ogeom_bail!(
            Construction,
            "a joint fit needs the same trace seen in every space"
        );
    }
    if !tolerance.is_finite() || tolerance <= 0.0 {
        ogeom_bail!(Construction, "a tolerance of {tolerance} is not a distance");
    }
    if degree == 0 || points.len() < 2 {
        ogeom_bail!(
            Construction,
            "a fit needs a degree of at least one and two points"
        );
    }
    let joined: Vec<[f64; 7]> = points
        .iter()
        .zip(on_a)
        .zip(on_b)
        .map(|((p, a), b)| [p.x, p.y, p.z, a.x, a.y, b.x, b.y])
        .collect();
    let parameters = chordal::<7>(&joined);
    if parameters.windows(2).any(|w| w[1] <= w[0]) {
        ogeom_bail!(
            Construction,
            "two consecutive samples coincide and leave no parameter between them"
        );
    }
    let kept: Vec<(f64, usize)> = stations
        .iter()
        .filter_map(|&i| parameters.get(i))
        .filter(|&&u| u > 0.0 && u < 1.0)
        .map(|&u| (u, 1))
        .collect();
    let (knots, mut controls) = fit_family::<7>(
        std::slice::from_ref(&joined),
        &parameters,
        degree,
        tolerance,
        closed,
        false,
        &kept,
        within,
    )?;
    let Some(control) = controls.pop() else {
        ogeom_bail!(NotDone, "the joint fit solved nothing");
    };
    let error = residuals::<7>(&knots, &control, &joined, &parameters)
        .iter()
        .fold(0.0_f64, |acc, e| acc.max(e.1));
    let curve = BSplineCurve::new(
        knots.clone(),
        control
            .iter()
            .map(|c| Point::new(c[0], c[1], c[2]))
            .collect(),
        tol,
    )?;
    let pa = BSpline2d::new(
        knots.clone(),
        control.iter().map(|c| Point2::new(c[3], c[4])).collect(),
        tol,
    )?;
    let pb = BSpline2d::new(
        knots,
        control.iter().map(|c| Point2::new(c[5], c[6])).collect(),
        tol,
    )?;
    let met = error <= tolerance;
    Ok((Fitted { curve, error, met }, pa, pb))
}

#[allow(clippy::type_complexity)]
fn fit_points_joint_inner(
    points: &[Point],
    on_a: &[Point2],
    on_b: &[Point2],
    degree: usize,
    tolerance: f64,
    smooth_loop: Option<usize>,
    tol: Tolerances,
) -> OgeomResult<(Fitted<BSplineCurve>, BSpline2d, BSpline2d)> {
    if points.len() != on_a.len() || points.len() != on_b.len() {
        ogeom_bail!(
            Construction,
            "a joint fit needs the same trace seen in every space"
        );
    }
    let joined: Vec<[f64; 7]> = points
        .iter()
        .zip(on_a)
        .zip(on_b)
        .map(|((p, a), b)| [p.x, p.y, p.z, a.x, a.y, b.x, b.y])
        .collect();
    // Centripetal first, as every free fit is, and by chord length where
    // that misses its target: the two disagree only where the samples are
    // spaced far from evenly, and there each is right for its own data
    // (the centripetal guess for a trace with a kink in it, chord length
    // for a smooth trace crowded at one end), so the closer of the two
    // stands.
    let first = fit_spaced::<7>(
        &joined,
        degree,
        tolerance,
        smooth_loop,
        Spacing::Centripetal,
        tol,
    )?;
    let (knots, control, error, met) = if first.3 {
        first
    } else {
        match fit_spaced::<7>(
            &joined,
            degree,
            tolerance,
            smooth_loop,
            Spacing::ChordLength,
            tol,
        ) {
            Ok(second) if second.2 < first.2 => second,
            _ => first,
        }
    };
    let curve = BSplineCurve::new(
        knots.clone(),
        control
            .iter()
            .map(|c| Point::new(c[0], c[1], c[2]))
            .collect(),
        tol,
    )?;
    let pa = BSpline2d::new(
        knots.clone(),
        control.iter().map(|c| Point2::new(c[3], c[4])).collect(),
        tol,
    )?;
    let pb = BSpline2d::new(
        knots,
        control.iter().map(|c| Point2::new(c[5], c[6])).collect(),
        tol,
    )?;
    Ok((Fitted { curve, error, met }, pa, pb))
}

/// Fit a pcurve at *fixed* parameters: the source curve's own.
///
/// The fixed-parameter twin of [`fit_points_2d`], and the difference is the
/// contract: parameter correction is what makes a free fit's residual honest,
/// and it is exactly what a *same-parameter* fit must never do, because the
/// parameters are not a guess to be improved. They are the 3D curve's own,
/// and drifting them is how a pcurve ends up evaluating away from the curve
/// it annotates. Here the parameters stay put, refinement adds knots where
/// the error says, and the reported error is the true same-parameter
/// deviation in the chart at the given points. Between them it is not
/// measured; [`fit_curve_2d_sampled`] measures there.
///
/// # Errors
///
/// As [`fit_points_2d`], and the parameters must be strictly increasing and
/// as many as the points.
pub fn fit_points_2d_at(
    parameters: &[f64],
    points: &[Point2],
    degree: usize,
    tolerance: f64,
    tol: Tolerances,
) -> OgeomResult<Fitted<BSpline2d>> {
    fit_points_2d_at_inner(parameters, points, degree, tolerance, false, tol)
}

/// Fit space points at *fixed* parameters: the caller's `t` values are the
/// curve's own, which is what keeps a replacement curve same-parameter with
/// every chart already speaking the old one. The error is measured at the
/// given points only; [`fit_curve_sampled`] measures between them too.
///
/// # Errors
///
/// As [`fit_points`], and the parameters must be strictly increasing and as
/// many as the points.
pub fn fit_points_at(
    parameters: &[f64],
    points: &[Point],
    degree: usize,
    tolerance: f64,
    tol: Tolerances,
) -> OgeomResult<Fitted<BSplineCurve>> {
    fit_points_at_inner(parameters, points, degree, tolerance, false, tol)
}

/// As [`fit_points_at`], with the loop's join made C1: the points close,
/// the last the first again.
///
/// # Errors
///
/// As [`fit_points_at`].
pub fn fit_points_at_closed(
    parameters: &[f64],
    points: &[Point],
    degree: usize,
    tolerance: f64,
    tol: Tolerances,
) -> OgeomResult<Fitted<BSplineCurve>> {
    fit_points_at_inner(parameters, points, degree, tolerance, true, tol)
}

fn fit_points_at_inner(
    parameters: &[f64],
    points: &[Point],
    degree: usize,
    tolerance: f64,
    closed: bool,
    tol: Tolerances,
) -> OgeomResult<Fitted<BSplineCurve>> {
    if parameters.len() != points.len() {
        ogeom_bail!(Construction, "one parameter per point, or the fit is a lie");
    }
    if parameters.windows(2).any(|w| w[1] <= w[0]) {
        ogeom_bail!(Construction, "fixed parameters must strictly increase");
    }
    if !tolerance.is_finite() || tolerance <= 0.0 {
        ogeom_bail!(Construction, "a tolerance of {tolerance} is not a distance");
    }
    if degree == 0 {
        ogeom_bail!(Construction, "a fit needs a degree of at least one");
    }
    let data: Vec<[f64; 3]> = points.iter().map(|p| [p.x, p.y, p.z]).collect();
    if data.len() < 2 {
        ogeom_bail!(Construction, "a fit needs at least two points");
    }
    let degree = degree.min(data.len() - 1);
    let (a, b) = (parameters[0], parameters[parameters.len() - 1]);
    let mut knots = single_span(degree, a, b)?;

    const ROUNDS: usize = 32;
    let mut best: Option<(KnotVector, Vec<[f64; 3]>, f64)> = None;
    for _ in 0..ROUNDS {
        let control = match least_squares::<3>(&knots, &data, parameters, closed) {
            Ok(control) => control,
            Err(e) => {
                if best.is_some() {
                    break;
                }
                return Err(e);
            }
        };
        let errors = residuals::<3>(&knots, &control, &data, parameters);
        let worst = errors.iter().fold(0.0_f64, |acc, e| acc.max(e.1));
        if best.as_ref().is_none_or(|(_, _, held)| worst < *held) {
            best = Some((knots.clone(), control.clone(), worst));
        }
        if worst <= tolerance {
            let curve = build_curve_3(knots, control, tol)?;
            return Ok(Fitted {
                curve,
                error: worst,
                met: true,
            });
        }
        if knots.control_point_count() >= data.len() {
            break;
        }
        let Some(refined) = refined_where_bad(&knots, &errors, tolerance, false)? else {
            break;
        };
        knots = refined_or_interpolating(refined, parameters, closed)?;
    }
    let Some((knots, control, worst)) = best else {
        ogeom_bail!(Construction, "the fit found no usable rounds");
    };
    let curve = build_curve_3(knots, control, tol)?;
    Ok(Fitted {
        curve,
        error: worst,
        met: false,
    })
}

/// The most spans an open sampled fit refines its samples to
/// ([`fit_curve_sampled`], [`fit_curve_2d_sampled`]).
pub const SAMPLED_SPANS: usize = 8192;

/// The most spans a closed sampled fit refines its samples to. A loop's
/// join couples its first and last control points, so its least-squares
/// system is solved whole rather than as a band, at a cost that grows
/// with the cube of the samples.
pub const SAMPLED_LOOP_SPANS: usize = 1024;

/// Fit a curve to the geometry `point` traces over its parameter, measured
/// between the samples as well as at them: the curve counterpart of
/// [`fit_surface_sampled`].
///
/// The fit is same-parameter with `point`, fitted at `ts` to start with:
/// the least-squares fit first, and where it misses, the spline through
/// every sample. Every span's eighth points are compared with `point`
/// there, and the widest misses climbed to the tops of their peaks between
/// them; a miss splits the span at its middle, and this repeats until
/// `tolerance` holds at every point measured or the samples would pass
/// [`SAMPLED_SPANS`] spans. A fit the eighth points pass is read again at
/// eight points of every knot span of its own, its peaks climbed, before
/// it is stated met: a peak at a corner of the trace can stand between
/// the eighth points. Where `closed`, the last of `ts` is the first
/// again, the fit crosses its join C1 and refines to at most
/// [`SAMPLED_LOOP_SPANS`] spans. The best fit is returned either
/// way, its error the worst distance measured.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the
/// tolerance is not a distance, the degree is zero or the parameters do not
/// strictly increase with at least two; whatever `point` refuses.
pub fn fit_curve_sampled(
    point: impl FnMut(f64) -> OgeomResult<Point>,
    ts: &[f64],
    closed: bool,
    degree: usize,
    tolerance: f64,
    tol: Tolerances,
) -> OgeomResult<Fitted<BSplineCurve>> {
    let (knots, control, error, met) =
        sampled::<3, Point>(point, |p| [p.x, p.y, p.z], ts, closed, degree, tolerance)?;
    Ok(Fitted {
        curve: build_curve_3(knots, control, tol)?,
        error,
        met,
    })
}

/// As [`fit_curve_sampled`], in the plane: a chart image fitted at the
/// parameters of the curve it annotates. The distances are in the trace's
/// own units.
///
/// # Errors
///
/// As [`fit_curve_sampled`].
pub fn fit_curve_2d_sampled(
    point: impl FnMut(f64) -> OgeomResult<Point2>,
    ts: &[f64],
    closed: bool,
    degree: usize,
    tolerance: f64,
    tol: Tolerances,
) -> OgeomResult<Fitted<BSpline2d>> {
    let (knots, control, error, met) =
        sampled::<2, Point2>(point, |p| [p.x, p.y], ts, closed, degree, tolerance)?;
    let curve = BSpline2d::new(
        knots,
        control
            .into_iter()
            .map(|c| Point2::new(c[0], c[1]))
            .collect(),
        tol,
    )?;
    Ok(Fitted { curve, error, met })
}

/// The best sampled fit: knots, control points, the worst distance
/// measured, and whether that is within the tolerance.
type SampledFit<const D: usize> = (KnotVector, Vec<[f64; D]>, f64, bool);

/// The refinement both sampled fits share, in `D` coordinates.
fn sampled<const D: usize, P>(
    mut trace: impl FnMut(f64) -> OgeomResult<P>,
    coords: impl Fn(&P) -> [f64; D],
    ts: &[f64],
    closed: bool,
    degree: usize,
    tolerance: f64,
) -> OgeomResult<SampledFit<D>> {
    if !(tolerance > 0.0 && tolerance.is_finite()) {
        ogeom_bail!(Construction, "a tolerance of {tolerance} is not a distance");
    }
    if degree == 0 {
        ogeom_bail!(Construction, "a fit needs a degree of at least one");
    }
    if ts.len() < 2 || ts.iter().any(|t| !t.is_finite()) || ts.windows(2).any(|w| w[1] <= w[0]) {
        ogeom_bail!(
            Construction,
            "a sampled fit needs strictly increasing parameters, at least two"
        );
    }
    let most = if closed {
        SAMPLED_LOOP_SPANS
    } else {
        SAMPLED_SPANS
    };
    let mut ts = ts.to_vec();
    // Every point evaluated, by parameter: a split span's new sample and
    // most of its new eighths were measured before.
    let mut seen: ogeom_core::FastMap<u64, [f64; D]> = ogeom_core::FastMap::default();
    let mut at = |t: f64| -> OgeomResult<[f64; D]> {
        if let Some(p) = seen.get(&t.to_bits()) {
            return Ok(*p);
        }
        let p = coords(&trace(t)?);
        seen.insert(t.to_bits(), p);
        Ok(p)
    };
    // The best fit so far, its error, and whether it was read densely.
    let mut best: Option<(KnotVector, Vec<[f64; D]>, f64, bool)> = None;
    loop {
        let mut data = ts.iter().map(|t| at(*t)).collect::<OgeomResult<Vec<_>>>()?;
        if closed {
            data[ts.len() - 1] = data[0];
        }
        // As the surface fit: the lean least-squares fit first, and where
        // it strays between the samples, the spline through every sample
        // decides where to refine.
        let mut split: Vec<bool> = Vec::new();
        for interpolate in [false, true] {
            let (knots, mut controls) = fit_family::<D>(
                std::slice::from_ref(&data),
                &ts,
                degree,
                tolerance * 0.5,
                closed,
                interpolate,
                &[],
                false,
            )?;
            let Some(control) = controls.pop() else {
                ogeom_bail!(NotDone, "the curve fit solved nothing");
            };
            let mut error = residuals::<D>(&knots, &control, &data, &ts)
                .iter()
                .fold(0.0_f64, |acc, e| acc.max(e.1));
            split = vec![false; ts.len() - 1];
            let mut worst = Vec::with_capacity(ts.len() - 1);
            for (i, pair) in ts.windows(2).enumerate() {
                let mut off = (0.0_f64, 0);
                for (k, t) in eighths(pair[0], pair[1]).into_iter().enumerate() {
                    let here = distance::<D>(&at(t)?, &position::<D>(&knots, &control, t));
                    if here > off.0 {
                        off = (here, k);
                    }
                }
                error = error.max(off.0);
                split[i] = off.0 > tolerance;
                worst.push(off);
            }
            let climbed = climbed_peaks(&ts, &worst, error, |t| {
                Ok(distance::<D>(&at(t)?, &position::<D>(&knots, &control, t)))
            });
            for (i, peak) in climbed {
                error = error.max(peak);
                split[i] |= peak > tolerance;
            }
            // A fit the eighth points pass is read densely before its error
            // is stated: a peak between them can stand over the target.
            let read = error <= tolerance;
            if read {
                let dense = read_densely(
                    &mut |t| {
                        Ok(at(t).map_or(0.0, |p| {
                            distance::<D>(&p, &position::<D>(&knots, &control, t))
                        }))
                    },
                    &knot_breaks(&knots, (ts[0], ts[ts.len() - 1])),
                    0,
                    error,
                    tolerance,
                )?;
                error = error.max(dense.widest.0);
                for t in dense.misses {
                    split[span_of(&ts, t)] = true;
                }
            }
            if error <= tolerance {
                return Ok((knots, control, error, true));
            }
            if best.as_ref().is_none_or(|b| error < b.2) {
                best = Some((knots, control, error, read));
            }
        }
        let mut next = Vec::with_capacity(ts.len() * 2);
        for (pair, &m) in ts.windows(2).zip(&split) {
            next.push(pair[0]);
            if m {
                next.push(f64::midpoint(pair[0], pair[1]));
            }
        }
        next.push(ts[ts.len() - 1]);
        if next.len() == ts.len() || next.len() > most + 1 {
            break;
        }
        ts = next;
    }
    let Some((knots, control, mut error, read)) = best else {
        ogeom_bail!(Construction, "the curve could not be sampled");
    };
    if !read {
        let dense = read_densely(
            &mut |t| {
                Ok(at(t).map_or(0.0, |p| {
                    distance::<D>(&p, &position::<D>(&knots, &control, t))
                }))
            },
            &knot_breaks(&knots, (ts[0], ts[ts.len() - 1])),
            0,
            error,
            tolerance,
        )?;
        error = error.max(dense.widest.0);
    }
    Ok((knots, control, error, false))
}

/// Fit a curve on a surface and its chart image together to the trace
/// `trace` gives over its parameter: at each parameter the point in space
/// and its chart position, `on_surface` the surface the chart belongs to.
///
/// One knot vector and one parameterization for both, fixed at `ts` to
/// start with, so the two are same-parameter with the trace and with each
/// other. Every sample and every span's eighth points are measured in
/// space: the curve and the surface at the image's chart position, each
/// against the trace's point and against each other, and the widest
/// misses climbed to the tops of their peaks between them, and once those
/// pass, every knot span of the fit at eight points, as
/// [`fit_curve_sampled`] reads it. A miss splits
/// the span, and this repeats until `tolerance` holds at every point measured or the samples would
/// pass 1024 spans; the error is the worst measured either way. The chart
/// coordinates are weighed in the fit by the surface's own speed along
/// them, so a residual in either space counts as the distance it is.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the
/// tolerance is not a distance, the degree is zero or the parameters do not
/// strictly increase with at least two; whatever `trace` refuses, and
/// whatever `on_surface` refuses other than a position off the surface's
/// domain, which a fitted image straying there misses by more than any
/// distance and is refined away from.
pub fn fit_trace_sampled(
    mut trace: impl FnMut(f64) -> OgeomResult<(Point, Point2)>,
    on_surface: impl Fn(Point2) -> OgeomResult<Point>,
    ts: &[f64],
    degree: usize,
    tolerance: f64,
    tol: Tolerances,
) -> OgeomResult<(Fitted<BSplineCurve>, BSpline2d)> {
    use crate::traits::Curve3d as _;
    const MOST: usize = 1024;
    if !(tolerance > 0.0 && tolerance.is_finite()) {
        ogeom_bail!(Construction, "a tolerance of {tolerance} is not a distance");
    }
    if degree == 0 {
        ogeom_bail!(Construction, "a fit needs a degree of at least one");
    }
    if ts.len() < 2 || ts.iter().any(|t| !t.is_finite()) || ts.windows(2).any(|w| w[1] <= w[0]) {
        ogeom_bail!(
            Construction,
            "a sampled fit needs strictly increasing parameters, at least two"
        );
    }
    let mut ts = ts.to_vec();
    let mut seen: ogeom_core::FastMap<u64, (Point, Point2)> = ogeom_core::FastMap::default();
    let mut at = |t: f64| -> OgeomResult<(Point, Point2)> {
        if let Some(p) = seen.get(&t.to_bits()) {
            return Ok(*p);
        }
        let p = trace(t)?;
        seen.insert(t.to_bits(), p);
        Ok(p)
    };
    // The surface's speed along each chart direction, at the first samples.
    let mut speed = [0.0_f64; 2];
    for t in &ts {
        let (_, uv) = at(*t)?;
        let here = on_surface(uv)?;
        for (axis, step) in [(0, Point2::new(1.0, 0.0)), (1, Point2::new(0.0, 1.0))] {
            let h = 1e-6 * (1.0 + uv.x.abs().max(uv.y.abs()));
            // Stepped back where the step forward leaves the surface's
            // domain, as at a sample on a sphere's pole.
            let mut moved = None;
            for sign in [1.0, -1.0] {
                let at = Point2::new(uv.x + sign * step.x * h, uv.y + sign * step.y * h);
                match on_surface(at) {
                    Ok(p) => {
                        moved = Some(p);
                        break;
                    }
                    Err(ogeom_core::OgeomError::Domain(_)) => {}
                    Err(e) => return Err(e),
                }
            }
            if let Some(p) = moved {
                speed[axis] = speed[axis].max(p.distance(here) / h);
            }
        }
    }
    let weight = speed.map(|s| if s > 0.0 && s.is_finite() { s } else { 1.0 });
    // The best fit so far, and whether it was read densely.
    let mut best: Option<((Fitted<BSplineCurve>, BSpline2d), bool)> = None;
    loop {
        let samples = ts.iter().map(|t| at(*t)).collect::<OgeomResult<Vec<_>>>()?;
        let data: Vec<[f64; 5]> = samples
            .iter()
            .map(|(p, uv)| [p.x, p.y, p.z, uv.x * weight[0], uv.y * weight[1]])
            .collect();
        let mut split: Vec<bool> = Vec::new();
        for interpolate in [false, true] {
            let (knots, mut controls) = fit_family::<5>(
                std::slice::from_ref(&data),
                &ts,
                degree,
                tolerance * 0.5,
                false,
                interpolate,
                &[],
                false,
            )?;
            let Some(control) = controls.pop() else {
                ogeom_bail!(NotDone, "the trace fit solved nothing");
            };
            let curve = build_curve_3(
                knots.clone(),
                control.iter().map(|c| [c[0], c[1], c[2]]).collect(),
                tol,
            )?;
            let image = BSpline2d::new(
                knots,
                control
                    .iter()
                    .map(|c| Point2::new(c[3] / weight[0], c[4] / weight[1]))
                    .collect(),
                tol,
            )?;
            let mut off_at = |t: f64| -> OgeomResult<f64> {
                let (p, _) = at(t)?;
                let c = curve.point_at(t, tol)?;
                let Some(s) = lifted(&on_surface, &image, t, tol)? else {
                    return Ok(f64::INFINITY);
                };
                Ok(c.distance(p).max(s.distance(p)).max(c.distance(s)))
            };
            let mut error = 0.0_f64;
            for t in &ts {
                error = error.max(off_at(*t)?);
            }
            split = vec![false; ts.len() - 1];
            let mut worst = Vec::with_capacity(ts.len() - 1);
            for (i, pair) in ts.windows(2).enumerate() {
                let mut off = (0.0_f64, 0);
                for (k, t) in eighths(pair[0], pair[1]).into_iter().enumerate() {
                    let here = off_at(t)?;
                    if here > off.0 {
                        off = (here, k);
                    }
                }
                error = error.max(off.0);
                split[i] = off.0 > tolerance;
                worst.push(off);
            }
            for (i, peak) in climbed_peaks(&ts, &worst, error, &mut off_at) {
                error = error.max(peak);
                split[i] |= peak > tolerance;
            }
            // As the curve fit: read densely before the error is stated.
            let read = error <= tolerance;
            if read {
                let dense = read_densely(
                    &mut |t| Ok(off_at(t).unwrap_or(0.0)),
                    &knot_breaks(curve.knots(), (ts[0], ts[ts.len() - 1])),
                    0,
                    error,
                    tolerance,
                )?;
                error = error.max(dense.widest.0);
                for t in dense.misses {
                    split[span_of(&ts, t)] = true;
                }
            }
            let met = error <= tolerance;
            let candidate = (Fitted { curve, error, met }, image);
            if met {
                return Ok(candidate);
            }
            if best.as_ref().is_none_or(|b| error < b.0.0.error) {
                best = Some((candidate, read));
            }
        }
        let mut next = Vec::with_capacity(ts.len() * 2);
        for (pair, &m) in ts.windows(2).zip(&split) {
            next.push(pair[0]);
            if m {
                next.push(f64::midpoint(pair[0], pair[1]));
            }
        }
        next.push(ts[ts.len() - 1]);
        if next.len() == ts.len() || next.len() > MOST + 1 {
            break;
        }
        ts = next;
    }
    let Some(((mut fitted, image), read)) = best else {
        ogeom_bail!(Construction, "the trace could not be sampled");
    };
    if !read {
        let mut off_at = |t: f64| -> OgeomResult<f64> {
            let (p, _) = at(t)?;
            let c = fitted.curve.point_at(t, tol)?;
            let Some(s) = lifted(&on_surface, &image, t, tol)? else {
                return Ok(f64::INFINITY);
            };
            Ok(c.distance(p).max(s.distance(p)).max(c.distance(s)))
        };
        let dense = read_densely(
            &mut |t| Ok(off_at(t).unwrap_or(0.0)),
            &knot_breaks(fitted.curve.knots(), (ts[0], ts[ts.len() - 1])),
            0,
            fitted.error,
            tolerance,
        )?;
        fitted.error = fitted.error.max(dense.widest.0);
    }
    Ok((fitted, image))
}

/// The surface at a fitted chart image's position at `t`; `None` where the
/// image strays off the surface's domain, a candidate that misses the
/// trace by more than any distance.
fn lifted(
    on_surface: &impl Fn(Point2) -> OgeomResult<Point>,
    image: &BSpline2d,
    t: f64,
    tol: Tolerances,
) -> OgeomResult<Option<Point>> {
    use crate::traits::Curve2d as _;
    match on_surface(image.point_at(t, tol)?) {
        Ok(p) => Ok(Some(p)),
        Err(ogeom_core::OgeomError::Domain(_)) => Ok(None),
        Err(e) => Err(e),
    }
}

/// The tops of the peaks of a curve fit's miss beside its widest eighth
/// point readings: for each span whose widest reading `worst` (the miss and
/// which eighth) stands near the widest of all, `error`, the miss `off`
/// climbed by golden sections between the eighths either side of it, at
/// most `PEAKS` spans, widest first. An eighth point stands up to a tenth
/// below the peak beside it where the miss swells between them. A point
/// `off` cannot read is passed over.
fn climbed_peaks(
    ts: &[f64],
    worst: &[(f64, usize)],
    error: f64,
    mut off: impl FnMut(f64) -> OgeomResult<f64>,
) -> Vec<(usize, f64)> {
    const NEAR_PEAK: f64 = 0.8;
    const PEAKS: usize = 16;
    let ratio = (5.0_f64.sqrt() - 1.0) / 2.0;
    let mut spans: Vec<usize> = (0..worst.len())
        .filter(|&i| worst[i].0 >= error * NEAR_PEAK && worst[i].0 > 0.0)
        .collect();
    spans.sort_by(|&a, &b| worst[b].0.total_cmp(&worst[a].0));
    let mut read = |t: f64| off(t).unwrap_or(0.0);
    spans
        .into_iter()
        .take(PEAKS)
        .map(|i| {
            let (lo, hi) = (ts[i], ts[i + 1]);
            let eighth = (hi - lo) / 8.0;
            // The k-th eighth point stands k + 1 eighths in.
            let at = lo + eighth * precise(worst[i].1 + 1);
            let (mut a, mut b) = (at - eighth, at + eighth);
            let mut top = worst[i].0;
            for _ in 0..30 {
                let (c, d) = (b - (b - a) * ratio, a + (b - a) * ratio);
                let (rc, rd) = (read(c), read(d));
                top = top.max(rc).max(rd);
                if rc > rd {
                    b = d;
                } else {
                    a = c;
                }
            }
            (i, top)
        })
        .collect()
}

/// The eighth points of a span, its middle among them: where a curve fit
/// is measured between two samples. A split makes the middle a sample and
/// these points the new spans' quarters and middles, so most are not
/// evaluated twice.
fn eighths(a: f64, b: f64) -> [f64; 7] {
    let mid = f64::midpoint(a, b);
    let (q1, q3) = (f64::midpoint(a, mid), f64::midpoint(mid, b));
    [
        f64::midpoint(a, q1),
        q1,
        f64::midpoint(q1, mid),
        mid,
        f64::midpoint(mid, q3),
        q3,
        f64::midpoint(q3, b),
    ]
}

fn build_curve_3(
    knots: KnotVector,
    control: Vec<[f64; 3]>,
    tol: Tolerances,
) -> OgeomResult<BSplineCurve> {
    BSplineCurve::new(
        knots,
        control
            .into_iter()
            .map(|c| Point::new(c[0], c[1], c[2]))
            .collect(),
        tol,
    )
}

/// As [`fit_points_2d_at`], with the loop's join made C1.
///
/// The chart image of a closed curve either returns to its first point or,
/// crossing its surface's seam, to that point one period over. Both are the
/// same loop, and the join constraint holds either way because it speaks
/// derivatives, not positions.
///
/// # Errors
///
/// As [`fit_points_2d_at`].
pub fn fit_points_2d_at_closed(
    parameters: &[f64],
    points: &[Point2],
    degree: usize,
    tolerance: f64,
    tol: Tolerances,
) -> OgeomResult<Fitted<BSpline2d>> {
    fit_points_2d_at_inner(parameters, points, degree, tolerance, true, tol)
}

fn fit_points_2d_at_inner(
    parameters: &[f64],
    points: &[Point2],
    degree: usize,
    tolerance: f64,
    closed: bool,
    tol: Tolerances,
) -> OgeomResult<Fitted<BSpline2d>> {
    if parameters.len() != points.len() {
        ogeom_bail!(Construction, "one parameter per point, or the fit is a lie");
    }
    if parameters.windows(2).any(|w| w[1] <= w[0]) {
        ogeom_bail!(Construction, "fixed parameters must strictly increase");
    }
    if !tolerance.is_finite() || tolerance <= 0.0 {
        ogeom_bail!(Construction, "a tolerance of {tolerance} is not a distance");
    }
    if degree == 0 {
        ogeom_bail!(Construction, "a fit needs a degree of at least one");
    }
    let data: Vec<[f64; 2]> = points.iter().map(|p| [p.x, p.y]).collect();
    if data.len() < 2 {
        ogeom_bail!(Construction, "a fit needs at least two points");
    }
    let degree = degree.min(data.len() - 1);
    let (a, b) = (parameters[0], parameters[parameters.len() - 1]);
    let mut knots = single_span(degree, a, b)?;

    const ROUNDS: usize = 32;
    let mut best: Option<(KnotVector, Vec<[f64; 2]>, f64)> = None;
    for _ in 0..ROUNDS {
        // A refinement can place a knot in a span the fixed parameters never
        // visit: clustered samples leave the system singular. That kills
        // the *round*, not the fit: the best earlier round still stands.
        let control = match least_squares::<2>(&knots, &data, parameters, closed) {
            Ok(control) => control,
            Err(e) => {
                if best.is_some() {
                    break;
                }
                return Err(e);
            }
        };
        let errors = residuals::<2>(&knots, &control, &data, parameters);
        let worst = errors.iter().fold(0.0_f64, |acc, e| acc.max(e.1));
        if best.as_ref().is_none_or(|(_, _, held)| worst < *held) {
            best = Some((knots.clone(), control.clone(), worst));
        }
        if worst <= tolerance {
            let curve = BSpline2d::new(
                knots,
                control
                    .into_iter()
                    .map(|c| Point2::new(c[0], c[1]))
                    .collect(),
                tol,
            )?;
            return Ok(Fitted {
                curve,
                error: worst,
                met: true,
            });
        }
        if knots.control_point_count() >= data.len() {
            break;
        }
        let Some(refined) = refined_where_bad(&knots, &errors, tolerance, false)? else {
            break;
        };
        knots = refined_or_interpolating(refined, parameters, closed)?;
    }
    let (knots, control, error) = best.ok_or_else(|| {
        ogeom_core::ogeom_err!(Construction, "the fixed-parameter fit never solved")
    })?;
    let curve = BSpline2d::new(
        knots,
        control
            .into_iter()
            .map(|c| Point2::new(c[0], c[1]))
            .collect(),
        tol,
    )?;
    Ok(Fitted {
        curve,
        error,
        met: false,
    })
}

/// How a free fit first assigns parameters to its samples, before the
/// correction rounds move them to the feet.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
enum Spacing {
    /// The square root of each chord: the general choice, which keeps a
    /// corner from pulling the parameterization through it.
    Centripetal,
    /// Each chord as it is: exact for a straight run however unevenly it
    /// was sampled, and right for anything traced smoothly. A marched
    /// section's walk halves its step to a micron at a window's rim and
    /// grows it back only twofold a point, so its samples crowd one end.
    /// There the centripetal guess sits so far from the feet that the
    /// correction rounds never reach them, and a straight section fits as
    /// a curve far off its own line.
    ChordLength,
}

/// The dimension-generic core.
///
/// Least squares is solved coordinate by coordinate: the collocation matrix
/// depends only on the parameters and the knots, so the expensive part is
/// shared and each coordinate is one more right-hand side.
fn fit<const D: usize>(
    points: &[[f64; D]],
    degree: usize,
    tolerance: f64,
    smooth_loop: bool,
    tol: Tolerances,
) -> OgeomResult<(KnotVector, Vec<[f64; D]>, f64, bool)> {
    fit_spaced::<D>(
        points,
        degree,
        tolerance,
        smooth_loop.then_some(D),
        Spacing::Centripetal,
        tol,
    )
}

fn fit_spaced<const D: usize>(
    points: &[[f64; D]],
    degree: usize,
    tolerance: f64,
    smooth_loop: Option<usize>,
    spacing: Spacing,
    tol: Tolerances,
) -> OgeomResult<(KnotVector, Vec<[f64; D]>, f64, bool)> {
    if !tolerance.is_finite() || tolerance <= 0.0 {
        ogeom_bail!(Construction, "a tolerance of {tolerance} is not a distance");
    }
    if degree == 0 {
        ogeom_bail!(Construction, "a fit needs a degree of at least one");
    }
    let points = collapse::<D>(points, tol);
    if points.len() < 2 {
        ogeom_bail!(
            Construction,
            "a fit needs at least two distinct points, got {}",
            points.len()
        );
    }
    let degree = degree.min(points.len() - 1);
    // A loop the caller asked to close smoothly: the first point repeated at
    // the end, and the join constrained to matching *tangents*. A surface
    // built over a C0 loop shows the crease. Opt-in, because the constraint
    // spends shape freedom at the join that an open fit keeps. The caller
    // names how many leading coordinates must meet: a section closes in
    // space while its chart images may end a period from where they began,
    // which the join holds just as well.
    let closed = smooth_loop.is_some_and(|meeting| {
        let (first, last) = (&points[0], &points[points.len() - 1]);
        first[..meeting]
            .iter()
            .zip(&last[..meeting])
            .map(|(x, y)| (x - y) * (x - y))
            .sum::<f64>()
            .sqrt()
            <= tol.confusion()
    });
    let parameters = match spacing {
        Spacing::Centripetal => centripetal::<D>(&points),
        Spacing::ChordLength => chordal::<D>(&points),
    };

    // Start with the fewest control points a clamped curve of this degree can
    // have: one Bézier span. Refinement adds knots only where the error says.
    let (a, b) = (parameters[0], parameters[parameters.len() - 1]);
    let mut knots = single_span(degree, a, b)?;

    // Each round may split every offending span, so the count grows by at most
    // a factor of two a round. A cap keeps a pathological input from spinning.
    const ROUNDS: usize = 32;
    let mut parameters = parameters;
    let mut best: Option<(KnotVector, Vec<[f64; D]>, f64)> = None;
    for _ in 0..ROUNDS {
        // Parameter correction in an earlier round can slide the data out of
        // a span that refinement checked *before* the slide, and the solve
        // reports itself singular. The rounds before it stand: keep the best
        // of them rather than promoting a bookkeeping casualty to an error.
        let control = match least_squares::<D>(&knots, &points, &parameters, closed) {
            Ok(control) => control,
            Err(e) => {
                if let Some((knots, control, worst)) = best {
                    return Ok((knots, control, worst, false));
                }
                return Err(e);
            }
        };
        // Parameter correction, and it is not a refinement. It is what makes
        // the residual mean anything. The residual is measured at each point's
        // assigned parameter, and the centripetal assignment is a guess: where
        // it drifts from the curve's own flow, a *perfect* curve still shows an
        // error at the assigned spot, the loop reads that as the curve's
        // fault, and refinement adds knots forever against a floor it can
        // never get under. Projecting each point onto the current curve
        // (Newton on the foot of the perpendicular) removes the
        // parameterization's share of the error and leaves the curve's.
        for _ in 0..2 {
            correct_parameters::<D>(&knots, &control, &points, &mut parameters, closed);
        }
        let mut errors = residuals::<D>(&knots, &control, &points, &parameters);
        wandering::<D>(&knots, &control, &points, &parameters, &mut errors);
        let worst = errors.iter().fold(0.0_f64, |acc, e| acc.max(e.1));
        if best.as_ref().is_none_or(|(_, _, held)| worst < *held) {
            best = Some((knots.clone(), control.clone(), worst));
        }
        if worst <= tolerance {
            return Ok((knots, control, worst, true));
        }

        // More control points than data points is interpolation wearing a
        // different name, and past that adding knots buys nothing. A closed
        // fit spent one control point on the C1 join, so its budget runs one
        // further.
        if knots.control_point_count() >= points.len() + usize::from(closed) {
            break;
        }
        let Some(refined) = refined_where_bad(&knots, &errors, tolerance, false)? else {
            break;
        };
        knots = refined_or_interpolating(refined, &parameters, closed)?;
    }

    #[allow(clippy::unwrap_used, reason = "at least one round always runs")]
    let (knots, control, worst) = best.unwrap();
    Ok((knots, control, worst, false))
}

/// Drop consecutive duplicates, which contribute a zero-length chord and make
/// the parameterization stall.
fn collapse<const D: usize>(points: &[[f64; D]], tol: Tolerances) -> Vec<[f64; D]> {
    let mut out: Vec<[f64; D]> = Vec::with_capacity(points.len());
    for p in points {
        if out
            .last()
            .is_some_and(|q| distance::<D>(p, q) <= tol.confusion() * 0.01)
        {
            continue;
        }
        out.push(*p);
    }
    out
}

fn distance<const D: usize>(a: &[f64; D], b: &[f64; D]) -> f64 {
    a.iter()
        .zip(b)
        .map(|(x, y)| (x - y) * (x - y))
        .sum::<f64>()
        .sqrt()
}

/// Fit a surface to *scattered* points: no grid required.
///
/// The points are parameterized by projection onto the cloud's own
/// principal plane (the two dominant directions of its covariance), and
/// fitted by tensor-product least squares with a bending penalty at
/// `smoothing`, which is what keeps the system solvable where the scatter
/// leaves basis functions unsupported. The reported error is the honest
/// maximum distance from any input to the surface at its assigned
/// parameters.
///
/// The cloud must be a *height field* over its principal plane: points
/// that fold over (a closed shell, a cliff) project onto each other, and
/// the fit answers with a large reported error rather than a lie.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if
/// there are too few points for the control budget, the budget cannot
/// carry the degree, or the weight is not finite and non-negative.
pub fn fit_surface_scattered(
    points: &[Point],
    degree: usize,
    controls: (usize, usize),
    smoothing: f64,
    tol: Tolerances,
) -> OgeomResult<Fitted<crate::BSplineSurface>> {
    use crate::traits::Surface as _;
    if !smoothing.is_finite() || smoothing < 0.0 {
        ogeom_bail!(
            Construction,
            "a smoothing weight of {smoothing} is not a weight"
        );
    }
    let (nu, nv) = controls;
    if nu < degree + 1 || nv < degree + 1 {
        ogeom_bail!(
            Construction,
            "{nu}x{nv} control points cannot carry degree {degree}"
        );
    }
    let m = points.len();
    if m < 4 {
        ogeom_bail!(Construction, "a scattered fit needs at least four points");
    }

    // The principal plane: centroid plus the covariance's two dominant
    // directions.
    let centroid = {
        let mut sum = ogeom_math::Vector::ZERO;
        for p in points {
            sum += p.to_vector();
        }
        sum / precise(m)
    };
    let mut covariance = nalgebra::Matrix3::<f64>::zeros();
    for p in points {
        let d = p.to_vector() - centroid;
        let v = nalgebra::Vector3::new(d.x, d.y, d.z);
        covariance += v * v.transpose();
    }
    let eigen = nalgebra::SymmetricEigen::new(covariance);
    let mut order: Vec<usize> = (0..3).collect();
    order.sort_by(|a, b| {
        eigen.eigenvalues[*b]
            .partial_cmp(&eigen.eigenvalues[*a])
            .unwrap_or(core::cmp::Ordering::Equal)
    });
    let axis = |i: usize| {
        let c = eigen.eigenvectors.column(order[i]);
        ogeom_math::Vector::new(c[0], c[1], c[2])
    };
    let (u_axis, v_axis) = (axis(0), axis(1));

    // Parameters over the knot domains, from the projected extents.
    let mut spans = Vec::with_capacity(m);
    let (mut ulo, mut uhi) = (f64::INFINITY, f64::NEG_INFINITY);
    let (mut vlo, mut vhi) = (f64::INFINITY, f64::NEG_INFINITY);
    for p in points {
        let d = p.to_vector() - centroid;
        let (pu, pv) = (d.dot(u_axis), d.dot(v_axis));
        ulo = ulo.min(pu);
        uhi = uhi.max(pu);
        vlo = vlo.min(pv);
        vhi = vhi.max(pv);
        spans.push((pu, pv));
    }
    if uhi - ulo <= tol.confusion() || vhi - vlo <= tol.confusion() {
        ogeom_bail!(
            Construction,
            "the cloud is flat in a principal direction; fit a curve"
        );
    }
    let u_knots = KnotVector::clamped_uniform(degree, nu)?;
    let v_knots = KnotVector::clamped_uniform(degree, nv)?;
    let (ka, kb) = (u_knots.domain_start(), u_knots.domain_end());
    let (la, lb) = (v_knots.domain_start(), v_knots.domain_end());
    let parameters: Vec<(f64, f64)> = spans
        .iter()
        .map(|(pu, pv)| {
            (
                ka + (kb - ka) * ((pu - ulo) / (uhi - ulo)),
                la + (lb - la) * ((pv - vlo) / (vhi - vlo)),
            )
        })
        .collect();

    // Tensor design matrix and the bending penalty in both directions.
    let n = nu * nv;
    let mut a = nalgebra::DMatrix::<f64>::zeros(m, n);
    for (k, &(pu, pv)) in parameters.iter().enumerate() {
        let uspan = u_knots.span_unchecked(pu);
        let vspan = v_knots.span_unchecked(pv);
        let ub = u_knots.basis(uspan, pu);
        let vb = v_knots.basis(vspan, pv);
        let ufirst = uspan - degree;
        let vfirst = vspan - degree;
        for (j, bv) in vb.iter().enumerate() {
            for (i, bu) in ub.iter().enumerate() {
                a[(k, (ufirst + i) * nv + (vfirst + j))] = bu * bv;
            }
        }
    }
    let mut k_full = a.transpose() * &a;
    let mut add_penalty = |along_u: bool| {
        // Second differences along one grid direction, in the grid's own
        // layout: flat index = i_u · nv + j_v.
        let (count_a, count_b) = if along_u { (nu, nv) } else { (nv, nu) };
        for jb in 0..count_b {
            for ia in 0..count_a.saturating_sub(2) {
                let base = |offset: usize| -> usize {
                    if along_u {
                        (ia + offset) * nv + jb
                    } else {
                        jb * nv + ia + offset
                    }
                };
                let idx = [base(0), base(1), base(2)];
                let w = [1.0, -2.0, 1.0];
                for x in 0..3 {
                    for y in 0..3 {
                        k_full[(idx[x], idx[y])] += smoothing * w[x] * w[y];
                    }
                }
            }
        }
    };
    add_penalty(true);
    add_penalty(false);

    let decomposition = k_full.clone().lu();
    let mut control = vec![Point::ORIGIN; n];
    for axis_i in 0..3 {
        let b = nalgebra::DVector::from_iterator(m, points.iter().map(|p| [p.x, p.y, p.z][axis_i]));
        let rhs = a.transpose() * &b;
        let Some(solved) = decomposition.solve(&rhs) else {
            ogeom_bail!(
                Numeric,
                "the scattered system is singular; raise the smoothing or lower the controls"
            );
        };
        for (slot, value) in solved.iter().enumerate() {
            match axis_i {
                0 => control[slot].x = *value,
                1 => control[slot].y = *value,
                _ => control[slot].z = *value,
            }
        }
    }

    let grid = ogeom_math::ControlGrid::new(control, nu, nv)?;
    let surface = crate::BSplineSurface::new(u_knots, v_knots, &grid, tol)?;
    let mut error = 0.0f64;
    for (p, &(pu, pv)) in points.iter().zip(&parameters) {
        error = error.max(surface.point_at(pu, pv, tol)?.distance(*p));
    }
    Ok(Fitted {
        curve: surface,
        error,
        met: true,
    })
}

/// Fit a surface to scattered points at parameters the caller chose, over
/// the caller's knot vectors, faired by the thin-plate energy.
///
/// Where [`fit_surface_scattered`] reads the parameters off the cloud's
/// principal plane, here they come from a chart of the caller's: one that
/// unfolds a region the principal plane would fold over. The control net
/// minimises the mean squared distance from each point to the surface at its
/// parameters plus `fairing` times the thin-plate energy, the integral of
/// `|S_uu|² + 2|S_uv|² + |S_vv|²` over the whole knot domain. The energy
/// settles the controls no point reaches (a span the points leave empty, a
/// margin past them) on the smoothest continuation of the rest. The
/// reported error is the largest distance from any point to the surface at
/// its parameters; `met` is always true, there being no target.
///
/// The normal equations are banded (a control couples only to those within
/// a degree of it in each direction), and are solved in that band.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the
/// points and parameters differ in number, there are fewer than four, a
/// parameter lies outside the knot domain, or the weight is not finite and
/// non-negative; [`OgeomError::Numeric`](ogeom_core::OgeomError::Numeric)
/// if the system is singular.
pub fn fit_surface_scattered_at(
    points: &[Point],
    parameters: &[(f64, f64)],
    u_knots: &KnotVector,
    v_knots: &KnotVector,
    fairing: f64,
    tol: Tolerances,
) -> OgeomResult<Fitted<crate::BSplineSurface>> {
    use crate::traits::Surface as _;
    if !fairing.is_finite() || fairing < 0.0 {
        ogeom_bail!(
            Construction,
            "a fairing weight of {fairing} is not a weight"
        );
    }
    let m = points.len();
    if m != parameters.len() {
        ogeom_bail!(
            Construction,
            "{m} points were given {} parameters",
            parameters.len()
        );
    }
    if m < 4 {
        ogeom_bail!(Construction, "a scattered fit needs at least four points");
    }
    let (ua, ub) = u_knots.domain();
    let (va, vb) = v_knots.domain();
    if let Some(&(u, v)) = parameters
        .iter()
        .find(|&&(u, v)| !(u >= ua && u <= ub && v >= va && v <= vb))
    {
        ogeom_bail!(
            Construction,
            "the parameter ({u}, {v}) is outside the domain [{ua}, {ub}] x [{va}, {vb}]"
        );
    }
    let (pu, pv) = (u_knots.degree(), v_knots.degree());
    let (nu, nv) = (u_knots.control_point_count(), v_knots.control_point_count());
    let n = nu * nv;
    // Two controls interact only within a degree of each other in both
    // directions, so with the controls numbered along the direction that has
    // more of them, row by row of the other, the system's band is a few rows
    // of the shorter direction wide.
    let rows_along_u = nu >= nv;
    let slot = |i: usize, j: usize| {
        if rows_along_u { i * nv + j } else { j * nu + i }
    };
    let band = if rows_along_u {
        pu * nv + pv
    } else {
        pv * nu + pu
    };
    let mut lower = vec![vec![0.0; band + 1]; n];
    let mut rhs = vec![[0.0; 3]; n];
    let add = |lower: &mut [Vec<f64>], r: usize, c: usize, value: f64| {
        let (r, c) = if r >= c { (r, c) } else { (c, r) };
        lower[r][r - c] += value;
    };

    let share = 1.0 / precise(m);
    let mut terms: Vec<(usize, f64)> = Vec::with_capacity((pu + 1) * (pv + 1));
    for (p, &(u, v)) in points.iter().zip(parameters) {
        let (us, vs) = (u_knots.span_unchecked(u), v_knots.span_unchecked(v));
        let (bu, bv) = (u_knots.basis(us, u), v_knots.basis(vs, v));
        terms.clear();
        for (a, x) in bu.iter().enumerate() {
            for (b, y) in bv.iter().enumerate() {
                terms.push((slot(us - pu + a, vs - pv + b), x * y));
            }
        }
        for (k, &(r, x)) in terms.iter().enumerate() {
            rhs[r][0] += share * x * p.x;
            rhs[r][1] += share * x * p.y;
            rhs[r][2] += share * x * p.z;
            for &(c, y) in &terms[..=k] {
                add(&mut lower, r, c, share * x * y);
            }
        }
    }

    // The thin-plate energy. A tensor product's second derivatives
    // separate, so the energy between controls (a, b) and (c, d) is
    // `U2[a][c] V0[b][d] + 2 U1[a][c] V1[b][d] + U0[a][c] V2[b][d]`, where
    // `Uk` holds the integrals of products of the `k`th derivatives of the
    // `u` basis functions, and `Vk` of the `v` ones.
    if fairing > 0.0 {
        let (gu, gv) = (derivative_grams(u_knots), derivative_grams(v_knots));
        for a in 0..nu {
            for b in 0..nv {
                let r = slot(a, b);
                for c in a.saturating_sub(pu)..(a + pu + 1).min(nu) {
                    for d in b.saturating_sub(pv)..(b + pv + 1).min(nv) {
                        let col = slot(c, d);
                        if col > r {
                            continue;
                        }
                        let (f, x) = (a.min(c), a.abs_diff(c));
                        let (e, y) = (b.min(d), b.abs_diff(d));
                        let energy = gu[2][f][x] * gv[0][e][y]
                            + 2.0 * gu[1][f][x] * gv[1][e][y]
                            + gu[0][f][x] * gv[2][e][y];
                        lower[r][r - col] += fairing * energy;
                    }
                }
            }
        }
    }

    let Some(factor) = banded_cholesky(&lower, band) else {
        ogeom_bail!(
            Numeric,
            "the scattered system is singular; raise the fairing or lower the controls"
        );
    };
    let solved: Vec<Vec<f64>> = (0..3)
        .map(|axis| {
            let b: Vec<f64> = rhs.iter().map(|r| r[axis]).collect();
            banded_substitute(&factor, band, &b)
        })
        .collect();
    // The net in the grid's own order, `u` outer.
    let mut control = Vec::with_capacity(n);
    for i in 0..nu {
        for j in 0..nv {
            let k = slot(i, j);
            control.push(Point::new(solved[0][k], solved[1][k], solved[2][k]));
        }
    }
    let grid = ogeom_math::ControlGrid::new(control, nu, nv)?;
    let surface = crate::BSplineSurface::new(u_knots.clone(), v_knots.clone(), &grid, tol)?;
    let mut error = 0.0_f64;
    for (p, &(u, v)) in points.iter().zip(parameters) {
        error = error.max(surface.point_at(u, v, tol)?.distance(*p));
    }
    Ok(Fitted {
        curve: surface,
        error,
        met: true,
    })
}

/// The integrals over a knot vector's domain of the products of its basis
/// functions' derivatives, for orders nought to two: `grams[k][i][d]` is
/// the integral of `N_i^(k) N_(i+d)^(k)`, for `d` up to the degree (the
/// functions further apart share no span).
fn derivative_grams(knots: &KnotVector) -> [Vec<Vec<f64>>; 3] {
    let p = knots.degree();
    let n = knots.control_point_count();
    let mut grams: [Vec<Vec<f64>>; 3] = core::array::from_fn(|_| vec![vec![0.0; p + 1]; n]);
    let distinct = knots.distinct();
    for w in distinct.windows(2) {
        let (a, b) = (w[0].0, w[1].0);
        for (t, weight) in ogeom_math::gauss_legendre_rule(a, b) {
            let span = knots.span_unchecked(t);
            let rows = knots.basis_derivatives(span, t, 2);
            for (k, gram) in grams.iter_mut().enumerate() {
                let Some(row) = rows.get(k) else {
                    continue;
                };
                for i in 0..=p {
                    for j in i..=p {
                        gram[span - p + i][j - i] += weight * row[i] * row[j];
                    }
                }
            }
        }
    }
    grams
}

/// Fill the region bounded by four curves with a fitted patch: the
/// transfinite Coons blend of the boundaries, sampled and fitted, its error
/// reported.
///
/// The curves must close corner to corner in the order given (`bottom`
/// runs with `u`, `top` above it, `left` and `right` with `v`), each
/// traversed over its own domain. The patch is fitted to the Coons surface
/// by [`fit_surface_sampled`], starting from `samples` (at least four) a
/// side, so the stated error is measured between the samples as well as at
/// them. The Coons surface itself interpolates the boundaries exactly, so
/// the fit error bounds the distance between the returned patch's borders
/// and the boundary it fills.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the
/// corners do not meet within `tol`, plus whatever the fit refuses.
pub fn fill_boundary(
    bottom: &crate::Curve,
    top: &crate::Curve,
    left: &crate::Curve,
    right: &crate::Curve,
    samples: usize,
    tolerance: f64,
    tol: Tolerances,
) -> OgeomResult<Fitted<crate::BSplineSurface>> {
    use crate::traits::Curve3d as _;
    let samples = samples.max(4);
    let at = |curve: &crate::Curve, t: f64| -> OgeomResult<Point> {
        let (lo, hi) = curve.domain();
        curve.point_at(lo + (hi - lo) * t, tol)
    };
    // The four corners, each named twice. They must agree.
    let c00 = at(bottom, 0.0)?;
    let c10 = at(bottom, 1.0)?;
    let c01 = at(top, 0.0)?;
    let c11 = at(top, 1.0)?;
    let slack = tol.confusion() * 1e3;
    for (name, a, b) in [
        ("bottom-left", c00, at(left, 0.0)?),
        ("top-left", c01, at(left, 1.0)?),
        ("bottom-right", c10, at(right, 0.0)?),
        ("top-right", c11, at(right, 1.0)?),
    ] {
        if a.distance(b) > slack {
            ogeom_bail!(
                Construction,
                "the {name} corner does not close: the boundaries miss by {}",
                a.distance(b)
            );
        }
    }

    // The bilinearly blended Coons point: ruled in each direction, the
    // doubly-ruled corner sheet subtracted once.
    let coons = |u: f64, v: f64| -> OgeomResult<Point> {
        let cu0 = at(bottom, u)?;
        let cu1 = at(top, u)?;
        let d0v = at(left, v)?;
        let d1v = at(right, v)?;
        let ruled_u = cu0.to_vector() * (1.0 - v) + cu1.to_vector() * v;
        let ruled_v = d0v.to_vector() * (1.0 - u) + d1v.to_vector() * u;
        let corners = c00.to_vector() * ((1.0 - u) * (1.0 - v))
            + c10.to_vector() * (u * (1.0 - v))
            + c01.to_vector() * ((1.0 - u) * v)
            + c11.to_vector() * (u * v);
        Ok(Point::from_vector(ruled_u + ruled_v - corners))
    };
    // Fitted at the Coons surface's own parameters, `samples` a side to
    // start with, and checked and refined between them.
    fit_surface_sampled(
        coons,
        &Sampling::even(((0.0, 1.0), (0.0, 1.0)), samples - 1),
        3,
        tolerance,
        tol,
    )
}

/// Fit a rectangular grid of points with a tensor-product B-spline surface.
///
/// A grid fit kept honest the way the curve fit is: rows first,
/// columns second, both passes at *fixed* parameters: a grid's
/// parameterization is shared property, and correcting it per row is how a
/// grid stops being one. Each pass adapts one shared knot vector against the
/// worst residual across its whole family, so every row rides the same
/// basis, which is what makes the second pass a fit over control points
/// rather than a guess. The reported error is measured at the end, surface
/// against every input point, and `met` does not round up.
///
/// The error is measured at the grid points only. Between them the surface
/// follows the spline through the grid, which for a smooth grid differs
/// from the sampled shape by an amount that shrinks with the spacing (the
/// parameters are assigned from the chords, so even a quadratic surface is
/// followed closely rather than reproduced), but it is not bounded here.
///
/// `rows[j][i]` runs `i` along `u` and `j` along `v`.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the grid is
/// not rectangular or is smaller than two by two, or the tolerance is not a
/// distance.
pub fn fit_surface_grid(
    rows: &[Vec<Point>],
    degree: usize,
    tolerance: f64,
    tol: Tolerances,
) -> OgeomResult<Fitted<crate::BSplineSurface>> {
    fit_surface_grid_inner(rows, degree, tolerance, false, false, tol)
}

/// As [`fit_surface_grid`], parameterized by chord length instead of the
/// centripetal assignment, the open counterpart of
/// [`fit_surface_grid_closed_v_chordal`], for a marched band that stops at
/// its run-out instead of closing on itself.
///
/// # Errors
///
/// As [`fit_surface_grid`].
pub fn fit_surface_grid_chordal(
    rows: &[Vec<Point>],
    degree: usize,
    tolerance: f64,
    tol: Tolerances,
) -> OgeomResult<Fitted<crate::BSplineSurface>> {
    fit_surface_grid_inner(rows, degree, tolerance, false, true, tol)
}

/// As [`fit_surface_grid`], with the rows placed across `v` by chord
/// length: sections lofted through, where the spacing between sections is
/// the skin's own and may change along it. Centripetal parameters squash a
/// change of spacing (a step twice another's gets a parameter step only
/// about one and a half times as long), and a cubic through sections so
/// placed sags between them. Each row keeps the centripetal assignment
/// along `u`.
///
/// # Errors
///
/// As [`fit_surface_grid`].
pub fn fit_surface_grid_sections(
    rows: &[Vec<Point>],
    degree: usize,
    tolerance: f64,
    tol: Tolerances,
) -> OgeomResult<Fitted<crate::BSplineSurface>> {
    fit_surface_grid_parameterized(rows, degree, tolerance, false, (false, true), tol)
}

/// As [`fit_surface_grid`], with the `v` direction closed into a smooth loop.
///
/// The rows must form a loop (the first row repeated at the end), and the
/// join across it is made C1 the way [`fit_points_closed`] makes a curve's:
/// the constraint is eliminated inside the shared solve, so the two border
/// control rows are equal and the loop crosses its own seam smoothly.
///
/// # Errors
///
/// As [`fit_surface_grid`], and additionally if the rows are not a loop.
pub fn fit_surface_grid_closed_v(
    rows: &[Vec<Point>],
    degree: usize,
    tolerance: f64,
    tol: Tolerances,
) -> OgeomResult<Fitted<crate::BSplineSurface>> {
    if rows.len() < 3 {
        ogeom_bail!(Construction, "a closed skin needs at least three rows");
    }
    let (first, last) = (&rows[0], &rows[rows.len() - 1]);
    if first.len() != last.len()
        || first
            .iter()
            .zip(last)
            .any(|(a, b)| a.distance(*b) > tol.confusion() * 100.0)
    {
        ogeom_bail!(
            Construction,
            "a closed skin needs a loop: the first row repeated at the end"
        );
    }
    fit_surface_grid_inner(rows, degree, tolerance, true, false, tol)
}

/// As [`fit_surface_grid_closed_v`], parameterized by chord length instead
/// of the centripetal assignment.
///
/// The choice matters for dense data: a marched grid samples smooth curves
/// finely but not evenly, and chord length keeps the parameter's speed
/// uniform across the spacing jumps a march's closure leaves behind.
///
/// # Errors
///
/// As [`fit_surface_grid_closed_v`].
pub fn fit_surface_grid_closed_v_chordal(
    rows: &[Vec<Point>],
    degree: usize,
    tolerance: f64,
    tol: Tolerances,
) -> OgeomResult<Fitted<crate::BSplineSurface>> {
    if rows.len() < 3 {
        ogeom_bail!(Construction, "a closed skin needs at least three rows");
    }
    let (first, last) = (&rows[0], &rows[rows.len() - 1]);
    if first.len() != last.len()
        || first
            .iter()
            .zip(last)
            .any(|(a, b)| a.distance(*b) > tol.confusion() * 100.0)
    {
        ogeom_bail!(
            Construction,
            "a closed skin needs a loop: the first row repeated at the end"
        );
    }
    fit_surface_grid_inner(rows, degree, tolerance, true, true, tol)
}

fn fit_surface_grid_inner(
    rows: &[Vec<Point>],
    degree: usize,
    tolerance: f64,
    closed_v: bool,
    by_chord: bool,
    tol: Tolerances,
) -> OgeomResult<Fitted<crate::BSplineSurface>> {
    fit_surface_grid_parameterized(rows, degree, tolerance, closed_v, (by_chord, by_chord), tol)
}

/// The grid fit, with chord-length or centripetal parameters chosen per
/// direction: `by_chord.0` along the rows (`u`), `by_chord.1` across them
/// (`v`).
fn fit_surface_grid_parameterized(
    rows: &[Vec<Point>],
    degree: usize,
    tolerance: f64,
    closed_v: bool,
    by_chord: (bool, bool),
    tol: Tolerances,
) -> OgeomResult<Fitted<crate::BSplineSurface>> {
    if !tolerance.is_finite() || tolerance <= 0.0 {
        ogeom_bail!(Construction, "a tolerance of {tolerance} is not a distance");
    }
    let nv = rows.len();
    if nv < 2 {
        ogeom_bail!(Construction, "a surface fit needs at least two rows");
    }
    let nu = rows[0].len();
    if nu < 2 || rows.iter().any(|r| r.len() != nu) {
        ogeom_bail!(Construction, "a surface fit needs a rectangular grid");
    }
    let raw: Vec<Vec<[f64; 3]>> = rows
        .iter()
        .map(|r| r.iter().map(|p| [p.x, p.y, p.z]).collect())
        .collect();

    let (u_params, v_params) = averaged_parameters(&raw, by_chord);
    fit_grid_on(
        rows,
        &raw,
        (&u_params, &v_params),
        degree,
        tolerance,
        (closed_v, false),
        (&[], &[]),
        tol,
    )
}

/// The parameters [`fit_surface_grid`] and its chord-length and section
/// variants assign a rectangular grid: one shared assignment per
/// direction, the average of every row's (`u`) and every column's (`v`),
/// chord length or centripetal as `by_chord` says per direction.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the
/// grid is not rectangular or is smaller than two by two.
pub fn grid_parameters(
    rows: &[Vec<Point>],
    by_chord: (bool, bool),
) -> OgeomResult<(Vec<f64>, Vec<f64>)> {
    let nu = rows.first().map_or(0, Vec::len);
    if rows.len() < 2 || nu < 2 || rows.iter().any(|r| r.len() != nu) {
        ogeom_bail!(Construction, "a surface fit needs a rectangular grid");
    }
    let raw: Vec<Vec<[f64; 3]>> = rows
        .iter()
        .map(|r| r.iter().map(|p| [p.x, p.y, p.z]).collect())
        .collect();
    Ok(averaged_parameters(&raw, by_chord))
}

/// Averaged parameters: one shared assignment per direction. A family
/// with no extent (a row collapsed to a single point, the apex a skin
/// narrows to) has no parameterization opinion: its zero chords assign
/// [0, …, 0, 1], and averaging that in compresses everyone else's
/// parameters toward the start and sends the fitted border curves on an
/// oscillating sprint over the tail. It rides the others' parameters
/// instead. A constant row fits exactly at any assignment.
fn averaged_parameters(raw: &[Vec<[f64; 3]>], by_chord: (bool, bool)) -> (Vec<f64>, Vec<f64>) {
    let average = |families: &[Vec<[f64; 3]>], by_chord: bool| -> Vec<f64> {
        let mut sums = vec![0.0; families[0].len()];
        let mut counted = 0.0_f64;
        for family in families {
            let extent: f64 = family
                .windows(2)
                .map(|pair| distance::<3>(&pair[0], &pair[1]))
                .sum();
            if extent <= 0.0 {
                continue;
            }
            counted += 1.0;
            let assigned = if by_chord {
                chordal::<3>(family)
            } else {
                centripetal::<3>(family)
            };
            for (s, p) in sums.iter_mut().zip(assigned) {
                *s += p;
            }
        }
        if counted == 0.0 {
            // Every family degenerate: uniform is the only honest assignment.
            let n = sums.len();
            return (0..n)
                .map(|i| {
                    #[allow(clippy::cast_precision_loss)]
                    let t = i as f64 / (n - 1).max(1) as f64;
                    t
                })
                .collect();
        }
        sums.iter().map(|s| s / counted).collect()
    };
    let nu = raw[0].len();
    let u_params = average(raw, by_chord.0);
    let columns: Vec<Vec<[f64; 3]>> = (0..nu)
        .map(|i| raw.iter().map(|r| r[i]).collect())
        .collect();
    let v_params = average(&columns, by_chord.1);
    (u_params, v_params)
}

/// Fit a rectangular grid at *fixed* parameters: `u_params[i]` is the
/// surface's own `u` at `rows[j][i]`, and `v_params[j]` its `v`.
///
/// The surface counterpart of [`fit_points_at`]: the caller's chart is
/// kept, which is what lets a replacement surface take over every pcurve
/// drawn on the old one. As [`fit_surface_grid`], rows first and columns
/// second on shared knots, and the reported error is the surface against
/// every input point at its parameters. Between the grid points the fit is
/// not measured here; a caller holding the exact surface measures it there.
///
/// # Errors
///
/// As [`fit_surface_grid`], and the parameters must strictly increase and
/// be one per row and one per column.
pub fn fit_surface_grid_at(
    u_params: &[f64],
    v_params: &[f64],
    rows: &[Vec<Point>],
    degree: usize,
    tolerance: f64,
    tol: Tolerances,
) -> OgeomResult<Fitted<crate::BSplineSurface>> {
    fit_surface_grid_at_with_knots(u_params, v_params, rows, degree, tolerance, (&[], &[]), tol)
}

/// Interior knots a fit keeps, `(parameter, multiplicity)`.
type Kept<'a> = &'a [(f64, usize)];

/// As [`fit_surface_grid_at`], each direction's knots starting from the
/// interior knots given for it, `(parameter, multiplicity)`, and keeping
/// them. Where the chart the points come from is only C0 or C1 along a
/// parameter line, a knot there of multiplicity `degree` or `degree - 1`
/// lets the fit turn as the chart does; without one, a fit smoother than
/// its data converges on such a line only as fast as its spans shrink.
/// A multiplicity past the degree counts as the degree, and a knot outside
/// the parameters' open range is not used.
///
/// # Errors
///
/// As [`fit_surface_grid_at`].
pub fn fit_surface_grid_at_with_knots(
    u_params: &[f64],
    v_params: &[f64],
    rows: &[Vec<Point>],
    degree: usize,
    tolerance: f64,
    knots: (Kept<'_>, Kept<'_>),
    tol: Tolerances,
) -> OgeomResult<Fitted<crate::BSplineSurface>> {
    if !tolerance.is_finite() || tolerance <= 0.0 {
        ogeom_bail!(Construction, "a tolerance of {tolerance} is not a distance");
    }
    if rows.len() < 2 || rows.len() != v_params.len() {
        ogeom_bail!(
            Construction,
            "a surface fit needs at least two rows, one v parameter each"
        );
    }
    if u_params.len() < 2 || rows.iter().any(|r| r.len() != u_params.len()) {
        ogeom_bail!(
            Construction,
            "a surface fit needs a rectangular grid, one u parameter per column"
        );
    }
    if u_params.windows(2).any(|w| w[1] <= w[0]) || v_params.windows(2).any(|w| w[1] <= w[0]) {
        ogeom_bail!(Construction, "fixed parameters must strictly increase");
    }
    let raw: Vec<Vec<[f64; 3]>> = rows
        .iter()
        .map(|r| r.iter().map(|p| [p.x, p.y, p.z]).collect())
        .collect();
    fit_grid_on(
        rows,
        &raw,
        (u_params, v_params),
        degree,
        tolerance,
        (false, false),
        knots,
        tol,
    )
}

/// Where [`fit_surface_sampled`] starts sampling and how far it may refine.
#[derive(Clone, Debug)]
pub struct Sampling {
    /// The first grid's column parameters (`u`), strictly increasing.
    pub us: Vec<f64>,
    /// The first grid's row parameters (`v`), strictly increasing.
    pub vs: Vec<f64>,
    /// Whether the geometry is known between the columns (`.0`) and
    /// between the rows (`.1`). Only a direction known between its samples
    /// is checked there and refined; the other is fitted through its
    /// samples, which are all there is of it.
    pub between: (bool, bool),
    /// Whether the fit closes smoothly across `v`: the last row is the
    /// first again, and the loop crosses its seam C1 as in
    /// [`fit_surface_grid_closed_v`].
    pub closed_v: bool,
    /// The most spans a direction is refined to.
    pub most: usize,
}

impl Sampling {
    /// `spans` even spans over each side of `domain`, both directions
    /// known between their samples, open, refined to at most 512 spans.
    #[must_use]
    pub fn even(domain: ((f64, f64), (f64, f64)), spans: usize) -> Self {
        let even = |(a, b): (f64, f64)| -> Vec<f64> {
            let n = spans.max(1);
            (0..=n)
                .map(|i| {
                    if i == n {
                        b
                    } else {
                        a + (b - a) * precise(i) / precise(n)
                    }
                })
                .collect()
        };
        Self {
            us: even(domain.0),
            vs: even(domain.1),
            between: (true, true),
            closed_v: false,
            most: 512,
        }
    }
}

/// Fit a surface to the geometry `point` traces over its parameters,
/// measured between the samples as well as at them.
///
/// The fit is same-parameter with `point`: it is fitted at the grid's own
/// parameters, and the fit at `(u, v)` is compared with `point(u, v)`,
/// which bounds the distance between the two from above. Each span of a
/// direction known between its samples is checked at its quarter points,
/// along the sample lines and across them; a miss splits the spans it lies
/// in, its middle becoming a sample, so the grid (and the fit's knots,
/// which follow its residuals) thickens only where the fit strays. Where
/// the least-squares fit strays between the samples, the spline through
/// every sample (on knots that follow the samples) is measured too and
/// decides where to refine, and the closer of the two stands. This repeats
/// until `tolerance` holds at every point measured or a direction would
/// pass `sampling.most` spans. A fit the quarter points pass is read again
/// along its sample lines (and midway between them), at eight points of
/// every knot span of its own and its peaks climbed, along them and
/// across, before it is stated met: a peak at a crease of the geometry can
/// stand between the quarter points. The best fit is returned either way,
/// read so, its error the worst measured and `met` whether that is within
/// `tolerance`.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the
/// tolerance is not a distance or the parameters do not strictly increase
/// with at least two a side; whatever `point` refuses.
pub fn fit_surface_sampled(
    mut point: impl FnMut(f64, f64) -> OgeomResult<Point>,
    sampling: &Sampling,
    degree: usize,
    tolerance: f64,
    tol: Tolerances,
) -> OgeomResult<Fitted<crate::BSplineSurface>> {
    if !(tolerance > 0.0 && tolerance.is_finite()) {
        ogeom_bail!(Construction, "a tolerance of {tolerance} is not a distance");
    }
    let increasing = |ps: &[f64]| {
        ps.len() >= 2 && ps.iter().all(|p| p.is_finite()) && ps.windows(2).all(|w| w[1] > w[0])
    };
    if !increasing(&sampling.us) || !increasing(&sampling.vs) {
        ogeom_bail!(
            Construction,
            "a sampled fit needs strictly increasing parameters, two a side"
        );
    }
    let (mut us, mut vs) = (sampling.us.clone(), sampling.vs.clone());
    // The best fit so far, and whether it was read densely.
    let mut best: Option<(Fitted<crate::BSplineSurface>, bool)> = None;
    // Every point evaluated, by parameters: a split span's new sample is
    // the midpoint already checked.
    let mut seen: ogeom_core::FastMap<(u64, u64), Point> = ogeom_core::FastMap::default();
    let mut at = |u: f64, v: f64| -> OgeomResult<Point> {
        if let Some(p) = seen.get(&(u.to_bits(), v.to_bits())) {
            return Ok(*p);
        }
        let p = point(u, v)?;
        seen.insert((u.to_bits(), v.to_bits()), p);
        Ok(p)
    };
    loop {
        let mut rows = vs
            .iter()
            .map(|v| us.iter().map(|u| at(*u, *v)).collect())
            .collect::<OgeomResult<Vec<Vec<Point>>>>()?;
        if sampling.closed_v {
            let first = rows[0].clone();
            if let Some(last) = rows.last_mut() {
                *last = first;
            }
        }
        let raw: Vec<Vec<[f64; 3]>> = rows
            .iter()
            .map(|r| r.iter().map(|p| [p.x, p.y, p.z]).collect())
            .collect();
        // Half the tolerance at the samples leaves the other half for
        // between them. The least-squares fit is the lean one, but where it
        // strays between the samples that may be its knots rather than the
        // sampling: the spline through every sample, on knots that follow
        // the samples, then measures the sampling alone and decides where
        // to refine.
        let mut splits: Option<(Vec<bool>, Vec<bool>)> = None;
        for interpolate in [false, true] {
            let fitted = fit_grid_on(
                &rows,
                &raw,
                (&us, &vs),
                degree,
                tolerance * 0.5,
                (sampling.closed_v, interpolate),
                (&[], &[]),
                tol,
            )?;
            let (between, mut split_u, mut split_v) = measured_between(
                &mut at,
                &fitted.curve,
                (&us, &vs),
                sampling.between,
                tolerance,
                tol,
            )?;
            let mut error = fitted.error.max(between);
            // A fit the quarter points pass is read densely before its
            // error is stated: a peak between them can stand over the
            // target.
            let read_densely = error <= tolerance;
            if read_densely {
                let along = measured_along_lines(
                    &mut at,
                    &fitted.curve,
                    (&us, &vs),
                    sampling.between,
                    tolerance,
                    error,
                    (&mut split_u, &mut split_v),
                    tol,
                )?;
                error = error.max(along);
            }
            let candidate = Fitted {
                curve: fitted.curve,
                error,
                met: error <= tolerance,
            };
            if candidate.met {
                return Ok(candidate);
            }
            if best.as_ref().is_none_or(|b| error < b.0.error) {
                best = Some((candidate, read_densely));
            }
            splits = Some((split_u, split_v));
        }
        let Some((split_u, split_v)) = splits else {
            break;
        };
        let split = |knots: &[f64], marked: &[bool]| -> Vec<f64> {
            let mut out = Vec::with_capacity(knots.len() * 2);
            for (pair, &m) in knots.windows(2).zip(marked) {
                out.push(pair[0]);
                if m {
                    out.push(f64::midpoint(pair[0], pair[1]));
                }
            }
            out.push(knots[knots.len() - 1]);
            out
        };
        let (next_u, next_v) = (split(&us, &split_u), split(&vs, &split_v));
        let grew = next_u.len() > us.len() || next_v.len() > vs.len();
        if !grew || next_u.len() > sampling.most + 1 || next_v.len() > sampling.most + 1 {
            break;
        }
        (us, vs) = (next_u, next_v);
    }
    let Some((mut best, read_densely)) = best else {
        ogeom_bail!(Construction, "the surface could not be sampled");
    };
    if !read_densely {
        let (mut split_u, mut split_v) = (vec![false; us.len() - 1], vec![false; vs.len() - 1]);
        let along = measured_along_lines(
            &mut at,
            &best.curve,
            (&us, &vs),
            sampling.between,
            tolerance,
            best.error,
            (&mut split_u, &mut split_v),
            tol,
        )?;
        best.error = best.error.max(along);
        best.met = best.error <= tolerance;
    }
    Ok(best)
}

/// How many points of each knot span a fit is read at once its samples'
/// own check passes ([`read_densely`]).
const DENSE_PER_SPAN: usize = 8;

/// The knot breaks of `knots` from `lo` to `hi`, both ends among them.
fn knot_breaks(knots: &KnotVector, (lo, hi): (f64, f64)) -> Vec<f64> {
    let mut breaks: Vec<f64> = knots
        .distinct()
        .into_iter()
        .map(|(k, _)| k)
        .filter(|k| (lo..=hi).contains(k))
        .collect();
    if breaks.first().is_none_or(|k| *k > lo) {
        breaks.insert(0, lo);
    }
    if breaks.last().is_none_or(|k| *k < hi) {
        breaks.push(hi);
    }
    breaks
}

/// A peak of a miss climbed along a line: its height, where it stands, and
/// the readings either side it was climbed between.
struct Peak {
    height: f64,
    at: f64,
    reach: (f64, f64),
}

/// What [`read_densely`] read along a line.
struct Dense {
    /// The widest miss read, and where.
    widest: (f64, f64),
    /// Every peak climbed.
    peaks: Vec<Peak>,
    /// Where a reading or a climbed peak stands over the tolerance.
    misses: Vec<f64>,
}

/// Golden sections of `off` between `reach`, from `top` (a height and
/// where it was read): the highest reading met, and where. Where `lean`,
/// each section reads one new point and keeps the other, `LEAN_STEPS`
/// times; otherwise it reads both afresh, `STEPS` times.
fn climb(
    off: &mut dyn FnMut(f64) -> OgeomResult<f64>,
    (mut a, mut b): (f64, f64),
    mut top: (f64, f64),
    lean: bool,
) -> OgeomResult<(f64, f64)> {
    const STEPS: usize = 40;
    const LEAN_STEPS: usize = 30;
    let ratio = (5.0_f64.sqrt() - 1.0) / 2.0;
    let mut keep = |o: f64, t: f64| {
        if o > top.0 {
            top = (o, t);
        }
    };
    if lean {
        let (mut c, mut d) = (b - (b - a) * ratio, a + (b - a) * ratio);
        let (mut oc, mut od) = (off(c)?, off(d)?);
        keep(oc, c);
        keep(od, d);
        for _ in 0..LEAN_STEPS {
            if oc > od {
                (b, d, od) = (d, c, oc);
                c = b - (b - a) * ratio;
                oc = off(c)?;
                keep(oc, c);
            } else {
                (a, c, oc) = (c, d, od);
                d = a + (b - a) * ratio;
                od = off(d)?;
                keep(od, d);
            }
        }
        return Ok(top);
    }
    for _ in 0..STEPS {
        let (c, d) = (b - (b - a) * ratio, a + (b - a) * ratio);
        let (oc, od) = (off(c)?, off(d)?);
        keep(oc, c);
        keep(od, d);
        if oc > od {
            b = d;
        } else {
            a = c;
        }
    }
    Ok(top)
}

/// A fit's miss `off` read densely along a line: at [`DENSE_PER_SPAN`]
/// points of every span between two of `breaks` (the fit's knots, which
/// crowd where the geometry turns sharply), and the readings that stand at
/// least as high as their neighbours climbed to the tops of their peaks
/// between the readings either side: the widest `always`, every other
/// within `NEAR_PEAK` of the reference (the widest reading, or `floor`,
/// the miss already measured elsewhere, if wider), and every other whose
/// neighbours' slopes could carry a crease between it and them that
/// near. Where the miss creases at a corner of the geometry the
/// readings either side of the corner can stand well below its peak, a
/// third below where the fit passes through the geometry at the knots
/// either side. `over` is the tolerance whose misses are reported.
fn read_densely(
    off: &mut dyn FnMut(f64) -> OgeomResult<f64>,
    breaks: &[f64],
    always: usize,
    floor: f64,
    over: f64,
) -> OgeomResult<Dense> {
    const NEAR_PEAK: f64 = 0.8;
    let mut ts: Vec<f64> = Vec::with_capacity(breaks.len() * DENSE_PER_SPAN);
    for w in breaks.windows(2) {
        for k in 0..DENSE_PER_SPAN {
            let f = precise(k) / precise(DENSE_PER_SPAN);
            ts.push(w[0] + (w[1] - w[0]) * f);
        }
    }
    if let Some(last) = breaks.last() {
        ts.push(*last);
    }
    let readings: Vec<f64> = ts.iter().map(|t| off(*t)).collect::<OgeomResult<_>>()?;
    let mut widest = (0.0_f64, ts.first().copied().unwrap_or(0.0));
    let mut misses = Vec::new();
    for (r, t) in readings.iter().zip(&ts) {
        if *r > widest.0 {
            widest = (*r, *t);
        }
        if *r > over {
            misses.push(*t);
        }
    }
    let mut standing: Vec<(f64, usize)> = (0..readings.len())
        .filter(|&i| {
            (i == 0 || readings[i - 1] <= readings[i])
                && readings.get(i + 1).is_none_or(|r| *r <= readings[i])
        })
        .map(|i| (readings[i], i))
        .collect();
    standing.sort_by(|a, b| b.0.total_cmp(&a.0));
    let reference = floor.max(widest.0);
    // The most a crease between two readings can rise over them, were the
    // miss as steep as anywhere beside them: half their gap times the
    // steepest slope read over the two gaps either side.
    let slope = |k: usize| (readings[k + 1] - readings[k]).abs() / (ts[k + 1] - ts[k]);
    let rise = |i: usize| -> f64 {
        let gaps = i.saturating_sub(2)..(i + 2).min(ts.len() - 1);
        let steepest = gaps.clone().map(slope).fold(0.0_f64, f64::max);
        let mut rise = 0.0_f64;
        for k in [i.checked_sub(1), Some(i)].into_iter().flatten() {
            if k + 1 < ts.len() {
                let crest = f64::midpoint(readings[k], readings[k + 1])
                    + steepest * (ts[k + 1] - ts[k]) / 2.0;
                rise = rise.max(crest);
            }
        }
        rise
    };
    let mut peaks = Vec::new();
    for (rank, &(reading, i)) in standing.iter().enumerate() {
        if rank >= always && reading.max(rise(i)) < NEAR_PEAK * reference {
            continue;
        }
        let reach = (ts[i.saturating_sub(1)], ts[(i + 1).min(ts.len() - 1)]);
        let (height, at) = climb(off, reach, (reading, ts[i]), rank >= always)?;
        if height > widest.0 {
            widest = (height, at);
        }
        if height > over {
            misses.push(at);
        }
        peaks.push(Peak { height, at, reach });
    }
    Ok(Dense {
        widest,
        peaks,
        misses,
    })
}

/// The span of `knots` that `x` lies in.
fn span_of(knots: &[f64], x: f64) -> usize {
    knots.partition_point(|k| *k <= x).clamp(1, knots.len() - 1) - 1
}

/// The worst distance between `surface` and `point` read densely along
/// lines of the sampling `params`, marking the spans that miss
/// `tolerance` in `splits`; `floor` is the miss already measured.
///
/// Along each direction known between its samples, every sample line
/// across it (and along `u`, where `v` is known between its samples too,
/// the line midway between two rows) is read as [`read_densely`] reads it,
/// over the knot spans of `surface` that way: a corner's peak can stand
/// between a span's quarter points, too narrow for them to see. The widest
/// peaks are then climbed across, between the lines either side, and along
/// again from there, and counted a hundredth over, as [`measured_between`]
/// counts its own: a ridge (a section's corner carried across a skin) can
/// run across the lines and peak between them. A miss on a sample line
/// splits its span along the line, any other the cell it lies in, as
/// [`measured_between`] does.
#[allow(
    clippy::too_many_arguments,
    reason = "the sampling, its directions, the spans to split and the miss so far"
)]
fn measured_along_lines(
    point: &mut impl FnMut(f64, f64) -> OgeomResult<Point>,
    surface: &crate::BSplineSurface,
    (us, vs): (&[f64], &[f64]),
    between: (bool, bool),
    tolerance: f64,
    floor: f64,
    (split_u, split_v): (&mut [bool], &mut [bool]),
    tol: Tolerances,
) -> OgeomResult<f64> {
    use crate::traits::Surface as _;
    const NEAR_PEAK: f64 = 0.8;
    const ACROSS_PEAKS: usize = 16;
    // Along `u` the widest peaks of every line are climbed whatever their
    // height, every section read afresh: a swept skin's rows are where its
    // sections lie, and are read the most thoroughly. Along `v` only the
    // peaks that could reach the widest are climbed.
    const ROW_PEAKS: usize = 8;
    const CLIMB_MARGIN: f64 = 0.01;
    // A point the geometry will not give is no reading.
    let mut off = |u: f64, v: f64| -> OgeomResult<f64> {
        Ok(match point(u, v) {
            Ok(p) => p.distance(surface.point_at(u, v, tol)?),
            Err(_) => 0.0,
        })
    };
    let mut error = 0.0_f64;
    let mut centres: Vec<(usize, usize)> = Vec::new();
    // Every peak climbed on a line the other direction is known across:
    // the peak, whether its line runs along `u`, the line, and the cell it
    // runs through where it is no sample line.
    let mut tops: Vec<(Peak, bool, f64, Option<usize>)> = Vec::new();
    for along_u in [true, false] {
        let (known_along, known_across) = if along_u {
            between
        } else {
            (between.1, between.0)
        };
        if !known_along {
            continue;
        }
        let (along, across) = if along_u { (us, vs) } else { (vs, us) };
        let knots = if along_u {
            surface.u_knots()
        } else {
            surface.v_knots()
        };
        let breaks = knot_breaks(knots, (along[0], along[along.len() - 1]));
        let mut lines: Vec<(f64, Option<usize>)> = across.iter().map(|x| (*x, None)).collect();
        if known_across && along_u {
            lines.extend(
                across
                    .windows(2)
                    .enumerate()
                    .map(|(j, w)| (f64::midpoint(w[0], w[1]), Some(j))),
            );
        }
        for (x, cell) in lines {
            let dense = read_densely(
                &mut |t| if along_u { off(t, x) } else { off(x, t) },
                &breaks,
                if along_u { ROW_PEAKS } else { 0 },
                floor.max(error),
                tolerance,
            )?;
            error = error.max(dense.widest.0);
            for t in dense.misses {
                let i = span_of(along, t);
                match (cell, along_u) {
                    (None, true) => split_u[i] = true,
                    (None, false) => split_v[i] = true,
                    (Some(j), true) => centres.push((i, j)),
                    (Some(j), false) => centres.push((j, i)),
                }
            }
            if known_across {
                tops.extend(dense.peaks.into_iter().map(|p| (p, along_u, x, cell)));
            }
        }
    }
    tops.sort_by(|a, b| b.0.height.total_cmp(&a.0.height));
    let near = NEAR_PEAK * floor.max(error);
    for (peak, along_u, x, cell) in tops.iter().take(ACROSS_PEAKS) {
        if peak.height < near {
            break;
        }
        let (along_u, x, t) = (*along_u, *x, peak.at);
        let across = if along_u { vs } else { us };
        // The lines either side of this one: the samples either side of a
        // midway line, the midway lines either side of a sample.
        let reach = match cell {
            Some(j) => (across[*j], across[j + 1]),
            None => {
                let j = across.partition_point(|y| *y < x);
                (
                    f64::midpoint(across[j.saturating_sub(1)], x),
                    f64::midpoint(x, across[(j + 1).min(across.len() - 1)]),
                )
            }
        };
        let (height, y) = climb(
            &mut |y| if along_u { off(t, y) } else { off(y, t) },
            reach,
            (peak.height, x),
            true,
        )?;
        let (height, s) = climb(
            &mut |s| if along_u { off(s, y) } else { off(y, s) },
            peak.reach,
            (height, t),
            true,
        )?;
        let height = height * (1.0 + CLIMB_MARGIN);
        error = error.max(height);
        if height > tolerance {
            let (u, v) = if along_u { (s, y) } else { (y, s) };
            centres.push((span_of(us, u), span_of(vs, v)));
        }
    }
    for (i, j) in centres {
        if !split_u[i] && !split_v[j] {
            split_u[i] = true;
            split_v[j] = true;
        }
    }
    Ok(error)
}

/// The worst distance between `surface` and `point` between the samples of
/// `params` (in the directions `between` says are known there), and which
/// spans miss `tolerance`.
///
/// Each span is read at its quarter points, along the sample lines and
/// across them. The widest of those readings that stand above their
/// neighbours are then each climbed to the top of their peak, within a
/// quarter span of the reading either way, and counted a hundredth over
/// the height climbed: a quarter point stands up to a tenth below the peak
/// beside it where the error swells between them, and the climb a few
/// thousandths below it at most.
/// A miss inside a cell is put down to whichever direction already misses
/// along its sample lines through that cell's spans, and to both only where
/// neither does: a skin that strays across its rows has no use for more
/// samples along them.
#[allow(clippy::type_complexity, reason = "the error and the spans to split")]
fn measured_between(
    point: &mut impl FnMut(f64, f64) -> OgeomResult<Point>,
    surface: &crate::BSplineSurface,
    (us, vs): (&[f64], &[f64]),
    between: (bool, bool),
    tolerance: f64,
    tol: Tolerances,
) -> OgeomResult<(f64, Vec<bool>, Vec<bool>)> {
    use crate::traits::Surface as _;
    // A reading this close to the widest may stand beside a peak wider
    // than every reading: at most `PEAKS` of them are climbed.
    const NEAR_PEAK: f64 = 0.8;
    const PEAKS: usize = 16;
    const CLIMB_MARGIN: f64 = 0.01;
    let checks = |knots: &[f64], between: bool| -> Vec<(f64, Option<usize>)> {
        let mut out = Vec::with_capacity(knots.len() * 2);
        for (i, pair) in knots.windows(2).enumerate() {
            out.push((pair[0], None));
            if between {
                let mid = f64::midpoint(pair[0], pair[1]);
                out.push((f64::midpoint(pair[0], mid), Some(i)));
                out.push((mid, Some(i)));
                out.push((f64::midpoint(mid, pair[1]), Some(i)));
            }
        }
        out.push((knots[knots.len() - 1], None));
        out
    };
    let u_checks = checks(us, between.0);
    let v_checks = checks(vs, between.1);
    let mut split_u = vec![false; us.len() - 1];
    let mut split_v = vec![false; vs.len() - 1];
    let mut centres: Vec<(usize, usize)> = Vec::new();
    let mut error = 0.0_f64;
    // Every reading, row by row of `v_checks`; a sample's own is not read
    // here (the fit reports it), and stands as `None`.
    let mut readings: Vec<Vec<Option<f64>>> = vec![vec![None; u_checks.len()]; v_checks.len()];
    for (row, &(v, v_span)) in v_checks.iter().enumerate() {
        for (column, &(u, u_span)) in u_checks.iter().enumerate() {
            if u_span.is_none() && v_span.is_none() {
                continue;
            }
            let off = point(u, v)?.distance(surface.point_at(u, v, tol)?);
            readings[row][column] = Some(off);
            error = error.max(off);
            if off > tolerance {
                match (u_span, v_span) {
                    (Some(i), Some(j)) => centres.push((i, j)),
                    (Some(i), None) => split_u[i] = true,
                    (None, Some(j)) => split_v[j] = true,
                    (None, None) => {}
                }
            }
        }
    }
    // The readings standing at least as high as their eight neighbours and
    // near the widest, widest first.
    let read = |row: usize, column: usize| readings[row][column].unwrap_or(0.0);
    let mut peaks: Vec<(f64, usize, usize)> = Vec::new();
    for (row, line) in readings.iter().enumerate() {
        for (column, reading) in line.iter().enumerate() {
            let Some(here) = *reading else {
                continue;
            };
            if here < error * NEAR_PEAK {
                continue;
            }
            let mut top = true;
            for (dr, dc) in [
                (-1, -1),
                (-1, 0),
                (-1, 1),
                (0, -1),
                (0, 1),
                (1, -1),
                (1, 0),
                (1, 1),
            ] {
                let (Some(r), Some(c)) =
                    (row.checked_add_signed(dr), column.checked_add_signed(dc))
                else {
                    continue;
                };
                if r < v_checks.len() && c < u_checks.len() && read(r, c) > here {
                    top = false;
                }
            }
            if top {
                peaks.push((here, row, column));
            }
        }
    }
    peaks.sort_by(|a, b| b.0.total_cmp(&a.0));
    // Climbed by golden sections along each known direction in turn, over
    // the reading's neighbourhood. A point the geometry will not give is
    // passed over: it is no reading.
    let ratio = (5.0_f64.sqrt() - 1.0) / 2.0;
    let mut off_at = |u: f64, v: f64| -> f64 {
        match (point(u, v), surface.point_at(u, v, tol)) {
            (Ok(p), Ok(q)) => p.distance(q),
            _ => 0.0,
        }
    };
    let reach = |checks: &[(f64, Option<usize>)], k: usize| -> (f64, f64) {
        (
            checks[k.saturating_sub(1)].0,
            checks[(k + 1).min(checks.len() - 1)].0,
        )
    };
    for &(_, row, column) in peaks.iter().take(PEAKS) {
        let (mut u, mut v) = (u_checks[column].0, v_checks[row].0);
        let mut widest = read(row, column);
        let (u_reach, v_reach) = (reach(&u_checks, column), reach(&v_checks, row));
        for _ in 0..3 {
            for along_u in [true, false] {
                if (along_u && !between.0) || (!along_u && !between.1) {
                    continue;
                }
                let (mut a, mut b) = if along_u { u_reach } else { v_reach };
                let mut at = |x: f64| if along_u { off_at(x, v) } else { off_at(u, x) };
                for _ in 0..24 {
                    let (c, d) = (b - (b - a) * ratio, a + (b - a) * ratio);
                    if at(c) > at(d) {
                        b = d;
                    } else {
                        a = c;
                    }
                }
                let x = f64::midpoint(a, b);
                let found = at(x);
                if found > widest {
                    widest = found;
                    if along_u {
                        u = x;
                    } else {
                        v = x;
                    }
                }
            }
        }
        // Climbed one direction at a time, a peak on a ridge running
        // across both is left up to a few thousandths under its top; a
        // hundredth over what was climbed covers it.
        let widest = widest * (1.0 + CLIMB_MARGIN);
        error = error.max(widest);
        if widest > tolerance {
            let span = |knots: &[f64], x: f64| {
                knots.partition_point(|k| *k <= x).clamp(1, knots.len() - 1) - 1
            };
            let (i, j) = (span(us, u), span(vs, v));
            match between {
                (true, true) => centres.push((i, j)),
                (true, false) => split_u[i] = true,
                (false, true) => split_v[j] = true,
                (false, false) => {}
            }
        }
    }
    for (i, j) in centres {
        if !split_u[i] && !split_v[j] {
            split_u[i] = true;
            split_v[j] = true;
        }
    }
    Ok((error, split_u, split_v))
}

/// The two-pass grid fit on parameters already assigned.
#[allow(
    clippy::too_many_arguments,
    reason = "the grid, its two parameterizations and the fit's own settings"
)]
fn fit_grid_on(
    rows: &[Vec<Point>],
    raw: &[Vec<[f64; 3]>],
    (u_params, v_params): (&[f64], &[f64]),
    degree: usize,
    tolerance: f64,
    (closed_v, interpolate): (bool, bool),
    (u_kept, v_kept): (Kept<'_>, Kept<'_>),
    tol: Tolerances,
) -> OgeomResult<Fitted<crate::BSplineSurface>> {
    use crate::traits::Surface as _;
    // Pass one: every row on one shared knot vector.
    let (u_knots, row_controls) = fit_family::<3>(
        raw,
        u_params,
        degree,
        tolerance * 0.5,
        false,
        interpolate,
        u_kept,
        false,
    )?;
    // Pass two: the columns of control points, against the v parameters.
    let k = u_knots.control_point_count();
    let control_columns: Vec<Vec<[f64; 3]>> = (0..k)
        .map(|i| row_controls.iter().map(|r| r[i]).collect())
        .collect();
    let (v_knots, column_controls) = fit_family::<3>(
        &control_columns,
        v_params,
        degree,
        tolerance * 0.5,
        closed_v,
        interpolate,
        v_kept,
        false,
    )?;
    let l = v_knots.control_point_count();

    // Assemble: `column_controls[i][j]` is the control at u-index i,
    // v-index j. The grid is row-major in u.
    let mut net: Vec<Point> = Vec::with_capacity(k * l);
    for column in &column_controls {
        for c in column {
            net.push(Point::new(c[0], c[1], c[2]));
        }
    }
    let grid = ogeom_math::ControlGrid::new(net, k, l)?;
    let surface = crate::BSplineSurface::new(u_knots, v_knots, &grid, tol)?;

    // The honest error: the surface against every input point, at the
    // grid's own parameters.
    let mut worst = 0.0_f64;
    for (j, row) in rows.iter().enumerate() {
        for (i, p) in row.iter().enumerate() {
            let at = surface.point_at(u_params[i], v_params[j], tol)?;
            worst = worst.max(at.distance(*p));
        }
    }
    Ok(Fitted {
        curve: surface,
        error: worst,
        met: worst <= tolerance,
    })
}

/// The best round a family fit reached: knots, one control row per member,
/// and the worst residual.
type FamilyRound<const D: usize> = (KnotVector, Vec<Vec<[f64; D]>>, f64);

/// Fit a family of point rows sharing parameters onto one knot vector,
/// refined against the worst residual across the whole family, parameters
/// held fixed. The `kept` interior knots, `(parameter, multiplicity)`,
/// are in the vector from the start and stay; they are not combined with
/// interpolation, whose knots follow the parameters alone. Where `within`,
/// a span is split only so each part holds a parameter strictly inside it
/// ([`refined_where_bad`]).
#[allow(
    clippy::too_many_arguments,
    reason = "the family, its parameters and the fit's terms"
)]
fn fit_family<const D: usize>(
    family: &[Vec<[f64; D]>],
    parameters: &[f64],
    degree: usize,
    tolerance: f64,
    closed: bool,
    interpolate: bool,
    kept: Kept<'_>,
    within: bool,
) -> OgeomResult<(KnotVector, Vec<Vec<[f64; D]>>)> {
    let degree = degree.min(parameters.len() - 1).max(1);
    let (a, b) = (parameters[0], parameters[parameters.len() - 1]);
    let mut knots = if interpolate {
        interpolating_knots(degree, parameters)?
    } else {
        let mut knots = single_span(degree, a, b)?;
        for &(at, multiplicity) in kept {
            if at > a && at < b {
                knots = knots.with_knot_inserted(at, multiplicity.min(degree))?;
            }
        }
        knots
    };
    let kept = !interpolate && kept.iter().any(|&(at, _)| at > a && at < b);
    const ROUNDS: usize = 24;
    let mut best: Option<FamilyRound<D>> = None;
    for _ in 0..ROUNDS {
        let mut controls = Vec::with_capacity(family.len());
        let mut merged: Vec<(f64, f64)> = parameters.iter().map(|u| (*u, 0.0)).collect();
        let mut solvable = true;
        for row in family {
            match least_squares::<D>(&knots, row, parameters, closed) {
                Ok(control) => {
                    for (slot, entry) in residuals::<D>(&knots, &control, row, parameters)
                        .iter()
                        .zip(merged.iter_mut())
                    {
                        entry.1 = entry.1.max(slot.1);
                    }
                    controls.push(control);
                }
                Err(_) => {
                    solvable = false;
                    break;
                }
            }
        }
        if !solvable {
            break;
        }
        let worst = merged.iter().fold(0.0_f64, |acc, e| acc.max(e.1));
        if best.as_ref().is_none_or(|(_, _, held)| worst < *held) {
            best = Some((knots.clone(), controls, worst));
        }
        if worst <= tolerance
            || knots.control_point_count() >= parameters.len() + usize::from(closed)
        {
            break;
        }
        let Some(refined) = refined_where_bad(&knots, &merged, tolerance, within)? else {
            break;
        };
        // A refinement past one control point per datum leaves spans with
        // no data and cannot be solved. Kept knots stay, which the averaged
        // vector of the interpolating spline would drop.
        if kept {
            if refined.control_point_count() > parameters.len() {
                break;
            }
            knots = refined;
        } else {
            knots = refined_or_interpolating(refined, parameters, closed)?;
        }
    }
    let Some((knots, controls, _)) = best else {
        ogeom_bail!(NotDone, "the family fit solved no round at all");
    };
    Ok((knots, controls))
}

/// The knots an open fit's next round stands on, given its refinement.
///
/// A refinement reaching one control point per datum is the interpolating
/// spline. On knots split at the residuals' medians that spline meets every
/// datum and can swing far from the data between them, a spline through a
/// smooth curve's samples straying by a sizeable fraction of their spacing;
/// on knots averaged from the parameters it follows the data. A closed fit
/// spends a control point on its join and keeps its refinement.
fn refined_or_interpolating(
    refined: KnotVector,
    parameters: &[f64],
    closed: bool,
) -> OgeomResult<KnotVector> {
    if closed || refined.control_point_count() < parameters.len() {
        return Ok(refined);
    }
    interpolating_knots(refined.degree(), parameters)
}

/// The knot vector a spline through every one of `parameters` stands on:
/// clamped at the ends, each interior knot the average of `degree`
/// consecutive parameters, so every span holds data.
fn interpolating_knots(degree: usize, parameters: &[f64]) -> OgeomResult<KnotVector> {
    let n = parameters.len();
    let (a, b) = (parameters[0], parameters[n - 1]);
    let mut knots = vec![a; degree + 1];
    for j in 1..n.saturating_sub(degree) {
        #[allow(clippy::cast_precision_loss, reason = "a small count")]
        let mean = parameters[j..j + degree].iter().sum::<f64>() / degree as f64;
        knots.push(mean);
    }
    knots.extend(std::iter::repeat_n(b, degree + 1));
    KnotVector::new(knots, degree)
}

/// Chord-length parameters, normalized to `[0, 1]`.
///
/// For densely sampled points of a smooth curve, chord length approximates
/// arc length whatever the sampling density does, so the parameter-to-point
/// map keeps a uniform speed across spacing jumps, exactly where the
/// centripetal assignment would fold a spacing jump into a speed kink.
fn chordal<const D: usize>(points: &[[f64; D]]) -> Vec<f64> {
    let mut out = Vec::with_capacity(points.len());
    out.push(0.0);
    let mut total = 0.0;
    for pair in points.windows(2) {
        total += distance::<D>(&pair[0], &pair[1]);
        out.push(total);
    }
    if total > 0.0 {
        for u in &mut out {
            *u /= total;
        }
    }
    if let Some(last) = out.last_mut() {
        *last = 1.0;
    }
    out
}

/// Centripetal parameters over the points, on `[0, 1]`.
fn centripetal<const D: usize>(points: &[[f64; D]]) -> Vec<f64> {
    let mut out = Vec::with_capacity(points.len());
    out.push(0.0);
    let mut total = 0.0;
    for pair in points.windows(2) {
        total += distance::<D>(&pair[0], &pair[1]).sqrt();
        out.push(total);
    }
    if total > 0.0 {
        for u in &mut out {
            *u /= total;
        }
    }
    // Exactness at the far end matters: the last parameter must be the domain
    // end, not a rounding neighbour of it.
    if let Some(last) = out.last_mut() {
        *last = 1.0;
    }
    out
}

/// A clamped knot vector with a single span.
fn single_span(degree: usize, a: f64, b: f64) -> OgeomResult<KnotVector> {
    let mut knots = Vec::with_capacity(2 * (degree + 1));
    knots.extend(core::iter::repeat_n(a, degree + 1));
    knots.extend(core::iter::repeat_n(b, degree + 1));
    KnotVector::new(knots, degree)
}

/// Least-squares control points for a fixed knot vector.
///
/// The ends are pinned to the first and last data points (they are where the
/// curve joins its neighbours), and the interior is solved. With no interior
/// there is nothing to solve and the pinned Bézier is the answer.
fn least_squares<const D: usize>(
    knots: &KnotVector,
    points: &[[f64; D]],
    parameters: &[f64],
    closed: bool,
) -> OgeomResult<Vec<[f64; D]>> {
    let n = knots.control_point_count();
    let m = points.len();
    let degree = knots.degree();

    let mut control = vec![[0.0; D]; n];
    control[0] = points[0];
    control[n - 1] = points[m - 1];
    if n <= 2 {
        return Ok(control);
    }

    // C1 across a closed loop's join: the clamped end derivatives are
    // degree/(t_{p+1}-t_1) * (P_1 - P_0) and
    // degree/(t_{n+p-1}-t_{n-1}) * (P_{n-1} - P_{n-2}). Setting them equal
    // makes P_{n-2} = P_{n-1} - r (P_1 - P_0): one unknown eliminated, the
    // constraint exact in the solve rather than patched on after. Both ends
    // are pinned, so this holds whether they coincide or sit one period
    // apart, as a chart image crossing its surface's seam does.
    let t = knots.knots();
    let ratio = if closed && n >= 4 {
        let d_start = t[degree + 1] - t[1];
        let d_end = t[n + degree - 1] - t[n - 1];
        (d_start > 0.0 && d_end > 0.0).then(|| d_end / d_start)
    } else {
        None
    };
    let eliminated = ratio.map(|_| n - 2);
    let unknown_count = match eliminated {
        Some(_) => n - 3,
        None => n - 2,
    };
    if unknown_count == 0 {
        if let (Some(r), Some(e)) = (ratio, eliminated) {
            // Only P_1 = P_{n-2} remains, and the constraint alone fixes it:
            // P_{n-2} = P_{n-1} - r (P_1 - P_0) with P_1 = P_{n-2} gives the
            // point dividing between the pinned ends.
            for d in 0..D {
                control[e][d] = r.mul_add(points[0][d], points[m - 1][d]) / (1.0 + r);
            }
            let _ = r;
        }
        return Ok(control);
    }

    // The collocation rows, with the pinned ends moved to the right-hand side.
    // An open curve's normal matrix is banded (each row touches `degree + 1`
    // consecutive unknowns), so it is held as its band. A loop's eliminated
    // column couples the far ends, and it is held whole.
    let mut normal = NormalMatrix::new(unknown_count, degree, ratio.is_some());
    let mut rhs = vec![vec![0.0; unknown_count]; D];

    let mut rows: Vec<(usize, Row)> = Vec::with_capacity(m);
    for (k, &u) in parameters.iter().enumerate() {
        let span = knots.span_unchecked(u);
        let basis = knots.basis(span, u);
        let first = span - degree;
        let mut row: Row = basis
            .iter()
            .enumerate()
            .map(|(j, b)| (first + j, *b))
            .collect();
        // Substitute the eliminated column: b_e P_e with
        // P_e = P_{n-1} + r P_0 - r P_1 becomes known shares on the pinned
        // ends and a coefficient of -r b_e on P_1.
        if let (Some(r), Some(e)) = (ratio, eliminated)
            && let Some(position) = row.iter().position(|(i, _)| *i == e)
        {
            let (_, b_e) = row.remove(position);
            row.push((usize::MAX, b_e * r));
            row.push((n - 1, b_e));
            row.push((1, -r * b_e));
        }
        rows.push((k, row));
    }

    for (k, row) in &rows {
        // The residual this row wants to explain, after the pinned ends. The
        // sentinel usize::MAX marks the eliminated column's share of P_0.
        let mut target = points[*k];
        for (index, b) in row {
            if *index == 0 || *index == usize::MAX {
                for d in 0..D {
                    target[d] -= b * points[0][d];
                }
            } else if *index == n - 1 {
                for d in 0..D {
                    target[d] -= b * points[m - 1][d];
                }
            }
        }
        let is_known = |i: usize| i == 0 || i == n - 1 || i == usize::MAX;
        for (i, bi) in row {
            if is_known(*i) {
                continue;
            }
            for (j, bj) in row {
                if is_known(*j) {
                    continue;
                }
                normal.add(i - 1, j - 1, bi * bj);
            }
            for d in 0..D {
                rhs[d][i - 1] += bi * target[d];
            }
        }
    }

    // The normal matrix can be singular when a span has no parameter in it:
    // a knot was placed where there is no data to say where the curve goes.
    let Some(solved) = normal.solve(&rhs) else {
        ogeom_bail!(
            NotDone,
            "the fitting system is singular: a knot span contains no data"
        );
    };
    for (d, column) in solved.iter().enumerate() {
        for i in 0..unknown_count {
            control[i + 1][d] = column[i];
        }
    }
    if let (Some(r), Some(e)) = (ratio, eliminated) {
        for d in 0..D {
            control[e][d] = points[m - 1][d] + r * (points[0][d] - control[1][d]);
        }
    }
    Ok(control)
}

/// A collocation row's columns and basis values: the `degree + 1` of the
/// span and, on a loop, two more for the eliminated column.
type Row = smallvec::SmallVec<[(usize, f64); 8]>;

/// A least-squares normal matrix, `BᵀB`, symmetric and positive definite
/// where the fit is well posed: held as its lower band where every row of
/// `B` touches only `band + 1` consecutive unknowns, and whole otherwise.
enum NormalMatrix {
    /// `lower[i][k]` is entry `(i, i - k)`, for `k` up to the band.
    Banded {
        lower: Vec<Vec<f64>>,
        band: usize,
    },
    Dense(nalgebra::DMatrix<f64>),
}

impl NormalMatrix {
    fn new(size: usize, band: usize, dense: bool) -> Self {
        if dense {
            Self::Dense(nalgebra::DMatrix::zeros(size, size))
        } else {
            Self::Banded {
                lower: vec![vec![0.0; band + 1]; size],
                band,
            }
        }
    }

    fn add(&mut self, i: usize, j: usize, value: f64) {
        match self {
            Self::Dense(m) => m[(i, j)] += value,
            Self::Banded { lower, band } => {
                // Symmetric: only the lower half is kept, each pair once.
                if j <= i && i - j <= *band {
                    lower[i][i - j] += value;
                }
            }
        }
    }

    /// The solution for each right-hand side, by Cholesky: a banded one in
    /// the band's own storage, costing the size times the band squared
    /// where an inverse cost the size cubed. `None` where the system is
    /// singular, or so near it that a pivot falls a trillion times below
    /// the largest diagonal entry: a knot span with (next to) no data in
    /// it, whose control points the data does not decide. Solved anyway,
    /// they come out anywhere, and the curve with them.
    fn solve(&self, rhs: &[Vec<f64>]) -> Option<Vec<Vec<f64>>> {
        match self {
            Self::Banded { lower, band } => {
                let factor = banded_cholesky(lower, *band)?;
                Some(
                    rhs.iter()
                        .map(|b| banded_substitute(&factor, *band, b))
                        .collect(),
                )
            }
            Self::Dense(m) => {
                let largest = m.diagonal().iter().fold(0.0_f64, |a, d| a.max(d.abs()));
                let factor = m.clone().cholesky()?;
                let floor = largest * SINGULAR;
                if factor.l_dirty().diagonal().iter().any(|d| d * d <= floor) {
                    return None;
                }
                Some(
                    rhs.iter()
                        .map(|b| {
                            factor
                                .solve(&nalgebra::DVector::from_column_slice(b))
                                .as_slice()
                                .to_vec()
                        })
                        .collect(),
                )
            }
        }
    }
}

/// How far below the largest diagonal entry a Cholesky pivot may fall
/// before the normal matrix counts as singular.
const SINGULAR: f64 = 1e-12;

/// The Cholesky factor of a banded symmetric matrix given by its lower
/// band, in the same storage; `None` at a non-positive pivot.
fn banded_cholesky(lower: &[Vec<f64>], band: usize) -> Option<Vec<Vec<f64>>> {
    let n = lower.len();
    let largest = lower.iter().fold(0.0_f64, |a, row| a.max(row[0].abs()));
    let floor = largest * SINGULAR;
    let mut factor = vec![vec![0.0; band + 1]; n];
    for i in 0..n {
        let first = i.saturating_sub(band);
        for j in first..=i {
            let mut sum = lower[i][i - j];
            for k in first.max(j.saturating_sub(band))..j {
                sum -= factor[i][i - k] * factor[j][j - k];
            }
            if i == j {
                if sum <= floor || !sum.is_finite() {
                    return None;
                }
                factor[i][0] = sum.sqrt();
            } else {
                factor[i][i - j] = sum / factor[j][0];
            }
        }
    }
    Some(factor)
}

/// Solve `L Lᵀ x = b` for a banded Cholesky factor `L`.
fn banded_substitute(factor: &[Vec<f64>], band: usize, b: &[f64]) -> Vec<f64> {
    let n = factor.len();
    let mut y = vec![0.0; n];
    for i in 0..n {
        let mut sum = b[i];
        for k in i.saturating_sub(band)..i {
            sum -= factor[i][i - k] * y[k];
        }
        y[i] = sum / factor[i][0];
    }
    let mut x = vec![0.0; n];
    for i in (0..n).rev() {
        let mut sum = y[i];
        for k in (i + 1)..n.min(i + band + 1) {
            sum -= factor[k][k - i] * x[k];
        }
        x[i] = sum / factor[i][0];
    }
    x
}

/// Move each parameter to the foot of the perpendicular from its point.
///
/// One Newton step per call on `g(u) = (C(u) - Q) · C'(u) = 0`. The ends stay
/// pinned: they are where the curve joins its neighbours, and letting them
/// slide would trade end accuracy for interior accuracy silently.
fn correct_parameters<const D: usize>(
    knots: &KnotVector,
    control: &[[f64; D]],
    points: &[[f64; D]],
    parameters: &mut [f64],
    looped: bool,
) {
    let (lo, hi) = knots.domain();
    let last = parameters.len() - 1;
    // On a loop, a trust region of a few sample spacings. Newton on a curve
    // that is still poor (the early rounds of a fit) can project a point
    // to the far side of the domain, and a loop's join makes that a cliff:
    // both ends of the domain are the same place in space, one wild foot
    // lands at the wrong end, and the monotonicity repair below drags every
    // parameter after it onto it. An open curve has no such identification,
    // and clamping its corrections only slows the endgame.
    #[allow(clippy::cast_precision_loss)]
    let max_step = if looped {
        (hi - lo) * 8.0 / parameters.len() as f64
    } else {
        hi - lo
    };
    for (k, u) in parameters.iter_mut().enumerate() {
        if k == 0 || k == last {
            continue;
        }
        let (at, d1, d2) = evaluate::<D>(knots, control, *u);
        let gap: [f64; D] = core::array::from_fn(|d| at[d] - points[k][d]);
        let dot = |a: &[f64; D], b: &[f64; D]| a.iter().zip(b).map(|(x, y)| x * y).sum::<f64>();
        let numerator = dot(&gap, &d1);
        let denominator = dot(&d1, &d1) + dot(&gap, &d2);
        if denominator.abs() <= f64::MIN_POSITIVE {
            continue;
        }
        let stepped = *u - (numerator / denominator).clamp(-max_step, max_step);
        if stepped.is_finite() {
            *u = stepped.clamp(lo, hi);
        }
    }
    // Projection can reorder neighbours near a tight turn. The fit assumes the
    // parameters walk forward with the points.
    for k in 1..parameters.len() {
        if parameters[k] < parameters[k - 1] {
            parameters[k] = parameters[k - 1];
        }
    }
}

/// A curve point and its first two derivatives, from raw knots and control.
fn evaluate<const D: usize>(
    knots: &KnotVector,
    control: &[[f64; D]],
    u: f64,
) -> ([f64; D], [f64; D], [f64; D]) {
    let degree = knots.degree();
    let span = knots.span_unchecked(u);
    let table = knots.basis_derivatives(span, u, 2);
    let first = span - degree;
    let mut out = [[0.0; D]; 3];
    for (order, row) in table.iter().enumerate().take(3) {
        for (j, b) in row.iter().enumerate() {
            for d in 0..D {
                out[order][d] += b * control[first + j][d];
            }
        }
    }
    (out[0], out[1], out[2])
}

/// The distance from each point to the curve at its parameter.
fn residuals<const D: usize>(
    knots: &KnotVector,
    control: &[[f64; D]],
    points: &[[f64; D]],
    parameters: &[f64],
) -> Vec<(f64, f64)> {
    parameters
        .iter()
        .zip(points)
        .map(|(&u, p)| (u, distance::<D>(&position::<D>(knots, control, u), p)))
        .collect()
}

/// The curve's point at `u`, from raw knots and control.
fn position<const D: usize>(knots: &KnotVector, control: &[[f64; D]], u: f64) -> [f64; D] {
    let span = knots.span_unchecked(u);
    let basis = knots.basis(span, u);
    let first = span - knots.degree();
    let mut at = [0.0; D];
    for (j, b) in basis.iter().enumerate() {
        for d in 0..D {
            at[d] += b * control[first + j][d];
        }
    }
    at
}

/// How far the curve strays between samples, charged to the samples.
///
/// Correction lets every point find its own foot, so a curve that loops
/// away between two samples and comes back for each still meets every
/// point, and a short trace can read within tolerance as a curve many
/// times its length. Midway between two samples' parameters the curve lies
/// within half their chord of the chord's middle, or the excess is the
/// fit's error at both samples, where refinement can act on it.
fn wandering<const D: usize>(
    knots: &KnotVector,
    control: &[[f64; D]],
    points: &[[f64; D]],
    parameters: &[f64],
    errors: &mut [(f64, f64)],
) {
    for k in 0..parameters.len().saturating_sub(1) {
        let middle = f64::midpoint(parameters[k], parameters[k + 1]);
        let (at, _, _) = evaluate::<D>(knots, control, middle);
        let chord: [f64; D] =
            core::array::from_fn(|d| f64::midpoint(points[k][d], points[k + 1][d]));
        let excess = distance::<D>(&at, &chord) - distance::<D>(&points[k], &points[k + 1]) / 2.0;
        if excess > 0.0 {
            errors[k].1 = errors[k].1.max(excess);
            errors[k + 1].1 = errors[k + 1].1.max(excess);
        }
    }
}

/// The knot vector with every offending span split.
///
/// Split at the *median parameter* inside the span, not its geometric middle.
/// The middle stalls where the data crowds into one half of a bad span: a
/// middle knot leaves the other half empty, an empty half makes the
/// least-squares system singular, and the span can never be refined again.
/// The fit then converges to just above the target and sticks there. The
/// median always leaves data on both sides.
///
/// `None` when no span can be split further: every bad span holds fewer than
/// two parameters, and a knot needs data on both sides to be supported.
/// Where `within`, a parameter on a span's own start knot is not counted,
/// so each span a split makes holds one strictly inside it.
///
/// Errors in strictly increasing parameter order, as a fit at fixed
/// parameters gives them, hold each span's samples as one run, found by
/// bisection; others are sifted span by span, in the order given.
fn refined_where_bad(
    knots: &KnotVector,
    errors: &[(f64, f64)],
    tolerance: f64,
    within: bool,
) -> OgeomResult<Option<KnotVector>> {
    let distinct = knots.distinct();
    let mut refined = knots.clone();
    let mut changed = false;
    let ordered = errors.windows(2).all(|w| w[0].0 < w[1].0);
    for window in distinct.windows(2) {
        let (lo, hi) = (window[0].0, window[1].0);
        let inside: std::borrow::Cow<'_, [(f64, f64)]> = if ordered {
            let from = errors.partition_point(|(u, _)| *u < lo);
            let to = errors.partition_point(|(u, _)| *u < hi);
            std::borrow::Cow::Borrowed(&errors[from..to])
        } else {
            errors
                .iter()
                .filter(|(u, _)| *u >= lo && *u < hi)
                .copied()
                .collect()
        };
        let bad = inside.iter().any(|(_, e)| *e > tolerance);
        let inside = if within {
            std::borrow::Cow::Owned(inside.iter().filter(|(u, _)| *u > lo).copied().collect())
        } else {
            inside
        };
        if !bad || inside.len() < 2 {
            continue;
        }
        // The knot between the two middle parameters, kept strictly interior.
        let at = f64::midpoint(inside[inside.len() / 2 - 1].0, inside[inside.len() / 2].0)
            .clamp(lo + (hi - lo) * 1e-6, hi - (hi - lo) * 1e-6);
        // Both sides must keep data, or the new span is unsupported.
        let left = inside.iter().any(|(u, _)| *u < at);
        let right = inside.iter().any(|(u, _)| *u >= at);
        if left && right {
            refined = refined.with_knot_inserted(at, 1)?;
            changed = true;
        }
    }
    Ok(if changed { Some(refined) } else { None })
}

#[cfg(test)]
#[allow(clippy::unwrap_used)]
mod grid_tests {
    use super::*;
    use crate::traits::Surface as _;
    use ogeom_core::Tolerances;

    const T: Tolerances = Tolerances::millimetres();

    /// A loop blended nine tenths of the way from a circle to a square,
    /// both run by length from their own starts, has a kink at each of the
    /// square's corners, and the fit's miss peaks between the eighth points
    /// it is read at. The error a sampled fit states is no less than its
    /// miss read two hundred thousand times round the loop.
    #[test]
    fn a_sampled_loop_fit_states_no_less_than_its_widest_miss() {
        use crate::traits::Curve3d as _;
        let square = |g: f64| -> Point {
            let s = g.rem_euclid(1.0) * 4.0;
            let r = s - s.floor();
            let (x, y) = match s.floor() {
                0.0 => (1.0, 2.0f64.mul_add(r, -1.0)),
                1.0 => (2.0f64.mul_add(-r, 1.0), 1.0),
                2.0 => (-1.0, 2.0f64.mul_add(-r, 1.0)),
                _ => (2.0f64.mul_add(r, -1.0), -1.0),
            };
            Point::new(x, y, 0.0)
        };
        let blend = |g: f64| {
            let (sin, cos) = g.mul_add(core::f64::consts::TAU, 0.3).sin_cos();
            let (c, q) = (Point::new(1.2 * cos, 1.2 * sin, 0.0), square(g));
            c + (q - c) * 0.9
        };
        let ts: Vec<f64> = (0..=64).map(|i| f64::from(i) / 64.0).collect();
        let fitted = fit_curve_sampled(|g| Ok(blend(g)), &ts, true, 3, 1e-4, T).unwrap();
        assert!(fitted.met);
        let n = 200_000;
        let mut worst = 0.0_f64;
        for i in 0..=n {
            let g = f64::from(i) / f64::from(n);
            worst = worst.max(blend(g).distance(fitted.curve.point_at(g, T).unwrap()));
        }
        assert!(
            worst <= fitted.error,
            "the fit misses by {worst:.4e} and states {:.4e}",
            fitted.error
        );
    }

    /// A closed polygon of 25 corners read by arc length: the fit strays
    /// most at the corners, in peaks too narrow for its eighth points to
    /// see and too many to climb them all. The error a sampled curve fit
    /// states is no less than its miss at every corner.
    #[test]
    fn a_sampled_curve_fit_states_no_less_than_its_miss_at_a_corner() {
        use crate::traits::Curve3d as _;
        let n = 25;
        let corners: Vec<Point> = (0..n)
            .map(|k| {
                let k = f64::from(k);
                let (sin, cos) = (k / f64::from(n) * core::f64::consts::TAU).sin_cos();
                let r = 3.0f64.mul_add((5.0 * k).sin(), 10.0);
                Point::new(r * cos, r * sin, (3.0 * k).cos())
            })
            .collect();
        let mut at_corner = vec![0.0_f64];
        for (k, c) in corners.iter().enumerate() {
            at_corner.push(at_corner[k] + c.distance(corners[(k + 1) % corners.len()]));
        }
        let round = at_corner[corners.len()];
        let polygon = |g: f64| -> Point {
            let s = round * g;
            let k = at_corner[1..]
                .partition_point(|l| *l < s)
                .min(corners.len() - 1);
            let f = (s - at_corner[k]) / (at_corner[k + 1] - at_corner[k]);
            let (a, b) = (corners[k], corners[(k + 1) % corners.len()]);
            a + (b - a) * f
        };
        let ts: Vec<f64> = (0..=16).map(|i| f64::from(i) / 16.0).collect();
        let fitted = fit_curve_sampled(|g| Ok(polygon(g)), &ts, true, 3, 1e-3, T).unwrap();
        assert!(fitted.met, "{}", fitted.error);
        let worst = at_corner[..corners.len()]
            .iter()
            .map(|l| {
                let g = l / round;
                polygon(g).distance(fitted.curve.point_at(g, T).unwrap())
            })
            .fold(0.0_f64, f64::max);
        assert!(
            worst <= fitted.error,
            "the fit misses by {worst:.4e} at a corner and states {:.4e}",
            fitted.error
        );
    }

    /// A triangle read by arc length, twisted as it rises: the surface's
    /// only creases are its corners, ridges running up the skin at fixed
    /// `u`, and the fit strays most along them, in peaks too narrow for a
    /// span's quarter points to see. The error a sampled surface fit
    /// states is no less than its miss read up every corner.
    #[test]
    fn a_sampled_surface_fit_states_no_less_than_its_miss_at_a_crease() {
        use crate::traits::Surface as _;
        let corners: [[f64; 2]; 3] = [[7.9, 8.2], [-11.2, 9.2], [-2.3, -10.0]];
        let mut at_corner = [0.0_f64; 4];
        for k in 0..3 {
            let (a, b) = (corners[k], corners[(k + 1) % 3]);
            at_corner[k + 1] = at_corner[k] + (a[0] - b[0]).hypot(a[1] - b[1]);
        }
        let round = at_corner[3];
        let twisted = |u: f64, v: f64| -> OgeomResult<Point> {
            let s = round * u;
            let k = at_corner[1..].partition_point(|l| *l < s).min(2);
            let f = (s - at_corner[k]) / (at_corner[k + 1] - at_corner[k]);
            let (a, b) = (corners[k], corners[(k + 1) % 3]);
            let (x, y) = (a[0] + (b[0] - a[0]) * f, a[1] + (b[1] - a[1]) * f);
            let (sin, cos) = (1.8 * v).sin_cos();
            Ok(Point::new(x * cos - y * sin, x * sin + y * cos, 31.0 * v))
        };
        let tolerance = 1e-3;
        let fitted = fit_surface_sampled(
            twisted,
            &Sampling::even(((0.0, 1.0), (0.0, 1.0)), 8),
            3,
            tolerance,
            T,
        )
        .unwrap();
        assert!(fitted.met, "{}", fitted.error);
        let mut worst = 0.0_f64;
        for corner in &at_corner[..3] {
            let u = corner / round;
            for j in 0..=1024 {
                let v = f64::from(j) / 1024.0;
                let off = twisted(u, v)
                    .unwrap()
                    .distance(fitted.curve.point_at(u, v, T).unwrap());
                worst = worst.max(off);
            }
        }
        assert!(
            worst <= fitted.error,
            "the fit misses by {worst:.4e} at a corner and states {:.4e}",
            fitted.error
        );
    }

    #[test]
    fn a_torus_patch_grid_fits_to_tolerance_on_and_off_the_grid() {
        let torus = ogeom_math::Torus::new(ogeom_math::Frame::WORLD, 2.0, 0.5, T).unwrap();
        let surface = crate::TorusSurface::new(torus);
        let (nu, nv) = (25, 17);
        let span_u = 1.2_f64;
        let span_v = 0.9_f64;
        let sample = |fu: f64, fv: f64| surface.point_at(span_u * fu, span_v * fv, T).unwrap();
        let mut rows = Vec::new();
        for j in 0..nv {
            let mut row = Vec::new();
            for i in 0..nu {
                row.push(sample(
                    f64::from(i) / f64::from(nu - 1),
                    f64::from(j) / f64::from(nv - 1),
                ));
            }
            rows.push(row);
        }
        let fitted = fit_surface_grid(&rows, 3, 1e-4, T).unwrap();
        assert!(fitted.met, "error {} above the target", fitted.error);

        // Off the grid too: the fit describes the surface, not just its
        // samples. The fitted chart and the torus's differ, so compare by
        // distance to the true surface rather than at matched parameters.
        let (ud, vd) = fitted.curve.domain();
        for i in 0..8 {
            for j in 0..8 {
                let u = ud.0 + (ud.1 - ud.0) * (0.07 + 0.9 * f64::from(i) / 7.0);
                let v = vd.0 + (vd.1 - vd.0) * (0.07 + 0.9 * f64::from(j) / 7.0);
                let p = fitted.curve.point_at(u, v, T).unwrap();
                let d = torus.distance_to(p);
                assert!(d < 5e-4, "off-grid deviation {d} at ({u}, {v})");
            }
        }
    }

    #[test]
    fn a_grid_fitted_at_its_own_parameters_keeps_the_chart() {
        // A torus patch sampled at its own (u, v), unevenly in u: the fit
        // must answer the torus's point at the torus's parameters, off the
        // grid as well as on it.
        let torus = ogeom_math::Torus::new(ogeom_math::Frame::WORLD, 2.0, 0.5, T).unwrap();
        let surface = crate::TorusSurface::new(torus);
        let us: Vec<f64> = (0..21)
            .map(|i| 0.3 + 1.2 * (f64::from(i) / 20.0).powi(2))
            .collect();
        let vs: Vec<f64> = (0..17).map(|j| -0.4 + 0.9 * f64::from(j) / 16.0).collect();
        let rows: Vec<Vec<Point>> = vs
            .iter()
            .map(|v| {
                us.iter()
                    .map(|u| surface.point_at(*u, *v, T).unwrap())
                    .collect()
            })
            .collect();
        let fitted = fit_surface_grid_at(&us, &vs, &rows, 3, 1e-5, T).unwrap();
        assert!(fitted.met, "error {} above the target", fitted.error);
        let ((u0, u1), (v0, v1)) = fitted.curve.domain();
        assert!((u0 - 0.3).abs() < 1e-12 && (u1 - 1.5).abs() < 1e-12);
        assert!((v0 + 0.4).abs() < 1e-12 && (v1 - 0.5).abs() < 1e-12);
        for i in 0..9 {
            for j in 0..9 {
                let u = 0.3 + 1.2 * (0.05 + 0.9 * f64::from(i) / 8.0);
                let v = -0.4 + 0.9 * (0.05 + 0.9 * f64::from(j) / 8.0);
                let gap = fitted
                    .curve
                    .point_at(u, v, T)
                    .unwrap()
                    .distance(surface.point_at(u, v, T).unwrap());
                assert!(gap < 1e-4, "chart drifts by {gap} at ({u}, {v})");
            }
        }
        assert!(fit_surface_grid_at(&vs, &us, &rows, 3, 1e-5, T).is_err());
    }

    #[test]
    fn a_grid_the_basis_can_represent_fits_to_rounding() {
        // Points from a bilinear patch: degree one in both directions.
        let corner = |x: f64, y: f64| Point::new(x, y, 0.3 * x - 0.2 * y);
        let mut rows = Vec::new();
        for j in 0..6 {
            let mut row = Vec::new();
            for i in 0..6 {
                row.push(corner(f64::from(i), 2.0 * f64::from(j)));
            }
            rows.push(row);
        }
        let fitted = fit_surface_grid(&rows, 1, 1e-9, T).unwrap();
        assert!(fitted.met, "error {} above rounding", fitted.error);
    }

    /// The largest height gap between `surface` and the graph of `height`,
    /// read on a 101 by 101 grid over the surface's domain: the samples and
    /// the points between them.
    fn height_gap(surface: &crate::BSplineSurface, height: impl Fn(f64, f64) -> f64) -> f64 {
        let ((u0, u1), (v0, v1)) = surface.domain();
        let mut worst = 0.0_f64;
        for i in 0..=100 {
            for j in 0..=100 {
                let u = (u1 - u0).mul_add(f64::from(i) / 100.0, u0);
                let v = (v1 - v0).mul_add(f64::from(j) / 100.0, v0);
                let p = surface.point_at(u, v, T).unwrap();
                worst = worst.max((p.z - height(p.x, p.y)).abs());
            }
        }
        worst
    }

    /// A grid sampled from a graph over `xs` by `ys`.
    fn graph_rows(xs: &[f64], ys: &[f64], height: impl Fn(f64, f64) -> f64) -> Vec<Vec<Point>> {
        ys.iter()
            .map(|&y| xs.iter().map(|&x| Point::new(x, y, height(x, y))).collect())
            .collect()
    }

    /// A saddle's 11 by 11 grid needs the spline through every point. Its
    /// rows curve, so their chord parameters are not even in `x`, and the
    /// fit follows the saddle between the points to the spacing's fourth
    /// power rather than reproducing it. The same spline on knots split at
    /// the residuals' medians meets every point and swings over a tenth of
    /// the spacing off the saddle near the corners.
    #[test]
    fn a_saddle_grid_fit_follows_the_saddle_between_its_points() {
        let saddle = |x: f64, y: f64| x.mul_add(x, -y * y) / 20.0;
        let side: Vec<f64> = (0..=10).map(|i| f64::from(2 * i - 10)).collect();
        let rows = graph_rows(&side, &side, saddle);
        for fit in [fit_surface_grid, fit_surface_grid_chordal] {
            let fitted = fit(&rows, 3, 1e-6, T).unwrap();
            assert!(fitted.met, "error {} at the points", fitted.error);
            let gap = height_gap(&fitted.curve, saddle);
            assert!(gap < 2e-3, "the fit leaves the saddle by {gap:.3e}");
        }
    }

    /// A polynomial of degree at most the fit's in each direction, sampled
    /// at its own parameters, is reproduced everywhere, on an even grid
    /// and an uneven one alike. A grid of straight rows evenly spaced gets
    /// even parameters from its chords, and the free grid fit reproduces
    /// its degree (1, 1) saddle too.
    #[test]
    fn a_polynomial_the_degree_holds_is_reproduced_between_the_samples() {
        let even: Vec<f64> = (0..=10).map(|i| f64::from(i) / 5.0 - 1.0).collect();
        let uneven: Vec<f64> = (0..=10)
            .map(|i| 2.0 * (f64::from(i) / 10.0).powf(1.6) - 1.0)
            .collect();
        let saddle = |x: f64, y: f64| x.mul_add(x, -y * y) / 20.0;
        let cubic = |x: f64, y: f64| (x * x).mul_add(x, -2.0 * x * y * y) + 0.5 * y;
        for (xs, ys) in [(&even, &even), (&uneven, &even), (&even, &uneven)] {
            for (degree, height) in [
                (2, &saddle as &dyn Fn(f64, f64) -> f64),
                (3, &saddle),
                (3, &cubic),
            ] {
                let rows = graph_rows(xs, ys, height);
                let fitted = fit_surface_grid_at(xs, ys, &rows, degree, 1e-9, T).unwrap();
                assert!(fitted.met, "error {} at the points", fitted.error);
                let gap = height_gap(&fitted.curve, height);
                assert!(gap < 1e-9, "degree {degree} misses by {gap:.3e}");
            }
        }
        let ruled = |x: f64, y: f64| x * y / 5.0;
        let rows = graph_rows(&even, &even, ruled);
        let fitted = fit_surface_grid(&rows, 3, 1e-9, T).unwrap();
        let gap = height_gap(&fitted.curve, ruled);
        assert!(gap < 1e-9, "the ruled saddle misses by {gap:.3e}");
    }

    #[test]
    fn a_ragged_grid_is_refused() {
        let rows = vec![
            vec![Point::ORIGIN, Point::new(1.0, 0.0, 0.0)],
            vec![Point::new(0.0, 1.0, 0.0)],
        ];
        assert!(fit_surface_grid(&rows, 2, 1e-6, T).is_err());
    }
}

#[cfg(test)]
#[allow(clippy::unwrap_used)]
mod tests {

    #[test]
    fn a_v_closed_grid_fits_a_ring_with_a_smooth_join() {
        use crate::traits::Surface as _;
        // A torus sampled as a loop of rows: each row a tube circle, the
        // rows running the major circle round and back to the start.
        let tau = core::f64::consts::TAU;
        let (major, minor) = (5.0_f64, 1.5_f64);
        let rows: Vec<Vec<Point>> = (0..=24)
            .map(|j| {
                let a = tau * f64::from(j) / 24.0;
                (0..=16)
                    .map(|i| {
                        let b = tau * f64::from(i) / 16.0;
                        let r = minor.mul_add(b.cos(), major);
                        Point::new(r * a.cos(), r * a.sin(), minor * b.sin())
                    })
                    .collect()
            })
            .collect();
        let fitted = fit_surface_grid_closed_v(&rows, 3, 5e-3, T).unwrap();
        assert!(fitted.met, "the ring fit should meet its tolerance");
        let surface = fitted.curve;

        // Exact closure: the two border control rows are equal, so the loop
        // crosses its seam with nothing to weld.
        let grid = surface.grid();
        let (k, l) = (grid.u_count(), grid.v_count());
        for i in 0..k {
            let first = grid.points()[i * l].point();
            let last = grid.points()[i * l + (l - 1)].point();
            assert!(
                first.distance(last) < 1e-12,
                "border control rows differ at u-index {i}"
            );
        }

        // And smoothly: the v-derivative agrees across the join.
        let (u_dom, v_dom) = surface.domain();
        for i in 0..5 {
            let u = u_dom.0 + (u_dom.1 - u_dom.0) * f64::from(i) / 4.0;
            let (_, dv0) = surface.d1_at(u, v_dom.0, T).unwrap();
            let (_, dv1) = surface.d1_at(u, v_dom.1, T).unwrap();
            let gap = (dv0 / dv0.magnitude() - dv1 / dv1.magnitude()).magnitude();
            assert!(gap < 1e-6, "the join kinks at u = {u}: {gap}");
        }
    }

    #[test]
    fn a_grid_that_is_not_a_loop_is_refused_by_the_closed_entry() {
        let rows: Vec<Vec<Point>> = (0..4)
            .map(|j| {
                (0..4)
                    .map(|i| Point::new(f64::from(i), f64::from(j), 0.0))
                    .collect()
            })
            .collect();
        assert!(fit_surface_grid_closed_v(&rows, 3, 1e-3, T).is_err());
    }

    use super::*;
    use crate::traits::{Curve2d as _, Curve3d as _};
    use core::f64::consts::TAU;

    /// A cubic through two samples at its ends that bulges ten units off
    /// the chord between them meets both samples exactly, and strays by
    /// its bulge less half the chord: charged to both.
    #[test]
    fn a_curve_straying_between_its_samples_is_charged_for_it() {
        let knots = KnotVector::new(vec![0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0], 3).unwrap();
        let control = [[0.0, 0.0], [0.0, 10.0], [1.0, 10.0], [1.0, 0.0]];
        let points = [[0.0, 0.0], [1.0, 0.0]];
        let parameters = [0.0, 1.0];
        let mut errors = residuals::<2>(&knots, &control, &points, &parameters);
        assert!(errors.iter().all(|e| e.1 < 1e-12));
        wandering::<2>(&knots, &control, &points, &parameters, &mut errors);
        for e in &errors {
            assert!((e.1 - 7.0).abs() < 1e-12, "{e:?}");
        }
    }

    /// Eleven points of a parabola, or of a sine at its own parameters in
    /// space or in a chart, need the spline through every point at these
    /// tolerances. The same
    /// spline on knots split at the residuals' medians meets the points and
    /// swings a fifth of the spacing (the sine: a quarter of its amplitude)
    /// away between them.
    #[test]
    fn a_spline_through_every_point_follows_the_curve_between_them() {
        let parabola: Vec<Point> = (0..=10)
            .map(|i| {
                let x = f64::from(i) - 5.0;
                Point::new(x, x * x / 20.0, 0.0)
            })
            .collect();
        let fitted = fit_points(&parabola, 3, 1e-9, T).unwrap();
        assert!(fitted.met, "error {}", fitted.error);
        let (a, b) = fitted.curve.domain();
        let mut worst = 0.0_f64;
        for i in 0..=2000 {
            let p = fitted
                .curve
                .point_at((b - a).mul_add(f64::from(i) / 2000.0, a), T)
                .unwrap();
            worst = worst.max((p.y - p.x * p.x / 20.0).abs());
        }
        assert!(worst < 5e-4, "the fit leaves the parabola by {worst:.3e}");

        let ts: Vec<f64> = (0..=10).map(|i| f64::from(i) / 10.0).collect();
        let sine = |t: f64| Point::new(t, (2.0 * t).sin(), 0.0);
        let points: Vec<Point> = ts.iter().map(|t| sine(*t)).collect();
        let fitted = fit_points_at(&ts, &points, 3, 1e-12, T).unwrap();
        assert!(fitted.met, "error {}", fitted.error);
        let mut worst = 0.0_f64;
        for i in 0..=2000 {
            let t = f64::from(i) / 2000.0;
            worst = worst.max(fitted.curve.point_at(t, T).unwrap().distance(sine(t)));
        }
        assert!(worst < 5e-4, "the fit leaves the sine by {worst:.3e}");

        let chart: Vec<Point2> = points.iter().map(|p| Point2::new(p.x, p.y)).collect();
        let fitted = fit_points_2d_at(&ts, &chart, 3, 1e-12, T).unwrap();
        assert!(fitted.met, "error {}", fitted.error);
        let mut worst = 0.0_f64;
        for i in 0..=2000 {
            let t = f64::from(i) / 2000.0;
            let p = crate::Curve2d::point_at(&fitted.curve, t, T).unwrap();
            let on = sine(t);
            worst = worst.max(p.distance(Point2::new(on.x, on.y)));
        }
        assert!(worst < 5e-4, "the chart fit leaves the sine by {worst:.3e}");
    }

    #[test]
    fn scattered_points_fit_without_a_grid() {
        // A paraboloid sampled at pseudo-random spots (a deterministic
        // congruential walk, no randomness in the test) fits within a
        // stated error.
        let mut state = 12345u64;
        let mut next = || {
            state = state
                .wrapping_mul(6364136223846793005)
                .wrapping_add(1442695040888963407);
            {
                #[allow(clippy::cast_precision_loss, reason = "31 bits fit exactly")]
                let r = (state >> 33) as f64;
                r / f64::from(u32::MAX) * 2.0 - 1.0
            }
        };
        let points: Vec<Point> = (0..200)
            .map(|_| {
                let (x, y) = (next() * 5.0, next() * 5.0);
                Point::new(x, y, 0.2 * (x * x + y * y))
            })
            .collect();
        let fitted = fit_surface_scattered(&points, 3, (8, 8), 1e-3, T).unwrap();
        assert!(
            fitted.error < 0.05,
            "the paraboloid fits to {}",
            fitted.error
        );
    }

    /// A bicubic patch sampled at its own parameters comes back itself over
    /// its own knots, the net numbered along either direction; points that
    /// cover only a corner of the domain still fit there, the fairing
    /// settling the rest, flat where the patch is a plane.
    #[test]
    fn scattered_points_at_their_parameters_fit_over_the_given_knots() {
        use crate::traits::Surface as _;
        let u_knots =
            KnotVector::new(vec![0.0, 0.0, 0.0, 0.0, 0.3, 0.6, 1.0, 1.0, 1.0, 1.0], 3).unwrap();
        let v_knots = KnotVector::clamped_uniform(3, 5).unwrap();
        let (nu, nv) = (6, 5);
        let mut net = Vec::new();
        for i in 0..nu {
            for j in 0..nv {
                let (x, y) = (precise(i), precise(j));
                net.push(Point::new(x, y, (x * 0.7).sin() + 0.1 * x * y));
            }
        }
        let grid = ogeom_math::ControlGrid::new(net, nu, nv).unwrap();
        let exact = crate::BSplineSurface::new(u_knots.clone(), v_knots.clone(), &grid, T).unwrap();
        let mut parameters = Vec::new();
        let mut points = Vec::new();
        for i in 0..=20 {
            for j in 0..=13 {
                let (u, v) = (precise(i) / 20.0, precise(j) / 13.0);
                parameters.push((u, v));
                points.push(exact.point_at(u, v, T).unwrap());
            }
        }
        let fitted =
            fit_surface_scattered_at(&points, &parameters, &u_knots, &v_knots, 0.0, T).unwrap();
        assert!(fitted.error < 1e-9, "{}", fitted.error);
        // Fewer `u` controls than `v` ones numbers the net the other way.
        let (v_wide, u_narrow) = (u_knots.clone(), v_knots.clone());
        let swapped: Vec<(f64, f64)> = parameters.iter().map(|&(u, v)| (v, u)).collect();
        let turned =
            fit_surface_scattered_at(&points, &swapped, &u_narrow, &v_wide, 0.0, T).unwrap();
        assert!(turned.error < 1e-9, "{}", turned.error);

        // A plane sampled over one corner of the domain: the fairing
        // carries the plane on over the rest.
        let corner: Vec<(f64, f64)> = (0..=8)
            .flat_map(|i| (0..=8).map(move |j| (precise(i) / 20.0, precise(j) / 20.0)))
            .collect();
        let plane: Vec<Point> = corner
            .iter()
            .map(|&(u, v)| Point::new(u, v, 0.5 * u - 0.25 * v))
            .collect();
        let fitted =
            fit_surface_scattered_at(&plane, &corner, &u_knots, &v_knots, 1e-6, T).unwrap();
        assert!(fitted.error < 1e-6, "{}", fitted.error);
        let far = fitted.curve.point_at(1.0, 1.0, T).unwrap();
        assert!((far.z - 0.25).abs() < 1e-3, "{far:?}");
        assert!(
            fit_surface_scattered_at(&plane, &corner[1..], &u_knots, &v_knots, 0.0, T).is_err()
        );
    }

    #[test]
    fn fairing_trades_closeness_for_straightness_and_says_the_price() {
        // Points on a line, kicked alternately off it: the pure fit chases
        // the noise, the faired one lays a batten through it.
        let points: Vec<Point> = (0..15)
            .map(|i| {
                let t = f64::from(i);
                let kick = if i % 2 == 0 { 0.1 } else { -0.1 };
                Point::new(t, t + kick, 0.0)
            })
            .collect();
        let chased = fit_points_faired(&points, 3, 10, 0.0, T).unwrap();
        let faired = fit_points_faired(&points, 3, 10, 50.0, T).unwrap();

        // Bending energy of the control polygon: fairing must lower it.
        let bend = |c: &BSplineCurve| -> f64 {
            let pts = c.control_points();
            (1..pts.len() - 1)
                .map(|i| {
                    let p0 = pts[i - 1].point();
                    let p1 = pts[i].point();
                    let p2 = pts[i + 1].point();
                    ((p2 - p1) - (p1 - p0)).magnitude().powi(2)
                })
                .sum()
        };
        assert!(
            bend(&faired.curve) < bend(&chased.curve) / 4.0,
            "fairing must straighten: {} vs {}",
            bend(&faired.curve),
            bend(&chased.curve)
        );
        // The price is stated: the faired error is larger (it includes the
        // tangential slip a strong weight causes), while the *geometric*
        // deviation from the underlying line stays inside the noise band.
        assert!(faired.error >= chased.error);
        use crate::traits::Curve3d as _;
        let (lo, hi) = faired.curve.domain();
        for i in 0..=32 {
            let u = lo + (hi - lo) * f64::from(i) / 32.0;
            let p = faired.curve.point_at(u, T).unwrap();
            let off_line = (p.y - p.x).abs() / core::f64::consts::SQRT_2;
            assert!(
                off_line < 0.12,
                "the batten stays in the noise band: {off_line}"
            );
        }
        // Ends interpolate exactly.
        assert!(faired.curve.point_at(lo, T).unwrap().distance(points[0]) < 1e-9);
        assert!(faired.curve.point_at(hi, T).unwrap().distance(points[14]) < 1e-9);
    }

    const T: Tolerances = Tolerances::millimetres();

    /// Distance from a point to a curve: scan, then ternary-refine.
    fn nearest(curve: &BSplineCurve, p: Point) -> f64 {
        let at = |u: f64| curve.point_at(u, T).map_or(f64::MAX, |q| p.distance(q));
        let mut best = (0.0, f64::MAX);
        for i in 0..=4000 {
            #[allow(clippy::cast_precision_loss)]
            let u = i as f64 / 4000.0;
            let d = at(u);
            if d < best.1 {
                best = (u, d);
            }
        }
        let (mut lo, mut hi) = ((best.0 - 5e-4).max(0.0), (best.0 + 5e-4).min(1.0));
        for _ in 0..100 {
            let one = lo + (hi - lo) / 3.0;
            let two = hi - (hi - lo) / 3.0;
            if at(one) < at(two) {
                hi = two;
            } else {
                lo = one;
            }
        }
        at(f64::midpoint(lo, hi)).min(best.1)
    }

    /// Points along a circle, the standard curve a spline cannot represent
    /// exactly and can approach as closely as asked.
    fn circle_points(n: usize, radius: f64) -> Vec<Point> {
        (0..=n)
            .map(|i| {
                #[allow(clippy::cast_precision_loss)]
                let a = TAU * i as f64 / n as f64;
                Point::new(radius * a.cos(), radius * a.sin(), 0.0)
            })
            .collect()
    }

    #[test]
    fn a_fit_meets_the_tolerance_it_was_asked_for_and_says_what_it_reached() {
        let points = circle_points(200, 5.0);
        for tolerance in [1e-2, 1e-4, 1e-6] {
            let fitted = fit_points(&points, 3, tolerance, T).unwrap();
            assert!(
                fitted.met,
                "target {tolerance:e} not met, got {:e}",
                fitted.error
            );
            assert!(
                fitted.error <= tolerance,
                "reported {:e} over target {tolerance:e}",
                fitted.error
            );
            // And the report is honest: measure independently, point to curve.
            // A bare scan reports its own step size: 2000 samples over a
            // circumference of thirty is an 8e-3 grid, and comparing that to a
            // 1e-6 fit measures the scan, so the scan brackets and a ternary
            // search finishes.
            let mut worst = 0.0_f64;
            for p in &points {
                worst = worst.max(nearest(&fitted.curve, *p));
            }
            assert!(
                worst <= tolerance * 1.5,
                "independent measurement found {worst:e} against {tolerance:e}"
            );
        }
    }

    #[test]
    fn a_tighter_tolerance_never_uses_fewer_control_points() {
        let points = circle_points(300, 3.0);
        let coarse = fit_points(&points, 3, 1e-2, T).unwrap();
        let fine = fit_points(&points, 3, 1e-6, T).unwrap();
        assert!(
            fine.curve.control_points().len() > coarse.curve.control_points().len(),
            "{} then {}",
            coarse.curve.control_points().len(),
            fine.curve.control_points().len()
        );
        // And the coarse one is genuinely coarse: far fewer control points
        // than input points, or the fit is interpolation in disguise.
        assert!(coarse.curve.control_points().len() < 30);
    }

    #[test]
    fn knots_go_where_the_error_is() {
        // A straight run with one tight corner. Uniform refinement would
        // spread knots evenly. Adaptive refinement must put them in the
        // corner.
        let mut points = Vec::new();
        for i in 0..=100 {
            points.push(Point::new(f64::from(i) * 0.1, 0.0, 0.0));
        }
        for i in 1..=50 {
            let a = f64::from(i) / 50.0 * core::f64::consts::FRAC_PI_2;
            points.push(Point::new(10.0 + a.sin() * 0.5, (1.0 - a.cos()) * 0.5, 0.0));
        }
        for i in 1..=100 {
            points.push(Point::new(10.5, 0.5 + f64::from(i) * 0.1, 0.0));
        }

        let fitted = fit_points(&points, 3, 1e-4, T).unwrap();
        assert!(fitted.met);

        // Knot parameters cluster around the corner, which sits at roughly
        // half way through the arc length.
        let distinct = fitted.curve.knots().distinct();
        let interior: Vec<f64> = distinct[1..distinct.len() - 1]
            .iter()
            .map(|(u, _)| *u)
            .collect();
        let near_corner = interior
            .iter()
            .filter(|u| (0.40..0.60).contains(*u))
            .count();
        assert!(
            near_corner * 2 > interior.len(),
            "only {near_corner} of {} interior knots are near the corner",
            interior.len()
        );
    }

    #[test]
    fn the_ends_are_honoured_exactly_and_a_closed_loop_stays_closed() {
        let points = circle_points(64, 2.0);
        let fitted = fit_points(&points, 3, 1e-3, T).unwrap();
        let (a, b) = fitted.curve.knots().domain();
        let start = fitted.curve.point_at(a, T).unwrap();
        let end = fitted.curve.point_at(b, T).unwrap();
        assert!(start.is_equal(points[0], T), "the start drifted");
        assert!(end.is_equal(*points.last().unwrap(), T), "the end drifted");
        assert!(start.is_equal(end, T), "the loop opened");
    }

    #[test]
    fn a_smooth_loop_closes_with_matching_tangents_at_the_join() {
        use crate::traits::Curve3d as _;
        let points = circle_points(128, 3.0);
        let fitted = fit_points_closed(&points, 3, 1e-4, T).unwrap();
        assert!(fitted.met, "target not met, reached {:e}", fitted.error);

        let curve: crate::curve::Curve = fitted.curve.clone().into();
        let (a, b) = fitted.curve.knots().domain();
        let start = curve.point_at(a, T).unwrap();
        let end = curve.point_at(b, T).unwrap();
        assert!(start.is_equal(end, T), "the loop opened");

        // The join is C1: the derivative leaving the seam equals the one
        // arriving, exactly. The constraint is solved, not approximated.
        let out = curve.d1_at(a, T).unwrap();
        let back = curve.d1_at(b, T).unwrap();
        assert!(
            (out - back).magnitude() <= 1e-9 * out.magnitude(),
            "the join creases: {out:?} vs {back:?}"
        );

        // And the open fit of the same data does *not* promise this, which is
        // why the closed fit is its own entry.
        let open = fit_points(&points, 3, 1e-4, T).unwrap();
        let ocurve: crate::curve::Curve = open.curve.into();
        let (oa, ob) = (a, b);
        let _ = (ocurve.d1_at(oa, T).unwrap(), ocurve.d1_at(ob, T).unwrap());
    }

    #[test]
    fn a_joint_closed_fit_is_smooth_in_all_seven_coordinates() {
        use crate::traits::{Curve2d as _, Curve3d as _};
        let n = 96;
        let mut points = Vec::new();
        let mut on_a = Vec::new();
        let mut on_b = Vec::new();
        for i in 0..=n {
            #[allow(clippy::cast_precision_loss)]
            let t = core::f64::consts::TAU * i as f64 / n as f64;
            points.push(Point::new(4.0 * t.cos(), 4.0 * t.sin(), 1.0));
            // Two different smooth closed chart images of the same loop.
            on_a.push(Point2::new(t.cos(), t.sin()));
            on_b.push(Point2::new(2.0 * t.sin(), t.cos() - 3.0));
        }
        let (space, pa, pb) = fit_points_joint_closed(&points, &on_a, &on_b, 3, 1e-4, T).unwrap();
        assert!(space.met, "reached {:e}", space.error);

        let curve: crate::curve::Curve = space.curve.into();
        let (lo, hi) = curve.domain();
        let out3 = curve.d1_at(lo, T).unwrap();
        let back3 = curve.d1_at(hi, T).unwrap();
        assert!((out3 - back3).magnitude() <= 1e-9 * out3.magnitude());
        for plane in [&pa, &pb] {
            let planar: crate::curve2d::PlanarCurve = plane.clone().into();
            let out2 = planar.d1_at(lo, T).unwrap();
            let back2 = planar.d1_at(hi, T).unwrap();
            assert!(
                (out2 - back2).magnitude() <= 1e-9 * out2.magnitude(),
                "a pcurve creases at its seam"
            );
        }
    }

    #[test]
    fn a_fit_that_is_not_a_loop_is_refused_by_the_closed_entry() {
        let mut points = circle_points(32, 1.0);
        points.pop();
        assert!(fit_points_closed(&points, 3, 1e-3, T).is_err());
    }

    #[test]
    fn an_impossible_target_is_reported_not_rounded_up_to_success() {
        // Five points cannot pin a curve to a picometre unless the curve
        // interpolates them, and past that adding knots buys nothing. The fit
        // must say it fell short and how far.
        let points = vec![
            Point::new(0.0, 0.0, 0.0),
            Point::new(1.0, 1.0, 0.0),
            Point::new(2.0, -1.0, 0.0),
            Point::new(3.0, 1.0, 0.0),
            Point::new(4.0, 0.0, 0.0),
        ];
        let fitted = fit_points(&points, 3, 1e-15, T).unwrap();
        // With as many control points as points it may interpolate and land
        // at rounding. Either way the flags must be consistent.
        assert_eq!(fitted.met, fitted.error <= 1e-15);
    }

    #[test]
    fn the_2d_fit_is_the_same_machinery() {
        let points: Vec<Point2> = (0..=100)
            .map(|i| {
                #[allow(clippy::cast_precision_loss)]
                let a = TAU * f64::from(i) / 100.0;
                Point2::new(3.0 * a.cos(), 3.0 * a.sin())
            })
            .collect();
        let fitted = fit_points_2d(&points, 3, 1e-4, T).unwrap();
        assert!(fitted.met);
        assert!(fitted.error <= 1e-4);
        // Spot-check a few points independently, with refinement past the
        // scan's own resolution.
        for p in points.iter().step_by(17) {
            let scan = |u: f64| {
                fitted
                    .curve
                    .point_at(u, T)
                    .map_or(f64::MAX, |q| p.distance(q))
            };
            let mut best = (0.0, f64::MAX);
            for i in 0..=2000 {
                #[allow(clippy::cast_precision_loss)]
                let u = i as f64 / 2000.0;
                let d = scan(u);
                if d < best.1 {
                    best = (u, d);
                }
            }
            let (mut lo, mut hi) = ((best.0 - 1e-3).max(0.0), (best.0 + 1e-3).min(1.0));
            for _ in 0..100 {
                let one = lo + (hi - lo) / 3.0;
                let two = hi - (hi - lo) / 3.0;
                if scan(one) < scan(two) {
                    hi = two;
                } else {
                    lo = one;
                }
            }
            let found = scan(f64::midpoint(lo, hi)).min(best.1);
            assert!(found < 2e-4, "a 2d point is {found:e} off the fit");
        }
    }

    #[test]
    fn inputs_that_describe_nothing_are_refused() {
        let p = Point::ORIGIN;
        assert!(fit_points(&[], 3, 1e-3, T).is_err());
        assert!(fit_points(&[p], 3, 1e-3, T).is_err());
        assert!(
            fit_points(&[p, p, p], 3, 1e-3, T).is_err(),
            "all duplicates"
        );
        let two = [p, Point::new(1.0, 0.0, 0.0)];
        assert!(fit_points(&two, 3, 0.0, T).is_err());
        assert!(fit_points(&two, 3, -1.0, T).is_err());
        assert!(fit_points(&two, 3, f64::NAN, T).is_err());
        assert!(fit_points(&two, 0, 1e-3, T).is_err());
        // Two points always fit: the segment between them.
        assert!(fit_points(&two, 3, 1e-9, T).unwrap().met);
    }

    /// A joint fit reproduces a straight run however unevenly it was walked.
    ///
    /// A marched section's walk starts on a window's rim with a step that
    /// halved to a micron and grows back twofold a point, so most of its
    /// samples crowd one end. Parameterised centripetally, a cubic through
    /// them is not the line they lie on, and the correction rounds do not
    /// reach the feet. Chord length makes the line exact.
    #[test]
    fn a_joint_fit_holds_a_straight_run_sampled_geometrically() {
        let mut points = Vec::new();
        let mut on_a = Vec::new();
        let mut on_b = Vec::new();
        let mut s = 0.0_f64;
        let mut step = 1e-5;
        while s < 16.0 {
            points.push(Point::new(-6.5, 9.0 - s, 8.0));
            on_a.push(Point2::new(-8.0 + s, -1.0));
            on_b.push(Point2::new(-0.2, 0.7 - s / 16.0));
            s += step;
            step = (step * 2.0).min(1.0);
        }
        points.push(Point::new(-6.5, -7.0, 8.0));
        on_a.push(Point2::new(8.0, -1.0));
        on_b.push(Point2::new(-0.2, -0.3));
        let (space, pa, pb) = fit_points_joint(&points, &on_a, &on_b, 3, 1e-6, T).unwrap();
        assert!(space.error < 1e-6, "the run is a line: {}", space.error);
        let (lo, hi) = space.curve.domain();
        for i in 0..=200 {
            let t = lo + (hi - lo) * f64::from(i) / 200.0;
            let p = space.curve.point_at(t, T).unwrap();
            assert!(
                (p.x + 6.5).abs() < 1e-6 && (p.z - 8.0).abs() < 1e-6,
                "off the line at {p:?}"
            );
            let a = pa.point_at(t, T).unwrap();
            assert!((a.y + 1.0).abs() < 1e-6 && (a.x - (-8.0 + (9.0 - p.y))).abs() < 1e-6);
            let b = pb.point_at(t, T).unwrap();
            assert!((b.x + 0.2).abs() < 1e-6);
        }
    }

    /// A planar loop weaving faster than its starting samples closes C1 and
    /// holds its tolerance everywhere, not only where it was measured.
    #[test]
    fn a_sampled_planar_loop_closes_and_holds_between_its_samples() {
        use crate::traits::Curve2d as _;
        let tau = core::f64::consts::TAU;
        let at = |t: f64| {
            let r = 0.1f64.mul_add((9.0 * tau * t).sin(), 1.0);
            Point2::new(r * (tau * t).cos(), r * (tau * t).sin())
        };
        let start: Vec<f64> = (0..=12).map(|i| f64::from(i) / 12.0).collect();
        let fitted = fit_curve_2d_sampled(|t| Ok(at(t)), &start, true, 3, 1e-6, T).unwrap();
        assert!(fitted.met, "{}", fitted.error);
        let curve = &fitted.curve;
        let (a, b) = (
            curve.point_at(0.0, T).unwrap(),
            curve.point_at(1.0, T).unwrap(),
        );
        assert!(
            a.distance(b) < 1e-12,
            "the loop is open by {}",
            a.distance(b)
        );
        let (da, db) = (
            curve.derivatives_at(0.0, 1, T).unwrap()[1],
            curve.derivatives_at(1.0, 1, T).unwrap()[1],
        );
        assert!(
            (da - db).magnitude() < 1e-6 * da.magnitude(),
            "{da:?} against {db:?}"
        );
        let mut worst = 0.0_f64;
        for i in 0..=7919 {
            let t = (f64::from(i) + 0.37) / 7920.0;
            worst = worst.max(curve.point_at(t, T).unwrap().distance(at(t)));
        }
        assert!(worst <= 1e-6, "{worst}");
    }
}