ogeom_intersect/

curve_surface.rs

//! Where a curve pierces a surface.
//!
//! Two consumers drive it: where an edge crosses a face in the boolean, and
//! the exact point-in-solid classifier, which is a ray/surface query per
//! face.
//!
//! # Well-posed
//!
//! `C(t) = S(u, v)` is three equations in three unknowns; unlike the
//! surface/surface system, nothing has to be pinned for Newton to converge to
//! a point. The analytic cases are still answered in closed form first: a line
//! against a plane or a quadric is a linear or quadratic equation, and solving
//! a quadratic by iteration would be slower and less exact than writing down
//! its roots.
//!
//! # A curve lying in the surface
//!
//! A line in a plane crosses it nowhere and everywhere. That is an overlap
//! (the parameter range of the curve that lies in the surface), and it is a
//! different answer from any list of points. Detected where the analytic
//! forms can see it, and on the general path as a run of piercings between
//! which the curve never leaves the surface.

use ogeom_core::{OgeomResult, Tolerances, ogeom_bail};
use ogeom_geom::{Curve, Curve3d, Surface, SurfaceGeometry};
use ogeom_math::{Point, Vector, solve};

use crate::march::{Cell, sample_by, segment_meets_triangle};

/// One piercing of a surface by a curve.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Piercing {
    /// The parameter on the curve.
    pub on_curve: f64,
    /// The parameters on the surface.
    pub on_surface: (f64, f64),
    /// Where, taken from the curve.
    pub point: Point,
    /// The distance between the two evaluations there.
    pub gap: f64,
}

/// What a curve does to a surface.
#[derive(Debug, Clone, PartialEq)]
pub struct CurveSurfaceIntersection {
    /// Isolated piercings, in order along the curve.
    pub crossings: Vec<Piercing>,
    /// Parameter ranges of the curve that lie *in* the surface.
    ///
    /// Detected for the analytic cases (a line in a plane) and, on the
    /// general path, as a run of piercings longer than the curve's sampling
    /// step between which the curve never leaves the surface.
    pub lying: Vec<(f64, f64)>,
}

impl CurveSurfaceIntersection {
    /// No contact found.
    #[must_use]
    pub fn is_empty(&self) -> bool {
        self.crossings.is_empty() && self.lying.is_empty()
    }

    const fn empty() -> Self {
        Self {
            crossings: Vec::new(),
            lying: Vec::new(),
        }
    }
}

/// How hard the general path looks.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct CurveSurfaceOptions {
    /// How many segments the curve is sampled into for seeding.
    pub samples: usize,
    /// How finely the surface is sampled, per direction.
    pub grid: usize,
    /// The widest gap that still counts as a piercing.
    pub gap: f64,
}

impl Default for CurveSurfaceOptions {
    fn default() -> Self {
        Self {
            samples: 128,
            grid: 24,
            gap: 1e-7,
        }
    }
}

/// Where a curve pierces a surface.
///
/// Analytic line/plane, line/sphere and line/cylinder are answered in closed
/// form; everything else is seeded polyhedrally and polished by Newton on the
/// well-posed three-by-three system.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the options
/// are unusable.
pub fn intersect_curve_surface(
    curve: &Curve,
    surface: &SurfaceGeometry,
    options: CurveSurfaceOptions,
    tol: Tolerances,
) -> OgeomResult<CurveSurfaceIntersection> {
    if let (Curve::Line(_), SurfaceGeometry::BSpline(spline)) = (curve, surface)
        && options.samples >= 2
        && options.gap.is_finite()
        && options.gap > 0.0
        && let Some(patches) = crate::bezier_clip::Patches::of(spline, tol)?
    {
        return clipped(curve, surface, &patches, options, tol);
    }
    intersect_with(curve, surface, options, tol, || {
        let cells = sample_by(surface, seeding(surface, options.grid), tol);
        let tree = CellTree::over(&cells, options.gap);
        (
            std::borrow::Cow::Owned(cells),
            std::borrow::Cow::Owned(tree),
        )
    })
}

/// Where a curve touches a surface without crossing it, solved from
/// `seed` on the curve: the touch, or `None` if there is none beside the
/// seed.
///
/// At a touch the gap between the curve and the surface falls to zero
/// without changing sign, so `C(t) = S(u, v)` has a double root there:
/// Newton on it slows to a halving of the distance per step and stops where
/// the gap, the square of the distance left, meets its tolerance, short of
/// the touch by the square root of that tolerance times the radius the two
/// part at. [`intersect_curve_surface`] lands anywhere in that valley. The
/// touch is solved here as one instead, from the seed's foot on the
/// surface: where that foot is the curve point's nearest (two rows) and the
/// curve runs in the surface's tangent plane (one row), rows that stay
/// regular where the two part at second order. It is kept where the curve
/// and the surface are within the confusion distance there.
#[must_use]
pub fn touch_curve_surface(
    curve: &Curve,
    surface: &SurfaceGeometry,
    seed: f64,
    tol: Tolerances,
) -> Option<Piercing> {
    let foot =
        ogeom_geom::project_on_surface(surface, curve.point_at(seed, tol).ok()?, 16, tol).ok()?;
    let along = |t: f64| {
        let d = curve.derivatives_at(t, 2, tol).ok()?;
        Some((Point::ORIGIN + *d.first()?, *d.get(1)?, *d.get(2)?))
    };
    let start = [seed, foot.parameters.0, foot.parameters.1];
    let [t, u, v] = touch_solved(along, surface, start, tol)?;
    let (lo, hi) = curve.domain();
    let ((ua, ub), (va, vb)) = surface.domain();
    let within = |x: f64, lo: f64, hi: f64, periodic: bool| {
        periodic || (x >= lo - tol.parametric() && x <= hi + tol.parametric())
    };
    if !within(t, lo, hi, curve.is_periodic())
        || !within(u, ua, ub, surface.is_periodic_u())
        || !within(v, va, vb, surface.is_periodic_v())
    {
        return None;
    }
    let point = curve.point_at(t, tol).ok()?;
    let gap = point.distance(surface.point_at(u, v, tol).ok()?);
    (gap <= tol.confusion()).then_some(Piercing {
        on_curve: t,
        on_surface: (u, v),
        point,
        gap,
    })
}

/// Where a curve touches a surface, solved from `start`, a parameter on the
/// curve and one on the surface beside the touch: the three there, or
/// `None` where the solve does not find a touch beside the start.
///
/// `along` gives the curve's point and first two derivatives at a
/// parameter. The curve touches the surface where the surface's foot of it
/// is its closest point (two rows), and the curve runs in the surface's
/// tangent plane (one row). The three rows are regular where the curve and
/// the surface part at second order. The gap along the normal is not a row,
/// as it touches zero without crossing it: the caller measures it.
///
/// The solve is kept only where its linear model holds over the move it
/// made, by Kantorovich's condition as [`crate::walk::settle`] asks it:
/// `‖J₀⁻¹·(J₁ − J₀)‖ ≤ ½` in the Frobenius norm, which is never less than
/// the operator norm. A start further from the touch than that, or one
/// beside a crossing, where the curve turns into the surface's tangent
/// plane only away from the crossing if at all, is not trusted to have
/// found a touch beside it.
pub(crate) fn touch_solved(
    along: impl Fn(f64) -> Option<(Point, Vector, Vector)>,
    surface: &SurfaceGeometry,
    start: [f64; 3],
    tol: Tolerances,
) -> Option<[f64; 3]> {
    let system = |x: &[f64; 3]| -> Option<([f64; 3], [[f64; 3]; 3])> {
        let (p, c, c2) = along(x[0])?;
        let (q, ou, ov) = surface.point_d1_at(x[1], x[2], tol).ok()?;
        let (ouu, ouv, ovv) = surface.d2_at(x[1], x[2], tol).ok()?;
        let gap = p - q;
        let normal = ou.cross(ov);
        let residual = [gap.dot(ou), gap.dot(ov), c.dot(normal)];
        let jacobian = [
            [
                c.dot(ou),
                gap.dot(ouu) - ou.dot(ou),
                gap.dot(ouv) - ov.dot(ou),
            ],
            [
                c.dot(ov),
                gap.dot(ouv) - ou.dot(ov),
                gap.dot(ovv) - ov.dot(ov),
            ],
            [
                c2.dot(normal),
                c.dot(ouu.cross(ov) + ou.cross(ouv)),
                c.dot(ouv.cross(ov) + ou.cross(ovv)),
            ],
        ];
        Some((residual, jacobian))
    };
    let (_, jacobian) = system(&start)?;
    // The rows are products of derivatives with no common scale, so the
    // solve is run until its steps stop moving the parameters, and what it
    // found is a root where the step from there would not move them either:
    // at the root the damped search cannot lower the residual further, and
    // reports running out.
    let criteria = solve::Criteria {
        residual: 0.0,
        step: tol.parametric(),
        max_iterations: 40,
    };
    let (found, _, _, _) = solve::newton_system_fixed(
        |x| system(x).unwrap_or(([f64::INFINITY; 3], [[0.0; 3]; 3])),
        start,
        criteria,
    )
    .ok()?;
    let (left, after) = system(&found)?;
    let step = solve::linear_fixed(after, left)?;
    if step.iter().any(|d| d.abs() > tol.parametric()) {
        return None;
    }
    let mut drift = 0.0;
    for column in 0..3 {
        let change: [f64; 3] =
            core::array::from_fn(|row| after[row][column] - jacobian[row][column]);
        let moved = solve::linear_fixed(jacobian, change)?;
        drift += moved.iter().map(|v| v * v).sum::<f64>();
    }
    (drift.sqrt() <= 0.5).then_some(found)
}

/// A surface made ready for many curve queries.
///
/// The general path seeds from a sampling of the surface into flat cells,
/// and that sampling depends only on the surface and the options. Taken
/// once, on the first query that needs it, and shared by every later one,
/// so a caller asking the same face about many curves (a ray per drawn
/// point, a ray per classified point) pays for it once. Every answer is the
/// one [`intersect_curve_surface`] gives for the same curve.
#[derive(Debug)]
pub struct PreparedSurface {
    surface: SurfaceGeometry,
    options: CurveSurfaceOptions,
    tol: Tolerances,
    cells: std::sync::OnceLock<(Vec<Cell>, CellTree)>,
    /// A spline's Bézier patches, for straight curves.
    patches: std::sync::OnceLock<Option<crate::bezier_clip::Patches>>,
}

impl PreparedSurface {
    /// Ready `surface` for queries at `options` and `tol`.
    #[must_use]
    pub const fn new(
        surface: SurfaceGeometry,
        options: CurveSurfaceOptions,
        tol: Tolerances,
    ) -> Self {
        Self {
            surface,
            options,
            tol,
            cells: std::sync::OnceLock::new(),
            patches: std::sync::OnceLock::new(),
        }
    }

    /// The surface queried.
    #[must_use]
    pub const fn surface(&self) -> &SurfaceGeometry {
        &self.surface
    }

    /// Where `curve` pierces the surface, as [`intersect_curve_surface`].
    ///
    /// # Errors
    ///
    /// As [`intersect_curve_surface`].
    pub fn intersect(&self, curve: &Curve) -> OgeomResult<CurveSurfaceIntersection> {
        let (surface, options, tol) = (&self.surface, self.options, self.tol);
        // A straight curve meets a spline where its patches' control hulls
        // say it may, and nowhere else.
        if let (Curve::Line(_), SurfaceGeometry::BSpline(spline)) = (curve, surface)
            && options.samples >= 2
            && options.gap.is_finite()
            && options.gap > 0.0
            && let Some(patches) = self
                .patches
                .get_or_init(|| crate::bezier_clip::Patches::of(spline, tol).ok().flatten())
        {
            return clipped(curve, surface, patches, options, tol);
        }
        intersect_with(curve, surface, options, tol, || {
            let (cells, tree) = self.cells.get_or_init(|| {
                let cells = sample_by(surface, seeding(surface, options.grid), tol);
                let tree = CellTree::over(&cells, options.gap);
                (cells, tree)
            });
            (
                std::borrow::Cow::Borrowed(cells.as_slice()),
                std::borrow::Cow::Borrowed(tree),
            )
        })
    }
}

fn intersect_with<'c>(
    curve: &Curve,
    surface: &SurfaceGeometry,
    options: CurveSurfaceOptions,
    tol: Tolerances,
    cells: impl FnOnce() -> (std::borrow::Cow<'c, [Cell]>, std::borrow::Cow<'c, CellTree>),
) -> OgeomResult<CurveSurfaceIntersection> {
    if options.samples < 2 || options.grid < 2 {
        ogeom_bail!(Construction, "seeding needs at least two steps each way");
    }
    if !options.gap.is_finite() || options.gap <= 0.0 {
        ogeom_bail!(Construction, "a gap of {} is not a distance", options.gap);
    }

    match (curve, surface) {
        (Curve::Line(line), SurfaceGeometry::Plane(p)) => {
            Ok(line_plane(line, p.plane(), curve, surface, tol))
        }
        (Curve::Line(line), SurfaceGeometry::Sphere(s)) => Ok(line_quadric(
            line,
            curve,
            surface,
            sphere_roots(line, s.sphere()),
            options,
            tol,
        )),
        (Curve::Line(line), SurfaceGeometry::Cylinder(c)) => Ok(line_quadric(
            line,
            curve,
            surface,
            cylinder_roots(line, c.cylinder()),
            options,
            tol,
        )),
        (Curve::Line(line), SurfaceGeometry::Cone(c)) => Ok(line_quadric(
            line,
            curve,
            surface,
            cone_roots(line, c.cone(), tol),
            options,
            tol,
        )),
        (Curve::Line(line), SurfaceGeometry::Torus(t)) => Ok(line_quadric(
            line,
            curve,
            surface,
            torus_roots(line, t.torus()),
            options,
            tol,
        )),
        _ => {
            let (cells, tree) = cells();
            general(curve, surface, &cells, &tree, options, tol)
        }
    }
}

// --- analytic ----------------------------------------------------------------

fn line_plane(
    line: &ogeom_geom::LineCurve,
    plane: ogeom_math::Plane,
    curve: &Curve,
    surface: &SurfaceGeometry,
    tol: Tolerances,
) -> CurveSurfaceIntersection {
    let axis = line.axis();
    let along = plane.normal().dot(axis.direction);
    let height = plane.signed_distance_to(axis.location);

    if along.abs() <= tol.angular() {
        // Parallel: in the plane, or never touching it.
        if height.abs() <= tol.confusion() {
            return CurveSurfaceIntersection {
                crossings: Vec::new(),
                lying: vec![line.domain()],
            };
        }
        return CurveSurfaceIntersection::empty();
    }

    let t = -height / along;
    let (lo, hi) = line.domain();
    if t < lo - tol.parametric() || t > hi + tol.parametric() {
        return CurveSurfaceIntersection::empty();
    }
    let point = axis.location + axis.direction.vector() * t;
    let Some(found) = invert(surface, point, curve, t, tol) else {
        return CurveSurfaceIntersection::empty();
    };
    if found.gap > tol.confusion() {
        // The crossing is real on the unbounded plane but outside this
        // surface's stated extents; the clamped polish says so as a gap.
        return CurveSurfaceIntersection::empty();
    }
    CurveSurfaceIntersection {
        crossings: vec![found],
        lying: Vec::new(),
    }
}

/// How far rounding can move a discriminant `p - q` from its true value.
fn rounding(p: f64, q: f64) -> f64 {
    8.0 * f64::EPSILON * p.abs().max(q.abs())
}

/// The line parameters at which a line meets a sphere.
fn sphere_roots(line: &ogeom_geom::LineCurve, sphere: ogeom_math::Sphere) -> Vec<f64> {
    let axis = line.axis();
    let d = axis.direction.vector();
    let m = axis.location - sphere.centre();
    // |m + t d|^2 = r^2, with |d| = 1.
    let b = m.dot(d);
    let c = sphere.radius().mul_add(-sphere.radius(), m.dot(m));
    let discriminant = b.mul_add(b, -c);
    // A tangent line has a discriminant of exactly zero, which rounding
    // leaves a few ulps either side; below zero by no more than that is the
    // tangency, and the polish that follows rejects a false one by its gap.
    if discriminant < -rounding(b * b, c) {
        return Vec::new();
    }
    let root = discriminant.max(0.0).sqrt();
    if root == 0.0 {
        vec![-b]
    } else {
        vec![-b - root, -b + root]
    }
}

/// The line parameters at which a line meets a cylinder.
/// A line's offset from a frame's origin and its direction, in the frame's
/// own axes.
fn in_frame(
    line: &ogeom_geom::LineCurve,
    frame: ogeom_math::Frame,
) -> (ogeom_math::Vector, ogeom_math::Vector) {
    let axis = line.axis();
    let local = |v: ogeom_math::Vector| {
        ogeom_math::Vector::new(
            v.dot(frame.x().vector()),
            v.dot(frame.y().vector()),
            v.dot(frame.z().vector()),
        )
    };
    (
        local(axis.location - frame.origin()),
        local(axis.direction.vector()),
    )
}

/// The line parameters at which a line meets a cone, either nappe: the
/// polish that follows keeps those on the cone's stated extent.
fn cone_roots(
    line: &ogeom_geom::LineCurve,
    cone: ogeom_math::Cone,
    tol: ogeom_core::Tolerances,
) -> Vec<f64> {
    let (m, d) = in_frame(line, cone.frame());
    // x^2 + y^2 = (r0 + k z)^2 along m + t d.
    let (r0, k) = (cone.reference_radius(), cone.half_angle().tan());
    let rim = k.mul_add(m.z, r0);
    let a = d.x.mul_add(d.x, d.y * d.y) - k * k * d.z * d.z;
    let b = 2.0 * (m.x.mul_add(d.x, m.y * d.y) - k * d.z * rim);
    let c = m.x.mul_add(m.x, m.y * m.y) - rim * rim;
    ogeom_math::solve::roots(&[c, b, a], tol.parametric()).unwrap_or_default()
}

/// The line parameters at which a line meets a torus: the real roots of
/// `(|w|^2 + R^2 - r^2)^2 = 4 R^2 (w_x^2 + w_y^2)` along `w = m + t d` in
/// the torus's frame.
fn torus_roots(line: &ogeom_geom::LineCurve, torus: ogeom_math::Torus) -> Vec<f64> {
    let (m, d) = in_frame(line, torus.frame());
    let (big, small) = (torus.major_radius(), torus.minor_radius());
    // Solved about the line's nearest approach to the centre, in units of
    // the torus's own size, so the roots stand near one: a torus hundreds
    // of millimetres down the line puts them there otherwise, and the
    // coefficients' spread costs them digits.
    let near = -m.dot(d) / d.dot(d);
    let size = big + small;
    let m = (m + d * near) / size;
    let (big, small) = (big / size, small / size);
    let a = d.dot(d);
    let b = 2.0 * m.dot(d);
    let c = m.dot(m) + big * big - small * small;
    let p = d.x.mul_add(d.x, d.y * d.y);
    let q = 2.0 * m.x.mul_add(d.x, m.y * d.y);
    let s = m.x.mul_add(m.x, m.y * m.y);
    let four = 4.0 * big * big;
    let coefficients = [
        c.mul_add(c, -four * s),
        2.0f64.mul_add(b * c, -four * q),
        b.mul_add(b, 2.0 * a * c) - four * p,
        2.0 * a * b,
        a * a,
    ];
    let (found, count) = quartic_roots(&coefficients);
    found[..count]
        .iter()
        .map(|sigma| near + sigma * size)
        .collect()
}

/// The real roots of a quartic whose roots stand near one, tangencies
/// included: a critical point where the quartic is within 1e-9 of its coefficients'
/// size of zero is a tangency, the coefficients being differences of terms
/// that size. Two roots a few roots of epsilon apart are merged into the
/// one tangency rounding split them from.
fn quartic_roots(c: &[f64; 5]) -> ([f64; 4], usize) {
    let mut found = [0.0; 4];
    let touch = 1e-9 * c.iter().map(|x| x.abs()).sum::<f64>();
    let count = ogeom_math::solve::real_roots(c, touch, &mut found);
    let mut merged = 0;
    for i in 0..count {
        let t = found[i];
        if merged > 0 && (t - found[merged - 1]).abs() <= 1e-6 * (1.0 + t.abs()) {
            found[merged - 1] = f64::midpoint(found[merged - 1], t);
        } else {
            found[merged] = t;
            merged += 1;
        }
    }
    (found, merged)
}

fn cylinder_roots(line: &ogeom_geom::LineCurve, cylinder: ogeom_math::Cylinder) -> Vec<f64> {
    let axis = line.axis();
    let w = cylinder.axis().direction.vector();
    // Strip the components along the cylinder's axis; what is left is a 2D
    // circle problem in the perpendicular plane.
    let d = axis.direction.vector();
    let m = axis.location - cylinder.axis().location;
    let d_perp = d - w * d.dot(w);
    let m_perp = m - w * m.dot(w);
    let a = d_perp.dot(d_perp);
    if a <= f64::MIN_POSITIVE {
        // The line runs along the axis direction: on the wall it would lie,
        // not pierce, and lying is not detected here.
        return Vec::new();
    }
    let b = d_perp.dot(m_perp);
    let c = cylinder
        .radius()
        .mul_add(-cylinder.radius(), m_perp.dot(m_perp));
    let discriminant = b.mul_add(b, -(a * c));
    if discriminant < -rounding(b * b, a * c) {
        return Vec::new();
    }
    let root = discriminant.max(0.0).sqrt();
    if root == 0.0 {
        vec![-b / a]
    } else {
        vec![(-b - root) / a, (-b + root) / a]
    }
}

/// Roots dressed as piercings, filtered by the line's own range and the
/// surface's extents.
fn line_quadric(
    line: &ogeom_geom::LineCurve,
    curve: &Curve,
    surface: &SurfaceGeometry,
    roots: Vec<f64>,
    options: CurveSurfaceOptions,
    tol: Tolerances,
) -> CurveSurfaceIntersection {
    let axis = line.axis();
    let (lo, hi) = line.domain();
    let mut crossings = Vec::new();
    for t in roots {
        if t < lo - tol.parametric() || t > hi + tol.parametric() {
            continue;
        }
        let point = axis.location + axis.direction.vector() * t;
        let Some(found) = invert(surface, point, curve, t, tol) else {
            continue;
        };
        // The extent check is the gap. The polish clamps the surface
        // parameters into the stated domain, so a root beyond the cylinder's
        // height converges to the rim with a gap of exactly how far past it
        // was: a piercing of the unbounded geometry, not of this surface.
        if found.gap > tol.confusion() {
            continue;
        }
        let _ = options;
        crossings.push(Piercing {
            on_curve: t,
            on_surface: found.on_surface,
            point,
            gap: found.gap,
        });
    }
    crossings.sort_by(|a, b| {
        a.on_curve
            .partial_cmp(&b.on_curve)
            .unwrap_or(core::cmp::Ordering::Equal)
    });
    CurveSurfaceIntersection {
        crossings,
        lying: Vec::new(),
    }
}

/// Surface parameters of a point known to lie on an analytic surface.
///
/// Closed-form inversion for the quadrics; refined by one Newton pass so the
/// reported parameters evaluate back onto the point to rounding.
fn invert(
    surface: &SurfaceGeometry,
    point: Point,
    curve: &Curve,
    on_curve: f64,
    tol: Tolerances,
) -> Option<Piercing> {
    let guess = match surface {
        SurfaceGeometry::Plane(p) => {
            let local = p.plane().frame().to_local(point);
            (local.x, local.y)
        }
        SurfaceGeometry::Sphere(s) => {
            let local = s.sphere().frame().to_local(point);
            let latitude = (local.z / s.sphere().radius()).clamp(-1.0, 1.0).asin();
            (
                local.y.atan2(local.x).rem_euclid(core::f64::consts::TAU),
                latitude,
            )
        }
        SurfaceGeometry::Cylinder(c) => {
            let local = c.cylinder().frame().to_local(point);
            (
                local.y.atan2(local.x).rem_euclid(core::f64::consts::TAU),
                local.z,
            )
        }
        SurfaceGeometry::Cone(c) => {
            ogeom_math::elementary::cone_parameters(&c.cone(), point, tol).ok()?
        }
        SurfaceGeometry::Torus(t) => {
            ogeom_math::elementary::torus_parameters(&t.torus(), point, tol).ok()?
        }
        _ => return None,
    };
    // One polish step against the curve point, so parameter rounding in the
    // inversion does not survive into the result, and the gap comes with it,
    // because the polish clamps into the surface's extents and the gap is
    // what says whether the clamped answer still touches the curve.
    polish(curve, surface, on_curve, guess, tol)
}

// --- general -----------------------------------------------------------------

fn general(
    curve: &Curve,
    surface: &SurfaceGeometry,
    cells: &[Cell],
    tree: &CellTree,
    options: CurveSurfaceOptions,
    tol: Tolerances,
) -> OgeomResult<CurveSurfaceIntersection> {
    let (lo, hi) = curve.domain();

    let mut points = Vec::with_capacity(options.samples + 1);
    for i in 0..=options.samples {
        #[allow(clippy::cast_precision_loss)]
        let t = lo + (hi - lo) * i as f64 / options.samples as f64;
        if let Ok(p) = curve.point_at(t, tol) {
            points.push((t, p));
        }
    }

    let segments: Vec<((f64, Point), (f64, Point))> =
        points.windows(2).map(|pair| (pair[0], pair[1])).collect();
    // Which segments come near which cells, by cell index.
    let mut pairs: Vec<(usize, usize)> = Vec::new();
    let meet = |i: usize, c: usize, pairs: &mut Vec<(usize, usize)>| {
        let cell = &cells[c];
        let ((_, p0), (_, p1)) = segments[i];
        if segment_near_cell(p0, p1, cell, options.gap.max(cell.sag)) {
            pairs.push((i, c));
        }
    };
    // Each segment asks the tree for the cells whose boxes, grown by their
    // own margins, its box meets.
    let mut found: Vec<usize> = Vec::new();
    for (i, &((_, p0), (_, p1))) in segments.iter().enumerate() {
        let low = Point::new(p0.x.min(p1.x), p0.y.min(p1.y), p0.z.min(p1.z));
        let high = Point::new(p0.x.max(p1.x), p0.y.max(p1.y), p0.z.max(p1.z));
        found.clear();
        tree.visit(low, high, &mut found);
        for &c in &found {
            meet(i, c, &mut pairs);
        }
    }
    pairs.sort_unstable();
    let mut crossings: Vec<Piercing> = Vec::new();
    for (segment, k) in pairs {
        let ((t0, p0), (t1, p1)) = segments[segment];
        let cell = &cells[k];
        {
            // Near the cell within the surface's own bow from it: a curve
            // crossing the surface in the gap between the flat cell and
            // the curved patch it stands for (a ray starting a few microns
            // from the wall it leaves by) meets no cell, and is seeded by
            // every segment passing within the cell's bow of its triangle.
            if segment_meets_triangle(p0, p1, cell.corners).is_none()
                && !(cell.sag > options.gap
                    && segment_near_cell(p0, p1, cell, cell.sag)
                    && segment_triangle_distance(p0, p1, cell.corners) <= cell.sag)
            {
                continue;
            }
            // Newton from where the segment meets the cell, not from the
            // cell's corner: on a wall bowing a few hundredths of a
            // millimetre over a cell, the corner can stand far enough off
            // that the first step leaves the chart, and clamped at its edge
            // the solve stalls there. The corner stays a second try.
            let (near_t, near_uv) = seed_in(cell, p0, p1, t0, t1);
            let Some(found) = [(near_t, near_uv), (f64::midpoint(t0, t1), cell.at)]
                .into_iter()
                .filter_map(|(t, uv)| polish(curve, surface, t, uv, tol))
                .find(|found| found.gap <= options.gap)
            else {
                continue;
            };
            let reach = tol.confusion() * 100.0;
            if !crossings
                .iter()
                .any(|c| c.point.distance(found.point) <= reach)
            {
                crossings.push(found);
            }
        }
    }
    crossings.sort_by(|a, b| {
        a.on_curve
            .partial_cmp(&b.on_curve)
            .unwrap_or(core::cmp::Ordering::Equal)
    });
    Ok(gathered(curve, surface, crossings, options, tol))
}

/// Where a straight curve pierces a spline surface, seeded where the
/// surface's patches' control hulls come within the gap of it.
fn clipped(
    curve: &Curve,
    surface: &SurfaceGeometry,
    patches: &crate::bezier_clip::Patches,
    options: CurveSurfaceOptions,
    tol: Tolerances,
) -> OgeomResult<CurveSurfaceIntersection> {
    let (lo, hi) = curve.domain();
    let (a, b) = (curve.point_at(lo, tol)?, curve.point_at(hi, tol)?);
    let reach = tol.confusion() * 100.0;
    let mut crossings: Vec<Piercing> = Vec::new();
    for (f, uv) in patches.seeds(a, b, options.gap) {
        let Some(found) = polish(curve, surface, (hi - lo).mul_add(f, lo), uv, tol) else {
            continue;
        };
        if found.gap > options.gap
            || crossings
                .iter()
                .any(|c| c.point.distance(found.point) <= reach)
        {
            continue;
        }
        crossings.push(found);
    }
    crossings.sort_by(|a, b| {
        a.on_curve
            .partial_cmp(&b.on_curve)
            .unwrap_or(core::cmp::Ordering::Equal)
    });
    Ok(gathered(curve, surface, crossings, options, tol))
}

/// Piercings gathered into the contacts they describe.
///
/// Newton lands a little apart each time it is started near a tangency,
/// where the gap touches zero without crossing it, and all along a stretch
/// of curve lying in the surface; a radius of rounding merges neither. So
/// neighbouring piercings between which the curve never leaves the surface
/// are one contact: a run shorter than the curve's sampling step is one
/// tangency, kept once (where the gap is least), and a longer one is a
/// stretch lying in the surface.
fn gathered(
    curve: &Curve,
    surface: &SurfaceGeometry,
    crossings: Vec<Piercing>,
    options: CurveSurfaceOptions,
    tol: Tolerances,
) -> CurveSurfaceIntersection {
    let stays_on = |x: &Piercing, y: &Piercing| -> bool {
        (1..=3).all(|k| {
            let t = x.on_curve + (y.on_curve - x.on_curve) * f64::from(k) / 4.0;
            let Ok(p) = curve.point_at(t, tol) else {
                return false;
            };
            // From either end: by a pole or an apex one side's chart
            // turns too fast for the other's seed.
            [x.on_surface, y.on_surface].into_iter().any(|seed| {
                crate::march::nearest_on(surface, seed, p, tol)
                    .is_some_and(|(_, q)| q.distance(p) <= options.gap)
            })
        })
    };
    let (lo, hi) = curve.domain();
    #[allow(clippy::cast_precision_loss)]
    let step = (hi - lo).abs() / options.samples.max(1) as f64;
    let mut runs: Vec<Vec<Piercing>> = Vec::new();
    for crossing in crossings {
        match runs.last_mut() {
            Some(run) if run.last().is_some_and(|last| stays_on(last, &crossing)) => {
                run.push(crossing);
            }
            _ => runs.push(vec![crossing]),
        }
    }
    let mut out = CurveSurfaceIntersection {
        crossings: Vec::new(),
        lying: Vec::new(),
    };
    for run in runs {
        let (Some(first), Some(last)) = (run.first(), run.last()) else {
            continue;
        };
        if run.len() > 1 && (last.on_curve - first.on_curve).abs() >= step {
            out.lying.push((first.on_curve, last.on_curve));
        } else if let Some(best) = run.iter().min_by(|a, b| {
            a.gap
                .partial_cmp(&b.gap)
                .unwrap_or(core::cmp::Ordering::Equal)
        }) {
            out.crossings.push(*best);
        }
    }
    out
}

/// How many seed cells to lay along each direction of a surface: the
/// asked grid, and for a spline at least two per knot span. A patch swept
/// several turns round an axis (a thread's flank) spans dozens of knots
/// along its length, and a grid of the asked size lays flat cells a turn's
/// fraction wide whose chords stand a tenth of a millimetre off the wall;
/// a curve crossing the wall inside that gap meets no cell and is missed.
fn seeding(surface: &SurfaceGeometry, grid: usize) -> (usize, usize) {
    const CAP: usize = 1024;
    let SurfaceGeometry::BSpline(spline) = surface else {
        return (grid, grid);
    };
    let spans = |knots: &ogeom_math::KnotVector| knots.distinct().len().saturating_sub(1);
    (
        grid.max(2 * spans(spline.u_knots())).min(CAP.max(grid)),
        grid.max(2 * spans(spline.v_knots())).min(CAP.max(grid)),
    )
}

/// Where to start Newton for a segment near a cell: the point of the cell
/// the segment passes through, or failing that the one nearest its middle,
/// with its parameters on the curve and on the surface read off the cell's
/// corners.
fn seed_in(cell: &Cell, p0: Point, p1: Point, t0: f64, t1: f64) -> (f64, (f64, f64)) {
    let [a, b, c] = cell.corners;
    let (t, at) = segment_meets_triangle(p0, p1, cell.corners).map_or_else(
        || (f64::midpoint(t0, t1), p0.midpoint(p1)),
        |x| {
            let length = p0.distance(p1);
            let f = if length > 0.0 {
                p0.distance(x) / length
            } else {
                0.5
            };
            (t0 + (t1 - t0) * f, x)
        },
    );
    // Barycentric weights of the point's foot in the cell's plane, pulled
    // back inside the cell.
    let (e1, e2, d) = (b - a, c - a, at - a);
    let (d11, d12, d22) = (e1.dot(e1), e1.dot(e2), e2.dot(e2));
    let (d1, d2) = (d.dot(e1), d.dot(e2));
    let det = d11 * d22 - d12 * d12;
    let (mut wb, mut wc) = if det > 0.0 {
        ((d22 * d1 - d12 * d2) / det, (d11 * d2 - d12 * d1) / det)
    } else {
        (1.0 / 3.0, 1.0 / 3.0)
    };
    wb = wb.clamp(0.0, 1.0);
    wc = wc.clamp(0.0, 1.0);
    if wb + wc > 1.0 {
        let sum = wb + wc;
        wb /= sum;
        wc /= sum;
    }
    let wa = 1.0 - wb - wc;
    let [pa, pb, pc] = cell.params;
    (
        t,
        (
            wa * pa.0 + wb * pb.0 + wc * pc.0,
            wa * pa.1 + wb * pb.1 + wc * pc.1,
        ),
    )
}

/// The distance between the segment `p..q` and the triangle `corners`:
/// nothing where it pierces the triangle, else the least of its ends'
/// distances to the triangle and its distances to the triangle's sides.
fn segment_triangle_distance(p: Point, q: Point, corners: [Point; 3]) -> f64 {
    if segment_meets_triangle(p, q, corners).is_some() {
        return 0.0;
    }
    let [a, b, c] = corners;
    let mut least = point_triangle_distance(p, corners).min(point_triangle_distance(q, corners));
    for (x, y) in [(a, b), (b, c), (c, a)] {
        least = least.min(segment_segment_distance(p, q, x, y));
    }
    least
}

/// The distance from `p` to the triangle `corners`.
fn point_triangle_distance(p: Point, [a, b, c]: [Point; 3]) -> f64 {
    let normal = (b - a).cross(c - a);
    let area = normal.magnitude();
    if area > 0.0 {
        let n = normal / area;
        let foot = p - n * (p - a).dot(n);
        let inside = [(a, b), (b, c), (c, a)]
            .iter()
            .all(|&(x, y)| (y - x).cross(foot - x).dot(n) >= 0.0);
        if inside {
            return (p - a).dot(n).abs();
        }
    }
    distance_to_segment(p, a, b)
        .min(distance_to_segment(p, b, c))
        .min(distance_to_segment(p, c, a))
}

/// The distance between the segments `p..q` and `r..s`.
fn segment_segment_distance(p: Point, q: Point, r: Point, s: Point) -> f64 {
    let (d1, d2, w) = (q - p, s - r, p - r);
    let (a, e, f) = (d1.dot(d1), d2.dot(d2), d2.dot(w));
    let (b, c) = (d1.dot(d2), d1.dot(w));
    let denominator = a * e - b * b;
    let mut t = if denominator > 1e-300 {
        ((b * f - c * e) / denominator).clamp(0.0, 1.0)
    } else {
        0.0
    };
    let mut u = if e > 0.0 { (b * t + f) / e } else { 0.0 };
    if u < 0.0 {
        u = 0.0;
        t = if a > 0.0 {
            (-c / a).clamp(0.0, 1.0)
        } else {
            0.0
        };
    } else if u > 1.0 {
        u = 1.0;
        t = if a > 0.0 {
            ((b - c) / a).clamp(0.0, 1.0)
        } else {
            0.0
        };
    }
    (p + d1 * t).distance(r + d2 * u)
}

/// The distance from `p` to the segment `a..b`.
fn distance_to_segment(p: Point, a: Point, b: Point) -> f64 {
    let ab = b - a;
    let length = ab.dot(ab);
    let f = if length > 0.0 {
        ((p - a).dot(ab) / length).clamp(0.0, 1.0)
    } else {
        0.0
    };
    (a + ab * f).distance(p)
}

/// The cells binned in a tree of boxes, each cell's box grown by the
/// farthest a curve may stand from it and still be seeded there: the gap,
/// or the cell's bow where that is wider. Built once for a surface and
/// asked once per segment of every curve.
#[derive(Debug, Clone)]
pub(crate) struct CellTree {
    nodes: Vec<TreeNode>,
    order: Vec<usize>,
}

#[derive(Debug, Clone)]
struct TreeNode {
    low: Point,
    high: Point,
    /// A leaf's range of `order`, or an inner node's two children.
    span: (usize, usize),
    leaf: bool,
}

impl CellTree {
    /// The most cells a leaf holds.
    const LEAF: usize = 8;

    pub(crate) fn over(cells: &[Cell], gap: f64) -> Self {
        let boxes: Vec<(Point, Point)> = cells
            .iter()
            .map(|c| {
                let m = gap.max(c.sag);
                (
                    Point::new(c.low.x - m, c.low.y - m, c.low.z - m),
                    Point::new(c.high.x + m, c.high.y + m, c.high.z + m),
                )
            })
            .collect();
        let mut tree = Self {
            nodes: Vec::new(),
            order: (0..cells.len()).collect(),
        };
        if !cells.is_empty() {
            tree.build(&boxes, 0, cells.len());
        }
        tree
    }

    fn build(&mut self, boxes: &[(Point, Point)], from: usize, to: usize) -> usize {
        let (mut low, mut high) = boxes[self.order[from]];
        for &i in &self.order[from..to] {
            let (l, h) = boxes[i];
            low = Point::new(low.x.min(l.x), low.y.min(l.y), low.z.min(l.z));
            high = Point::new(high.x.max(h.x), high.y.max(h.y), high.z.max(h.z));
        }
        let at = self.nodes.len();
        self.nodes.push(TreeNode {
            low,
            high,
            span: (from, to),
            leaf: true,
        });
        if to - from <= Self::LEAF {
            return at;
        }
        // Split at the median along the box's longest side.
        let extent = high - low;
        let axis = if extent.x >= extent.y && extent.x >= extent.z {
            0
        } else if extent.y >= extent.z {
            1
        } else {
            2
        };
        let key = |i: usize| {
            let (l, h) = boxes[i];
            match axis {
                0 => l.x + h.x,
                1 => l.y + h.y,
                _ => l.z + h.z,
            }
        };
        let middle = from + (to - from) / 2;
        self.order[from..to]
            .select_nth_unstable_by(middle - from, |&a, &b| key(a).total_cmp(&key(b)));
        let left = self.build(boxes, from, middle);
        let right = self.build(boxes, middle, to);
        self.nodes[at].span = (left, right);
        self.nodes[at].leaf = false;
        at
    }

    /// The cells whose grown boxes meet the box `low..high`, added to `out`.
    pub(crate) fn visit(&self, low: Point, high: Point, out: &mut Vec<usize>) {
        if self.nodes.is_empty() {
            return;
        }
        let mut stack = vec![0_usize];
        while let Some(n) = stack.pop() {
            let node = &self.nodes[n];
            if node.low.x > high.x
                || low.x > node.high.x
                || node.low.y > high.y
                || low.y > node.high.y
                || node.low.z > high.z
                || low.z > node.high.z
            {
                continue;
            }
            if node.leaf {
                out.extend_from_slice(&self.order[node.span.0..node.span.1]);
            } else {
                stack.push(node.span.0);
                stack.push(node.span.1);
            }
        }
    }
}

/// Whether a segment's box comes near a cell's.
fn segment_near_cell(a: Point, b: Point, cell: &Cell, margin: f64) -> bool {
    let low = Point::new(a.x.min(b.x), a.y.min(b.y), a.z.min(b.z));
    let high = Point::new(a.x.max(b.x), a.y.max(b.y), a.z.max(b.z));
    low.x <= cell.high.x + margin
        && cell.low.x <= high.x + margin
        && low.y <= cell.high.y + margin
        && cell.low.y <= high.y + margin
        && low.z <= cell.high.z + margin
        && cell.low.z <= high.z + margin
}

/// Newton on the well-posed system `C(t) = S(u, v)`.
fn polish(
    curve: &Curve,
    surface: &SurfaceGeometry,
    seed_t: f64,
    seed_uv: (f64, f64),
    tol: Tolerances,
) -> Option<Piercing> {
    let clamp_t = |t: f64| {
        let (lo, hi) = curve.domain();
        if curve.is_periodic() {
            let span = hi - lo;
            if span > 0.0 {
                return lo + (t - lo).rem_euclid(span);
            }
        }
        t.clamp(lo, hi)
    };
    let clamp_uv = |u: f64, v: f64| {
        let ((ua, ub), (va, vb)) = surface.domain();
        let fold = |x: f64, lo: f64, hi: f64, periodic: bool| {
            if periodic {
                let span = hi - lo;
                if span > 0.0 {
                    return lo + (x - lo).rem_euclid(span);
                }
            }
            x.clamp(lo, hi)
        };
        (
            fold(u, ua, ub, surface.is_periodic_u()),
            fold(v, va, vb, surface.is_periodic_v()),
        )
    };

    // The gap and its Jacobian apart: the damped step measures many trial
    // points and steps onto few, and a spline's derivatives cost more than
    // its point.
    let gap_at = |x: &[f64; 3]| {
        let t = clamp_t(x[0]);
        let (u, v) = clamp_uv(x[1], x[2]);
        let (Ok(pc), Ok(ps)) = (curve.point_at(t, tol), surface.point_at(u, v, tol)) else {
            // Nowhere to measure from: infinite, so the damped step backs off.
            return [f64::INFINITY; 3];
        };
        let gap = pc - ps;
        [gap.x, gap.y, gap.z]
    };
    let jacobian_at = |x: &[f64; 3]| {
        let t = clamp_t(x[0]);
        let (u, v) = clamp_uv(x[1], x[2]);
        let (Ok(dc), Ok((du, dv))) = (curve.d1_at(t, tol), surface.d1_at(u, v, tol)) else {
            return None;
        };
        Some([
            [dc.x, -du.x, -dv.x],
            [dc.y, -du.y, -dv.y],
            [dc.z, -du.z, -dv.z],
        ])
    };
    let criteria = solve::Criteria {
        residual: tol.confusion() * 0.01,
        step: tol.parametric(),
        max_iterations: 40,
    };
    let found = solve::newton_system_fixed_lazy(
        gap_at,
        jacobian_at,
        [seed_t, seed_uv.0, seed_uv.1],
        criteria,
    )
    .ok()?;
    let t = clamp_t(found.0[0]);
    let (u, v) = clamp_uv(found.0[1], found.0[2]);
    let pc = curve.point_at(t, tol).ok()?;
    let ps = surface.point_at(u, v, tol).ok()?;
    Some(Piercing {
        on_curve: t,
        on_surface: (u, v),
        point: pc,
        gap: pc.distance(ps),
    })
}

#[cfg(test)]
#[allow(clippy::unwrap_used)]
mod tests {
    use super::*;
    use ogeom_geom::{
        BSplineCurve, CircleCurve, CylinderSurface, LineCurve, PlaneSurface, SphereSurface,
    };
    use ogeom_math::{Circle, Cylinder, Direction, Frame, KnotVector, Plane, 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, height: (f64, f64)) -> SurfaceGeometry {
        CylinderSurface::new(Cylinder::new(Frame::WORLD, radius, T).unwrap(), height)
            .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 segment(from: Point, to: Point) -> Curve {
        LineCurve::segment(from, to, T).unwrap().into()
    }

    #[test]
    fn a_line_through_a_sphere_pierces_it_where_the_quadratic_says() {
        let ball = sphere(2.0);
        let ray = segment(Point::new(-5.0, 0.0, 0.0), Point::new(5.0, 0.0, 0.0));
        let found =
            intersect_curve_surface(&ray, &ball, CurveSurfaceOptions::default(), T).unwrap();
        assert_eq!(found.crossings.len(), 2);
        assert!(
            found.crossings[0]
                .point
                .is_equal(Point::new(-2.0, 0.0, 0.0), T)
        );
        assert!(
            found.crossings[1]
                .point
                .is_equal(Point::new(2.0, 0.0, 0.0), T)
        );
        for hit in &found.crossings {
            assert!(hit.gap < 1e-12);
            // The surface parameters evaluate back onto the point.
            let lifted = ball
                .point_at(hit.on_surface.0, hit.on_surface.1, T)
                .unwrap();
            assert!(lifted.is_equal(hit.point, T));
        }

        // Tangent: one root. Missing: none.
        let grazing = segment(Point::new(-5.0, 0.0, 2.0), Point::new(5.0, 0.0, 2.0));
        assert_eq!(
            intersect_curve_surface(&grazing, &ball, CurveSurfaceOptions::default(), T)
                .unwrap()
                .crossings
                .len(),
            1
        );
        let missing = segment(Point::new(-5.0, 0.0, 3.0), Point::new(5.0, 0.0, 3.0));
        assert!(
            intersect_curve_surface(&missing, &ball, CurveSurfaceOptions::default(), T)
                .unwrap()
                .is_empty()
        );
    }

    #[test]
    fn a_line_through_a_cylinder_respects_its_height() {
        let drum = cylinder(2.0, (-1.0, 1.0));
        // Crosses the infinite cylinder at z = 0: inside the height, two hits.
        let level = segment(Point::new(-5.0, 0.0, 0.0), Point::new(5.0, 0.0, 0.0));
        assert_eq!(
            intersect_curve_surface(&level, &drum, CurveSurfaceOptions::default(), T)
                .unwrap()
                .crossings
                .len(),
            2
        );
        // Crosses at z = 3: the unbounded geometry meets it, this surface
        // does not reach there.
        let high = segment(Point::new(-5.0, 0.0, 3.0), Point::new(5.0, 0.0, 3.0));
        assert!(
            intersect_curve_surface(&high, &drum, CurveSurfaceOptions::default(), T)
                .unwrap()
                .is_empty()
        );
    }

    #[test]
    fn a_line_lying_in_a_plane_is_an_overlap_not_a_crossing_list() {
        let ground = plane(Point::ORIGIN, Vector::Z);
        let lying = segment(Point::new(-3.0, 1.0, 0.0), Point::new(3.0, 1.0, 0.0));
        let found =
            intersect_curve_surface(&lying, &ground, CurveSurfaceOptions::default(), T).unwrap();
        assert!(found.crossings.is_empty());
        assert_eq!(found.lying.len(), 1);

        let crossing = segment(Point::new(0.0, 0.0, -1.0), Point::new(0.0, 0.0, 1.0));
        let found =
            intersect_curve_surface(&crossing, &ground, CurveSurfaceOptions::default(), T).unwrap();
        assert_eq!(found.crossings.len(), 1);
        assert!(found.crossings[0].point.is_equal(Point::ORIGIN, T));

        let parallel = segment(Point::new(-3.0, 0.0, 1.0), Point::new(3.0, 0.0, 1.0));
        assert!(
            intersect_curve_surface(&parallel, &ground, CurveSurfaceOptions::default(), T)
                .unwrap()
                .is_empty()
        );
    }

    #[test]
    fn a_circle_pierces_a_plane_twice_through_the_general_path() {
        // A circle in the xz plane against the ground: no analytic case
        // handles circle/plane here, so this is the seeded Newton path, and
        // the answer is known exactly anyway.
        let ring: Curve = CircleCurve::new(
            Circle::new(
                Frame::new(Point::new(0.0, 0.0, 0.0), -Direction::Y, Direction::X, T).unwrap(),
                2.0,
                T,
            )
            .unwrap(),
        )
        .into();
        let ground = plane(Point::ORIGIN, Vector::Z);
        let found =
            intersect_curve_surface(&ring, &ground, CurveSurfaceOptions::default(), T).unwrap();
        assert_eq!(found.crossings.len(), 2);
        for hit in &found.crossings {
            assert!(hit.gap < 1e-9);
            assert!(hit.point.z.abs() < 1e-9);
            assert!((hit.point.to_vector().magnitude() - 2.0).abs() < 1e-9);
        }
    }

    #[test]
    fn a_prepared_surface_answers_as_the_one_shot_query() {
        // A half-pipe spline strip, which only the seeded path can answer,
        // asked about many lines in turn: the stored sampling must give
        // every one the answer a fresh sampling gives, piercings and all.
        use ogeom_geom::BSplineSurface;
        use ogeom_math::ControlGrid;
        let mut points = Vec::new();
        for i in 0..7 {
            let a = core::f64::consts::PI * f64::from(i) / 6.0;
            for j in 0..2 {
                points.push(Point::new(2.0 * a.cos(), 3.0 * f64::from(j), 2.0 * a.sin()));
            }
        }
        let strip: SurfaceGeometry = BSplineSurface::new(
            KnotVector::clamped_uniform(3, 7).unwrap(),
            KnotVector::clamped_uniform(1, 2).unwrap(),
            &ControlGrid::new(points, 7, 2).unwrap(),
            T,
        )
        .unwrap()
        .into();
        let options = CurveSurfaceOptions::default();
        let prepared = PreparedSurface::new(strip.clone(), options, T);
        let mut pierced = 0;
        for k in 0..12 {
            let x = -2.5 + 5.0 * f64::from(k) / 11.0;
            let lines: [Curve; 2] = [
                LineCurve::new(ogeom_math::Axis::new(
                    Point::new(x, 1.5, 0.0),
                    Direction::new(Vector::new(0.1, 0.2, 1.0), T).unwrap(),
                ))
                .into(),
                segment(Point::new(x, -1.0, -1.0), Point::new(-x, 4.0, 3.0)),
            ];
            for line in &lines {
                let once = intersect_curve_surface(line, &strip, options, T).unwrap();
                pierced += once.crossings.len();
                assert_eq!(prepared.intersect(line).unwrap(), once, "line {k}");
            }
        }
        assert!(pierced > 0, "some lines pass through the strip");
    }

    #[test]
    fn a_spline_through_a_sphere_is_found_and_polished() {
        // A spline wandering through the ball: piercings with no closed form
        // anywhere, verified implicitly: each reported point is on the
        // sphere to the gap it claims.
        let wander: Curve = BSplineCurve::new(
            KnotVector::new(vec![0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0], 3).unwrap(),
            vec![
                Point::new(-4.0, -1.0, -1.0),
                Point::new(-1.0, 2.0, 1.0),
                Point::new(1.0, -2.0, -1.0),
                Point::new(4.0, 1.0, 1.0),
            ],
            T,
        )
        .unwrap()
        .into();
        let ball = sphere(2.0);
        let found =
            intersect_curve_surface(&wander, &ball, CurveSurfaceOptions::default(), T).unwrap();
        assert!(!found.crossings.is_empty(), "the spline passes through");
        for hit in &found.crossings {
            assert!(hit.gap < 1e-9);
            let SurfaceGeometry::Sphere(s) = &ball else {
                unreachable!()
            };
            assert!(s.sphere().distance_to(hit.point).abs() < 1e-9);
        }
    }

    #[test]
    fn unusable_options_are_refused() {
        let ball = sphere(1.0);
        let ray = segment(Point::new(-5.0, 0.0, 0.0), Point::new(5.0, 0.0, 0.0));
        for options in [
            CurveSurfaceOptions {
                samples: 1,
                ..CurveSurfaceOptions::default()
            },
            CurveSurfaceOptions {
                grid: 1,
                ..CurveSurfaceOptions::default()
            },
            CurveSurfaceOptions {
                gap: 0.0,
                ..CurveSurfaceOptions::default()
            },
        ] {
            assert!(intersect_curve_surface(&ray, &ball, options, T).is_err());
        }
    }

    /// The general path (a torus has no closed form against a line) answers
    /// each tangency once, and a curve lying in the surface as the stretch
    /// it lies along rather than as hundreds of piercings.
    /// A torus far down a line is met where it is: its quartic, solved where
    /// the line starts, has roots hundreds of units out, and the eigenvalues
    /// lost them.
    #[test]
    fn a_torus_far_down_a_line_is_still_met() {
        use ogeom_geom::TorusSurface;
        use ogeom_math::Torus;
        let options = CurveSurfaceOptions::default();
        let tilted = Frame::new(
            Point::new(0.0, 785.0, -140.0),
            Direction::new(Vector::new(0.0, -0.999_390_827, 0.034_899_497), T).unwrap(),
            Direction::X,
            T,
        )
        .unwrap();
        let torus: SurfaceGeometry =
            TorusSurface::new(Torus::new(tilted, 120.0, 12.0, T).unwrap()).into();
        let w = Vector::new(1.0, 1.0, 1.0) / 3.0_f64.sqrt();
        let from = Point::new(-160.0, 531.0, -320.0);
        let line = segment(from, from + w * 800.0);
        let found = intersect_curve_surface(&line, &torus, options, T).unwrap();
        let cells = sample_by(&torus, seeding(&torus, options.grid), T);
        let tree = CellTree::over(&cells, options.gap);
        let general = general(&line, &torus, &cells, &tree, options, T).unwrap();
        assert!(!general.crossings.is_empty());
        assert_eq!(found.crossings.len(), general.crossings.len());
        for (a, b) in found.crossings.iter().zip(&general.crossings) {
            assert!(
                a.point.distance(b.point) < 1e-6,
                "{:?} against {:?}",
                a.point,
                b.point
            );
        }
    }

    /// A line through a torus's middle crosses the tube four times, and one
    /// through a cone's axis crosses its wall twice, at the closed forms'
    /// points.
    #[test]
    fn lines_cross_a_torus_and_a_cone_where_the_closed_forms_say() {
        use ogeom_geom::{ConeSurface, TorusSurface};
        use ogeom_math::{Cone, Torus};
        let options = CurveSurfaceOptions::default();
        let torus: SurfaceGeometry =
            TorusSurface::new(Torus::new(Frame::WORLD, 60.0, 20.0, T).unwrap()).into();
        let across = segment(Point::new(-100.0, 0.0, 0.0), Point::new(100.0, 0.0, 0.0));
        let found = intersect_curve_surface(&across, &torus, options, T).unwrap();
        let xs: Vec<f64> = found.crossings.iter().map(|c| c.point.x).collect();
        assert_eq!(xs.len(), 4, "{xs:?}");
        for (got, want) in xs.iter().zip([-80.0, -40.0, 40.0, 80.0]) {
            assert!((got - want).abs() < 1e-9, "{xs:?}");
        }
        let cone: SurfaceGeometry = ConeSurface::new(
            Cone::new(Frame::WORLD, 10.0, core::f64::consts::FRAC_PI_4, T).unwrap(),
            (0.0, 20.0),
        )
        .unwrap()
        .into();
        let level = segment(Point::new(-50.0, 0.0, 5.0), Point::new(50.0, 0.0, 5.0));
        let found = intersect_curve_surface(&level, &cone, options, T).unwrap();
        let xs: Vec<f64> = found.crossings.iter().map(|c| c.point.x).collect();
        assert_eq!(xs.len(), 2, "{xs:?}");
        for (got, want) in xs.iter().zip([-15.0, 15.0]) {
            assert!((got - want).abs() < 1e-9, "{xs:?}");
        }
    }

    /// A line anchored far from a cone meets it where one anchored beside
    /// it does: the quadratic's constant term grows with the anchor's
    /// distance, and its leading term must not be dropped for that.
    #[test]
    fn a_line_anchored_far_away_still_crosses_a_cone() {
        use ogeom_geom::ConeSurface;
        use ogeom_math::Cone;
        let options = CurveSurfaceOptions::default();
        let cone: SurfaceGeometry = ConeSurface::new(
            Cone::new(Frame::WORLD, 10.0, core::f64::consts::FRAC_PI_4, T).unwrap(),
            (0.0, 20.0),
        )
        .unwrap()
        .into();
        let far = segment(Point::new(-1e5, 0.0, 5.0), Point::new(50.0, 0.0, 5.0));
        let found = intersect_curve_surface(&far, &cone, options, T).unwrap();
        let xs: Vec<f64> = found.crossings.iter().map(|c| c.point.x).collect();
        assert_eq!(xs.len(), 2, "{xs:?}");
        for (got, want) in xs.iter().zip([-15.0, 15.0]) {
            assert!((got - want).abs() < 1e-6, "{xs:?}");
        }
    }

    /// Lines that miss a torus through its hole or along its axis meet
    /// nothing; one grazing the outer equator touches it once, at the
    /// tangent point; a segment ending on the tube keeps its end crossings.
    #[test]
    fn lines_through_the_hole_along_the_axis_and_grazing_a_torus() {
        use ogeom_geom::TorusSurface;
        use ogeom_math::Torus;
        let options = CurveSurfaceOptions::default();
        let tilted = Frame::new(
            Point::new(3.0, -2.0, 1.0),
            Direction::new(Vector::new(0.1, 0.2, 1.0), T).unwrap(),
            Direction::new(Vector::new(1.0, 0.0, -0.1), T).unwrap(),
            T,
        )
        .unwrap();
        let torus: SurfaceGeometry =
            TorusSurface::new(Torus::new(tilted, 60.0, 20.0, T).unwrap()).into();
        let local = |x: f64, y: f64, z: f64| tilted.to_world(Point::new(x, y, z));
        let meets = |from: Point, to: Point| {
            intersect_curve_surface(&segment(from, to), &torus, options, T)
                .unwrap()
                .crossings
        };
        assert!(meets(local(0.0, 0.0, -100.0), local(0.0, 0.0, 100.0)).is_empty());
        assert!(meets(local(5.0, 3.0, -100.0), local(-15.0, 3.0, 100.0)).is_empty());
        let grazing = meets(local(-100.0, 80.0, 0.0), local(100.0, 80.0, 0.0));
        assert_eq!(grazing.len(), 1, "{grazing:?}");
        assert!(grazing[0].point.distance(local(0.0, 80.0, 0.0)) < 1e-6);
        let ends = meets(local(-80.0, 0.0, 0.0), local(80.0, 0.0, 0.0));
        let xs: Vec<f64> = ends.iter().map(|c| tilted.to_local(c.point).x).collect();
        assert_eq!(xs.len(), 4, "{xs:?}");
        for (got, want) in xs.iter().zip([-80.0, -40.0, 40.0, 80.0]) {
            assert!((got - want).abs() < 1e-9, "{xs:?}");
        }
    }

    #[test]
    fn tangencies_come_back_once_and_lying_curves_as_stretches() {
        use ogeom_geom::TorusSurface;
        use ogeom_math::Torus;
        let torus: SurfaceGeometry =
            TorusSurface::new(Torus::new(Frame::WORLD, 60.0, 20.0, T).unwrap()).into();
        let options = CurveSurfaceOptions::default();
        // Along the tube's top: tangent at two points.
        let top = segment(Point::new(-100.0, 0.0, 20.0), Point::new(100.0, 0.0, 20.0));
        let found = intersect_curve_surface(&top, &torus, options, T).unwrap();
        assert_eq!(found.crossings.len(), 2, "{:?}", found.crossings);
        assert!(found.lying.is_empty());
        // Touching the inner equator and crossing the outer twice.
        let inner = segment(Point::new(-100.0, 40.0, 0.0), Point::new(100.0, 40.0, 0.0));
        let found = intersect_curve_surface(&inner, &torus, options, T).unwrap();
        assert_eq!(found.crossings.len(), 3, "{:?}", found.crossings);
        // A parallel of the torus lies in it the whole way round.
        let parallel: Curve = CircleCurve::new(
            Circle::new(
                Frame::new(Point::new(0.0, 0.0, 20.0), Direction::Z, Direction::X, T).unwrap(),
                60.0,
                T,
            )
            .unwrap(),
        )
        .into();
        let found = intersect_curve_surface(&parallel, &torus, options, T).unwrap();
        assert!(found.crossings.is_empty(), "{:?}", found.crossings);
        assert_eq!(found.lying.len(), 1);
        let (from, to) = found.lying[0];
        assert!((to - from).abs() > 6.0, "{from} to {to}");
    }

    /// A wall leaning off a round's torus band, its foot back, cuts the
    /// band along a section that leaves the wall's face tangent to two of
    /// its edges: the top edge, where the band is tangent to the top, and
    /// the wall's seam with the drum, where the band is tangent to the
    /// drum. Each edge touches the torus there, the gap falling to zero
    /// without changing sign, and the touch comes back on its closed-form
    /// point from a seed a millimetre's thousandth along the edge, where a
    /// solve for the gap stops in the valley. A line crossing the band is
    /// not a touch.
    #[test]
    fn an_edge_tangent_to_a_band_touches_it_where_the_band_meets_its_face() {
        use ogeom_geom::TorusSurface;
        use ogeom_math::Torus;
        let band: SurfaceGeometry = TorusSurface::new(
            Torus::new(
                Frame::new(Point::new(4.0, 4.0, 1.5), Direction::Z, -Direction::X, T).unwrap(),
                2.0,
                0.5,
                T,
            )
            .unwrap(),
        )
        .into();
        let lean = 0.1_f64;
        let wall = plane(Point::new(4.0, 4.0, 2.0), Vector::new(0.0, 1.0, -lean));
        let drum = CylinderSurface::new(
            Cylinder::new(
                Frame::new(Point::new(4.0, 4.0, 0.0), Direction::Z, Direction::X, T).unwrap(),
                2.5,
                T,
            )
            .unwrap(),
            (0.0, 2.0),
        )
        .unwrap()
        .into();
        let crate::section::SurfaceIntersection::Along(seams) = crate::section::intersect_surfaces(
            &wall,
            &drum,
            crate::section::IntersectOptions::default(),
            T,
        )
        .unwrap() else {
            panic!("the wall meets the drum");
        };
        let foot = Point::new(4.0 + 2.5_f64.mul_add(2.5, -0.05 * 0.05).sqrt(), 3.95, 1.5);
        // The nearest of a fine sampling: a seed, not an answer.
        let nearest = |c: &Curve, p: Point| {
            let (lo, hi) = c.domain();
            (0..=20_000)
                .map(|k| lo + (hi - lo) * f64::from(k) / 20_000.0)
                .map(|t| (t, c.point_at(t, T).unwrap().distance(p)))
                .min_by(|a, b| a.1.total_cmp(&b.1))
                .unwrap()
        };
        let seam = seams
            .iter()
            .map(|s| &s.curve)
            .find(|c| nearest(c, foot).1 < 1e-2)
            .unwrap();
        let top = segment(Point::new(1.5, 4.0, 2.0), Point::new(6.5, 4.0, 2.0));
        for (edge, want) in [(&top, Point::new(6.0, 4.0, 2.0)), (seam, foot)] {
            let at = nearest(edge, want).0;
            let step = 1e-3 / edge.d1_at(at, T).unwrap().magnitude();
            for seed in [at - step, at + step] {
                let touch = touch_curve_surface(edge, &band, seed, T).unwrap();
                assert!(
                    touch.point.distance(want) < 1e-9 && touch.gap < 1e-9,
                    "{touch:?} against {want:?}"
                );
            }
        }
        let across = segment(Point::new(6.2, 4.0, 3.0), Point::new(6.2, 4.0, 0.0));
        let found =
            intersect_curve_surface(&across, &band, CurveSurfaceOptions::default(), T).unwrap();
        assert_eq!(found.crossings.len(), 2, "{:?}", found.crossings);
        for crossing in found.crossings {
            assert!(touch_curve_surface(&across, &band, crossing.on_curve, T).is_none());
        }
    }
}