ogeom_geom/

fit.rs

//! Fitting a B-spline to points, to a stated error target.
//!
//! The missing half of fitting. Interpolation and fixed-count approximation
//! exist in `ogeom-algo`; what they cannot do is *choose*: the caller names a
//! control-point count and hopes. This module is the loop that closes that:
//! fit, measure where the fit is worst, refine the knots exactly there, and
//! repeat until the error target is met.
//!
//! `docs/PLAN.md` carries this with a warning worth repeating: 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.
//!
//! # 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;

#[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.
    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. The boolean
/// found that: pcurves claiming 1e-7 evaluated millimetres from their 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, false, 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, true, tol)
}

#[allow(clippy::type_complexity)]
fn fit_points_joint_inner(
    points: &[Point],
    on_a: &[Point2],
    on_b: &[Point2],
    degree: usize,
    tolerance: f64,
    smooth_loop: 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"
        );
    }
    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.
///
/// # 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.
///
/// # 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>> {
    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, false) {
            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)? else {
            break;
        };
        knots = refined;
    }
    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,
    })
}

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)? else {
            break;
        };
        knots = refined;
    }
    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,
    })
}

/// 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.
#[allow(clippy::type_complexity)]
/// 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. Marched
    /// sections are sampled by a walk whose step halves to a micron at a
    /// window's rim and grows back only twofold a point; the centripetal
    /// guess sits so far from the feet there that two rounds of correction
    /// never reach them, and a straight section came back twelve hundred
    /// millimetres off its own line.
    ChordLength,
}

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,
        Spacing::Centripetal,
        tol,
    )
}

fn fit_spaced<const D: usize>(
    points: &[[f64; D]],
    degree: usize,
    tolerance: f64,
    smooth_loop: bool,
    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.
    let closed =
        smooth_loop && distance::<D>(&points[0], &points[points.len() - 1]) <= tol.confusion();
    let parameters = match spacing {
        Spacing::Centripetal => centripetal::<D>(&points),
        Spacing::ChordLength => chord_length::<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)? else {
            break;
        };
        knots = refined;
    }

    #[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,
    })
}

/// 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 *interpolates the Coons
/// surface's samples* to the stated tolerance; the Coons surface itself
/// interpolates the boundaries exactly, so the fit error is the whole
/// distance between the returned patch 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)
            );
        }
    }

    let mut rows: Vec<Vec<Point>> = Vec::with_capacity(samples);
    for j in 0..samples {
        let v = precise(j) / precise(samples - 1);
        let mut row = Vec::with_capacity(samples);
        for i in 0..samples {
            let u = precise(i) / precise(samples - 1);
            // The bilinearly blended Coons point: ruled in each direction,
            // the doubly-ruled corner sheet subtracted once.
            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);
            row.push(Point::from_vector(ruled_u + ruled_v - corners));
        }
        rows.push(row);
    }
    fit_surface_grid(&rows, 3, tolerance, tol)
}

/// Centripetal parameters over the points, on `[0, 1]`.
/// Fit a rectangular grid of points with a tensor-product B-spline surface.
///
/// The deferred 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.
///
/// `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 `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>> {
    use crate::traits::Surface as _;
    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();

    // 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.
    let average = |families: &[Vec<[f64; 3]>]| -> 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 u_params = average(&raw);
    let columns: Vec<Vec<[f64; 3]>> = (0..nu)
        .map(|i| raw.iter().map(|r| r[i]).collect())
        .collect();
    let v_params = average(&columns);

    // 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)?;
    // 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,
    )?;
    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,
    })
}

/// 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 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);

fn fit_family<const D: usize>(
    family: &[Vec<[f64; D]>],
    parameters: &[f64],
    degree: usize,
    tolerance: f64,
    closed: 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 = single_span(degree, a, 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)? else {
            break;
        };
        knots = refined;
    }
    let Some((knots, controls, _)) = best else {
        ogeom_bail!(NotDone, "the family fit solved no round at all");
    };
    Ok((knots, controls))
}

/// 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
}

/// Parameters proportional to the running chord length, ending exactly at
/// one.
fn chord_length<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
}

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.
    let mut normal = nalgebra::DMatrix::<f64>::zeros(unknown_count, unknown_count);
    let mut rhs = vec![nalgebra::DVector::<f64>::zeros(unknown_count); D];

    let mut rows: Vec<(usize, Vec<(usize, f64)>)> = 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: Vec<(usize, f64)> = 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[(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(inverted) = normal.clone().try_inverse() else {
        ogeom_bail!(
            NotDone,
            "the fitting system is singular: a knot span contains no data"
        );
    };
    for d in 0..D {
        let solved = &inverted * &rhs[d];
        for i in 0..unknown_count {
            control[i + 1][d] = solved[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)
}

/// 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)> {
    let degree = knots.degree();
    parameters
        .iter()
        .zip(points)
        .map(|(&u, p)| {
            let span = knots.span_unchecked(u);
            let basis = knots.basis(span, u);
            let first = span - 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];
                }
            }
            (u, distance::<D>(&at, p))
        })
        .collect()
}

/// 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: a trace of fifty millimetres read within tolerance as a curve of
/// three metres. 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 is where a textbook puts it and it stalls in practice: when 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.
fn refined_where_bad(
    knots: &KnotVector,
    errors: &[(f64, f64)],
    tolerance: f64,
) -> OgeomResult<Option<KnotVector>> {
    let distinct = knots.distinct();
    let mut refined = knots.clone();
    let mut changed = false;
    for window in distinct.windows(2) {
        let (lo, hi) = (window[0].0, window[1].0);
        let inside: Vec<f64> = errors
            .iter()
            .filter(|(u, _)| *u >= lo && *u < hi)
            .map(|(u, _)| *u)
            .collect();
        let bad = errors
            .iter()
            .any(|(u, e)| *u >= lo && *u < hi && *e > tolerance);
        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], inside[inside.len() / 2])
            .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();

    #[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_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);
    }

    #[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:?}");
        }
    }

    #[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
        );
    }

    #[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 could not be the line they lie on, and the correction rounds
    /// never reached the feet: a straight section came back twelve hundred
    /// millimetres off its own line. 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);
        }
    }
}