axiolid_brep_boolean/

section.rs

//! Section edges: where two exact B-reps' faces cross (ADR 0075, step 2).
//!
//! For every pair of faces, one from each operand:
//!
//! 1. The support surfaces' exact intersection curves come from
//!    `exact_surface_intersection` (#119): lines and conics, ruled and
//!    torus sections (ADR 0076), traced sections (ADR 0077). A ruled or
//!    torus section piece is read as a traced section over its span, so
//!    every curve past this point is a line, a conic or an
//!    `ImplicitSection3`.
//! 2. Each curve is cut where it crosses a boundary edge of either face.
//!    Two curves on one surface meet only up to the rounding of their
//!    constructed doubles, so the crossing is found transversally instead:
//!    where the curve crosses the ADJACENT face's surface
//!    (`exact_curve_surface_intersection`), kept when the hit lies on the
//!    edge within tolerance and inside its span. Across a seam, where the
//!    adjacent face shares the surface, the cutter is the plane through the
//!    ruling and the surface normal.
//! 3. Each piece between consecutive cuts lies wholly inside or wholly
//!    outside each face, so its midpoint decides it, through the face's
//!    certified domain (`FaceDomain`). A piece inside both faces is a section
//!    edge.
//!
//! Faces that meet without crossing are handled, not refused:
//!
//! - **Tangent supports** touch in a point or along a curve without
//!   crossing, so they add no section: neither face changes sides there.
//! - **A section along an existing edge** (the curve lies in the adjacent
//!   face's surface too) is cut at the edge's ends and marked as lying
//!   along that face's boundary: it splits the other face, not this one. A
//!   traced piece along a line or conic edge is that edge's own curve over
//!   its stretch, so one stretch has one curve on every face (#382: a
//!   revolved torus wall's cap edge against the plane through its axis).
//! - **Coincident supports** share a patch of surface. Each face's boundary
//!   edges are imprinted on the other where they run inside it, so the
//!   shared patch becomes a region of both faces; classification then sees
//!   it on the other solid's boundary (see `crate::boolean`). An imprinted
//!   edge's cuts are merged as a section's, and a cut within tolerance of
//!   the edge's own end is that end (#276, `imprint_cuts`).
//!
//! Work that cannot matter is skipped through sound boxes (`crate::bounds`,
//! enlarged by the tolerance): face pairs whose boxes are apart are not
//! sectioned, and a line or conic section is cut only against edges whose
//! boxes meet the two faces' common box (widened by `2^-20` of its
//! diagonal, so its planes stay off the faces' extremes), plus where it
//! leaves that box; pieces outside it lie off one of the faces.
//!
//! Operands placed by independent rigid motions agree only up to rounding,
//! and the exact predicates above see that residue (#228). Within
//! tolerance, then:
//!
//! - a plane parallel or perpendicular to a cylinder's axis cuts rulings
//!   or a circle, not a degenerate ellipse;
//! - a curve running along a boundary edge is read as running along it;
//! - a curve that touches the edge's adjacent surface (a double root,
//!   which rounding splits in two or loses) is cut where it meets the edge
//!   itself, and cuts closer together than tolerance are one cut;
//! - a plane touching a cylinder is read once for the two supports, and a
//!   curve on either one is cut by the other where it meets their contact
//!   ruling, in every face pair alike (`crate::contact`, #243).
//!
//! Along a smooth edge (a fillet running tangent into its neighbour) a
//! section is cut where the edge crosses the section's other surface, and
//! a trace that cannot be decided in a face's box of a periodic carrier is
//! taken over its whole turn (`crate::tangency`, #249).

use axiolid_brep::ExactBRep;
use axiolid_core::{Interval, Point2, Point3, Scalar, Tolerance, Vec3};
use axiolid_curve::Curve3;
use axiolid_evaluate::surface::{locate, normal};
use axiolid_evaluate::{curve::locate3, evaluate3};
use axiolid_measure::FaceDomain;
use axiolid_nurbs::{
    exact_curve_curve_intersection3, exact_curve_surface_intersection, exact_surface_intersection,
    spline_pair_intersection, ExactCurveIntersection, ExactIntersectionRefusal,
};
use axiolid_surface::{Plane, Surface};
use axiolid_topology::FaceId;
use core::f64::consts::TAU;

use crate::bounds::{edge_box, face_box, Aabb};
use crate::contact::{Contact, Contacts};
use crate::predicate;
use crate::report::{self, Reading, ToleranceDecisionKind};
use crate::support::{cleaned, same_support, window};
use crate::BooleanError;

/// Where along a piece its inside/outside sample is taken: off-centre, so
/// the midpoint of a symmetric section (a meridian's pole) is never it.
const SAMPLE: Scalar = 0.414_213_562_373_095;

/// How far a section's window reaches past the two faces' common box,
/// relative to its diagonal: `2^-20`. Widening is sound (the window only
/// skips work outside it), and keeps the window's planes off the faces'
/// extremes, where a tangent section would leave it next to the touching
/// point (#243).
const WINDOW_PAD: Scalar = 1.0 / 1_048_576.0;

/// A piece of an intersection curve lying on a face of each operand.
#[derive(Debug, Clone, PartialEq)]
pub struct SectionEdge {
    /// The face of the first operand the edge lies on.
    pub face_a: FaceId,
    /// The face of the second operand the edge lies on.
    pub face_b: FaceId,
    /// The exact curve.
    pub curve: Curve3,
    /// The curve's parameter span; for a circle or ellipse it may run past
    /// `2 pi`.
    pub span: Interval,
    /// The curve's point at `span.start`.
    pub start: Point3,
    /// The curve's point at `span.end`.
    pub end: Point3,
    /// Whether the edge runs along a boundary edge `face_a` already has, so
    /// it does not split `face_a`.
    pub along_a: bool,
    /// Whether the edge runs along a boundary edge `face_b` already has.
    pub along_b: bool,
    /// The surface that meets `face_a`'s surface along this edge: `face_b`'s
    /// for a crossing, or for an edge imprinted from a coincident face the
    /// surface that bounds it there. Its equation on `face_a` is the edge's
    /// pcurve.
    pub other_a: Surface,
    /// The surface that meets `face_b`'s surface along this edge.
    pub other_b: Surface,
}

/// Every section edge between the faces of `a` and the faces of `b`.
///
/// # Errors
///
/// A face pair whose section is not a line, circle or ellipse, or a piece
/// too close to a face boundary to classify.
pub fn section_edges(
    a: &ExactBRep,
    b: &ExactBRep,
    tolerance: Tolerance,
) -> Result<Vec<SectionEdge>, BooleanError> {
    let session = report::open(&[a, b], tolerance);
    let (edges, _) = sections_and_contacts(a, b, tolerance)?;
    session.finish()?;
    Ok(edges)
}

/// [`section_edges`], with the plane/cylinder contacts read within
/// tolerance on the way (`crate::contact`), which classification must
/// respect too.
pub(crate) fn sections_and_contacts(
    a: &ExactBRep,
    b: &ExactBRep,
    tolerance: Tolerance,
) -> Result<(Vec<SectionEdge>, Contacts), BooleanError> {
    let side_a = Side::new(a, tolerance)?;
    let side_b = Side::new(b, tolerance)?;
    let mut out = Vec::new();
    // Each plane/cylinder pair of supports is read once, over the largest
    // window any pair of their faces needs, so every face pair of the two
    // agrees (#243).
    let (readings, contacts) = plane_cylinder_readings(&side_a, &side_b, tolerance)?;
    // Surfaces are shared by many faces: each pair's closed form once.
    #[allow(clippy::type_complexity)]
    let mut closed_forms: Vec<(
        Surface,
        Surface,
        Result<axiolid_nurbs::ExactIntersectionCurve, ExactIntersectionRefusal>,
    )> = Vec::new();
    for fa in 0..side_a.faces.len() {
        for fb in 0..side_b.faces.len() {
            // Faces whose boxes, enlarged by the tolerance, are apart cannot
            // cross, touch or share a patch (#228).
            // Where both faces have boxes, the section matters only in their
            // common part: the window. It is widened a little past the
            // faces' boxes, whose planes run through the faces' extremes: a
            // section touching one there (a hole tangent to a face) would
            // leave the window a hair from the touching point, and the
            // piece between could not be told inside or outside (#243).
            let common = match (&side_a.boxes[fa], &side_b.boxes[fb]) {
                (Some(a), Some(b)) => match a.intersection(b) {
                    Some(common) => Some(common.widened(WINDOW_PAD * common.diagonal())),
                    None => continue,
                },
                _ => None,
            };
            let (sa, sb) = (side_a.surface(fa)?, side_b.surface(fb)?);
            if same_support(sa, sb, tolerance) {
                imprint(
                    &side_a, fa, &side_b, fb, true, &contacts, tolerance, &mut out,
                )?;
                imprint(
                    &side_b, fb, &side_a, fa, false, &contacts, tolerance, &mut out,
                )?;
                continue;
            }
            // Lines and conics in closed form; every other section traced in
            // one face's parameter box (ADR 0077), which holds every part
            // of it that can matter.
            let splines = matches!((sa, sb), (Surface::BSpline(_), Surface::BSpline(_)));
            // A plane parallel or perpendicular to a cylinder's axis to
            // within rounding, as read for the two supports.
            let rulings = readings
                .iter()
                .find(|(x, y, _)| x == sa && y == sb)
                .and_then(|(_, _, reading)| reading.as_ref())
                .map(|reading| reading.curves.clone());
            if rulings.as_ref().is_some_and(Vec::is_empty) {
                continue;
            }
            // Two B-splines are traced below, in the faces' own boxes.
            let closed_form = if splines || rulings.is_some() {
                None
            } else {
                let known = closed_forms
                    .iter()
                    .find(|(a, b, _)| a == sa && b == sb)
                    .map(|(_, _, r)| r.clone());
                let result = known.unwrap_or_else(|| {
                    let r = exact_surface_intersection(&cleaned(sa), &cleaned(sb));
                    closed_forms.push((sa.clone(), sb.clone(), r.clone()));
                    r
                });
                match result {
                    Ok(curve) => {
                        let conic = curve.branches.iter().zip(&curve.spans).all(|(b, s)| {
                            matches!(
                                (b, s),
                                (Curve3::Line(_), _)
                                    | (Curve3::Circle(_) | Curve3::Ellipse(_), None)
                            )
                        });
                        conic.then_some(curve)
                    }
                    // Apart, or touching without crossing: no section.
                    Err(
                        ExactIntersectionRefusal::Disjoint
                        | ExactIntersectionRefusal::NotRegularCurve,
                    ) => continue,
                    Err(_) => None,
                }
            };
            let branches: Vec<(Curve3, Option<Interval>)> = match (rulings, closed_form) {
                (Some(lines), _) => lines.into_iter().map(|line| (line, None)).collect(),
                (None, Some(curve)) => curve.branches.into_iter().zip(curve.spans).collect(),
                (None, None) => {
                    // Only a B-spline can carry a section with a B-spline
                    // (it has no equation to read elsewhere); otherwise the
                    // compact and low-degree surfaces carry best.
                    let rank = |s: &Surface| match s {
                        Surface::BSpline(_) => 4,
                        Surface::Torus(_) => 3,
                        Surface::Sphere(_) => 2,
                        Surface::Plane(_) => 0,
                        _ => 1,
                    };
                    // Two B-splines: the section carried on both, traced in
                    // both faces' boxes (ADR 0077).
                    if let (Surface::BSpline(ba), Surface::BSpline(bb)) = (sa, sb) {
                        let (la, ha) = side_a.domains[fa].bounds();
                        let (lb, hb) = side_b.domains[fb].bounds();
                        match spline_pair_intersection(
                            ba,
                            bb,
                            Some((window(sa, la, ha), window(sb, lb, hb))),
                        ) {
                            Ok(sections) => sections
                                .into_iter()
                                .map(|s| {
                                    let end = s.end();
                                    (Curve3::PairSection(s), Some(Interval::new(0.0, end)))
                                })
                                .collect(),
                            Err(ExactIntersectionRefusal::Disjoint) => continue,
                            Err(_) => return Err(BooleanError::UnsupportedSection),
                        }
                    } else {
                        let (carrier, other, side, face) = if rank(sa) >= rank(sb) {
                            (sa, sb, &side_a, fa)
                        } else {
                            (sb, sa, &side_b, fb)
                        };
                        let (lo, hi) = side.domains[face].bounds();
                        match crate::tangency::trace(carrier, other, lo, hi) {
                            Ok(sections) => sections
                                .into_iter()
                                .map(|s| {
                                    let end = s.curve.end();
                                    (Curve3::ImplicitSection(s), Some(Interval::new(0.0, end)))
                                })
                                .collect(),
                            // Apart, or touching only at isolated points: no
                            // section, as for the closed forms.
                            Err(
                                ExactIntersectionRefusal::Disjoint
                                | ExactIntersectionRefusal::NotRegularCurve,
                            ) => continue,
                            Err(_) => return Err(BooleanError::UnsupportedSection),
                        }
                    }
                }
            };
            for (index, (branch, span)) in branches.iter().enumerate() {
                // A ray (a cone's ruling from its apex) is a line bounded
                // at its finite ends.
                let bounds = match (branch, span) {
                    (Curve3::Line(_), Some(span)) => Some(*span),
                    _ => None,
                };
                let branch = branch.clone();
                let branch = &branch;
                // The window cuts lines and conics in closed form; a traced
                // section is cut against every edge, as before.
                let common = common.filter(|_| {
                    matches!(
                        branch,
                        Curve3::Line(_) | Curve3::Circle(_) | Curve3::Ellipse(_)
                    )
                });
                let (mut cuts, along_a) =
                    side_a.cuts(fa, branch, sb, common.as_ref(), &contacts, tolerance)?;
                let (more, along_b) =
                    side_b.cuts(fb, branch, sa, common.as_ref(), &contacts, tolerance)?;
                cuts.extend(more);
                // Branches of one section meet only where the surfaces touch
                // (the two ellipses of a Steinmetz pair): such a point splits
                // both, so the faces' graphs get a vertex there.
                for (other, _) in branches
                    .iter()
                    .skip(index + 1)
                    .chain(branches.iter().take(index))
                {
                    if let Ok(ExactCurveIntersection::Points(hits)) =
                        exact_curve_curve_intersection3(branch, other)
                    {
                        cuts.extend(hits.iter().map(|h| h.parameter.approx()));
                    }
                }
                if let Some(b) = bounds {
                    let (lo, hi) = (b.start.min(b.end), b.start.max(b.end));
                    cuts.retain(|&c| c >= lo && c <= hi);
                    cuts.extend([lo, hi].into_iter().filter(|x| x.is_finite()));
                }
                let mut cuts = merge_close(branch, cuts, tolerance)?;
                // The window's own crossings come after the merge: they are
                // where the window ends, not features to be merged.
                if let Some(w) = &common {
                    let ends = window_ends(&cuts, w.crossings(branch));
                    cuts.extend(ends.into_iter().filter(|&c| {
                        bounds.is_none_or(|b| c >= b.start.min(b.end) && c <= b.start.max(b.end))
                    }));
                }
                for (piece_curve, span) in pieces(branch, cuts)? {
                    // Off-centre, so a symmetric section's pole or seam
                    // crossing never becomes the sample.
                    let mid = sample(&piece_curve, span, [sa, sb], &contacts)?;
                    // Tangent contact adds no section. A traced section is
                    // never mere touching (the trace drops touching points
                    // and curves), so where its surfaces are tangent they
                    // cross: it stays.
                    let traced = matches!(
                        piece_curve,
                        Curve3::ImplicitSection(_) | Curve3::PairSection(_)
                    );
                    // Outside the window the piece is off one of the faces.
                    if common.as_ref().is_some_and(|w| !w.contains(mid)) {
                        continue;
                    }
                    if !traced && touching(sa, sb, mid, tolerance)? {
                        continue;
                    }
                    let at_a = side_a.locate(fa, mid, &along_a, tolerance)?;
                    if at_a == Place::Outside {
                        continue;
                    }
                    let at_b = side_b.locate(fb, mid, &along_b, tolerance)?;
                    if at_b == Place::Outside {
                        continue;
                    }
                    // A traced piece along a line or conic edge is that edge's
                    // own curve over its stretch (#382).
                    let edges: Vec<&(Curve3, Interval)> = along_a
                        .iter()
                        .filter(|_| at_a == Place::Boundary)
                        .chain(along_b.iter().filter(|_| at_b == Place::Boundary))
                        .collect();
                    let (piece_curve, span) =
                        along_edge_curve(piece_curve, span, mid, &edges, tolerance)?;
                    out.push(SectionEdge {
                        face_a: side_a.faces[fa],
                        face_b: side_b.faces[fb],
                        start: evaluate3(&piece_curve, span.start)
                            .map_err(|_| BooleanError::Evaluation)?,
                        end: evaluate3(&piece_curve, span.end)
                            .map_err(|_| BooleanError::Evaluation)?,
                        curve: piece_curve,
                        span,
                        along_a: at_a == Place::Boundary,
                        along_b: at_b == Place::Boundary,
                        other_a: sb.clone(),
                        other_b: sa.clone(),
                    });
                }
            }
        }
    }
    Ok((out, contacts))
}

/// One reading per pair of supports, a plane of one operand and a cylinder
/// of the other ([`plane_cylinder_within_rounding`]), over the largest
/// common box of any pair of their faces (none when a pair has no box),
/// with the contacts it read within tolerance.
#[allow(clippy::type_complexity)]
fn plane_cylinder_readings(
    side_a: &Side<'_>,
    side_b: &Side<'_>,
    tolerance: Tolerance,
) -> Result<(Vec<(Surface, Surface, Option<PlaneCylinder>)>, Contacts), BooleanError> {
    let mut extents: Vec<(Surface, Surface, Option<Scalar>)> = Vec::new();
    for fa in 0..side_a.faces.len() {
        for fb in 0..side_b.faces.len() {
            let (sa, sb) = (side_a.surface(fa)?, side_b.surface(fb)?);
            if !matches!(
                (sa, sb),
                (Surface::Plane(_), Surface::Cylinder(_))
                    | (Surface::Cylinder(_), Surface::Plane(_))
            ) {
                continue;
            }
            let extent = match (&side_a.boxes[fa], &side_b.boxes[fb]) {
                (Some(a), Some(b)) => match a.intersection(b) {
                    Some(common) => Some(common.diagonal()),
                    None => continue,
                },
                _ => None,
            };
            match extents.iter_mut().find(|(x, y, _)| x == sa && y == sb) {
                Some((_, _, known)) => *known = known.zip(extent).map(|(p, q)| p.max(q)),
                None => extents.push((sa.clone(), sb.clone(), extent)),
            }
        }
    }
    let mut contacts = Contacts::default();
    let mut readings = Vec::with_capacity(extents.len());
    for (sa, sb, extent) in extents {
        let reading = plane_cylinder_within_rounding(&sa, &sb, extent, tolerance);
        if let Some(PlaneCylinder {
            curves,
            contact: true,
        }) = &reading
        {
            if let [Curve3::Line(ruling)] = curves.as_slice() {
                let (plane, cylinder) = match sa {
                    Surface::Plane(_) => (sa.clone(), sb.clone()),
                    _ => (sb.clone(), sa.clone()),
                };
                contacts.push(Contact {
                    plane,
                    cylinder,
                    ruling: *ruling,
                });
            }
        }
        readings.push((sa, sb, reading));
    }
    Ok((readings, contacts))
}

/// The point of a piece that decides whether it lies in a face: off-centre
/// at [`SAMPLE`], or, where that point lies where a contact one of the
/// piece's two `surfaces` is read in is disputed (`crate::contact`), the
/// first point further along that is not. There the given operands and the
/// moved ones the reading stands for put the point on different sides of
/// the plane, and the piece was cut for the moved ones (#243).
fn sample(
    curve: &Curve3,
    span: Interval,
    surfaces: [&Surface; 2],
    contacts: &Contacts,
) -> Result<Point3, BooleanError> {
    let at = |f: Scalar| {
        evaluate3(curve, span.start + f * (span.end - span.start))
            .map_err(|_| BooleanError::Evaluation)
    };
    let disputed = |p: Point3| surfaces.iter().any(|s| contacts.disputed(s, p));
    let first = at(SAMPLE)?;
    if !disputed(first) {
        return Ok(first);
    }
    for f in [0.75, 0.9, 0.25, 0.1] {
        let p = at(f)?;
        if !disputed(p) {
            return Ok(p);
        }
    }
    Ok(first)
}

/// The section of a plane and a cylinder whose axis the plane is parallel
/// or perpendicular to within the linear tolerance over where the section
/// can matter, or `None` for any other pair.
///
/// The exact closed form decides parallel and perpendicular exactly, on the
/// numbers given. Two operands placed by independent rigid motions -- an
/// arched opening's jamb plane against its own soffit cylinder, or a wall's
/// face against it, after a general rotation -- leave the normal and axis
/// `1e-17` off, and the exact section is then an ellipse `1e16` long, or
/// tilted by the residue, that nothing can evaluate. Read within rounding,
/// the section is what the operands describe:
///
/// - parallel: two rulings, one where the plane touches the cylinder
///   (within linear tolerance), or none (an empty list);
/// - perpendicular: the circle where the axis pierces the plane.
///
/// Each curve lies on the cylinder exactly. It lies on a copy of the plane
/// turned (and, for a touching ruling, moved) by at most the caller's
/// linear tolerance `eps` over the region it can matter in, so the section
/// is the exact section of operands perturbed by at most `eps`:
///
/// - perpendicular when `r sin(theta) <= eps`, `theta` the angle between
///   the normal and the axis: the circle leaves the plane by at most
///   `r sin(theta)`;
/// - parallel when `|n . a| * extent <= eps`, `extent` the diagonal of the
///   largest common box of a pair of the two supports' faces (`None`, a
///   pair with no box: never read parallel; one reading per pair of
///   supports, so all their face pairs agree, #243): a
///   ruling leaves the plane by at most that over the box;
/// - one ruling when the plane is within `eps` of touching the cylinder.
///
/// Outside those bounds the exact closed form decides (#228). Exactly
/// parallel or perpendicular (`crate::predicate`), it decides too, unless
/// the plane is read as touching (#236).
fn plane_cylinder_within_rounding(
    a: &Surface,
    b: &Surface,
    extent: Option<Scalar>,
    tolerance: Tolerance,
) -> Option<PlaneCylinder> {
    let (cylinder, plane) = match (a, b) {
        (Surface::Cylinder(c), Surface::Plane(p)) | (Surface::Plane(p), Surface::Cylinder(c)) => {
            (c, p)
        }
        _ => return None,
    };
    let axis = cylinder.frame.z.normalize_or_zero();
    let normal = plane.frame.z.normalize_or_zero();
    if axis == Vec3::ZERO || normal == Vec3::ZERO {
        return None;
    }
    let along = axis.dot(normal);
    // Exactly perpendicular or parallel is the closed form's to decide; only
    // a reading the exact predicates do not confirm is taken here, and
    // reported (#236). These readings bound how far the plane moves, not how
    // far it turns, so only the linear measure gates them.
    let linear = Tolerance::new(tolerance.linear(), Scalar::MAX).unwrap_or(tolerance);
    match report::support(
        ToleranceDecisionKind::PlanePerpendicularToAxis,
        cylinder.radius * axis.cross(normal).length(),
        0.0,
        linear,
        || predicate::plane_perpendicular(plane, cylinder),
    ) {
        Reading::Exact => return None,
        Reading::Apart => {}
        Reading::Within => {
            // Perpendicular to the axis within rounding: the circle where the
            // axis pierces the plane, not an ellipse tilted by the residue.
            let centre = cylinder.frame.origin
                + axis * (normal.dot(plane.frame.origin - cylinder.frame.origin) / along);
            let x = (cylinder.frame.x - axis * axis.dot(cylinder.frame.x)).normalize_or_zero();
            if x == Vec3::ZERO || !centre.is_finite() {
                return None;
            }
            return Some(PlaneCylinder {
                curves: vec![Curve3::Circle(axiolid_curve::Circle3 {
                    frame: axiolid_core::Frame3 {
                        origin: centre,
                        x,
                        y: axis.cross(x),
                        z: axis,
                    },
                    radius: cylinder.radius,
                })],
                contact: false,
            });
        }
    }
    // No common box: never read parallel.
    let extent = extent?;
    let parallel = report::support(
        ToleranceDecisionKind::PlaneParallelToAxis,
        along.abs() * extent,
        0.0,
        linear,
        || predicate::plane_parallel(plane, cylinder),
    );
    if parallel == Reading::Apart {
        return None;
    }
    // The normal with its rounding along the axis removed.
    let normal = (normal - axis * along).normalize_or_zero();
    if normal == Vec3::ZERO {
        return None;
    }
    let distance = normal.dot(cylinder.frame.origin - plane.frame.origin);
    let foot = cylinder.frame.origin - normal * distance;
    let across = axis.cross(normal);
    let r = cylinder.radius;
    let touching = report::support(
        ToleranceDecisionKind::PlaneTouchesCylinder,
        (distance.abs() - r).abs(),
        0.0,
        linear,
        || predicate::plane_touches(plane, cylinder),
    );
    // Exactly parallel, the closed form decides, unless the plane is read
    // as touching within tolerance.
    if parallel == Reading::Exact && touching != Reading::Within {
        return None;
    }
    let offsets: Vec<Scalar> = if touching.holds() {
        vec![0.0]
    } else if distance.abs() > r {
        Vec::new()
    } else {
        let half = (r * r - distance * distance).sqrt();
        vec![half, -half]
    };
    Some(PlaneCylinder {
        curves: offsets
            .into_iter()
            .map(|offset| {
                Curve3::Line(axiolid_curve::Line3 {
                    origin: foot + across * offset,
                    direction: axis,
                })
            })
            .collect(),
        contact: touching == Reading::Within,
    })
}

/// A plane/cylinder section read within rounding or tolerance
/// ([`plane_cylinder_within_rounding`]).
#[derive(Debug, Clone)]
struct PlaneCylinder {
    /// The rulings, or the circle, of the section; empty when the plane
    /// misses the cylinder.
    curves: Vec<Curve3>,
    /// Whether the plane was read as touching the cylinder within
    /// tolerance (not exactly): `curves` is then the contact ruling, which
    /// every other crossing of the two surfaces must agree with
    /// (`crate::contact`, #243).
    contact: bool,
}

/// The points where a surface's angle parameter has no value: a sphere's
/// poles, a cone's apex.
fn poles(surface: &Surface) -> Vec<Point3> {
    match surface {
        Surface::Sphere(s) => {
            let z = s.frame.z.normalize() * s.radius;
            vec![s.frame.origin + z, s.frame.origin - z]
        }
        Surface::Cone(c) => {
            let slope = c.semi_angle.tan();
            if slope == 0.0 {
                Vec::new()
            } else {
                vec![c.frame.origin - c.frame.z.normalize() * (c.radius / slope)]
            }
        }
        _ => Vec::new(),
    }
}

/// Whether two surfaces through `point` share their tangent plane there: a
/// line where a cylinder rests on a plane, which it touches without
/// crossing. Analytic surfaces tangent along a whole curve lie on one side
/// of each other, so such a stretch splits neither face.
fn touching(
    a: &Surface,
    b: &Surface,
    point: Point3,
    tolerance: Tolerance,
) -> Result<bool, BooleanError> {
    let at = |s: &Surface| -> Result<axiolid_core::Vec3, BooleanError> {
        let (u, v) =
            locate(s, point, report::floored(tolerance)).map_err(|_| BooleanError::Evaluation)?;
        Ok(normal(s, u, v)
            .map_err(|_| BooleanError::Evaluation)?
            .normalize())
    };
    Ok(report::near_angle(
        ToleranceDecisionKind::Contact,
        at(a)?.cross(at(b)?).length(),
        tolerance,
    ))
}

/// The boundary edges of `from`'s face `face` that run inside `onto`'s face
/// `other`, on the same surface, as section edges (`first` when `from` is
/// the first operand).
#[allow(clippy::too_many_arguments)]
fn imprint(
    from: &Side<'_>,
    face: usize,
    onto: &Side<'_>,
    other: usize,
    first: bool,
    contacts: &Contacts,
    tolerance: Tolerance,
    out: &mut Vec<SectionEdge>,
) -> Result<(), BooleanError> {
    for (edge, curve, span) in from.edges_of(face)? {
        // The surface that bounds the imprinted edge in its own operand:
        // together with the shared surface it defines the edge's curve.
        let bounding = from.cutter(face, edge, &curve, span, tolerance, false)?;
        let (cuts, along) = onto.cuts(other, &curve, &bounding, None, contacts, tolerance)?;
        let (cuts, span) = imprint_cuts(&curve, span, cuts, tolerance)?;
        let shared = onto.surface(other)?;
        for piece in pieces_within(&curve, span, cuts) {
            let mid = sample(&curve, piece, [shared, &bounding], contacts)?;
            let at = onto.locate(other, mid, &along, tolerance)?;
            if at == Place::Outside {
                continue;
            }
            let (face_a, face_b, along_a, along_b) = if first {
                (
                    from.faces[face],
                    onto.faces[other],
                    true,
                    at == Place::Boundary,
                )
            } else {
                (
                    onto.faces[other],
                    from.faces[face],
                    at == Place::Boundary,
                    true,
                )
            };
            out.push(SectionEdge {
                face_a,
                face_b,
                curve: curve.clone(),
                span: piece,
                start: evaluate3(&curve, piece.start).map_err(|_| BooleanError::Evaluation)?,
                end: evaluate3(&curve, piece.end).map_err(|_| BooleanError::Evaluation)?,
                along_a,
                along_b,
                other_a: bounding.clone(),
                other_b: bounding.clone(),
            });
        }
    }
    Ok(())
}

/// Where a point on a face's support lies relative to the face.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum Place {
    Inside,
    Boundary,
    Outside,
}

/// Where `curve` crosses `edge` within tolerance, as a parameter on
/// `curve`: two lines at their closest approach, or a line and a circle or
/// ellipse where the line pierces the conic's plane. `None` when they stay
/// further apart than tolerance, run parallel (a line in the conic's
/// plane, parallel lines), or are another pair of families.
///
/// Both pairs are well conditioned where they cross transversally in
/// space, which is the case a tangent contact with the edge's adjacent
/// surface leaves (#228).
fn near_crossing(
    curve: &Curve3,
    edge: &Curve3,
    tolerance: Tolerance,
) -> Result<Option<Scalar>, BooleanError> {
    let conic_plane = |c: &Curve3| match c {
        Curve3::Circle(c) => Some(c.frame),
        Curve3::Ellipse(e) => Some(e.frame),
        _ => None,
    };
    let point = match (curve, edge) {
        (Curve3::Line(a), Curve3::Line(b)) => {
            let (d1, d2) = (a.direction, b.direction);
            let r = a.origin - b.origin;
            let (aa, bb, ab) = (d1.dot(d1), d2.dot(d2), d1.dot(d2));
            let denominator = aa * bb - ab * ab;
            if denominator <= 1e-12 * aa * bb {
                return Ok(None);
            }
            let (c, f) = (d1.dot(r), d2.dot(r));
            let s = (ab * f - c * bb) / denominator;
            let t = (aa * f - ab * c) / denominator;
            let (p, q) = (a.origin + d1 * s, b.origin + d2 * t);
            if !report::near(
                ToleranceDecisionKind::TangentCrossing,
                (p - q).length(),
                tolerance,
            ) {
                return Ok(None);
            }
            return Ok(Some(s));
        }
        (Curve3::Line(line), conic) | (conic, Curve3::Line(line)) => {
            let Some(frame) = conic_plane(conic) else {
                return Ok(None);
            };
            let n = frame.z.normalize();
            let along = line.direction.dot(n);
            if along.abs() <= 1e-9 * line.direction.length() {
                return Ok(None);
            }
            let s = (frame.origin - line.origin).dot(n) / along;
            line.origin + line.direction * s
        }
        _ => return Ok(None),
    };
    // The pierce point must lie on the conic, within tolerance.
    let conic = if matches!(curve, Curve3::Line(_)) {
        edge
    } else {
        curve
    };
    let wide = report::floored(tolerance);
    let Ok(t) = locate3(conic, point, wide) else {
        return Ok(None);
    };
    let on = evaluate3(conic, t).map_err(|_| BooleanError::Evaluation)?;
    if !report::near(
        ToleranceDecisionKind::TangentCrossing,
        (on - point).length(),
        tolerance,
    ) {
        return Ok(None);
    }
    match locate3(curve, point, wide) {
        Ok(t) => Ok(Some(t)),
        Err(_) => Ok(None),
    }
}

/// Cuts on a line, conic or traced section that lie within tolerance of
/// each other in space, merged into one at their mean parameter.
///
/// A curve crossing a face's edge where it touches the edge's other
/// surface -- a wall's face meeting an arched opening's jamb where the jamb
/// runs tangent into the arch -- has a double root there, and the exact
/// intersection returns it as two roots about `sqrt(eps)` apart. Kept
/// apart, they leave a sliver piece `1e-8` long whose ends no other piece
/// meets, and the face does not split (#228). A traced section does the same
/// where it touches a face's edge: a plane at a ring's radius touches the
/// latitude circle where a revolved circle's quarter faces meet, and the
/// trace's double root splits into roots `2e-8` apart (#381; the trace's
/// spurious roots are TODO(#397)). A traced loop is not merged across its
/// start. Other curves are returned unchanged.
///
/// Each merged cut moves by at most the caller's linear tolerance. A chain
/// of cuts each within tolerance of the next but spanning more than it
/// would move one of them further, so it is refused
/// ([`BooleanError::NearCoincidence`]).
fn merge_close(
    curve: &Curve3,
    cuts: Vec<Scalar>,
    tolerance: Tolerance,
) -> Result<Vec<Scalar>, BooleanError> {
    let periodic = matches!(curve, Curve3::Circle(_) | Curve3::Ellipse(_));
    if !(periodic || matches!(curve, Curve3::Line(_) | Curve3::ImplicitSection(_))) {
        return Ok(cuts);
    }
    // Cuts within rounding of each other are one cut and need no tolerance;
    // further apart, merging them is a reading within tolerance (#236).
    let eps = tolerance.linear().max(report::rounding());
    let merged = |gap: Scalar| report::near(ToleranceDecisionKind::MergedPoints, gap, tolerance);
    let at = |t: Scalar| evaluate3(curve, t).map_err(|_| BooleanError::Evaluation);
    let mut placed: Vec<(Scalar, Point3)> = Vec::with_capacity(cuts.len());
    for cut in cuts {
        if !cut.is_finite() {
            continue;
        }
        placed.push((cut, at(cut)?));
    }
    placed.sort_by(|a, b| a.0.total_cmp(&b.0));
    // (parameter sum, count, first parameter, first point, last point) per
    // cluster.
    let mut clusters: Vec<(Scalar, Scalar, Scalar, Point3, Point3)> = Vec::new();
    for (cut, point) in placed {
        match clusters.last_mut() {
            Some((sum, count, start, first, last)) if (point - *last).length() <= eps => {
                if !merged((point - *first).length()) {
                    return Err(BooleanError::NearCoincidence);
                }
                // On a closed conic one point has parameters a whole turn
                // apart (a seam cut at `0` and at `2 pi`): averaged as they
                // are, they would name the opposite side (#243).
                let cut = if periodic {
                    cut + TAU * ((*start - cut) / TAU).round()
                } else {
                    cut
                };
                *sum += cut;
                *count += 1.0;
                *last = point;
            }
            _ => clusters.push((cut, 1.0, cut, point, point)),
        }
    }
    let mut out: Vec<Scalar> = clusters
        .iter()
        .map(|(sum, count, _, _, _)| sum / count)
        .collect();
    // A closed conic's cuts either side of its parameter origin are one
    // point.
    if periodic && clusters.len() >= 2 {
        let (first, last) = (at(out[0])?, at(out[out.len() - 1])?);
        if merged((first - last).length()) {
            out.pop();
        }
    }
    Ok(out)
}

/// The cuts that split an imprinted edge (`curve` over `span`), as a
/// section's cuts are settled ([`merge_close`]), and with every cut within
/// the tolerance of one of the edge's own ends read as that end (#276).
///
/// The edge's ends are vertices of its own operand, where the faces that
/// bound it meet; a cut there is the other operand's boundary crossing the
/// edge where that operand's face is read as coincident with, or as
/// running along, the face the edge ends on. Kept apart, the cut would
/// leave a piece a fraction of the tolerance long that lies on the other
/// face's boundary, which no point can classify (a door's jamb on a wall's
/// floor face, with the door a rounding error longer than the wall is
/// thick). Merging it moves the other operand's boundary onto the end by
/// at most the tolerance, the reading its coincident faces already stand
/// for; recorded as [`ToleranceDecisionKind::MergedPoints`] above the
/// rounding floor, and at [`Tolerance::ZERO`] only a rounding residue
/// merges.
///
/// The wider merge can itself leave faces that do not sew (a door a
/// micrometre to half a millimetre past a wall's face, under a millimetre
/// tolerance); the boolean is then cut again with ends merged only within
/// the rounding floor, as before #276 (#291, `report::imprint_end_tolerance`,
/// [`crate::boolean_with_report`]).
///
/// A line edge that ends short of the other face's boundary, by a cut
/// within tolerance outside its span, is where the door past the face has
/// its edge cross the boundary: in the last retry (#409,
/// `report::extends_ends`) the span reaches that cut, so the piece ends on
/// the boundary and welds there, not the splitter's fixed slack away.
fn imprint_cuts(
    curve: &Curve3,
    span: Interval,
    cuts: Vec<Scalar>,
    tolerance: Tolerance,
) -> Result<(Vec<Scalar>, Interval), BooleanError> {
    let cuts = merge_close(curve, cuts, tolerance)?;
    let at = |t: Scalar| evaluate3(curve, t).map_err(|_| BooleanError::Evaluation);
    let ends = [at(span.start)?, at(span.end)?];
    let (mut lo, mut hi) = (span.start.min(span.end), span.start.max(span.end));
    let mut out = Vec::with_capacity(cuts.len());
    for cut in cuts {
        let point = at(cut)?;
        let gap = ends
            .iter()
            .map(|end| (point - *end).length())
            .fold(Scalar::INFINITY, Scalar::min);
        let ends = report::imprint_end_tolerance(tolerance);
        if !report::near(ToleranceDecisionKind::MergedPoints, gap, ends) {
            out.push(cut);
        } else if gap > report::rounding() {
            report::merged_an_end();
            // A line's end short of the other face's boundary: the edge
            // reaches the crossing, so its piece ends on the boundary and
            // not a fraction of the tolerance before it (#409).
            if matches!(curve, Curve3::Line(_)) && report::extends_ends() {
                lo = lo.min(cut);
                hi = hi.max(cut);
            }
        }
    }
    let span = if span.start <= span.end {
        Interval::new(lo, hi)
    } else {
        Interval::new(hi, lo)
    };
    Ok((out, span))
}

/// The pieces of a bounded edge span between the cuts inside it.
fn pieces_within(curve: &Curve3, span: Interval, cuts: Vec<Scalar>) -> Vec<Interval> {
    let (lo, hi) = (span.start.min(span.end), span.start.max(span.end));
    let slack = 1e-12 * (1.0 + lo.abs().max(hi.abs()));
    let periodic = matches!(curve, Curve3::Circle(_) | Curve3::Ellipse(_));
    let mut inside = vec![lo, hi];
    for cut in cuts {
        let candidates: &[Scalar] = if periodic {
            &[cut - TAU, cut, cut + TAU, cut + 2.0 * TAU]
        } else {
            &[cut]
        };
        for &c in candidates {
            if c > lo + slack && c < hi - slack {
                inside.push(c);
            }
        }
    }
    inside.sort_by(Scalar::total_cmp);
    inside.dedup_by(|x, y| (*x - *y).abs() <= 1e-12 * (1.0 + x.abs()));
    inside
        .windows(2)
        .map(|pair| Interval::new(pair[0], pair[1]))
        .collect()
}

/// Whether two parameters are one cut: [`pieces`]' deduplication slack, relative
/// to the parameter. In world-charted parameters it is not small: at
/// |t| = 5.6e6 it is 5.6e-6, wider than the 1.2e-6 by which a window's
/// crossing can pass a face (#445).
fn same_cut(x: Scalar, y: Scalar) -> bool {
    (x - y).abs() <= 1e-12 * (1.0 + x.abs().max(y.abs()))
}

/// The crossings of a window that may end a piece: those that are not
/// [`same_cut`] with a feature cut. A window end within the slack of a
/// feature cut would otherwise sort first and replace it in [`pieces`]'
/// deduplication, moving the section's end off the feature's exact vertex
/// (#445). The cut ends the piece.
fn window_ends(cuts: &[Scalar], crossings: Vec<Scalar>) -> Vec<Scalar> {
    crossings
        .into_iter()
        .filter(|&c| !cuts.iter().any(|&f| same_cut(c, f)))
        .collect()
}

/// The stretches between consecutive cuts, each with the curve it lies on:
/// finite stretches of a line (its unbounded ends leave every bounded
/// face), the cyclic arcs of a conic, and for a traced section the
/// stretches of its span. A traced loop's stretch across the loop's start
/// becomes a curve of its own.
fn pieces(curve: &Curve3, mut cuts: Vec<Scalar>) -> Result<Vec<(Curve3, Interval)>, BooleanError> {
    cuts.sort_by(Scalar::total_cmp);
    cuts.dedup_by(|x, y| same_cut(*x, *y));
    let plain = |spans: Vec<Interval>| spans.into_iter().map(|s| (curve.clone(), s)).collect();
    Ok(match curve {
        Curve3::Line(_) => plain(
            cuts.windows(2)
                .map(|pair| Interval::new(pair[0], pair[1]))
                .collect(),
        ),
        Curve3::ImplicitSection(section) => {
            let n = section.curve.end();
            let (pu, pv) = section.carrier.periodic();
            let closure = section.curve.closure(pu, pv);
            let slack = 1e-9 * (1.0 + n);
            // A cut at a loop's own start (a chain closing through a vertex
            // where branches cross, #249) is a cut like any other: the
            // loop's stretches then run between its start and the cuts.
            let at_start = cuts.iter().any(|&c| c <= slack || c >= n - slack);
            let mut inner: Vec<Scalar> = cuts
                .into_iter()
                .filter(|&c| c > slack && c < n - slack)
                .collect();
            inner.dedup_by(|x, y| (*x - *y).abs() <= slack);
            let closure = closure.filter(|_| !at_start);
            let own = |a: Scalar, b: Scalar| -> Result<(Curve3, Interval), BooleanError> {
                // One cut on a loop: the whole loop, starting there.
                let sub = if (a - b).abs() <= slack {
                    closure.and_then(|c| section.curve.rotated(a, c))
                } else {
                    section.curve.sub(a, b, closure)
                };
                let sub = sub.ok_or(BooleanError::Evaluation)?;
                let end = sub.end();
                Ok((
                    Curve3::ImplicitSection(axiolid_curve::ImplicitSection3 {
                        carrier: section.carrier.clone(),
                        curve: sub,
                    }),
                    Interval::new(0.0, end),
                ))
            };
            if closure.is_some() {
                if inner.is_empty() {
                    return Ok(vec![(curve.clone(), Interval::new(0.0, n))]);
                }
                let mut out = Vec::new();
                for pair in inner.windows(2) {
                    out.push(own(pair[0], pair[1])?);
                }
                // Across the loop's own start.
                out.push(own(inner[inner.len() - 1], inner[0])?);
                out
            } else {
                let mut ends = vec![0.0];
                ends.extend(inner);
                ends.push(n);
                let mut out = Vec::new();
                for pair in ends.windows(2) {
                    out.push(own(pair[0], pair[1])?);
                }
                out
            }
        }
        Curve3::PairSection(section) => {
            let n = section.end();
            let slack = 1e-9 * (1.0 + n);
            let mut inner: Vec<Scalar> = cuts
                .into_iter()
                .filter(|&c| c > slack && c < n - slack)
                .collect();
            inner.dedup_by(|x, y| (*x - *y).abs() <= slack);
            let (first, last) = (
                section.nodes[0].point,
                section.nodes[section.nodes.len() - 1].point,
            );
            let closed = first == last;
            let own = |a: Scalar, b: Scalar| -> Result<(Curve3, Interval), BooleanError> {
                let sub = if a < b {
                    section.sub(a, b)
                } else {
                    // Across a loop's own start.
                    section
                        .sub(a, n)
                        .zip(section.sub(0.0, b))
                        .map(|(mut x, y)| {
                            x.nodes.extend(y.nodes.into_iter().skip(1));
                            x
                        })
                };
                let sub = sub.ok_or(BooleanError::Evaluation)?;
                let end = sub.end();
                Ok((Curve3::PairSection(sub), Interval::new(0.0, end)))
            };
            let mut out = Vec::new();
            if closed {
                if inner.is_empty() {
                    return Ok(vec![(curve.clone(), Interval::new(0.0, n))]);
                }
                for pair in inner.windows(2) {
                    out.push(own(pair[0], pair[1])?);
                }
                out.push(own(inner[inner.len() - 1], inner[0])?);
            } else {
                let mut ends = vec![0.0];
                ends.extend(inner);
                ends.push(n);
                for pair in ends.windows(2) {
                    out.push(own(pair[0], pair[1])?);
                }
            }
            out
        }
        _ => {
            if cuts.is_empty() {
                return Ok(vec![(curve.clone(), Interval::new(0.0, TAU))]);
            }
            let mut out: Vec<Interval> = cuts
                .windows(2)
                .map(|pair| Interval::new(pair[0], pair[1]))
                .collect();
            out.push(Interval::new(cuts[cuts.len() - 1], cuts[0] + TAU));
            plain(out)
        }
    })
}

/// One operand, with each face's certified domain built once.
struct Side<'a> {
    brep: &'a ExactBRep,
    faces: Vec<FaceId>,
    domains: Vec<FaceDomain<'a>>,
    /// A sound box around each face, enlarged by the tolerance.
    boxes: Vec<Option<Aabb>>,
    /// A sound box around each edge, enlarged by the tolerance.
    edge_boxes: Vec<Option<Aabb>>,
    /// For each edge, the faces whose loops use it.
    edge_faces: Vec<Vec<usize>>,
}

impl<'a> Side<'a> {
    fn new(brep: &'a ExactBRep, tolerance: Tolerance) -> Result<Self, BooleanError> {
        let topology = brep.topology();
        let mut faces = Vec::with_capacity(topology.faces().len());
        let mut domains = Vec::with_capacity(topology.faces().len());
        for index in 0..topology.faces().len() {
            let face = topology
                .face_id_at(index)
                .ok_or(BooleanError::DanglingReference)?;
            let domain = FaceDomain::new(brep, face, tolerance)
                .map_err(BooleanError::Measure)?
                .ok_or(BooleanError::UnsupportedTrim)?;
            faces.push(face);
            domains.push(domain);
        }
        let mut edge_faces = vec![Vec::new(); topology.edges().len()];
        for (index, face) in topology.faces().iter().enumerate() {
            for bound in &face.bounds {
                let wire = topology
                    .loops()
                    .get(bound.loop_id.index())
                    .ok_or(BooleanError::DanglingReference)?;
                for use_ in &wire.edges {
                    edge_faces[use_.edge.index()].push(index);
                }
            }
        }
        let mut boxes = Vec::with_capacity(domains.len());
        for (index, domain) in domains.iter().enumerate() {
            let surface = topology.faces()[index]
                .surface
                .and_then(|id| brep.surfaces().get(id.index()))
                .ok_or(BooleanError::DanglingReference)?;
            let (lo, hi) = domain.bounds();
            boxes.push(face_box(surface, lo, hi, tolerance));
        }
        let mut edge_boxes = Vec::with_capacity(topology.edges().len());
        for (index, edge) in topology.edges().iter().enumerate() {
            let curve = edge.curve.and_then(|id| brep.curves3().get(id.index()));
            let span = topology
                .edge_id_at(index)
                .and_then(|id| brep.edge_interval(id));
            edge_boxes.push(match (curve, span) {
                (Some(curve), Some(span)) => edge_box(curve, span, tolerance),
                _ => None,
            });
        }
        Ok(Self {
            brep,
            faces,
            domains,
            boxes,
            edge_boxes,
            edge_faces,
        })
    }

    fn surface(&self, face: usize) -> Result<&'a Surface, BooleanError> {
        self.brep.topology().faces()[face]
            .surface
            .and_then(|id| self.brep.surfaces().get(id.index()))
            .ok_or(BooleanError::DanglingReference)
    }

    /// Where `point`, on the face's support, lies: on one of the edges the
    /// curve through it runs along (`along`), or inside or outside the face.
    fn locate(
        &self,
        face: usize,
        point: Point3,
        along: &[(Curve3, Interval)],
        tolerance: Tolerance,
    ) -> Result<Place, BooleanError> {
        for (curve, span) in along {
            if on_edge(curve, *span, point, tolerance)? {
                return Ok(Place::Boundary);
            }
        }
        let (u, v) = locate(self.surface(face)?, point, report::floored(tolerance))
            .map_err(|_| BooleanError::Evaluation)?;
        match self.domains[face]
            .contains(Point2::new(u, v))
            .map_err(BooleanError::Measure)?
        {
            Some(true) => Ok(Place::Inside),
            Some(false) => Ok(Place::Outside),
            None => Err(BooleanError::Undecided),
        }
    }

    /// The face's boundary edges, each once, with their curves and spans.
    fn edges_of(&self, face: usize) -> Result<Vec<(usize, Curve3, Interval)>, BooleanError> {
        let topology = self.brep.topology();
        let mut out = Vec::new();
        let mut seen = Vec::new();
        for bound in &topology.faces()[face].bounds {
            let wire = topology
                .loops()
                .get(bound.loop_id.index())
                .ok_or(BooleanError::DanglingReference)?;
            for use_ in &wire.edges {
                if seen.contains(&use_.edge) {
                    continue;
                }
                seen.push(use_.edge);
                let curve = topology.edges()[use_.edge.index()]
                    .curve
                    .and_then(|id| self.brep.curves3().get(id.index()))
                    .ok_or(BooleanError::DanglingReference)?;
                let span = self
                    .brep
                    .edge_interval(use_.edge)
                    .ok_or(BooleanError::DanglingReference)?;
                out.push((use_.edge.index(), curve.clone(), span));
            }
        }
        Ok(out)
    }

    /// Parameters on `curve` where it crosses a boundary edge of the face,
    /// and the edges it runs along (lying in the adjacent face's surface
    /// too), which cut it at their ends.
    #[allow(clippy::type_complexity)]
    ///
    /// `meets` is the other surface `curve` lies on. Where the curve cannot
    /// be intersected with the adjacent face's surface (a B-spline), the
    /// edge is intersected with `meets` instead: the curve crosses the edge
    /// exactly where the edge crosses `meets`.
    ///
    /// With a `window` (a box around everything the section can matter in),
    /// an edge whose own box misses it is skipped: any crossing with it lies
    /// outside the window, where the caller has cut the curve at the
    /// window's boundary and drops the pieces (#228).
    ///
    /// Where the edge's adjacent surface and one of the two the curve lies
    /// on are a plane and a cylinder read as touching (`contacts`), the
    /// curve is cut where it meets their contact ruling, as for every other
    /// face pair of the two (`crate::contact`, #243).
    fn cuts(
        &self,
        face: usize,
        curve: &Curve3,
        meets: &Surface,
        window: Option<&Aabb>,
        contacts: &Contacts,
        tolerance: Tolerance,
    ) -> Result<(Vec<Scalar>, Vec<(Curve3, Interval)>), BooleanError> {
        let topology = self.brep.topology();
        let mut out = Vec::new();
        let mut along = Vec::new();
        let mut seen = Vec::new();
        // A pole or apex the curve passes through cuts it: its pcurve jumps
        // round the angle there.
        for pole in poles(self.surface(face)?) {
            if let Ok(t) = locate3(curve, pole, report::floored(tolerance)) {
                let on = evaluate3(curve, t).map_err(|_| BooleanError::Evaluation)?;
                if report::near(
                    ToleranceDecisionKind::IncidentPoint,
                    (on - pole).length(),
                    tolerance,
                ) {
                    out.push(t);
                }
            }
        }
        for bound in &topology.faces()[face].bounds {
            let wire = topology
                .loops()
                .get(bound.loop_id.index())
                .ok_or(BooleanError::DanglingReference)?;
            for use_ in &wire.edges {
                if seen.contains(&use_.edge) {
                    continue;
                }
                seen.push(use_.edge);
                if let (Some(w), Some(b)) = (window, &self.edge_boxes[use_.edge.index()]) {
                    if !w.overlaps(b) {
                        continue;
                    }
                }
                let edge = &topology.edges()[use_.edge.index()];
                let edge_curve = edge
                    .curve
                    .and_then(|id| self.brep.curves3().get(id.index()))
                    .ok_or(BooleanError::DanglingReference)?;
                let span = self
                    .brep
                    .edge_interval(use_.edge)
                    .ok_or(BooleanError::DanglingReference)?;
                let mut cutter =
                    self.cutter(face, use_.edge.index(), edge_curve, span, tolerance, false)?;
                let mut result = exact_curve_surface_intersection(curve, &cutter);
                if matches!(result, Ok(ExactCurveIntersection::Contained)) {
                    // The adjacent face may lie on the same surface under a
                    // different frame (a column's half-walls): cut across
                    // the seam instead.
                    if let Ok(seam) =
                        self.cutter(face, use_.edge.index(), edge_curve, span, tolerance, true)
                    {
                        cutter = seam;
                        result = exact_curve_surface_intersection(curve, &cutter);
                    }
                }
                // Two operands placed independently put one plane at two
                // roundings: the curve then misses the adjacent face's
                // support by 1e-17 and the exact test calls it apart or
                // crossing, though it runs along the edge. Within tolerance
                // along the whole edge, it is contained (#228).
                if !matches!(result, Ok(ExactCurveIntersection::Contained))
                    && runs_along(curve, edge_curve, span, tolerance)?
                {
                    result = Ok(ExactCurveIntersection::Contained);
                }
                if !matches!(result, Ok(ExactCurveIntersection::Contained)) {
                    let own = self.surface(face)?;
                    // Along a smooth edge (a fillet running tangent into
                    // the adjacent face) the curve touches the cutter: it
                    // is cut where the edge crosses `meets` (#249).
                    if let Some(hits) = crate::tangency::smooth_edge_cuts(
                        own, &cutter, meets, edge_curve, span, curve, contacts, tolerance,
                    )? {
                        out.extend(hits);
                        continue;
                    }
                    if let Some(hits) =
                        contacts.crossing([own, meets], &cutter, curve, tolerance)?
                    {
                        for t in hits {
                            let point =
                                evaluate3(curve, t).map_err(|_| BooleanError::Evaluation)?;
                            if on_edge(edge_curve, span, point, tolerance)? {
                                out.push(t);
                            }
                        }
                        continue;
                    }
                }
                match result {
                    Ok(ExactCurveIntersection::Points(hits)) => {
                        let before = out.len();
                        for hit in hits {
                            if on_edge(edge_curve, span, hit.point, tolerance)? {
                                out.push(hit.parameter.approx());
                            }
                        }
                        // Where the curve touches the adjacent surface at
                        // the edge (a tangent crossing), rounding can leave
                        // the double root with no real root at all: the
                        // curve still meets the edge itself, transversally,
                        // and that crossing is read from the two curves.
                        if out.len() == before {
                            if let Some(t) = near_crossing(curve, edge_curve, tolerance)? {
                                let point =
                                    evaluate3(curve, t).map_err(|_| BooleanError::Evaluation)?;
                                if on_edge(edge_curve, span, point, tolerance)? {
                                    out.push(t);
                                }
                            }
                        }
                    }
                    // The curve lies in the adjacent face's surface as well
                    // as this one's: where it meets the edge it runs along
                    // it, so it is cut at the edge's ends.
                    Ok(ExactCurveIntersection::Contained) => {
                        for t in [span.start, span.end] {
                            let end =
                                evaluate3(edge_curve, t).map_err(|_| BooleanError::Evaluation)?;
                            if let Ok(s) = locate3(curve, end, report::floored(tolerance)) {
                                let on =
                                    evaluate3(curve, s).map_err(|_| BooleanError::Evaluation)?;
                                if report::near(
                                    ToleranceDecisionKind::IncidentPoint,
                                    (on - end).length(),
                                    tolerance,
                                ) {
                                    out.push(s);
                                }
                            }
                        }
                        along.push((edge_curve.clone(), span));
                    }
                    Ok(_) => return Err(BooleanError::UnsupportedSection),
                    Err(_) => {
                        let hits = match exact_curve_surface_intersection(edge_curve, meets) {
                            Ok(ExactCurveIntersection::Points(hits)) => hits,
                            _ => return Err(BooleanError::UnsupportedTrim),
                        };
                        for hit in hits {
                            if !on_span(edge_curve, span, hit.point, tolerance)? {
                                continue;
                            }
                            if let Ok(s) = locate3(curve, hit.point, report::floored(tolerance)) {
                                let on =
                                    evaluate3(curve, s).map_err(|_| BooleanError::Evaluation)?;
                                if report::near(
                                    ToleranceDecisionKind::IncidentPoint,
                                    (on - hit.point).length(),
                                    tolerance,
                                ) {
                                    out.push(s);
                                }
                            }
                        }
                    }
                }
            }
        }
        Ok((out, along))
    }
}

impl Side<'_> {
    /// A surface the section crosses exactly where it crosses the edge: the
    /// adjacent face's support, or across a seam (the same support on both
    /// sides) the plane through the ruling and the surface normal.
    fn cutter(
        &self,
        face: usize,
        edge: usize,
        edge_curve: &Curve3,
        span: Interval,
        tolerance: Tolerance,
        across_seam: bool,
    ) -> Result<Surface, BooleanError> {
        let own = self.surface(face)?;
        let adjacent = self.edge_faces[edge]
            .iter()
            .copied()
            .find(|&other| other != face);
        if let (Some(other), false) = (adjacent, across_seam) {
            let theirs = self.surface(other)?;
            if !same_support(own, theirs, tolerance) {
                return Ok(theirs.clone());
            }
        }
        // A seam: the same support on both sides (or a free edge). A seam
        // circle (a sphere's meridian or latitude, a torus's tube or ring
        // circle) is cut by the surface its normals sweep: coaxial with the
        // circle, a plane, cylinder or cone, which crosses the face
        // transversally along the whole circle.
        if let Curve3::Circle(circle) = edge_curve {
            return normal_sweep(own, circle, tolerance);
        }
        let Curve3::Line(line) = edge_curve else {
            return Err(BooleanError::UnsupportedTrim);
        };
        let mid = evaluate3(edge_curve, 0.5 * (span.start + span.end))
            .map_err(|_| BooleanError::Evaluation)?;
        let (u, v) =
            locate(own, mid, report::floored(tolerance)).map_err(|_| BooleanError::Evaluation)?;
        let n = normal(own, u, v).map_err(|_| BooleanError::Evaluation)?;
        let across = line.direction.cross(n);
        let length = across.length();
        if length == 0.0 || !length.is_finite() {
            return Err(BooleanError::UnsupportedTrim);
        }
        let z = across / length;
        let x = line.direction.normalize();
        Ok(Surface::Plane(Plane {
            frame: axiolid_core::Frame3 {
                origin: mid,
                x,
                y: z.cross(x),
                z,
            },
        }))
    }
}

/// The surface swept by `surface`'s normal lines along a circle on it,
/// when the surface is one of revolution about the circle's axis (or the
/// circle is a meridian, whose normals stay in its plane): a plane, a
/// cylinder or a cone coaxial with the circle.
fn normal_sweep(
    surface: &Surface,
    circle: &axiolid_curve::Circle3,
    tolerance: Tolerance,
) -> Result<Surface, BooleanError> {
    let f = circle.frame;
    let axis = f.x.cross(f.y).normalize();
    let x = f.x.normalize();
    let y = axis.cross(x);
    let p = f.origin + x * circle.radius;
    let (u, v) =
        locate(surface, p, report::floored(tolerance)).map_err(|_| BooleanError::Evaluation)?;
    let n = normal(surface, u, v)
        .map_err(|_| BooleanError::Evaluation)?
        .normalize();
    let frame = axiolid_core::Frame3 {
        origin: f.origin,
        x,
        y,
        z: axis,
    };
    let along = n.dot(axis);
    let radial = n.dot(x);
    let plane = || {
        Surface::Plane(Plane {
            frame: axiolid_core::Frame3 {
                origin: f.origin,
                x,
                y,
                z: axis,
            },
        })
    };
    // Normals within the circle's plane: the plane itself.
    if along.abs() <= 1e-12 {
        return Ok(plane());
    }
    // Normals along the axis: the cylinder through the circle.
    if radial.abs() <= 1e-12 {
        return Ok(Surface::Cylinder(axiolid_surface::Cylinder {
            frame,
            radius: circle.radius,
        }));
    }
    // Otherwise the cone of normal lines: at height `h` along the axis the
    // line is `radius + h * radial / along` from it.
    Ok(Surface::Cone(axiolid_surface::Cone {
        frame,
        radius: circle.radius,
        semi_angle: (radial / along).atan(),
    }))
}

/// Whether every point of the edge (`edge` over `span`) lies within the
/// caller's linear tolerance of `curve`, proved in closed form:
///
/// - two lines: the distance to a line is convex along a segment, so the
///   segment's two ends bound it;
/// - two circles: a point of one is within `|c1 - c2| + |r1 - r2| +
///   r2 * theta` of the other, `theta` the angle between their planes
///   (bounded by `2 |n1 x n2|`).
///
/// Any other pair is never read as running along (#228).
fn runs_along(
    curve: &Curve3,
    edge: &Curve3,
    span: Interval,
    tolerance: Tolerance,
) -> Result<bool, BooleanError> {
    let along =
        |residue: Scalar| report::near(ToleranceDecisionKind::SectionAlongEdge, residue, tolerance);
    match (curve, edge) {
        (Curve3::Line(line), Curve3::Line(_)) => {
            let d = line.direction.normalize_or_zero();
            if d == Vec3::ZERO {
                return Ok(false);
            }
            let mut residue: Scalar = 0.0;
            for t in [span.start, span.end] {
                let p = evaluate3(edge, t).map_err(|_| BooleanError::Evaluation)?;
                let off = p - line.origin;
                residue = residue.max((off - d * off.dot(d)).length());
            }
            Ok(along(residue))
        }
        (Curve3::Circle(a), Curve3::Circle(b)) => {
            let (na, nb) = (a.frame.z.normalize_or_zero(), b.frame.z.normalize_or_zero());
            if na == Vec3::ZERO || nb == Vec3::ZERO {
                return Ok(false);
            }
            let tilt = 2.0 * na.cross(nb).length();
            Ok(along(
                (a.frame.origin - b.frame.origin).length()
                    + (a.radius - b.radius).abs()
                    + a.radius * tilt,
            ))
        }
        _ => Ok(false),
    }
}

/// A traced piece (`curve` over `span`, its `sample` inside) that runs along
/// one of `edges` that is a line, circle or ellipse, as that edge's own
/// curve over the stretch between the piece's ends; any other piece
/// unchanged (#382).
///
/// A plane through a torus's axis meets the torus in the circles its faces
/// start from, found as traced sections; the coincident cap on that plane
/// imprints the same circles in closed form. A parametric midpoint decides
/// whether two copies of a stretch are one (`crate::split`), and a traced
/// curve's parameter is not the circle's angle, so the copies would be
/// imprinted twice and the faces would not sew. One stretch has one curve.
///
/// The ends move onto the edge by their residue, adopted only within the
/// tolerance (a decision [`ToleranceDecisionKind::SectionAlongEdge`], read
/// as for a section along an edge); a piece cut at its edge's ends, as
/// every piece along an edge is, does not move further than the incident
/// tolerance.
fn along_edge_curve(
    curve: Curve3,
    span: Interval,
    sample: Point3,
    edges: &[&(Curve3, Interval)],
    tolerance: Tolerance,
) -> Result<(Curve3, Interval), BooleanError> {
    if !matches!(curve, Curve3::ImplicitSection(_) | Curve3::PairSection(_)) {
        return Ok((curve, span));
    }
    let at = |c: &Curve3, t: Scalar| evaluate3(c, t).map_err(|_| BooleanError::Evaluation);
    let (start, end) = (at(&curve, span.start)?, at(&curve, span.end)?);
    for (edge, edge_span) in edges.iter().copied() {
        let periodic = matches!(edge, Curve3::Circle(_) | Curve3::Ellipse(_));
        if !(periodic || matches!(edge, Curve3::Line(_)))
            || !on_edge(edge, *edge_span, sample, tolerance)?
        {
            continue;
        }
        let locate = |p: Point3| {
            locate3(edge, p, report::floored(tolerance)).map_err(|_| BooleanError::Evaluation)
        };
        let (ts, te, tm) = (locate(start)?, locate(end)?, locate(sample)?);
        let adopted = if periodic {
            // Round the stretch from `ts` through `tm` to `te`, whichever
            // way it runs: `te` and `tm` in `(ts, ts + 2 pi]`.
            let up = |t: Scalar| t - TAU * ((t - ts) / TAU).ceil() + TAU;
            let (te_up, tm_up) = (up(te), up(tm));
            if tm_up < te_up {
                Interval::new(ts, te_up)
            } else {
                Interval::new(te_up - TAU, ts)
            }
        } else {
            Interval::new(ts.min(te), ts.max(te))
        };
        let residue = (at(edge, ts)? - start)
            .length()
            .max((at(edge, te)? - end).length());
        if !report::near(ToleranceDecisionKind::SectionAlongEdge, residue, tolerance) {
            return Ok((curve, span));
        }
        return Ok((edge.clone(), adopted));
    }
    Ok((curve, span))
}

/// Whether `point` lies on the edge within tolerance and inside its span.
fn on_edge(
    curve: &Curve3,
    span: Interval,
    point: Point3,
    tolerance: Tolerance,
) -> Result<bool, BooleanError> {
    let Ok(t) = locate3(curve, point, report::floored(tolerance)) else {
        return Ok(false);
    };
    let on = evaluate3(curve, t).map_err(|_| BooleanError::Evaluation)?;
    if !on_span(curve, span, point, tolerance)? {
        return Ok(false);
    }
    Ok(report::near(
        ToleranceDecisionKind::IncidentPoint,
        (on - point).length(),
        tolerance,
    ))
}

/// Whether `point`, on `curve`, lies within the edge's span (a closed
/// conic's span may start anywhere and run past a full turn).
pub(crate) fn on_span(
    curve: &Curve3,
    span: Interval,
    point: Point3,
    tolerance: Tolerance,
) -> Result<bool, BooleanError> {
    let t =
        locate3(curve, point, report::floored(tolerance)).map_err(|_| BooleanError::Evaluation)?;
    let (lo, hi) = (span.start.min(span.end), span.start.max(span.end));
    let slack = 1e-9 * (1.0 + lo.abs().max(hi.abs()));
    let periodic = matches!(curve, Curve3::Circle(_) | Curve3::Ellipse(_));
    let candidates: &[Scalar] = if periodic {
        &[t - TAU, t, t + TAU, t + 2.0 * TAU]
    } else {
        &[t]
    };
    Ok(candidates
        .iter()
        .any(|c| *c >= lo - slack && *c <= hi + slack))
}

#[cfg(test)]
mod tests {
    use super::*;
    use axiolid_core::Frame3;
    use axiolid_curve::{Circle3, Line3};

    const EPS: Scalar = 1e-6;

    #[test]
    fn a_window_end_never_replaces_a_cut() {
        // #445: the corbel's section line at world y = 5.6e6. The window's
        // crossings pass the feature cuts by 1.2e-6, inside the dedup slack.
        let line = Curve3::Line(Line3 {
            origin: Point3::new(600000.2, 0.0, 203.5),
            direction: -Vec3::Y,
        });
        let cuts = vec![-5600000.1, -5599999.9];
        let ends = window_ends(&cuts, vec![-5599999.8999988, -5600000.1000012]);
        assert!(ends.is_empty());
        let spans: Vec<Interval> = pieces(&line, cuts)
            .expect("pieces")
            .into_iter()
            .map(|p| p.1)
            .collect();
        assert_eq!(spans, vec![Interval::new(-5600000.1, -5599999.9)]);
        // Near the origin the slack is 1e-12 and both ends stay.
        let kept = window_ends(
            &[-0.1, 0.1],
            vec![0.10000119073677866, -0.10000119073677866],
        );
        assert_eq!(kept.len(), 2);
    }

    fn tol() -> Tolerance {
        Tolerance::METRE
    }

    fn x_axis() -> Curve3 {
        Curve3::Line(Line3 {
            origin: Point3::ZERO,
            direction: Vec3::X,
        })
    }

    #[test]
    fn cuts_within_tolerance_merge_and_a_longer_chain_is_refused() {
        let line = x_axis();
        let merged = merge_close(&line, vec![1.0, 1.0 + 0.5 * EPS, 2.0], tol()).unwrap();
        assert_eq!(merged.len(), 2);
        assert!((merged[0] - (1.0 + 0.25 * EPS)).abs() < 1e-15);
        // Ten tolerances apart: two cuts, decided exactly.
        let kept = merge_close(&line, vec![1.0, 1.0 + 10.0 * EPS], tol()).unwrap();
        assert_eq!(kept.len(), 2);
        // Each within tolerance of the next, 1.2 tolerances end to end.
        let chain = merge_close(&line, vec![1.0, 1.0 + 0.6 * EPS, 1.0 + 1.2 * EPS], tol());
        assert_eq!(chain, Err(BooleanError::NearCoincidence));
    }

    #[test]
    fn cuts_at_a_circles_seam_from_either_side_merge_at_the_seam() {
        // One point named at `0` and twice at `2 pi` (a seam cut, and the
        // ends of edges along the seam): one cut there, not their mean
        // `4 pi / 3` on the other side of the circle (#243).
        let circle = Curve3::Circle(Circle3 {
            frame: Frame3 {
                origin: Point3::ZERO,
                x: Vec3::X,
                y: Vec3::Y,
                z: Vec3::Z,
            },
            radius: 1.0,
        });
        let merged = merge_close(&circle, vec![0.0, TAU, TAU], tol()).unwrap();
        assert_eq!(merged.len(), 1);
        let at = evaluate3(&circle, merged[0]).unwrap();
        assert!((at - Point3::X).length() < 1e-12, "{at:?}");
        // Close either side of the seam: their mean is at the seam.
        let merged = merge_close(&circle, vec![TAU - 0.2 * EPS, 0.2 * EPS], tol()).unwrap();
        assert_eq!(merged.len(), 1);
        let at = evaluate3(&circle, merged[0]).unwrap();
        assert!((at - Point3::X).length() < 0.25 * EPS, "{at:?}");
    }

    #[test]
    fn an_edge_runs_along_a_curve_only_when_proved_within_tolerance() {
        let line = x_axis();
        let edge = |dy: Scalar| {
            Curve3::Line(Line3 {
                origin: Point3::new(0.0, dy, 0.0),
                direction: Vec3::new(1.0, dy, 0.0),
            })
        };
        let span = Interval::new(0.0, 1.0);
        // From 0.5 eps to 1.0 eps off the axis: within.
        assert!(runs_along(&line, &edge(0.5 * EPS), span, tol()).unwrap());
        // Its far end ten tolerances off: not.
        assert!(!runs_along(&line, &edge(5.0 * EPS), span, tol()).unwrap());
        let circle = |r: Scalar, tilt: Scalar| {
            Curve3::Circle(Circle3 {
                frame: Frame3 {
                    origin: Point3::ZERO,
                    x: Vec3::X,
                    y: Vec3::new(0.0, tilt.cos(), tilt.sin()),
                    z: Vec3::new(0.0, -tilt.sin(), tilt.cos()),
                },
                radius: r,
            })
        };
        let whole = Interval::new(0.0, core::f64::consts::TAU);
        assert!(runs_along(
            &circle(1.0, 0.0),
            &circle(1.0 + 0.5 * EPS, 0.0),
            whole,
            tol()
        )
        .unwrap());
        assert!(!runs_along(&circle(1.0, 0.0), &circle(1.0, 10.0 * EPS), whole, tol()).unwrap());
        // Other families are never read as running along.
        assert!(!runs_along(&line, &circle(1.0, 0.0), whole, tol()).unwrap());
    }

    #[test]
    fn a_plane_reads_parallel_to_an_axis_only_within_tolerance_over_its_extent() {
        let cylinder = Surface::Cylinder(axiolid_surface::Cylinder {
            frame: Frame3 {
                origin: Point3::ZERO,
                x: Vec3::X,
                y: Vec3::Y,
                z: Vec3::Z,
            },
            radius: 1.0,
        });
        let plane = |tilt: Scalar, offset: Scalar| {
            let z = Vec3::new(1.0, 0.0, tilt).normalize();
            Surface::Plane(Plane {
                frame: Frame3 {
                    origin: Point3::new(offset, 0.0, 0.0),
                    x: Vec3::Y,
                    y: z.cross(Vec3::Y),
                    z,
                },
            })
        };
        // Tilted 1e-9 over a 10 m extent: 1e-8 off, read parallel; the
        // plane through the axis cuts two rulings.
        let rulings =
            plane_cylinder_within_rounding(&plane(1e-9, 0.0), &cylinder, Some(10.0), tol());
        assert_eq!(
            rulings.map(|r| (r.curves.len(), r.contact)),
            Some((2, false))
        );
        // Tilted 1e-6 over 10 m: 1e-5 off, left to the exact closed form.
        assert!(
            plane_cylinder_within_rounding(&plane(1e-6, 0.0), &cylinder, Some(10.0), tol())
                .is_none()
        );
        // No extent: never read parallel.
        assert!(
            plane_cylinder_within_rounding(&plane(1e-9, 0.0), &cylinder, None, tol()).is_none()
        );
        // Within tolerance of touching: one ruling.
        let touching = plane_cylinder_within_rounding(
            &plane(0.0, 1.0 + 0.5 * EPS),
            &cylinder,
            Some(10.0),
            tol(),
        );
        assert_eq!(
            touching.map(|r| (r.curves.len(), r.contact)),
            Some((1, true))
        );
        // Exactly parallel and ten tolerances clear, or exactly parallel
        // through the axis: no reading, the exact closed form decides
        // (#236).
        for offset in [1.0 + 10.0 * EPS, 0.0] {
            assert!(plane_cylinder_within_rounding(
                &plane(0.0, offset),
                &cylinder,
                Some(10.0),
                tol()
            )
            .is_none());
        }
    }

    #[test]
    fn only_readings_the_exact_predicates_do_not_confirm_are_reported() {
        let cylinder = Surface::Cylinder(axiolid_surface::Cylinder {
            frame: Frame3 {
                origin: Point3::ZERO,
                x: Vec3::X,
                y: Vec3::Y,
                z: Vec3::Z,
            },
            radius: 1.0,
        });
        let plane = |z: Vec3, x: Vec3| {
            Surface::Plane(Plane {
                frame: Frame3 {
                    origin: Point3::new(0.0, 0.0, 0.5),
                    x,
                    y: z.cross(x),
                    z,
                },
            })
        };
        // Exactly perpendicular: the closed form's circle, nothing read.
        let session = report::open(&[], tol());
        assert!(plane_cylinder_within_rounding(
            &plane(Vec3::Z, Vec3::X),
            &cylinder,
            Some(10.0),
            tol()
        )
        .is_none());
        assert!(session.finish().expect("exact").is_exact());
        // Tilted by a rotation's rounding residue: read perpendicular, and
        // reported.
        let session = report::open(&[], tol());
        let tilted = Vec3::new(0.0, (core::f64::consts::FRAC_PI_2).cos(), 1.0);
        let circle =
            plane_cylinder_within_rounding(&plane(tilted, Vec3::X), &cylinder, Some(10.0), tol());
        assert_eq!(
            circle.map(|c| (c.curves.len(), c.contact)),
            Some((1, false))
        );
        let report = session.finish().expect("within tolerance");
        assert!(report.contains(ToleranceDecisionKind::PlanePerpendicularToAxis));
        assert!(report.linear() > 0.0 && report.linear() <= EPS);
    }
}