ogeom_intersect/

approx.rs

//! The approximation stage: a traced branch becomes curves.
//!
//! A traced branch is a polyline with a stated chord tolerance: honest, and
//! not what anything downstream wants to hold. An edge wants a curve in space;
//! a face wants that curve in its *own parameter space*, because splitting a
//! face happens there and a curve the face cannot express is a curve it cannot
//! be split along (`docs/DATA_MODEL.md` §6).
//!
//! So one branch becomes three fits sharing one tolerance: the 3D curve, and
//! one pcurve per surface, each fitted from the samples the tracer already
//! recorded. The tracer kept the parameters on both surfaces at every point
//! precisely for this moment; re-deriving them here would be a projection per
//! point, solving again what the marcher already solved.
//!
//! # The tolerance story, stated once
//!
//! The result's tolerance is a *sum of stated parts*, not a hope: the trace
//! sits within its chord tolerance of the true intersection, and the fit sits
//! within its own reported error of the trace. Both numbers are carried, and
//! the total is what an edge built on this curve must widen its tolerance to.
//! Nothing here rounds a miss up to a hit; a fit that could not reach its
//! target says so, and the caller decides whether the looser curve is usable.
//!
//! # Seams
//!
//! A branch crossing a periodic surface's seam has parameter samples that jump
//! by a period: the pcurve polyline tears even though the curve in space is
//! smooth. The samples are unwrapped before fitting: each step is folded to
//! the nearest image, so the pcurve runs continuously past the seam and may
//! legitimately leave `[0, 2π)`. That is what a pcurve on a periodic surface
//! is; folding it back would re-tear it.

use ogeom_core::{OgeomResult, Tolerances, ogeom_bail};
use ogeom_geom::{BSpline2d, BSplineCurve, Surface, SurfaceGeometry};
use ogeom_math::Point2;

use crate::march::Traced;

/// A branch of an intersection, as curves.
#[derive(Debug, Clone, PartialEq)]
pub struct IntersectionCurve {
    /// The curve in space.
    pub curve: BSplineCurve,
    /// The same curve in the first surface's parameter space.
    pub on_a: BSpline2d,
    /// And in the second's.
    pub on_b: BSpline2d,
    /// How far the *fits* may sit from the traced polyline, in millimetres.
    ///
    /// Measured: the curve against the trace's samples, each pcurve lifted
    /// against the curve, the curve against the surfaces' crossing between
    /// the samples, and, where the surfaces meet at a shallow angle, how far
    /// along them their crossing may sit from the curve. The
    /// distance to the true intersection adds the trace's own chord
    /// tolerance on top; both are stated so an edge built on this knows
    /// what to carry.
    pub fit_error: f64,
    /// Whether every fit met the tolerance it was asked for.
    pub met: bool,
    /// Whether the branch is a closed loop.
    pub closed: bool,
}

/// Fit one traced branch to curves, within `tolerance`.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the branch has
/// fewer than two points or the tolerance is not a positive distance.
pub fn approximate_branch(
    a: &SurfaceGeometry,
    b: &SurfaceGeometry,
    branch: &Traced,
    tolerance: f64,
    tol: Tolerances,
) -> OgeomResult<IntersectionCurve> {
    if branch.points.len() < 2 {
        ogeom_bail!(
            Construction,
            "a branch of {} points is not a curve",
            branch.points.len()
        );
    }

    // Marching correction can leave consecutive samples closer than the
    // rounding it converged within, and two samples at one chord-length
    // parameter are a knot span with no data in it: the fitting system
    // reports itself singular where the real defect is the duplicate. Thin
    // them here, where the trace's own step says what "too close" means.
    let mut points: Vec<ogeom_math::Point> = Vec::with_capacity(branch.points.len());
    let mut kept_a = Vec::with_capacity(branch.on_a.len());
    let mut kept_b = Vec::with_capacity(branch.on_b.len());
    // A sample is one point seen three ways, and where the three disagree
    // it is not data: through a point where the surfaces touch, the tracer
    // can report a step's position with its neighbour's parameters, and
    // the joint fit, asked to pass through both descriptions at once,
    // stalls a thousand times above its budget at that one sample.
    let agrees = |i: usize, p: &ogeom_math::Point| -> bool {
        let limit = tolerance.max(tol.confusion());
        let (ua, va) = branch.on_a[i];
        let (ub, vb) = branch.on_b[i];
        a.point_at(ua, va, tol)
            .is_ok_and(|q| q.distance(*p) <= limit)
            && b.point_at(ub, vb, tol)
                .is_ok_and(|q| q.distance(*p) <= limit)
    };
    for (i, p) in branch.points.iter().enumerate() {
        let end = i == 0 || i + 1 == branch.points.len();
        if let Some(last) = points.last()
            && last.distance(*p) <= tol.confusion() * 10.0
            && i + 1 != branch.points.len()
        {
            continue;
        }
        if !end && !agrees(i, p) {
            continue;
        }
        points.push(*p);
        kept_a.push(branch.on_a[i]);
        kept_b.push(branch.on_b[i]);
    }
    if points.len() < 2 {
        ogeom_bail!(Construction, "a branch of coincident points is not a curve");
    }

    // One fit in seven dimensions: the curve and both parameter images
    // together. Fitted separately, each fit's parameter correction drifts
    // its parameterization independently and the three results silently stop
    // being same-parameter: a pcurve claiming 1e-7 can evaluate millimetres
    // from its own curve. Jointly, one
    // parameterization and one knot vector serve all three, and the reported
    // error bounds every coordinate.
    let unwrapped_a = unwrap_periodic(a, &kept_a, tol);
    let unwrapped_b = unwrap_periodic(b, &kept_b, tol);
    // A closed branch takes the loop-smoothing fit: the join's tangents are
    // constrained to agree in all seven coordinates, so the section curve and
    // both pcurves cross their own seam without a crease. A loop winding
    // once round a periodic surface ends its image there a period from where
    // it began, and that fit, closing only where every image does, takes it
    // as an open trace whose ends lie close together, which can stall far
    // from the trace. Where it does, the loop is fitted closed in space with
    // the images' join C1 across the seam, right where the chart runs at one
    // speed across it, as a periodic surface's does (a patch that merely
    // meets itself at its seam need not), and the closer of the two stands.
    let winds_periodically = |surface: &SurfaceGeometry, image: &[Point2]| {
        let (first, last) = (image[0], image[image.len() - 1]);
        ((last.x - first.x).abs() <= tol.parametric() || surface.is_periodic_u())
            && ((last.y - first.y).abs() <= tol.parametric() || surface.is_periodic_v())
    };
    let free = || -> OgeomResult<Joint> {
        if branch.closed() {
            let closed = ogeom_geom::fit::fit_points_joint_closed(
                &points,
                &unwrapped_a,
                &unwrapped_b,
                3,
                tolerance,
                tol,
            )?;
            if !closed.0.met
                && winds_periodically(a, &unwrapped_a)
                && winds_periodically(b, &unwrapped_b)
            {
                let winding = ogeom_geom::fit::fit_points_joint_winding(
                    &points,
                    &unwrapped_a,
                    &unwrapped_b,
                    3,
                    tolerance,
                    tol,
                )?;
                if winding.0.error < closed.0.error {
                    return Ok(winding);
                }
            }
            Ok(closed)
        } else {
            ogeom_geom::fit::fit_points_joint(
                &points,
                &unwrapped_a,
                &unwrapped_b,
                3,
                tolerance,
                tol,
            )
        }
    };
    // The walk steps by how far its chord sags, which for a cubic is far
    // finer than the fit needs: between samples a step `h` apart on a curve
    // turning at `k`, a cubic misses by about `h^4 k^3 / 384` where the
    // chord sags by `h^2 k / 8`. So the knots start where that estimate
    // places a cubic's spans, every sample measures the fit, and only a
    // fit that misses one is made the free way as well, the closer of the
    // two standing.
    let stations = cubic_stations(&points, &unwrapped_a, &unwrapped_b, tolerance);
    let seeded = |closed: bool| {
        ogeom_geom::fit::fit_points_joint_from(
            &points,
            &unwrapped_a,
            &unwrapped_b,
            &stations,
            closed,
            3,
            tolerance,
            tol,
        )
    };
    let held = |closed: bool| {
        ogeom_geom::fit::fit_points_joint_from_held(
            &points,
            &unwrapped_a,
            &unwrapped_b,
            &stations,
            closed,
            3,
            tolerance,
            tol,
        )
    };
    let meets_itself = {
        let n = points.len();
        let (pa, pb) = (
            unwrapped_a[n - 1] - unwrapped_a[0],
            unwrapped_b[n - 1] - unwrapped_b[0],
        );
        let gap = (points[n - 1] - points[0]).square_magnitude()
            + pa.square_magnitude()
            + pb.square_magnitude();
        gap.sqrt() <= tol.confusion()
    };
    let mut fitted = seeded(branch.closed() && meets_itself).ok();
    if branch.closed()
        && !meets_itself
        && fitted.as_ref().is_none_or(|f| f.0.error > tolerance)
        && winds_periodically(a, &unwrapped_a)
        && winds_periodically(b, &unwrapped_b)
        && let Ok(winding) = seeded(true)
        && fitted.as_ref().is_none_or(|f| winding.0.error < f.0.error)
    {
        fitted = Some(winding);
    }
    let mut tried_free = true;
    let mut chosen = match fitted {
        Some(fitted) if fitted.0.error <= tolerance => {
            tried_free = false;
            fitted
        }
        Some(fitted) => {
            let other = free()?;
            if other.0.error < fitted.0.error {
                other
            } else {
                fitted
            }
        }
        None => free()?,
    };

    let measured =
        |(space, on_a, on_b): &Joint| measured_error(a, b, (space, on_a, on_b), &points, tol);
    let mut fit_error = measured(&chosen);
    // A fit that misses is fitted again from the stations, splitting only
    // spans that hold samples inside them
    // ([`ogeom_geom::fit::fit_points_joint_from_held`]): where the first
    // swung between two samples, this one stays near them. Measured the
    // same way, the closer stands.
    if fit_error > tolerance
        && let Ok(held) = held(branch.closed() && meets_itself)
    {
        let again = measured(&held);
        if again < fit_error {
            (chosen, fit_error) = (held, again);
        }
    }
    // A seeded fit within its tolerance at every sample stands without the
    // free one being made, but it can pass every sample and still swing
    // between two of them. Where it misses, measured, by more than a caller
    // takes whole (as beside a pole, below), the free fit is made and
    // measured the same way, and the closer stands.
    if fit_error > tolerance * 1e2
        && !tried_free
        && let Ok(other) = free()
    {
        let again = measured(&other);
        if again < fit_error {
            (chosen, fit_error) = (other, again);
        }
    }
    // Beside a pole of either chart the joint fit can miss by the chart's
    // swing: there the branch is fitted in space alone, where it is smooth,
    // and each image read off that curve ([`projected`]). Measured the same
    // way, the closer stands. A miss within a hundred times the tolerance
    // is one a caller takes whole (`fitted_in_pieces` in `section.rs`),
    // and stands.
    if fit_error > tolerance * 1e2
        && (!crate::march::chart_poles(a, tol).is_empty()
            || !crate::march::chart_poles(b, tol).is_empty())
        && let Some(projected) = projected(
            (a, b),
            &points,
            (&unwrapped_a, &unwrapped_b),
            branch.closed(),
            tolerance,
            tol,
        )
    {
        let again = measured(&projected.0).max(projected.1);
        if again < fit_error {
            (chosen, fit_error) = (projected.0, again);
        }
    }
    // An image fitted up to a pole of its chart can overshoot the pole's
    // bound between its samples, where the chart has no point. Its control
    // points are held to the bound, so the image is (a spline lies in its
    // control points' hull), and the fit measured again as it then stands.
    if let (Some(on_a), Some(on_b)) = (
        held_off_poles(a, &chosen.1, tol),
        held_off_poles(b, &chosen.2, tol),
    ) && (on_a != chosen.1 || on_b != chosen.2)
    {
        chosen = (chosen.0, on_a, on_b);
        fit_error = measured(&chosen);
    }
    let (space, on_a, on_b) = chosen;
    Ok(IntersectionCurve {
        fit_error,
        met: space.error <= tolerance,
        curve: space.curve,
        on_a,
        on_b,
        closed: branch.closed(),
    })
}

/// The branch through `points` fitted in space alone, at their chord
/// lengths, and each surface's image fitted to the curve's own points
/// inverted on it, at the curve's parameters, with how far either image
/// lifted stands from the curve where its fit read it; `None` where
/// either fit or an inversion fails.
///
/// Beside a pole of a chart the section turns gently in space while its
/// image swings: half a turn of longitude within a few times its distance
/// from the pole, which may be a few samples of the trace or none. Fitted
/// together, the image's swing holds the curve's knots and the curve the
/// image's, and the fit can stall far above its tolerance. Fitted apart,
/// the curve needs only the trace, and the image is read where the fit
/// asks for it ([`ogeom_geom::fit::fit_curve_2d_sampled`]), between the
/// samples as well as at them, each point inverted from the nearest one
/// already read and unfolded against it. The image is held, in its chart,
/// to `tolerance` over the surface's speed along the trace, and lifted it
/// is measured against the curve at every parameter the fit read. It is
/// fitted open, its ends pinned where the curve's are inverted: a loop's
/// image returns to its start or, winding round the chart, a period on.
fn projected(
    (a, b): (&SurfaceGeometry, &SurfaceGeometry),
    points: &[ogeom_math::Point],
    (image_a, image_b): (&[Point2], &[Point2]),
    closed: bool,
    tolerance: f64,
    tol: Tolerances,
) -> Option<(Joint, f64)> {
    use ogeom_geom::{Curve2d as _, Curve3d as _};
    // Halving the step past the mantissa's bits leaves it where it started.
    const HALVINGS: u32 = f64::MANTISSA_DIGITS;
    // As `lift_error` climbs its peaks.
    const PEAKS: usize = 16;
    const CLIMB: usize = 16;
    let mut along = Vec::with_capacity(points.len());
    let mut length = 0.0;
    along.push(length);
    for pair in points.windows(2) {
        let step = pair[0].distance(pair[1]);
        if step <= 0.0 {
            return None;
        }
        length += step;
        along.push(length);
    }
    let space = if closed {
        ogeom_geom::fit::fit_points_at_closed(&along, points, 3, tolerance, tol)
    } else {
        ogeom_geom::fit::fit_points_at(&along, points, 3, tolerance, tol)
    }
    .ok()?;
    let image = |surface: &SurfaceGeometry, seeds: &[Point2]| -> Option<(BSpline2d, f64)> {
        let ((ua, ub), (va, vb)) = surface.domain();
        let periods = (
            (surface.is_periodic_u() || surface.is_closed_u(tol)).then_some(ub - ua),
            (surface.is_periodic_v() || surface.is_closed_v(tol)).then_some(vb - va),
        );
        let fold = |near: f64, x: f64, period: Option<f64>| {
            period.map_or(x, |period| x - ((x - near) / period).round() * period)
        };
        let speed = seeds
            .iter()
            .filter_map(|uv| surface.d1_at(uv.x, uv.y, tol).ok())
            .map(|(du, dv)| du.magnitude() + dv.magnitude())
            .fold(0.0_f64, f64::max);
        if !(speed.is_finite() && speed > 0.0) {
            return None;
        }
        // Every image read, by parameter: the parameters run from zero, so
        // their bits sort as they do.
        let mut read: std::collections::BTreeMap<u64, Point2> = along
            .iter()
            .map(|t| t.to_bits())
            .zip(seeds.iter().copied())
            .collect();
        // A foot farther off than the curve may stand from the surface is
        // Newton's step thrown past where the chart turns (onto a pole's
        // bound, its longitude stuck): the point is reached in halves from
        // the one read before it instead.
        let limit = tolerance.max(space.error);
        let foot_from = |t: f64, near: (f64, Point2)| -> OgeomResult<Point2> {
            let mut from = near;
            let mut to = t;
            let mut halvings = 0;
            loop {
                let p = space.curve.point_at(to, tol)?;
                let foot =
                    ogeom_geom::project_on_surface_from(surface, p, (from.1.x, from.1.y), tol)?;
                if foot.distance <= limit {
                    let uv = Point2::new(
                        fold(from.1.x, foot.parameters.0, periods.0),
                        fold(from.1.y, foot.parameters.1, periods.1),
                    );
                    if to == t {
                        return Ok(uv);
                    }
                    from = (to, uv);
                    to = t;
                    continue;
                }
                halvings += 1;
                if halvings > HALVINGS {
                    ogeom_bail!(NotDone, "a point of the curve inverts on no foot near it");
                }
                to = f64::midpoint(from.0, to);
            }
        };
        let invert = |t: f64| -> OgeomResult<Point2> {
            let key = t.to_bits();
            let before = read
                .range(..key)
                .next_back()
                .map(|(s, uv)| (f64::from_bits(*s), *uv));
            let after = read
                .range(key..)
                .next()
                .map(|(s, uv)| (f64::from_bits(*s), *uv));
            let near = match (before, after) {
                (Some(before), Some(after)) => {
                    if t - before.0 < after.0 - t {
                        before
                    } else {
                        after
                    }
                }
                (Some(near), None) | (None, Some(near)) => near,
                (None, None) => ogeom_bail!(Construction, "no image read to start from"),
            };
            let uv = foot_from(t, near)?;
            read.entry(key).or_insert(uv);
            Ok(uv)
        };
        let fitted =
            ogeom_geom::fit::fit_curve_2d_sampled(invert, &along, false, 3, tolerance / speed, tol)
                .ok()?;
        // The lifted image against the curve at every parameter read and
        // midway between each two, and the widest of those misses climbed
        // to the top of its peak between its neighbours.
        let gap = |t: f64| -> Option<f64> {
            let uv = fitted.curve.point_at(t, tol).ok()?;
            let lifted = surface.point_at(uv.x, uv.y, tol).ok()?;
            Some(lifted.distance(space.curve.point_at(t, tol).ok()?))
        };
        let reads: Vec<f64> = read.keys().map(|&s| f64::from_bits(s)).collect();
        let mut stations = Vec::with_capacity(2 * reads.len());
        for pair in reads.windows(2) {
            stations.push(pair[0]);
            stations.push(f64::midpoint(pair[0], pair[1]));
        }
        stations.extend(reads.last());
        let misses = stations
            .iter()
            .map(|&t| gap(t))
            .collect::<Option<Vec<f64>>>()?;
        let mut off = misses.iter().fold(0.0_f64, |m, &x| m.max(x));
        let last = stations.len() - 1;
        let mut peaks: Vec<usize> = (0..=last)
            .filter(|&k| {
                (k == 0 || misses[k] >= misses[k - 1]) && (k == last || misses[k] >= misses[k + 1])
            })
            .collect();
        peaks.sort_by(|&i, &j| misses[j].total_cmp(&misses[i]));
        let ratio = (5.0_f64.sqrt() - 1.0) / 2.0;
        for k in peaks.into_iter().take(PEAKS) {
            let (mut left, mut right) =
                (stations[k.saturating_sub(1)], stations[(k + 1).min(last)]);
            for _ in 0..CLIMB {
                let (c, d) = (
                    right - (right - left) * ratio,
                    left + (right - left) * ratio,
                );
                let (rc, rd) = (gap(c)?, gap(d)?);
                off = off.max(rc).max(rd);
                if rc > rd {
                    right = d;
                } else {
                    left = c;
                }
            }
        }
        Some((fitted.curve, off))
    };
    let (on_a, off_a) = image(a, image_a)?;
    let (on_b, off_b) = image(b, image_b)?;
    Some(((space, on_a, on_b), off_a.max(off_b)))
}

/// `image` with each control point held inside the bound of every pole of
/// `surface`'s chart (see [`crate::march::chart_poles`]) it stands past;
/// `None` where a held control point is no control point.
fn held_off_poles(
    surface: &SurfaceGeometry,
    image: &BSpline2d,
    tol: Tolerances,
) -> Option<BSpline2d> {
    let poles = crate::march::chart_poles(surface, tol);
    if poles.is_empty() {
        return Some(image.clone());
    }
    let ((ua, ub), (va, vb)) = surface.domain();
    let mut changed = false;
    let mut control = Vec::with_capacity(image.control_points().len());
    for held in image.control_points() {
        let mut point = held.point();
        for &(k, bound, _) in &poles {
            let (coordinate, lo, hi) = if k == 0 {
                (&mut point.x, ua, ub)
            } else {
                (&mut point.y, va, vb)
            };
            let past = if (bound - hi).abs() <= (bound - lo).abs() {
                *coordinate > bound
            } else {
                *coordinate < bound
            };
            if past {
                *coordinate = bound;
                changed = true;
            }
        }
        control.push(ogeom_math::Weighted::new(point, held.weight, tol).ok()?);
    }
    if !changed {
        return Some(image.clone());
    }
    if image.is_rational() {
        BSpline2d::rational(image.knots().clone(), control).ok()
    } else {
        let points = control
            .into_iter()
            .map(ogeom_math::Weighted::point)
            .collect();
        BSpline2d::new(image.knots().clone(), points, tol).ok()
    }
}

/// How far a joint fit stands from the section it follows, in
/// millimetres: what [`IntersectionCurve::fit_error`] states.
///
/// The joint fit's own residual mixes space and chart coordinates and
/// says little about either alone, so it serves only as a bound on what
/// is measured here.
fn measured_error(
    a: &SurfaceGeometry,
    b: &SurfaceGeometry,
    (space, on_a, on_b): (
        &ogeom_geom::fit::Fitted<BSplineCurve>,
        &BSpline2d,
        &BSpline2d,
    ),
    points: &[ogeom_math::Point],
    tol: Tolerances,
) -> f64 {
    let lifted = lift_error(a, b, on_a, on_b, &space.curve, tol);
    // Where the surfaces meet at a shallow angle, a curve microns off both
    // can still sit far along them from where they cross. That distance is
    // the surfaces' gap over the sine of their angle, and it is no more than
    // how far the lifted pcurves stand from the crossing itself, found at
    // the same stations. Where a station's crossing could not be found, the
    // fit's chart residual carried through each surface's stretch bounds it
    // at the samples. The smaller of the two stands.
    let crossing = if lifted.unsettled {
        let charted =
            space_error(a, on_a, space.error, tol).max(space_error(b, on_b, space.error, tol));
        lifted.lifted_off.max(charted.max(space.error))
    } else {
        lifted.lifted_off
    };
    let lifted = lifted
        .along
        .max(lifted.curve_off)
        .max(lifted.across.min(crossing));
    // Each sample's distance from the curve, which the fit's residual
    // already bounds: read only where that bound would decide the error.
    if space.error <= lifted {
        lifted
    } else {
        lifted.max(trace_error(&space.curve, points, space.error, tol))
    }
}

/// A joint fit: the curve, and its images on the two surfaces.
type Joint = (ogeom_geom::fit::Fitted<BSplineCurve>, BSpline2d, BSpline2d);

/// The samples a cubic fit's knots start at, by index, the first and last
/// among them: each stretch between two holds a cubic's estimated miss,
/// `L Θ³ / 384` for a stretch `L` long turning through `Θ`, to an eighth
/// of `tolerance`, in space and in either chart, the turn read off the
/// trace's own polylines.
fn cubic_stations(
    points: &[ogeom_math::Point],
    image_a: &[Point2],
    image_b: &[Point2],
    tolerance: f64,
) -> Vec<usize> {
    let n = points.len();
    let target = tolerance / 8.0;
    let traces: [Vec<[f64; 3]>; 3] = [
        points.iter().map(|p| [p.x, p.y, p.z]).collect(),
        image_a.iter().map(|p| [p.x, p.y, 0.0]).collect(),
        image_b.iter().map(|p| [p.x, p.y, 0.0]).collect(),
    ];
    // Per trace: each segment's length, and the turn at each sample
    // between two segments.
    let shape: Vec<(Vec<f64>, Vec<f64>)> = traces
        .iter()
        .map(|trace| {
            let segment =
                |k: usize| -> [f64; 3] { core::array::from_fn(|d| trace[k + 1][d] - trace[k][d]) };
            let dot = |u: [f64; 3], v: [f64; 3]| u[0] * v[0] + u[1] * v[1] + u[2] * v[2];
            let lengths: Vec<f64> = (0..n - 1)
                .map(|k| dot(segment(k), segment(k)).sqrt())
                .collect();
            let mut turns = vec![0.0; n];
            for k in 1..n - 1 {
                let scale = lengths[k - 1] * lengths[k];
                if scale > 0.0 {
                    turns[k] = (dot(segment(k - 1), segment(k)) / scale)
                        .clamp(-1.0, 1.0)
                        .acos();
                }
            }
            (lengths, turns)
        })
        .collect();
    let mut out = vec![0];
    let mut from = 0;
    while from + 1 < n {
        let mut sums: Vec<(f64, f64)> = shape
            .iter()
            .map(|(lengths, _)| (lengths[from], 0.0))
            .collect();
        let mut to = from + 1;
        while to + 1 < n {
            let grown: Vec<(f64, f64)> = sums
                .iter()
                .zip(&shape)
                .map(|(&(length, turn), (lengths, turns))| (length + lengths[to], turn + turns[to]))
                .collect();
            if grown
                .iter()
                .any(|&(length, turn)| length * turn.powi(3) / 384.0 > target)
            {
                break;
            }
            sums = grown;
            to += 1;
        }
        out.push(to);
        from = to;
    }
    out
}

/// What lifting the pcurves onto their surfaces shows.
struct Lifted {
    /// How far either pcurve, lifted, stands from the curve at the same
    /// parameter.
    along: f64,
    /// That gap over the sine of the surfaces' angle there: how far the
    /// surfaces' true crossing may sit from the curve.
    across: f64,
    /// How far the curve stands from the surfaces' crossing in its normal
    /// plane.
    curve_off: f64,
    /// How far either pcurve, lifted, stands from that crossing.
    lifted_off: f64,
    /// Whether the crossing went unfound at some station.
    unsettled: bool,
}

/// Each pcurve lifted through its surface against the curve, at four
/// stations in every span of the curve's knots and at least two hundred
/// along it, between the trace's samples as well as at them. Near a cone's
/// apex or a sphere's pole the chart turns fast, and a fit that holds at
/// every sample can wander between them.
///
/// At each station the surfaces' crossing in the curve's normal plane is
/// found by Newton's method from the pcurves' points, and the curve and
/// both lifted pcurves are measured against it: what the fit misses
/// between the samples of the trace, read on the geometry the trace
/// follows. That miss runs smoothly along each span of the fit, so the
/// stations near its widest are climbed to the top of the peak beside
/// them.
fn lift_error(
    a: &SurfaceGeometry,
    b: &SurfaceGeometry,
    on_a: &BSpline2d,
    on_b: &BSpline2d,
    curve: &BSplineCurve,
    tol: Tolerances,
) -> Lifted {
    // A station stands up to a fifth below the peak beside it where the
    // miss swells between two of them; each such peak, widest first, is
    // climbed by golden sections between its neighbours.
    const NEAR_PEAK: f64 = 0.8;
    const PEAKS: usize = 16;
    const CLIMB: usize = 16;
    let (lo, hi) = curve.knots().domain();
    let spans = [curve.knots(), on_a.knots(), on_b.knots()]
        .into_iter()
        .map(|knots| knots.distinct().len().saturating_sub(1))
        .max()
        .unwrap_or(0)
        .max(1);
    let stations = (4 * spans).max(200);
    #[allow(clippy::cast_precision_loss)]
    let step = (hi - lo) / stations as f64;
    let mut out = Lifted {
        along: 0.0,
        across: 0.0,
        curve_off: 0.0,
        lifted_off: 0.0,
        unsettled: false,
    };
    #[allow(clippy::cast_precision_loss)]
    let at = |k: usize| lo + (hi - lo) * k as f64 / stations as f64;
    let read = |t: f64| station(a, b, on_a, on_b, curve, t, step, tol);
    // Each station's miss against the crossing, zero where none was read.
    let mut misses = Vec::with_capacity(stations + 1);
    for k in 0..=stations {
        let Some(here) = read(at(k)) else {
            misses.push(0.0);
            continue;
        };
        out.along = out.along.max(here.gap);
        out.across = out.across.max(here.across);
        if let Some((curve_off, lifted_off)) = here.off {
            out.curve_off = out.curve_off.max(curve_off);
            out.lifted_off = out.lifted_off.max(lifted_off);
            misses.push(curve_off.max(lifted_off));
        } else {
            out.unsettled = true;
            misses.push(0.0);
        }
    }
    let widest = misses.iter().fold(0.0_f64, |m, &x| m.max(x));
    let mut peaks: Vec<usize> = (0..=stations)
        .filter(|&k| {
            misses[k] > 0.0
                && misses[k] >= widest * NEAR_PEAK
                && (k == 0 || misses[k] >= misses[k - 1])
                && (k == stations || misses[k] >= misses[k + 1])
        })
        .collect();
    peaks.sort_by(|&i, &j| misses[j].total_cmp(&misses[i]));
    let ratio = (5.0_f64.sqrt() - 1.0) / 2.0;
    let mut climb = |t: f64| {
        read(t)
            .and_then(|s| s.off)
            .map_or(0.0, |(curve_off, lifted_off)| {
                out.curve_off = out.curve_off.max(curve_off);
                out.lifted_off = out.lifted_off.max(lifted_off);
                curve_off.max(lifted_off)
            })
    };
    for k in peaks.into_iter().take(PEAKS) {
        let (mut left, mut right) = (at(k.saturating_sub(1)), at((k + 1).min(stations)));
        let mut c = right - (right - left) * ratio;
        let mut d = left + (right - left) * ratio;
        let (mut rc, mut rd) = (climb(c), climb(d));
        for _ in 0..CLIMB {
            if rc > rd {
                right = d;
                (d, rd) = (c, rc);
                c = right - (right - left) * ratio;
                rc = climb(c);
            } else {
                left = c;
                (c, rc) = (d, rd);
                d = left + (right - left) * ratio;
                rd = climb(d);
            }
        }
    }
    out
}

/// What the curve reads at one parameter.
struct Station {
    /// How far either pcurve, lifted, stands from the curve.
    gap: f64,
    /// That gap over the sine of the surfaces' angle, zero where either
    /// normal is not determined.
    across: f64,
    /// How far the curve and either lifted pcurve stand from the surfaces'
    /// crossing in the curve's normal plane, `None` where it is not found.
    off: Option<(f64, f64)>,
}

/// A surface's chart point, its point there and first derivatives.
type Jet = (
    Point2,
    ogeom_math::Point,
    ogeom_math::Vector,
    ogeom_math::Vector,
);

/// The curve and both pcurves at `t`, measured against each other and
/// against the surfaces' crossing in the curve's normal plane there.
/// `None` where the curve cannot be evaluated. A crossing farther off
/// than `step` along the curve is another stretch of the crossing met in
/// this normal plane, not the one the curve follows, and is not read.
#[allow(clippy::too_many_arguments)]
fn station(
    a: &SurfaceGeometry,
    b: &SurfaceGeometry,
    on_a: &BSpline2d,
    on_b: &BSpline2d,
    curve: &BSplineCurve,
    t: f64,
    step: f64,
    tol: Tolerances,
) -> Option<Station> {
    use ogeom_geom::{Curve2d as _, Curve3d as _};
    let on = curve.point_at(t, tol).ok()?;
    let mut gap = 0.0_f64;
    let mut normals = [None; 2];
    let mut jets: [Option<Jet>; 2] = [None; 2];
    for (k, (surface, pcurve)) in [(a, on_a), (b, on_b)].into_iter().enumerate() {
        let Ok(chart) = pcurve.point_at(t, tol) else {
            continue;
        };
        let Ok((lifted, du, dv)) = surface.point_d1_at(chart.x, chart.y, tol) else {
            continue;
        };
        gap = gap.max(lifted.distance(on));
        jets[k] = Some((chart, lifted, du, dv));
        // The normal as the surface states it: none where the tangents
        // collapse against their own length.
        let scale = du.magnitude().max(dv.magnitude());
        let cross = du.cross(dv);
        let length = cross.magnitude();
        if length > tol.angular() * scale * scale {
            normals[k] = Some(cross * (1.0 / length));
        }
    }
    let across = match normals {
        [Some(na), Some(nb)] => gap / na.cross(nb).magnitude().max(tol.angular()),
        _ => 0.0,
    };
    let off = match (jets, curve.d1_at(t, tol)) {
        ([Some(ja), Some(jb)], Ok(tangent)) => {
            let reach = tangent.magnitude() * step;
            crossing(a, b, ja, jb, on, tangent, reach, tol)
                .filter(|found| found.distance(on) <= reach)
                .map(|found| {
                    (
                        on.distance(found),
                        ja.1.distance(found).max(jb.1.distance(found)),
                    )
                })
        }
        _ => None,
    };
    Some(Station { gap, across, off })
}

/// Where the surfaces cross in the plane through `anchor` normal to
/// `along`, by damped Newton steps from the two jets given, each step
/// first cut to carry neither surface's point past `reach` and then
/// halved until the surfaces draw closer. From the curve's own station
/// the crossing is a few steps off where there is one; where a step finds
/// no closer point within `HALVINGS`, or the steps do not settle within
/// `MOST`, no crossing is near and the answer is `None`.
#[allow(clippy::too_many_arguments)]
fn crossing(
    a: &SurfaceGeometry,
    b: &SurfaceGeometry,
    on_a: Jet,
    on_b: Jet,
    anchor: ogeom_math::Point,
    along: ogeom_math::Vector,
    reach: f64,
    tol: Tolerances,
) -> Option<ogeom_math::Point> {
    const MOST: usize = 8;
    const HALVINGS: usize = 8;
    let settled = tol.confusion() * 0.01;
    let residual_of = |pa: ogeom_math::Point, pb: ogeom_math::Point| {
        let gap = pa - pb;
        [gap.x, gap.y, gap.z, (pa - anchor).dot(along)]
    };
    let jacobian_of = |(au, av): (ogeom_math::Vector, ogeom_math::Vector),
                       (bu, bv): (ogeom_math::Vector, ogeom_math::Vector)| {
        [
            [au.x, av.x, -bu.x, -bv.x],
            [au.y, av.y, -bu.y, -bv.y],
            [au.z, av.z, -bu.z, -bv.z],
            [au.dot(along), av.dot(along), 0.0, 0.0],
        ]
    };
    let norm = |r: &[f64; 4]| r.iter().map(|v| v * v).sum::<f64>().sqrt();
    let mut x = [on_a.0.x, on_a.0.y, on_b.0.x, on_b.0.y];
    let mut point = on_a.1;
    let mut residual = residual_of(on_a.1, on_b.1);
    let mut jacobian = jacobian_of((on_a.2, on_a.3), (on_b.2, on_b.3));
    let mut size = norm(&residual);
    for _ in 0..MOST {
        if size <= settled {
            break;
        }
        let delta = solve4(jacobian, residual)?;
        // How far the step would carry each surface's point, to first
        // order. A crossing past `reach` is not the one sought, so the
        // halving starts at the first scale that stays within it.
        let carried = (0..3)
            .map(|row| {
                let on_a = jacobian[row][0] * delta[0] + jacobian[row][1] * delta[1];
                let on_b = jacobian[row][2] * delta[2] + jacobian[row][3] * delta[3];
                (on_a * on_a, on_b * on_b)
            })
            .fold((0.0, 0.0), |(sa, sb), (a, b)| (sa + a, sb + b));
        let carried = carried.0.max(carried.1).sqrt();
        let mut scale = 1.0;
        while scale * carried > reach {
            scale *= 0.5;
        }
        let mut accepted = None;
        for _ in 0..HALVINGS {
            let trial: [f64; 4] = core::array::from_fn(|i| x[i] - delta[i] * scale);
            let (ua, va) = crate::march::clamp(a, trial[0], trial[1]);
            let (ub, vb) = crate::march::clamp(b, trial[2], trial[3]);
            if let (Ok(pa), Ok(pb)) = (a.point_at(ua, va, tol), b.point_at(ub, vb, tol)) {
                let r = residual_of(pa, pb);
                let trial_size = norm(&r);
                if (trial_size < size || trial_size <= settled)
                    && let (Ok(da), Ok(db)) = (a.d1_at(ua, va, tol), b.d1_at(ub, vb, tol))
                {
                    accepted = Some((trial, pa, r, jacobian_of(da, db), trial_size));
                    break;
                }
            }
            scale *= 0.5;
        }
        let Some((next, at, r, j, next_size)) = accepted else {
            break;
        };
        let step = norm(&core::array::from_fn(|i| next[i] - x[i]));
        (x, point, residual, jacobian, size) = (next, at, r, j, next_size);
        if step <= tol.parametric() {
            break;
        }
    }
    (size <= tol.confusion()).then_some(point)
}

/// `m x = r` for four unknowns, by elimination with partial pivoting;
/// `None` at a zero pivot.
fn solve4(mut m: [[f64; 4]; 4], mut r: [f64; 4]) -> Option<[f64; 4]> {
    for col in 0..4 {
        let pivot = (col..4).max_by(|&i, &j| m[i][col].abs().total_cmp(&m[j][col].abs()))?;
        if m[pivot][col] == 0.0 {
            return None;
        }
        m.swap(col, pivot);
        r.swap(col, pivot);
        let head = m[col];
        for row in col + 1..4 {
            let factor = m[row][col] / head[col];
            for (entry, above) in m[row].iter_mut().zip(&head).skip(col) {
                *entry -= factor * above;
            }
            r[row] -= factor * r[col];
        }
    }
    for col in (0..4).rev() {
        r[col] /= m[col][col];
        for row in 0..col {
            r[row] -= m[row][col] * r[col];
        }
    }
    Some(r)
}

/// How far the curve stands from the trace it was fitted to: each sample's
/// distance to its nearest point on the curve.
///
/// The samples run in order along the curve, so each one's foot is found by
/// Newton steps from the last one's. Where that does not settle within
/// `bound`, the fit's own residual, which already bounds each sample's
/// distance from the curve at the fit's parameter, the nearest of the
/// curve's stations seeds the search instead, and the result never exceeds
/// `bound`.
fn trace_error(
    curve: &BSplineCurve,
    samples: &[ogeom_math::Point],
    bound: f64,
    tol: Tolerances,
) -> f64 {
    use ogeom_geom::Curve3d as _;
    let (lo, hi) = curve.knots().domain();
    let spans = curve.knots().distinct().len().saturating_sub(1).max(1);
    let count = (4 * spans).max(2 * samples.len()).max(200);
    #[allow(clippy::cast_precision_loss)]
    let at = |k: usize| lo + (hi - lo) * k as f64 / count as f64;
    let mut stations: Option<Vec<(f64, ogeom_math::Point)>> = None;
    let foot = |p: ogeom_math::Point, mut t: f64| -> (f64, f64) {
        let mut best = (t, f64::INFINITY);
        for _ in 0..8 {
            let (Ok(q), Ok(d)) = (curve.point_at(t, tol), curve.d1_at(t, tol)) else {
                break;
            };
            let gap = q.distance(p);
            if gap < best.1 {
                best = (t, gap);
            }
            let speed = d.dot(d);
            if speed <= f64::MIN_POSITIVE {
                break;
            }
            let next = (t + (p - q).dot(d) / speed).clamp(lo, hi);
            if (next - t).abs() <= (hi - lo) * 1e-12 {
                break;
            }
            t = next;
        }
        if let Ok(q) = curve.point_at(t, tol)
            && q.distance(p) < best.1
        {
            best = (t, q.distance(p));
        }
        best
    };
    let mut t = lo;
    let mut worst = 0.0_f64;
    for p in samples {
        let mut found = foot(*p, t);
        if found.1 > bound {
            let stations = stations.get_or_insert_with(|| {
                (0..=count)
                    .filter_map(|k| curve.point_at(at(k), tol).ok().map(|q| (at(k), q)))
                    .collect()
            });
            if let Some(k) = (0..stations.len()).min_by(|&x, &y| {
                stations[x]
                    .1
                    .distance(*p)
                    .total_cmp(&stations[y].1.distance(*p))
            }) {
                // Scanned finely over the stations either side, then
                // narrowed by golden section: where the curve all but stops
                // in space (its chart image swinging round a pole) or kinks
                // between two samples, Newton's steps are blind and the
                // distance has more than one dip between stations.
                let gap = |u: f64| {
                    curve
                        .point_at(u, tol)
                        .map_or(f64::INFINITY, |q| q.distance(*p))
                };
                let (from, to) = (
                    stations[k.saturating_sub(2)].0,
                    stations[(k + 2).min(stations.len() - 1)].0,
                );
                const FINE: u32 = 256;
                let h = (to - from) / f64::from(FINE);
                let start = (0..=FINE)
                    .map(|j| from + h * f64::from(j))
                    .min_by(|&x, &y| gap(x).total_cmp(&gap(y)))
                    .unwrap_or(from);
                let (mut a, mut b) = ((start - h).max(lo), (start + h).min(hi));
                let ratio = 0.5 * (5.0_f64.sqrt() - 1.0);
                for _ in 0..60 {
                    let (x, y) = (b - ratio * (b - a), a + ratio * (b - a));
                    if gap(x) <= gap(y) {
                        b = y;
                    } else {
                        a = x;
                    }
                }
                let again = foot(*p, 0.5 * (a + b));
                if again.1 < found.1 {
                    found = again;
                }
            }
        }
        t = found.0;
        worst = worst.max(found.1.min(bound));
    }
    worst
}

/// The fit's chart residual carried into space through the surface's
/// stretch along the pcurve: a bound on how far the lifted pcurve stands
/// from the trace, loose by the stretch where the residual is mostly in the
/// other coordinates, and so used only to cap an estimate, never stated.
fn space_error(
    surface: &SurfaceGeometry,
    pcurve: &BSpline2d,
    parameter_error: f64,
    tol: Tolerances,
) -> f64 {
    use ogeom_geom::Curve2d;
    // Convert the parameter-space error back through the surface's local
    // stretch at a few places; take the worst.
    let (lo, hi) = pcurve.domain();
    let mut worst = 0.0_f64;
    for i in 0..=16 {
        #[allow(clippy::cast_precision_loss)]
        let u = lo + (hi - lo) * f64::from(i) / 16.0;
        let Ok(at) = pcurve.point_at(u, tol) else {
            continue;
        };
        let Ok((du, dv)) = surface.d1_at(at.x, at.y, tol) else {
            continue;
        };
        let stretch = du.magnitude().max(dv.magnitude());
        worst = worst.max(parameter_error * stretch);
    }
    worst
}

/// Unfold parameter samples across a periodic surface's seam.
///
/// Each step is folded to the nearest image of the next sample, so a branch
/// crossing `u = 0` continues to `-0.1` rather than tearing to `2π - 0.1`. The
/// result may leave the surface's stated domain, which is what a pcurve
/// crossing a seam *is*.
fn unwrap_periodic(
    surface: &SurfaceGeometry,
    samples: &[(f64, f64)],
    tol: Tolerances,
) -> Vec<Point2> {
    let ((ua, ub), (va, vb)) = surface.domain();
    // Closure as well as periodicity: a converted drum is a clamped patch
    // that meets itself at its seam, and a loop walked round it lands on
    // either side of that seam by the walk's own rounding. Folded by the
    // chart's span like a period, the trace is the continuous curve it is.
    // Left as sampled, it jumps a whole span at the seam and the closed
    // fit chases the jump far from the trace.
    let u_period = if surface.is_periodic_u() || surface.is_closed_u(tol) {
        Some(ub - ua)
    } else {
        None
    };
    let v_period = if surface.is_periodic_v() || surface.is_closed_v(tol) {
        Some(vb - va)
    } else {
        None
    };
    let fold = |previous: f64, next: f64, period: Option<f64>| match period {
        None => next,
        Some(period) => {
            let mut candidate = next;
            while candidate - previous > period * 0.5 {
                candidate -= period;
            }
            while previous - candidate > period * 0.5 {
                candidate += period;
            }
            candidate
        }
    };

    let mut out = Vec::with_capacity(samples.len());
    let mut at = Point2::new(samples[0].0, samples[0].1);
    out.push(at);
    for sample in &samples[1..] {
        at = Point2::new(
            fold(at.x, sample.0, u_period),
            fold(at.y, sample.1, v_period),
        );
        out.push(at);
    }
    out
}

#[cfg(test)]
#[allow(clippy::unwrap_used)]
mod tests {
    use super::*;
    use crate::march::{Marching, branches};
    use ogeom_geom::{Curve2d, Curve3d, CylinderSurface, PlaneSurface, SphereSurface};
    use ogeom_math::{Cylinder, Direction, Frame, Plane, Point, Sphere, Vector};

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

    fn sphere(radius: f64) -> SurfaceGeometry {
        SphereSurface::new(Sphere::centred(Point::ORIGIN, radius, T).unwrap()).into()
    }

    fn cylinder(radius: f64) -> SurfaceGeometry {
        CylinderSurface::new(Cylinder::new(Frame::WORLD, radius, T).unwrap(), (-4.0, 4.0))
            .unwrap()
            .into()
    }

    fn plane(origin: Point, normal: Vector) -> SurfaceGeometry {
        PlaneSurface::over(
            Plane::through(origin, Direction::new(normal, T).unwrap()),
            (-6.0, 6.0),
            (-6.0, 6.0),
        )
        .unwrap()
        .into()
    }

    fn options() -> Marching {
        Marching {
            chord: 1e-5,
            ..Marching::default()
        }
    }

    /// The distance of a fitted curve from both surfaces, sampled densely.
    ///
    /// This is the measure the whole stage exists for: the *fit* (not the
    /// polyline it came from) is what downstream code holds, so the fit is
    /// what must lie on both surfaces.
    fn fitted_deviation(a: &SurfaceGeometry, b: &SurfaceGeometry, curve: &BSplineCurve) -> f64 {
        let off = |surface: &SurfaceGeometry, p: Point| match surface {
            SurfaceGeometry::Plane(x) => x.plane().distance_to(p),
            SurfaceGeometry::Sphere(x) => x.sphere().distance_to(p),
            SurfaceGeometry::Cylinder(x) => x.cylinder().distance_to(p),
            _ => 0.0,
        };
        let (lo, hi) = curve.knots().domain();
        let mut worst = 0.0_f64;
        for i in 0..=800 {
            #[allow(clippy::cast_precision_loss)]
            let u = lo + (hi - lo) * f64::from(i) / 800.0;
            if let Ok(p) = curve.point_at(u, T) {
                worst = worst.max(off(a, p).abs().max(off(b, p).abs()));
            }
        }
        worst
    }

    /// The error stated covers the fit between the trace's samples where
    /// the surfaces cross at a shallow angle. A bowl a hundred millimetres
    /// round at its foot dips two hundredths below a plane, both charts
    /// running evenly with `x` and `y`, and nine samples of half the circle
    /// they cross in are fitted through every one. Between the samples the
    /// fit and both pcurves leave the circle alike, within the plane, where
    /// the bowl stands only a fiftieth of that above it; the miss is read
    /// against the crossing itself.
    #[test]
    fn a_shallow_crossing_states_the_fits_miss_between_its_samples() {
        use crate::march::Stopped;
        use ogeom_geom::BSplineSurface;
        use ogeom_math::{ControlGrid, KnotVector};
        let (radius, dip, half) = (100.0, 0.02, 3.0);
        let a = plane(Point::ORIGIN, Vector::Z);
        // z = (x^2 + y^2) / (2 radius) - dip, exactly, as a quadratic patch
        // over [-half, half] in both.
        let square = [half * half, -half * half, half * half];
        let side = [-half, 0.0, half];
        let mut control = Vec::new();
        for i in 0..3 {
            for j in 0..3 {
                let z = (square[i] + square[j]) / (2.0 * radius) - dip;
                control.push(Point::new(side[i], side[j], z));
            }
        }
        let knots = || KnotVector::new(vec![0.0, 0.0, 0.0, 1.0, 1.0, 1.0], 2).unwrap();
        let b: SurfaceGeometry = BSplineSurface::new(
            knots(),
            knots(),
            &ControlGrid::new(control, 3, 3).unwrap(),
            T,
        )
        .unwrap()
        .into();
        let circle = (2.0 * radius * dip).sqrt();
        let SurfaceGeometry::Plane(flat) = &a else {
            unreachable!()
        };
        let mut branch = Traced {
            points: Vec::new(),
            on_a: Vec::new(),
            on_b: Vec::new(),
            stopped: Stopped::LeftTheDomain,
        };
        for k in 0..=8 {
            let angle = core::f64::consts::PI * f64::from(k) / 8.0;
            let p = Point::new(circle * angle.cos(), circle * angle.sin(), 0.0);
            let local = flat.plane().frame().to_local(p);
            let on_b = ((p.x + half) / (2.0 * half), (p.y + half) / (2.0 * half));
            assert!(b.point_at(on_b.0, on_b.1, T).unwrap().distance(p) < 1e-9);
            branch.points.push(p);
            branch.on_a.push((local.x, local.y));
            branch.on_b.push(on_b);
        }
        let fitted = approximate_branch(&a, &b, &branch, 1e-9, T).unwrap();
        let (lo, hi) = fitted.curve.knots().domain();
        let mut worst = 0.0_f64;
        for i in 0..=4000 {
            let p = fitted
                .curve
                .point_at(lo + (hi - lo) * f64::from(i) / 4000.0, T)
                .unwrap();
            let flat = Point::new(p.x, p.y, 0.0).to_vector().magnitude();
            worst = worst.max((flat - circle).abs().max(p.z.abs()));
        }
        assert!(
            worst > 1e-6,
            "the fit misses the circle between samples: {worst:.3e}"
        );
        assert!(
            fitted.fit_error >= 0.9 * worst,
            "the miss is stated: {:.3e} for {worst:.3e}",
            fitted.fit_error
        );
    }

    #[test]
    fn a_fitted_branch_lies_on_both_surfaces_to_the_stated_total() {
        // The tolerance story end to end: trace within 1e-5, fit within 1e-4,
        // so the fitted curve is within the sum of the two of the true
        // intersection, measured against the surfaces, not the polyline.
        let a = sphere(3.0);
        let b = cylinder(1.5);
        let found = branches(&a, &b, options(), T).unwrap();
        assert_eq!(found.len(), 2);

        for branch in &found {
            let fitted = approximate_branch(&a, &b, branch, 1e-4, T).unwrap();
            assert!(fitted.met, "fit error {:e}", fitted.fit_error);
            assert!(fitted.closed);
            let off = fitted_deviation(&a, &b, &fitted.curve);
            assert!(
                off <= 1e-4 + 1e-5,
                "the fitted curve is {off:e} off the surfaces"
            );
            // And it is compact: a curve, not a decorated polyline.
            assert!(
                fitted.curve.control_points().len() * 4 < branch.points.len(),
                "{} control points for {} samples",
                fitted.curve.control_points().len(),
                branch.points.len()
            );
        }
    }

    #[test]
    fn the_pcurves_lift_back_onto_the_curve() {
        // A pcurve is only worth having if evaluating it and lifting through
        // its surface lands on the intersection. Checked through both
        // surfaces at matched ends and sampled interiors.
        let a = sphere(3.0);
        let b = cylinder(1.5);
        let found = branches(&a, &b, options(), T).unwrap();
        let branch = &found[0];
        let fitted = approximate_branch(&a, &b, branch, 1e-4, T).unwrap();

        for (surface, pcurve) in [(&a, &fitted.on_a), (&b, &fitted.on_b)] {
            let (lo, hi) = pcurve.domain();
            for i in 0..=200 {
                #[allow(clippy::cast_precision_loss)]
                let u = lo + (hi - lo) * f64::from(i) / 200.0;
                let at = pcurve.point_at(u, T).unwrap();
                let lifted = surface.point_at(at.x, at.y, T).unwrap();
                // The lifted point is on its own surface by construction; what
                // matters is that it is on the *other* one too, i.e. on the
                // intersection.
                let off = match (surface as &SurfaceGeometry, &a, &b) {
                    _ if core::ptr::eq(surface, &a) => match &b {
                        SurfaceGeometry::Cylinder(c) => c.cylinder().distance_to(lifted),
                        _ => 0.0,
                    },
                    _ => match &a {
                        SurfaceGeometry::Sphere(s) => s.sphere().distance_to(lifted),
                        _ => 0.0,
                    },
                };
                assert!(
                    off.abs() < 5e-4,
                    "a lifted pcurve point is {off:e} off the intersection"
                );
            }
        }
    }

    #[test]
    fn a_branch_across_the_seam_gets_a_continuous_pcurve() {
        // A plane through a cylinder's axis at an angle produces an ellipse
        // whose pcurve crosses the cylinder's u = 0 seam. Folded naively the
        // pcurve tears by 2π; unwrapped it runs smoothly and leaves the stated
        // domain, which is what crossing a seam means.
        let a = cylinder(2.0);
        let b = plane(Point::ORIGIN, Vector::new(0.0, 0.4, 1.0));
        let found = branches(&a, &b, options(), T).unwrap();
        assert_eq!(found.len(), 1, "an oblique plane cuts one ellipse");
        let fitted = approximate_branch(&a, &b, &found[0], 1e-4, T).unwrap();

        // Continuity: no two adjacent samples of the fitted pcurve jump by
        // anything near a period.
        let (lo, hi) = fitted.on_a.domain();
        let mut previous = fitted.on_a.point_at(lo, T).unwrap();
        for i in 1..=400 {
            #[allow(clippy::cast_precision_loss)]
            let u = lo + (hi - lo) * f64::from(i) / 400.0;
            let at = fitted.on_a.point_at(u, T).unwrap();
            assert!(
                (at.x - previous.x).abs() < 1.0,
                "the pcurve tears at the seam: {} to {}",
                previous.x,
                at.x
            );
            previous = at;
        }
    }

    /// A loop walked round a converted drum is closed, seam or no seam.
    ///
    /// A cylinder converted to a patch is clamped, not periodic: it meets
    /// itself at its seam. A plane across it cuts a circle the walk reaches
    /// the seam on from both sides, each half stopping a fraction of a step
    /// short of it, and the joined branch has coincident ends. Left flagged
    /// as having left the domain, the arrangement downstream would hold a
    /// circle with two ends at one point. It is closed, and fitted as a loop
    /// whose chart image runs continuously across the seam.
    #[test]
    fn a_loop_cut_at_a_converted_drum_s_seam_is_closed() {
        let drum: SurfaceGeometry = cylinder(2.0).to_bspline(T).unwrap().into();
        assert!(matches!(drum, SurfaceGeometry::BSpline(_)));
        let cut = plane(Point::new(0.0, 0.0, 1.0), Vector::new(0.0, 0.2, 1.0));
        let found = branches(&drum, &cut, options(), T).unwrap();
        assert_eq!(found.len(), 1, "an oblique plane cuts one loop");
        assert!(found[0].closed(), "the loop closes on the seam");
        let fitted = approximate_branch(&drum, &cut, &found[0], 1e-4, T).unwrap();
        assert!(fitted.closed);
        assert!(
            fitted.fit_error < 1e-3,
            "the loop fits as one: {}",
            fitted.fit_error
        );
        let (lo, hi) = fitted.on_a.domain();
        let mut previous = fitted.on_a.point_at(lo, T).unwrap();
        for i in 1..=400 {
            let u = lo + (hi - lo) * f64::from(i) / 400.0;
            let at = fitted.on_a.point_at(u, T).unwrap();
            assert!(
                (at.x - previous.x).abs() < 0.5,
                "the chart image tears at the seam: {} to {}",
                previous.x,
                at.x
            );
            previous = at;
        }
    }

    /// A bore across a converted drum states the error its curves have.
    ///
    /// The joint fit's residual mixes millimetres with the drum chart's
    /// coordinates, and carried through the drum's stretch it reads hundreds
    /// of times larger than the distance the fitted curves stand from the
    /// two surfaces. The surfaces cross steeply here, so that distance is
    /// what is stated, within a small multiple.
    #[test]
    fn a_section_across_a_wide_drum_states_its_measured_error() {
        let radius = 23.6;
        let wall = CylinderSurface::new(
            Cylinder::new(Frame::WORLD, radius, T).unwrap(),
            (-30.0, 30.0),
        )
        .unwrap();
        let drum: SurfaceGeometry = SurfaceGeometry::from(wall).to_bspline(T).unwrap().into();
        let bore = Cylinder::new(
            Frame::new(
                Point::new(0.0, 0.0, 3.0),
                Direction::new(Vector::X, T).unwrap(),
                Direction::new(Vector::Y, T).unwrap(),
                T,
            )
            .unwrap(),
            9.0,
            T,
        )
        .unwrap();
        let drill: SurfaceGeometry = CylinderSurface::new(bore, (-40.0, 40.0)).unwrap().into();
        let marching = Marching {
            chord: 1e-4,
            ..Marching::default()
        };
        let found = branches(&drum, &drill, marching, T).unwrap();
        assert!(!found.is_empty());
        for branch in &found {
            let fitted = approximate_branch(&drum, &drill, branch, 1e-4, T).unwrap();
            let (lo, hi) = fitted.curve.knots().domain();
            let mut off = 0.0_f64;
            for i in 0..=2000 {
                let t = lo + (hi - lo) * f64::from(i) / 2000.0;
                let p = fitted.curve.point_at(t, T).unwrap();
                off = off.max(
                    wall.cylinder()
                        .distance_to(p)
                        .abs()
                        .max(bore.distance_to(p).abs()),
                );
            }
            // Honest both ways: no less than the curve stands off, and not
            // hundreds of times more.
            assert!(
                fitted.fit_error + marching.chord >= off,
                "states {:e} for a curve {off:e} off its surfaces",
                fitted.fit_error
            );
            assert!(
                fitted.fit_error <= 10.0 * off.max(marching.chord),
                "states {:e} for a curve {off:e} off its surfaces",
                fitted.fit_error
            );
        }
    }

    #[test]
    fn what_cannot_be_fitted_is_refused() {
        let a = sphere(1.0);
        let b = plane(Point::ORIGIN, Vector::Z);
        let found = branches(&a, &b, options(), T).unwrap();
        assert!(approximate_branch(&a, &b, &found[0], 0.0, T).is_err());
        assert!(approximate_branch(&a, &b, &found[0], -1.0, T).is_err());

        let empty = Traced {
            points: vec![],
            on_a: vec![],
            on_b: vec![],
            stopped: crate::march::Stopped::Stalled,
        };
        assert!(approximate_branch(&a, &b, &empty, 1e-4, T).is_err());
    }

    /// A drum whose wall holds a ball's axis meets the ball in two loops,
    /// each through a pole of its chart, where the chart's longitude turns
    /// half a turn and the curve does not. Each loop is cut at its pole and
    /// fitted open, pole to pole, its image on the ball ending on the pole
    /// row each side and inside the chart throughout; the curve and both
    /// images lifted lie on both surfaces to the stated error. Fitted
    /// whole, a loop closing on its pole joined the two ends half a turn
    /// apart, and its fit ran centimetres off both surfaces.
    #[test]
    fn loops_through_a_balls_poles_are_cut_there() {
        use core::f64::consts::FRAC_PI_2;
        let ball = sphere(10.0);
        for (cx, cy, r) in [(0.0, 1.0, 1.0), (0.2, 0.0, 0.2), (0.0, -0.3, 0.3)] {
            let frame =
                Frame::new(Point::new(cx, cy, -15.0), Direction::Z, Direction::Y, T).unwrap();
            let drum: SurfaceGeometry =
                CylinderSurface::new(Cylinder::new(frame, r, T).unwrap(), (0.0, 30.0))
                    .unwrap()
                    .into();
            let traced = branches(&ball, &drum, options(), T).unwrap();
            assert_eq!(traced.len(), 2, "drum {r} at ({cx}, {cy}): two loops");
            for branch in &traced {
                assert!(
                    !branch.closed(),
                    "drum {r} at ({cx}, {cy}): cut at its pole"
                );
                let fitted = approximate_branch(&ball, &drum, branch, 1e-6, T).unwrap();
                let (lo, hi) = fitted.curve.knots().domain();
                for t in [lo, hi] {
                    let p = fitted.curve.point_at(t, T).unwrap();
                    assert!(
                        (p.z.abs() - 10.0).abs() <= 1e-9,
                        "an end {p:?} off the pole"
                    );
                    let uv = fitted.on_a.point_at(t, T).unwrap();
                    assert!(
                        (uv.y.abs() - FRAC_PI_2).abs() <= 1e-12,
                        "an image end off the row"
                    );
                }
                let budget = options().chord + fitted.fit_error;
                assert!(budget < 1e-4, "drum {r} at ({cx}, {cy}): stated {budget}");
                for k in 0..=2000 {
                    let t = lo + (hi - lo) * f64::from(k) / 2000.0;
                    let p = fitted.curve.point_at(t, T).unwrap();
                    let on_ball = (p - Point::ORIGIN).magnitude() - 10.0;
                    let on_drum = (p.x - cx).hypot(p.y - cy) - r;
                    assert!(
                        on_ball.abs().max(on_drum.abs()) <= budget,
                        "drum {r} at ({cx}, {cy}): {p:?} off the surfaces"
                    );
                    let uv = fitted.on_a.point_at(t, T).unwrap();
                    assert!(uv.y.abs() <= FRAC_PI_2, "the image leaves the chart");
                    let lifted = ball.point_at(uv.x, uv.y, T).unwrap();
                    assert!(
                        lifted.distance(p) <= budget,
                        "the image stands off the curve"
                    );
                }
            }
        }
    }

    /// A drum whose wall passes a hair beside a ball's poles, farther than
    /// the trace's chord, meets the ball in two loops whose images on the
    /// ball swing half a turn of longitude within a few times that
    /// distance. Each fits within its tolerance, and the curve and both
    /// images lifted lie on both surfaces within the stated error, read
    /// densely. Fitted jointly alone, the loops missed by up to 0.35.
    #[test]
    fn loops_beside_a_balls_poles_fit_within_their_tolerance() {
        let ball = sphere(10.0);
        for (r, off) in [(1.0, 1e-3), (1.0, 1e-4), (2.5, 1e-3), (4.0, 1e-2)] {
            let wall = r * (1.0 + off);
            let frame =
                Frame::new(Point::new(0.0, -wall, -15.0), Direction::Z, Direction::Y, T).unwrap();
            let drum: SurfaceGeometry =
                CylinderSurface::new(Cylinder::new(frame, r, T).unwrap(), (0.0, 30.0))
                    .unwrap()
                    .into();
            let marching = Marching::default();
            let traced = branches(&ball, &drum, marching, T).unwrap();
            assert_eq!(traced.len(), 2, "drum {r} {off} off: two loops");
            for branch in &traced {
                let fitted = approximate_branch(&ball, &drum, branch, 1e-6, T).unwrap();
                assert!(
                    fitted.fit_error <= 2e-5,
                    "drum {r} {off} off: fit misses by {}",
                    fitted.fit_error
                );
                let (lo, hi) = fitted.curve.knots().domain();
                for k in 0..=20_000 {
                    let t = lo + (hi - lo) * f64::from(k) / 20_000.0;
                    let p = fitted.curve.point_at(t, T).unwrap();
                    let on_ball = (p - Point::ORIGIN).magnitude() - 10.0;
                    let on_drum = p.x.hypot(p.y + wall) - r;
                    assert!(
                        on_ball.abs().max(on_drum.abs()) <= marching.chord + fitted.fit_error,
                        "drum {r} {off} off: {p:?} off the surfaces"
                    );
                    for (surface, image) in [(&ball, &fitted.on_a), (&drum, &fitted.on_b)] {
                        let uv = image.point_at(t, T).unwrap();
                        let lifted = surface.point_at(uv.x, uv.y, T).unwrap();
                        assert!(
                            lifted.distance(p) <= fitted.fit_error,
                            "drum {r} {off} off: an image {} off the curve, stating {}",
                            lifted.distance(p),
                            fitted.fit_error
                        );
                    }
                }
            }
        }
    }
}