ogeom_mesh/

triangulate.rs

//! Triangulating a face.
//!
//! The second half of tessellation. A face is a trimmed region of a surface, so
//! the triangulation is built in the surface's `(u, v)` parameter space (where
//! the region is an ordinary polygon with holes) and then lifted back into
//! space by evaluating the surface at each vertex.
//!
//! # Why parameter space
//!
//! Triangulating in 3D would mean deciding which side of a curved boundary a
//! point falls on, in space, which is the point-in-solid problem. In parameter
//! space the boundary is a closed 2D polygon and the question is a winding
//! count. The surface does the rest.
//!
//! The cost is that parameter space is distorted: equal steps in `(u, v)` cover
//! very different distances near a sphere's pole than near its equator. So the
//! interior points are chosen by measuring deflection *in space* and the
//! triangulation is done in parameter space: measuring where the answer
//! matters, connecting where it is easy.
//!
//! # Watertightness
//!
//! A face's boundary points come from discretizing the *edge's* 3D curve and
//! evaluating the pcurve at those same parameters. Two faces sharing an edge
//! therefore place their boundary vertices at identical spatial positions, and
//! the join has no gap. Discretizing each face's pcurve independently would
//! give each face its own idea of where the edge runs, and the seams would show.

use ogeom_core::clock::Instant;
use ogeom_core::{Exact, OgeomResult, Predicates, Tolerances, ogeom_bail};
use ogeom_geom::Curve3d as _;
use ogeom_geom::{Curve2d, Surface, SurfaceGeometry};
use ogeom_math::{Direction, Point, Point2, Vector};
use ogeom_topo::{EdgeRepr, Model, NodeData, Orientation, Shape, ShapeType, Triangulation};
use spade::handles::FixedVertexHandle;
use spade::{
    ConstrainedDelaunayTriangulation, Point2 as SpadePoint, Triangulation as _, mitigate_underflow,
};

use crate::discretize::{Deflection, discretize};

/// Triangulate one face.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if `face` is not a
/// face or its geometry is missing;
/// [`OgeomError::Dangling`](ogeom_core::OgeomError::Dangling) if a handle fails to
/// resolve; [`OgeomError::NotDone`](ogeom_core::OgeomError::NotDone) if the boundary
/// cannot be triangulated.
pub fn triangulate_face(
    model: &Model,
    face: &Shape,
    deflection: Deflection,
    tol: Tolerances,
) -> OgeomResult<Triangulation> {
    // One pass at the caller's chord, and a second only where the first
    // came back short, the same shape as the whole-shape path, so a caller
    // meshing face by face pays for one triangulation per face, not two.
    let nothing = EdgeChords::new();
    let (mesh, verdict) = triangulate_reporting(model, face, deflection, Some(&nothing), tol)?;
    if verdict != Verdict::Short {
        return Ok(mesh);
    }
    triangulate_with(model, face, deflection, None, None, tol)
}

/// As [`triangulate_face`], taking the face from `kept` where it drew this
/// face at this deflection before, and keeping what it draws there.
///
/// # Errors
///
/// As [`triangulate_face`].
pub fn triangulate_face_kept(
    model: &Model,
    face: &Shape,
    deflection: Deflection,
    kept: &FaceMeshCache,
    tol: Tolerances,
) -> OgeomResult<Triangulation> {
    let nothing = EdgeChords::new();
    let (mesh, verdict) = reporting_kept(
        model,
        face,
        deflection,
        Some(&nothing),
        None,
        Some(kept),
        Pass::Agreeing,
        tol,
    )?;
    if verdict != Verdict::Short {
        return Ok(mesh);
    }
    triangulate_with(model, face, deflection, None, Some(kept), tol)
}

/// Face meshes kept by what each face is, for a caller that meshes a run of
/// models in which most faces recur unchanged, as a construction retried
/// after a few of its faces changed does.
///
/// A face is known by everything its mesh reads: its content as
/// [`ogeom_topo::face_content`] writes it out, then the deflection, the
/// tolerances, and the chord each of its edges was asked to be drawn to.
/// A face matches a kept one only where every one of these is the same to
/// the bit, and the kept mesh is the one drawing it again gives. A face
/// whose content cannot be read is drawn and not kept.
#[derive(Debug, Default)]
pub struct FaceMeshCache {
    held: std::sync::Mutex<ogeom_core::FastMap<Vec<u8>, (Triangulation, Verdict)>>,
    /// The chords a face that came up short asked its edges to be drawn
    /// at, under the same key as its mesh: each edge's in the order the
    /// face's wires walk them, `None` for one it did not name.
    chords: std::sync::Mutex<ogeom_core::FastMap<Vec<u8>, Vec<Option<f64>>>>,
}

impl FaceMeshCache {
    /// An empty cache.
    #[must_use]
    pub fn new() -> Self {
        Self::default()
    }

    fn get(&self, key: &[u8]) -> Option<(Triangulation, Verdict)> {
        self.held
            .lock()
            .unwrap_or_else(std::sync::PoisonError::into_inner)
            .get(key)
            .cloned()
    }

    fn put(&self, key: Vec<u8>, found: (Triangulation, Verdict)) {
        self.held
            .lock()
            .unwrap_or_else(std::sync::PoisonError::into_inner)
            .insert(key, found);
    }

    fn chords(&self, key: &[u8]) -> Option<Vec<Option<f64>>> {
        self.chords
            .lock()
            .unwrap_or_else(std::sync::PoisonError::into_inner)
            .get(key)
            .cloned()
    }

    fn keep_chords(&self, key: Vec<u8>, found: Vec<Option<f64>>) {
        self.chords
            .lock()
            .unwrap_or_else(std::sync::PoisonError::into_inner)
            .insert(key, found);
    }
}

/// The key a face's mesh is kept under: see [`FaceMeshCache`]. `finer` is
/// the shared chord map the face is drawn to, or `None` where it works out
/// its own from its content alone; `pass` is whether the face is drawn
/// again if it comes up short.
fn mesh_key(
    model: &Model,
    face: &Shape,
    deflection: Deflection,
    finer: Option<&EdgeChords>,
    pass: Pass,
    tol: Tolerances,
) -> Option<Vec<u8>> {
    let mut key = ogeom_topo::face_content(model, face)?;
    key.extend_from_slice(format!("{deflection:?}{tol:?}{pass:?}").as_bytes());
    match finer {
        None => key.push(b'o'),
        Some(finer) => {
            key.push(b's');
            for wire in model.ordered_children_of(face).ok()? {
                for edge in model.ordered_children_of(&wire).ok()? {
                    match finer.get(&edge.node().index()) {
                        Some(chord) => {
                            key.push(1);
                            key.extend_from_slice(&chord.to_bits().to_le_bytes());
                        }
                        None => key.push(0),
                    }
                }
            }
        }
    }
    Some(key)
}

/// [`triangulate_reporting_from`], through `kept` where there is one.
#[allow(
    clippy::too_many_arguments,
    reason = "the pass's settings, each passed through once"
)]
fn reporting_kept(
    model: &Model,
    face: &Shape,
    deflection: Deflection,
    finer: Option<&EdgeChords>,
    prepared: Option<Trimming>,
    kept: Option<&FaceMeshCache>,
    pass: Pass,
    tol: Tolerances,
) -> OgeomResult<(Triangulation, Verdict)> {
    let Some(kept) = kept else {
        return triangulate_reporting_from(model, face, deflection, finer, prepared, pass, tol);
    };
    let key = mesh_key(model, face, deflection, finer, pass, tol);
    if let Some(key) = &key
        && let Some(found) = kept.get(key)
    {
        return Ok(found);
    }
    let found = triangulate_reporting_from(model, face, deflection, finer, prepared, pass, tol)?;
    if let Some(key) = key {
        kept.put(key, found.clone());
    }
    Ok(found)
}

/// What one pass over a face found.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
enum Verdict {
    /// The boundary enclosed a region and it was drawn.
    Whole,
    /// The boundary crossed itself: what came back is fragments, or nothing.
    /// The edges want drawing finer.
    Short,
}

/// Whether a pass over a face is the last it gets.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
enum Pass {
    /// Drawn again at finer chords if it comes up short: a boundary drawn
    /// coarsely enough to cross itself is drawn finer, not taken apart.
    Agreeing,
    /// Not drawn again: a ring still crossing itself sheds its tails (see
    /// [`crate::tails`]) if that leaves it whole.
    Last,
}

/// One face, with the finer edge chords the whole shape agreed on.
///
/// `None` means this face is being meshed on its own and may work out its
/// own: alone it has no neighbour to disagree with.
fn triangulate_with(
    model: &Model,
    face: &Shape,
    deflection: Deflection,
    finer: Option<&EdgeChords>,
    kept: Option<&FaceMeshCache>,
    tol: Tolerances,
) -> OgeomResult<Triangulation> {
    let (mesh, verdict) =
        reporting_kept(model, face, deflection, finer, None, kept, Pass::Last, tol)?;
    if verdict == Verdict::Short && mesh.triangles.is_empty() {
        ogeom_bail!(
            NotDone,
            "the face's boundary enclosed no triangulable region"
        );
    }
    Ok(mesh)
}

/// The surface's unit normal at `(u, v)`, or where the surface is
/// degenerate there (a cone's apex, a sphere's pole, a patch's collapsed
/// corner), the normal a step inside along the vertex's own column.
///
/// A degenerate point has no normal of its own, and `None` for it would
/// shade the vertex black and drag every normal it is welded with towards
/// nothing. It has a *limit* normal along any line approaching
/// it: the apex of a cone seen up one ruling is that ruling's normal, and
/// a mesh vertex at the apex carries the `u` of the ruling it closes. The
/// step is a millionth of the domain towards its middle, first along
/// `v`, then `u`, then both; `None` only where all three are degenerate.
fn limit_normal(surface: &SurfaceGeometry, u: f64, v: f64, tol: Tolerances) -> Option<Vector> {
    if let Ok(n) = surface.normal_at(u, v, tol) {
        return Some(n.vector());
    }
    let ((ua, ub), (va, vb)) = surface.domain();
    let su = if u < f64::midpoint(ua, ub) { 1.0 } else { -1.0 };
    let sv = if v < f64::midpoint(va, vb) { 1.0 } else { -1.0 };
    // Widening steps: a patch whose corner collapses with its tangent
    // (the control rows drawn together and the next row too) is degenerate
    // to first order for a stretch, and a millionth of the domain is still
    // inside it; on a sphere a millimetre across the tangents a millionth
    // in from the pole cross to less than a direction resolves. A
    // hundredth of the domain is past both, and on a patch that small the
    // normal there is the corner's for every purpose.
    for scale in [1e-6, 1e-4, 1e-2] {
        let du = (ub - ua).abs().max(f64::EPSILON) * scale * su;
        let dv = (vb - va).abs().max(f64::EPSILON) * scale * sv;
        let found = [(u, v + dv), (u + du, v), (u + du, v + dv)]
            .into_iter()
            .find_map(|(nu, nv)| surface.normal_at(nu, nv, tol).ok().map(|n| n.vector()));
        if found.is_some() {
            return found;
        }
    }
    None
}

/// One face, and whether the rings it was built from crossed themselves.
///
/// The flag is how the shape-wide pass learns which faces need their edges
/// drawn finer without building every face's rings twice: the rings are
/// already in hand here, and the sweep over them is the only extra cost a
/// face that does not cross ever pays.
fn triangulate_reporting(
    model: &Model,
    face: &Shape,
    deflection: Deflection,
    finer: Option<&EdgeChords>,
    tol: Tolerances,
) -> OgeomResult<(Triangulation, Verdict)> {
    triangulate_reporting_from(model, face, deflection, finer, None, Pass::Agreeing, tol)
}

/// A face's rings walked ahead of drawing, and the chord its edges want
/// if it is narrower than a few of the caller's.
type Walked = (Trimming, Option<f64>);

/// How a drum or a cone is drawn between its rims, by what its rims say
/// of the face. Its two rims alone where every edge keeps to the face
/// within a hundredth of the chord and meets its neighbour across a
/// crease. A row of the face's own through its middle where an edge
/// strays further, as a fitted intersection in the file does, or the face
/// runs on smoothly into another curved one: drawn from its rims alone, a
/// face shows the rims' slop instead of itself, and at a smooth seam no
/// crease in the mesh says where it ends. Every fourth column of that row
/// where it runs on smoothly into a flat face only.
fn rim_holds(
    model: &Model,
    face: &Shape,
    deflection: Deflection,
    tol: Tolerances,
) -> OgeomResult<Rows> {
    use ogeom_geom::{Curve2d as _, Surface as _};
    let mut rows = Rows::Ends;
    for edge in ogeom_topo::explore(model, face, ogeom_topo::Filter::OfType(ShapeType::Edge))? {
        let Some(data) = model.node(&edge).and_then(|n| n.data().as_edge()) else {
            continue;
        };
        if data.tolerance.get() > deflection.chord * RIM_SLOP {
            return Ok(Rows::Full);
        }
        // The two faces' normals at the edge's middle, each read through
        // its own pcurve: whether the face runs on smoothly across it, and
        // into a flat face or a curved one.
        let mut normals = Vec::new();
        for repr in &data.representations {
            if let ogeom_topo::EdgeRepr::PCurve {
                curve,
                surface,
                range,
                ..
            } = repr
                && let (Some(c), Some(s)) = (
                    model.geometry().pcurve(*curve),
                    model.geometry().surface(*surface),
                )
                && let Ok(uv) = c.point_at(f64::midpoint(range.0, range.1), tol)
                && let Ok(n) = s.normal_at(uv.x, uv.y, tol)
            {
                normals.push((
                    *surface,
                    n,
                    matches!(s.kind(), ogeom_geom::SurfaceKind::Plane),
                ));
            }
        }
        for (i, (sa, na, pa)) in normals.iter().enumerate() {
            for (sb, nb, pb) in &normals[i + 1..] {
                if sa != sb && na.dot(*nb).abs() > SMOOTH_RIM {
                    if *pa || *pb {
                        rows = Rows::Sparse;
                    } else {
                        return Ok(Rows::Full);
                    }
                }
            }
        }
    }
    Ok(rows)
}

/// The widest tolerance of a face's edges and vertices: how far apart two
/// stretches of its boundary may stand and be one.
fn widest_tolerance(model: &Model, face: &Shape) -> OgeomResult<f64> {
    let mut widest = 0.0_f64;
    for kind in [ShapeType::Edge, ShapeType::Vertex] {
        for shape in ogeom_topo::explore_unique(model, face, kind)? {
            widest = widest.max(model.tolerance_of(&shape)?.map_or(0.0, |t| t.get()));
        }
    }
    Ok(widest)
}

/// `planar` with each tail's strip added, wound to meet the region across
/// the chord at the tail's mouth; `None` where the region does not hold
/// that chord.
fn with_strips(mut planar: PlanarMesh, strips: &[crate::tails::Strip]) -> Option<PlanarMesh> {
    let key = |p: Point2| (p.x.to_bits(), p.y.to_bits());
    let mut index: ogeom_core::FastMap<(u64, u64), u32> = ogeom_core::FastMap::default();
    for (i, &(u, v)) in planar.parameters.iter().enumerate() {
        index.insert(key(Point2::new(u, v)), u32::try_from(i).ok()?);
    }
    let mut sides: ogeom_core::FastSet<(u32, u32)> = ogeom_core::FastSet::default();
    for t in &planar.triangles {
        for k in 0..3 {
            sides.insert((t[k], t[(k + 1) % 3]));
        }
    }
    for strip in strips {
        let (from, to) = (
            *index.get(&key(strip.mouth.0))?,
            *index.get(&key(strip.mouth.1))?,
        );
        // The strip walks its mouth from `to` back to `from`; the region
        // must walk it the other way.
        let turned = if sides.contains(&(from, to)) {
            false
        } else if sides.contains(&(to, from)) {
            true
        } else {
            return None;
        };
        for triangle in &strip.triangles {
            let mut corners = [0_u32; 3];
            for (corner, p) in corners.iter_mut().zip(triangle) {
                *corner = match index.get(&key(*p)) {
                    Some(&i) => i,
                    None => {
                        let i = u32::try_from(planar.parameters.len()).ok()?;
                        planar.parameters.push((p.x, p.y));
                        index.insert(key(*p), i);
                        i
                    }
                };
            }
            planar.triangles.push(if turned {
                [corners[0], corners[2], corners[1]]
            } else {
                corners
            });
        }
    }
    Some(planar)
}

/// How a drum or a cone is drawn between its rims.
#[derive(Clone, Copy, PartialEq, Eq, Debug)]
enum Rows {
    /// A row through its middle, a point on every column.
    Full,
    /// That row with a point on every fourth column.
    Sparse,
    /// Its two rims alone.
    Ends,
}

/// The cosine of the turn under which two faces meet smoothly at an edge:
/// thirty degrees, the converter's crease.
const SMOOTH_RIM: f64 = 0.866_025_403_784_438_6;

/// How far an edge may stray from its face, against the chord, and the
/// face still be drawn from its rims alone: a hundredth.
const RIM_SLOP: f64 = 1e-2;

/// The deflection a face is drawn to in its own frame: the caller's chord
/// is a distance where the face stands, and a placement that scales the
/// face scales every chord drawn before it with it.
fn in_own_frame(model: &Model, face: &Shape, deflection: Deflection) -> Deflection {
    let scale = face
        .transform(model.datums())
        .map_or(1.0, |t| t.scale_factor().abs());
    if scale > 0.0 && scale.is_finite() && (scale - 1.0).abs() > f64::EPSILON {
        Deflection {
            chord: deflection.chord / scale,
            ..deflection
        }
    } else {
        deflection
    }
}

/// [`triangulate_reporting`] with the rings already walked, where the
/// caller has them and none of the face's edges were told to draw finer
/// since; the rings depend on nothing else.
fn triangulate_reporting_from(
    model: &Model,
    face: &Shape,
    deflection: Deflection,
    finer: Option<&EdgeChords>,
    prepared: Option<Trimming>,
    pass: Pass,
    tol: Tolerances,
) -> OgeomResult<(Triangulation, Verdict)> {
    deflection.validate()?;
    if model.kind_of(face)? != ShapeType::Face {
        ogeom_bail!(Construction, "expected a face");
    }
    let deflection = in_own_frame(model, face, deflection);
    let Some(node) = model.node(face) else {
        ogeom_bail!(Dangling, "face is not in this model");
    };
    let NodeData::Face(data) = node.data() else {
        ogeom_bail!(Construction, "face node holds no face data");
    };
    let Some(surface) = model.geometry().surface(data.surface) else {
        ogeom_bail!(Dangling, "face refers to a surface not in this model");
    };
    let placement = face.transform(model.datums())?;

    let own;
    let finer = match finer {
        Some(shared) => shared,
        None => {
            own = face_chords(
                model,
                face,
                data.surface,
                surface,
                deflection,
                &EdgeChords::new(),
                tol,
            )?;
            &own
        }
    };
    let phase = Instant::now();
    let Trimming {
        rings: uv,
        anchors,
        met,
        walked,
    } = match prepared {
        Some(trim) => trim,
        None => trimming_rings(
            model,
            face,
            data.surface,
            surface,
            deflection,
            finer,
            false,
            tol,
        )?,
    };
    let rings_ms = phase.elapsed().as_secs_f64() * 1e3;
    if uv.is_empty() {
        // Bounded, and bounding next to nothing: a sliver whose sides are
        // one line in the chart. It is drawn as a fan across each ring its
        // edges walked, at the points its neighbours share along those
        // edges, so the whole mesh still closes over it.
        let reversed = face.orientation() == Orientation::Reversed;
        return Ok((
            sliver_fan(&walked, surface, &placement, reversed, tol),
            Verdict::Whole,
        ));
    }
    let phase = Instant::now();
    let rows = rim_holds(model, face, deflection, tol)?;
    let mut planar = triangulate_region(&uv, surface, deflection, rows, tol)?;
    // A ring that crosses itself where it runs out and back along its own
    // way bounds nothing over that stretch: the tail comes off and is drawn
    // as a strip at its edges' points, where the rest triangulates whole.
    let mut drawn = None;
    if planar.crossed
        && pass == Pass::Last
        && let Some(peeled) = crate::tails::peel(
            &uv,
            surface,
            &placement,
            widest_tolerance(model, face)?,
            tol,
        )?
    {
        let again = triangulate_region(&peeled.rings, surface, deflection, rows, tol)?;
        if !again.crossed
            && let Some(joined) = with_strips(again, &peeled.strips)
        {
            planar = joined;
            drawn = Some(peeled.rings);
        }
    }
    let drawn = drawn.as_ref().unwrap_or(&uv);
    let region_ms = phase.elapsed().as_secs_f64() * 1e3;
    let phase = Instant::now();

    // Boundary-only triangulation checks the exact count before interior
    // points can hide missing triangles. Its minimum includes shared junctions
    // between touching lobes, which reduce the simple polygon's count.
    // Filtering degenerate triangles must still leave that minimum.
    //
    // And a face narrower than a few chords is short too, whatever its
    // count: drawn at the caller's chord its boundary sags by more than
    // the face is wide, and every triangle across the width stands off the
    // surface by that sag. Its edges want a chord under the width, which
    // is [`face_chords`]' first move; it is asked here so the whole-shape
    // pass tells the neighbours to draw those edges finer too.
    let boundary: usize = drawn.iter().map(Vec::len).sum();
    let narrow = narrow_chord(surface, drawn, deflection, tol)
        .is_some_and(|want| !edges_already_at(model, face, finer, want));
    let crossed = narrow || planar.crossed || planar.triangles.len() < planar.boundary_triangles;
    if *MESH_DEBUG_REFINE {
        eprintln!(
            "SHORT {} triangles against {boundary} boundary points in {} rings: crossed {crossed} (boundary pass {})",
            planar.triangles.len(),
            drawn.len(),
            planar.crossed
        );
    }

    // Boundary vertices take their positions from their edges' own curves
    // (the shared authority), keyed by their exact parameter-space bits.
    let mut anchored: ogeom_core::FastMap<(u64, u64), Point> = ogeom_core::FastMap::default();
    for (ring, ring_anchor) in uv.iter().zip(&anchors) {
        for (p, a) in ring.iter().zip(ring_anchor) {
            if let Some(point) = a {
                anchored.insert((p.x.to_bits(), p.y.to_bits()), *point);
            }
        }
    }

    // Lift into space. The normal follows the face's orientation, not the
    // surface's: a reversed face presents the other side, and a renderer or a
    // volume computation that ignored that would have the solid inside out.
    // A reflecting placement carries the surface's normal to the image's
    // as it is, but turns the chart over: triangles wound counter-clockwise
    // in the chart wind the other way in space, so the winding turns once
    // more and the normal does not.
    let reversed = face.orientation() == Orientation::Reversed;
    let flip = reversed != !face.location().preserves_handedness(model.datums())?;
    let mut mesh = Triangulation::new();
    mesh.deflection_met = met;
    let mut hits = 0_usize;
    for (u, v) in planar.parameters {
        // Anchors are already in world coordinates (their edges' own
        // placements applied), where surface lifts still need the face's.
        let point = match anchored.get(&(u.to_bits(), v.to_bits())) {
            Some(anchor) => {
                hits += 1;
                *anchor
            }
            None => placement.apply(surface.point_at(u, v, tol)?),
        };
        let normal =
            limit_normal(surface, u, v, tol).map_or(Vector::ZERO, |n| placement.apply_vector(n));
        mesh.positions.push(point);
        mesh.normals.push(if reversed { -normal } else { normal });
        mesh.parameters.push((u, v));
    }
    if *MESH_DEBUG_REFINE && (rings_ms + region_ms) > 50.0 {
        eprintln!(
            "PHASE rings {rings_ms:.0}ms region {region_ms:.0}ms lift {:.0}ms  ({} ring points, {} tris)",
            phase.elapsed().as_secs_f64() * 1e3,
            uv.iter().map(Vec::len).sum::<usize>(),
            planar.triangles.len()
        );
    }
    if *MESH_DEBUG {
        eprintln!(
            "DBG anchors map={} hits={} verts={}",
            anchored.len(),
            hits,
            mesh.positions.len()
        );
    }
    mesh.triangles = planar
        .triangles
        .into_iter()
        .map(|t| if flip { [t[0], t[2], t[1]] } else { t })
        .collect();
    split_flat_triangles(&mut mesh, tol);
    // The chord is met where the mesh says it is. The grid holds each row
    // and column to it, and a cell's middle can stand off by both sides'
    // sags together: each triangle's middle, lifted onto the surface, is
    // measured against the triangle's plane, where the face stands, and
    // against the triangle's own middle along the surface's normal there.
    // The plane alone is blind to a sliver laid along a line the surface
    // bends under: its plane holds the curve its corners lie on, and its
    // middle stands inside the bend by most of the bend's sag.
    if mesh.deflection_met && !matches!(surface.kind(), ogeom_geom::SurfaceKind::Plane) {
        let chord = deflection.chord * placement.scale_factor().abs();
        for t in &mesh.triangles {
            let [a, b, c] = t.map(|i| i as usize);
            let (pa, pb, pc) = (mesh.positions[a], mesh.positions[b], mesh.positions[c]);
            let Ok(normal) = (pb - pa).cross(pc - pa).normalized(tol) else {
                continue;
            };
            let (ua, ub, uc) = (mesh.parameters[a], mesh.parameters[b], mesh.parameters[c]);
            let middle = ((ua.0 + ub.0 + uc.0) / 3.0, (ua.1 + ub.1 + uc.1) / 3.0);
            let Ok(on) = surface.point_at(middle.0, middle.1, tol) else {
                continue;
            };
            let on = placement.apply(on);
            let mut off = (on - pa).dot(normal).abs();
            if let Ok(across) = surface.normal_at(middle.0, middle.1, tol)
                && let Ok(across) = placement.apply_vector(across.vector()).normalized(tol)
            {
                let centre =
                    Point::from_vector((pa.to_vector() + pb.to_vector() + pc.to_vector()) / 3.0);
                off = off.max((on - centre).dot(across).abs());
            }
            if off > chord + tol.confusion() {
                mesh.deflection_met = false;
                break;
            }
        }
    }
    Ok((
        mesh,
        if crossed {
            Verdict::Short
        } else {
            Verdict::Whole
        },
    ))
}

/// Triangulate every face below a shape, welded into one mesh.
///
/// # Errors
///
/// As [`triangulate_face`].
pub fn triangulate(
    model: &Model,
    shape: &Shape,
    deflection: Deflection,
    tol: Tolerances,
) -> OgeomResult<Triangulation> {
    Ok(triangulate_with_chords(model, shape, deflection, tol)?.0)
}

/// As [`triangulate`], with the edge chords the faces agreed on: what
/// [`edge_chords_for`] answers for the same shape and deflection, without
/// drawing every face a second time to find it.
///
/// # Errors
///
/// As [`triangulate_face`].
pub fn triangulate_with_chords(
    model: &Model,
    shape: &Shape,
    deflection: Deflection,
    tol: Tolerances,
) -> OgeomResult<(Triangulation, EdgeChords)> {
    let (mesh, chords, _) = triangulate_with_face_meshes(model, shape, deflection, None, tol)?;
    Ok((mesh, chords))
}

/// As [`triangulate_with_chords`], with every face's own mesh beside the
/// welded whole: each face below `shape` in exploration order, and what
/// [`triangulate_face_with`] answers for it with the returned chords. The
/// whole is welded from those same meshes, so they are drawn once. A face
/// `kept` has drawn before to the same chords is taken from it.
///
/// # Errors
///
/// As [`triangulate_face`].
pub fn triangulate_with_face_meshes(
    model: &Model,
    shape: &Shape,
    deflection: Deflection,
    kept: Option<&FaceMeshCache>,
    tol: Tolerances,
) -> OgeomResult<(Triangulation, EdgeChords, Vec<FaceMesh>)> {
    // Faces in two phases, as `tessellate` does: each face is meshed from a
    // model nothing is writing to, in parallel; the pieces are then appended
    // sequentially in face order. The split is what keeps the answer
    // bit-identical at any thread count: scheduling decides only who does
    // which face, never where its triangles land.
    let (faces, chords) = face_meshes_with_chords(model, shape, deflection, kept, tol)?;

    let mut mesh = Triangulation::new();
    let mut pieces: Vec<(usize, usize)> = Vec::with_capacity(faces.len());
    for (_, piece) in &faces {
        let piece = piece.as_ref().map_err(Clone::clone)?;
        let t0 = mesh.triangles.len();
        mesh.append(piece);
        pieces.push((t0, mesh.triangles.len()));
    }
    orient_pieces(&mut mesh, &pieces);
    let mesh = mesh.welded(tol);

    // A second, border-only pass at the tolerance the model itself recorded:
    // imported edges carry the file's slop in their widened tolerances, and
    // two faces lifting the same edge through disagreeing geometry land that
    // far apart. Interior edges are already manifold and are not touched.
    let mut reach = 0.0_f64;
    for kind in [ShapeType::Edge, ShapeType::Vertex] {
        for shape in ogeom_topo::explore(model, shape, ogeom_topo::Filter::OfType(kind))? {
            let recorded = model.node(&shape).map_or(0.0, |n| match n.data() {
                NodeData::Edge(d) => d.tolerance.get(),
                NodeData::Vertex(d) => d.tolerance.get(),
                _ => 0.0,
            });
            reach = reach.max(recorded);
        }
    }
    // Floored at a tenth of the chord the caller asked for: a border gap
    // smaller than that is below the resolution of the mesh they accepted,
    // whether or not the model recorded the slop that caused it.
    let reach = reach.max(deflection.chord * 0.1);
    let mesh = if reach > tol.confusion() {
        let reach = reach + tol.confusion();
        mesh.border_welded(reach).border_stitched(reach)
    } else {
        mesh
    };
    // What is left open narrower than the chord asked for is two faces
    // sampling a shared corner differently, below the mesh's own
    // resolution, and is sealed; a wider opening is the shape's.
    let faces = faces
        .into_iter()
        .map(|(face, piece)| {
            let piece =
                piece.and_then(|piece| own_answer(model, &face, piece, deflection, &chords, tol));
            (face, piece)
        })
        .collect();
    Ok((mesh.sealed(deflection.chord), chords, faces))
}

/// A face and its mesh, or why it would not draw.
pub type FaceMesh = (Shape, OgeomResult<Triangulation>);

/// Every face below `shape` in exploration order, each with what
/// [`triangulate_face_with`] answers for it with the chords the faces agree
/// on (what [`edge_chords_for`] answers), and those chords: the agreement
/// and the meshes drawn in one pass. A face `kept` has drawn before to the
/// same chords is taken from it.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the
/// deflection is unusable; a face that will not draw carries its own error.
pub fn face_meshes_for(
    model: &Model,
    shape: &Shape,
    deflection: Deflection,
    kept: Option<&FaceMeshCache>,
    tol: Tolerances,
) -> OgeomResult<(Vec<FaceMesh>, EdgeChords)> {
    deflection.validate()?;
    let (faces, chords) = face_meshes_with_chords(model, shape, deflection, kept, tol)?;
    let faces = faces
        .into_iter()
        .map(|(face, piece)| {
            let piece =
                piece.and_then(|piece| own_answer(model, &face, piece, deflection, &chords, tol));
            (face, piece)
        })
        .collect();
    Ok((faces, chords))
}

/// The faces below `shape` paired with [`face_meshes`]' answers for them.
fn face_meshes_with_chords(
    model: &Model,
    shape: &Shape,
    deflection: Deflection,
    kept: Option<&FaceMeshCache>,
    tol: Tolerances,
) -> OgeomResult<(Vec<FaceMesh>, EdgeChords)> {
    let faces: Vec<Shape> =
        ogeom_topo::explore(model, shape, ogeom_topo::Filter::OfType(ShapeType::Face))?;
    let (computed, chords) = face_meshes(model, &faces, deflection, kept, tol)?;
    Ok((faces.into_iter().zip(computed).collect(), chords))
}

/// A face's mesh from the shape-wide pass as [`triangulate_face_with`]
/// answers it: the pass keeps a face that drew nothing, where the face
/// asked alone is refused, so such a face is asked again.
fn own_answer(
    model: &Model,
    face: &Shape,
    piece: Triangulation,
    deflection: Deflection,
    chords: &EdgeChords,
    tol: Tolerances,
) -> OgeomResult<Triangulation> {
    if piece.triangles.is_empty() {
        triangulate_face_with(model, face, deflection, chords, tol)
    } else {
        Ok(piece)
    }
}

/// Every face below a shape drawn to the chords the faces agree on, in face
/// order, and those chords.
///
/// The first pass draws each face at the caller's chord on top of the
/// narrow faces' asks; a face whose boundary crossed itself folds finer
/// chords into the agreement, and only the faces touching an edge whose
/// chord changed are drawn again. A face still short after that is drawn
/// once more at the same chords as its last pass ([`Pass::Last`]). The
/// whole-shape mesh welds these, and the stored tessellation keeps them
/// face by face: the same meshes, drawn once.
pub(crate) fn face_meshes(
    read_model: &Model,
    faces: &[Shape],
    deflection: Deflection,
    kept: Option<&FaceMeshCache>,
    tol: Tolerances,
) -> OgeomResult<(Vec<OgeomResult<Triangulation>>, EdgeChords)> {
    let FirstPass {
        mut finer,
        mut computed,
        mut changed,
        short,
        mut drawn,
    } = first_pass(read_model, faces, deflection, kept, tol)?;
    let mut still = vec![false; faces.len()];
    for index in short {
        still[index] = true;
    }
    // Every face across a changed edge is drawn again at the new chords,
    // and asked again whether its rings crossed: a narrow strip drawn whole
    // with both its rims coarse can cross once one rim is drawn finer by
    // the face across it, the coarse rim's chords sagging over the fine
    // one. Its edges are then drawn finer in turn, and so on until no face
    // comes up short, or the rounds run out; the last round's meshes stand.
    for _ in 0..SETTLING_ROUNDS {
        if changed.is_empty() {
            break;
        }
        let again: Vec<usize> = (0..faces.len())
            .filter(|&i| {
                ogeom_topo::explore(
                    read_model,
                    &faces[i],
                    ogeom_topo::Filter::OfType(ShapeType::Edge),
                )
                .is_ok_and(|es| es.iter().any(|e| changed.contains(&e.node().index())))
            })
            .collect();
        let redone: Vec<OgeomResult<(Triangulation, Verdict)>> =
            ogeom_core::parallel::map_ordered(&again, |_, &index| {
                ogeom_core::progress::checkpoint()?;
                reporting_kept(
                    read_model,
                    &faces[index],
                    deflection,
                    Some(&finer),
                    None,
                    kept,
                    Pass::Agreeing,
                    tol,
                )
            });
        let mut short: Vec<usize> = Vec::new();
        for (index, one) in again.into_iter().zip(redone) {
            still[index] = matches!(one, Ok((_, Verdict::Short)));
            computed[index] = match one {
                Ok((mesh, Verdict::Short)) if mesh.triangles.is_empty() => {
                    Err(ogeom_core::ogeom_err!(
                        NotDone,
                        "the face's boundary enclosed no triangulable region"
                    ))
                }
                Ok((mesh, verdict)) => {
                    if verdict == Verdict::Short {
                        short.push(index);
                    }
                    Ok(mesh)
                }
                Err(e) => Err(e),
            };
        }
        drawn = finer.clone();
        changed = refine_short(read_model, faces, &short, &mut finer, deflection, kept, tol)?;
    }
    // A face still short is drawn once more at the chords it was drawn to,
    // as its last pass. The chords stay as they are, and so do its
    // neighbours. Its mesh is taken only where that pass comes out whole.
    let last: Vec<usize> = (0..faces.len()).filter(|&i| still[i]).collect();
    let redone: Vec<OgeomResult<(Triangulation, Verdict)>> =
        ogeom_core::parallel::map_ordered(&last, |_, &index| {
            ogeom_core::progress::checkpoint()?;
            reporting_kept(
                read_model,
                &faces[index],
                deflection,
                Some(&drawn),
                None,
                kept,
                Pass::Last,
                tol,
            )
        });
    for (index, one) in last.into_iter().zip(redone) {
        if let Ok((mesh, Verdict::Whole)) = one {
            computed[index] = Ok(mesh);
        }
    }
    Ok((computed, finer))
}

/// How many times the faces across changed edges are drawn again.
const SETTLING_ROUNDS: usize = 4;

/// The chords each face of `short` (indices into `faces`) wants its edges
/// drawn at, merged into `finer` where finer than held: the edges whose
/// chord changed.
fn refine_short(
    read_model: &Model,
    faces: &[Shape],
    short: &[usize],
    finer: &mut EdgeChords,
    deflection: Deflection,
    kept: Option<&FaceMeshCache>,
    tol: Tolerances,
) -> OgeomResult<ogeom_core::FastSet<u32>> {
    let mut changed: ogeom_core::FastSet<u32> = ogeom_core::FastSet::default();
    // Each face asked against the chords held before any of them answers,
    // so the answers do not hang on the order the faces are taken in.
    // That also lets them be asked side by side; the answers are merged
    // in the faces' order, the first error standing as it would one by one.
    let held = finer.clone();
    let asked: Vec<OgeomResult<Option<EdgeChords>>> =
        ogeom_core::parallel::map_ordered(short, |_, &index| {
            ogeom_core::progress::checkpoint()?;
            let face = &faces[index];
            let Some(node) = read_model.node(face) else {
                return Ok(None);
            };
            let NodeData::Face(data) = node.data() else {
                return Ok(None);
            };
            let Some(surface) = read_model.geometry().surface(data.surface) else {
                return Ok(None);
            };
            face_chords_kept(
                read_model,
                face,
                data.surface,
                surface,
                deflection,
                &held,
                kept,
                tol,
            )
            .map(Some)
        });
    for answer in asked {
        let Some(chords) = answer? else {
            continue;
        };
        for (edge, chord) in chords {
            let held = finer.entry(edge).or_insert(f64::INFINITY);
            if chord < *held {
                *held = chord;
                changed.insert(edge);
            }
        }
    }
    Ok(changed)
}

/// What the first pass over a shape's faces produced: the chords its faces
/// agreed to draw their shared edges to, the meshes drawn at the caller's
/// chord, and the edges whose chord the crossed faces changed after those
/// meshes were drawn; the faces touching one of them are drawn again.
struct FirstPass {
    finer: EdgeChords,
    computed: Vec<OgeomResult<Triangulation>>,
    changed: ogeom_core::FastSet<u32>,
    /// The faces that came up short.
    short: Vec<usize>,
    /// The chords the meshes were drawn to, before the short faces'.
    drawn: EdgeChords,
}

/// Every face drawn once at the caller's deflection, each saying whether
/// the rings it was drawn from crossed themselves, on top of the chords
/// the narrow faces asked for. Nearly none cross, and those faces are
/// finished. A face whose boundary crossed itself needs its edges drawn
/// finer, and so does every face that shares one of them, or the two
/// sides of that edge arrive with a different number of points, which is
/// a worse crack than the sliver the refinement was for. Their chords
/// are folded into the map and the edges that changed are named.
fn first_pass(
    read_model: &Model,
    faces: &[Shape],
    deflection: Deflection,
    cache: Option<&FaceMeshCache>,
    tol: Tolerances,
) -> OgeomResult<FirstPass> {
    // `Some(&finer)`, never `None`: a face left to itself refines its own
    // edges, which is right when it is meshed alone and wrong here, where
    // its neighbours must be told to refine the same ones. Every face is
    // drawn at exactly what the caller asked and reports what crossed.
    //
    // Before that, every face's rings are walked at the caller's chord and
    // its width read off them: a face narrower than a few chords wants its
    // edges drawn finer, and so do the faces across those edges. Known
    // before anything is drawn, those chords go into this pass, and the
    // big faces round a thousand small fillets are drawn once rather than
    // once and again. The rings are kept for the faces they still
    // describe: every face none of whose edges the map names.
    let (mut finer, kept) = walk_for_chords(read_model, faces, deflection, tol);
    let touched = |face: &Shape| -> bool {
        ogeom_topo::explore(
            read_model,
            face,
            ogeom_topo::Filter::OfType(ShapeType::Edge),
        )
        .is_ok_and(|es| es.iter().any(|e| finer.contains_key(&e.node().index())))
    };
    // The rings are handed over by the job that draws the face; a shared
    // slice cannot give them away, so each sits behind a lock it is taken
    // from once.
    let jobs: Vec<(&Shape, std::sync::Mutex<Option<Trimming>>)> = faces
        .iter()
        .zip(kept)
        .map(|(face, prep)| {
            let trim = prep.filter(|_| !touched(face));
            (face, std::sync::Mutex::new(trim))
        })
        .collect();
    let first: Vec<OgeomResult<(Triangulation, Verdict)>> =
        ogeom_core::parallel::map_ordered(&jobs, |_, (face, slot)| {
            ogeom_core::progress::checkpoint()?;
            let trim = slot.lock().ok().and_then(|mut held| held.take());
            reporting_kept(
                read_model,
                face,
                deflection,
                Some(&finer),
                trim,
                cache,
                Pass::Agreeing,
                tol,
            )
        });
    let mut computed: Vec<OgeomResult<Triangulation>> = Vec::with_capacity(faces.len());
    let mut crossed: Vec<usize> = Vec::new();
    for (index, one) in first.into_iter().enumerate() {
        match one {
            Ok((mesh, Verdict::Short)) => {
                crossed.push(index);
                computed.push(Ok(mesh));
            }
            Ok((mesh, Verdict::Whole)) => computed.push(Ok(mesh)),
            Err(e) => computed.push(Err(e)),
        }
    }
    if *MESH_DEBUG_REFINE && !crossed.is_empty() {
        eprintln!(
            "REFINE {} of {} faces came up short",
            crossed.len(),
            faces.len()
        );
    }
    // On top of the walk's map, not instead of it: a neighbour drawn again
    // must still draw the edges the walk held finer at that chord, or the
    // two sides of one of them disagree. Only an edge whose chord
    // *changed* sends its faces back.
    let drawn = finer.clone();
    let changed = refine_short(
        read_model, faces, &crossed, &mut finer, deflection, cache, tol,
    )?;
    Ok(FirstPass {
        finer,
        computed,
        changed,
        short: crossed,
        drawn,
    })
}

/// Every face's rings walked at the caller's chord, and the chords the
/// narrow faces want their edges drawn to, agreed across the shape: the
/// finest any face asks of an edge is what every face draws it to.
///
/// The rings are kept for the faces they still describe (every face none
/// of whose edges the map names), so the whole-shape pass draws those once.
fn walk_for_chords(
    read_model: &Model,
    faces: &[Shape],
    deflection: Deflection,
    tol: Tolerances,
) -> (EdgeChords, Vec<Option<Trimming>>) {
    let nothing = EdgeChords::new();
    let prepared: Vec<OgeomResult<Option<Walked>>> =
        ogeom_core::parallel::map_ordered(faces, |_, face| {
            ogeom_core::progress::checkpoint()?;
            let Some(node) = read_model.node(face) else {
                return Ok(None);
            };
            let NodeData::Face(data) = node.data() else {
                return Ok(None);
            };
            let Some(surface) = read_model.geometry().surface(data.surface) else {
                return Ok(None);
            };
            let deflection = in_own_frame(read_model, face, deflection);
            let trim = trimming_rings(
                read_model,
                face,
                data.surface,
                surface,
                deflection,
                &nothing,
                false,
                tol,
            )?;
            let narrow = narrow_chord(surface, &trim.rings, deflection, tol);
            Ok(Some((trim, narrow)))
        });
    let mut finer = EdgeChords::new();
    let mut kept: Vec<Option<Trimming>> = Vec::with_capacity(faces.len());
    for (face, one) in faces.iter().zip(prepared) {
        match one {
            Ok(Some((prep, narrow))) => {
                if let Some(chord) = narrow
                    && let Ok(edges) = ogeom_topo::explore(
                        read_model,
                        face,
                        ogeom_topo::Filter::OfType(ShapeType::Edge),
                    )
                {
                    for edge in edges {
                        let held = finer.entry(edge.node().index()).or_insert(chord);
                        *held = held.min(chord);
                    }
                }
                kept.push(Some(prep));
            }
            Ok(None) | Err(_) => kept.push(None),
        }
    }
    (finer, kept)
}

/// The chords the faces below `shape` agree to draw their shared edges to.
///
/// A face narrower than a few chords draws its edges finer than the
/// caller asked, and so does a face whose boundary crosses itself at that
/// chord; the face across each of those edges must draw it the same or
/// the two meshes disagree along it: a crack in every mesh assembled
/// face by face. [`triangulate`] agrees this for itself; a caller meshing
/// face by face (a viewer keeping one mesh per face) asks here once per
/// shape and hands the answer to [`triangulate_face_with`] for every face.
///
/// The answer costs a draw of every face at the caller's chord, since a
/// boundary is only known to cross itself once it is triangulated; a
/// caller meshing face by face pays that draw twice over, and one that
/// wants only the welded whole should ask [`triangulate`] instead.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the
/// deflection is unusable; a face that will not draw asks for nothing
/// rather than failing the shape.
pub fn edge_chords_for(
    model: &Model,
    shape: &Shape,
    deflection: Deflection,
    tol: Tolerances,
) -> OgeomResult<EdgeChords> {
    deflection.validate()?;
    let faces: Vec<Shape> =
        ogeom_topo::explore(model, shape, ogeom_topo::Filter::OfType(ShapeType::Face))?;
    Ok(first_pass(model, &faces, deflection, None, tol)?.finer)
}

/// As [`triangulate_face`], with the edge chords the shape agreed on (the
/// answer of [`edge_chords_for`]), so the face's boundary matches its
/// neighbours' point for point along every edge the map names.
///
/// The agreement already holds every refinement the face's own boundary
/// asked for, so the face is drawn to it and left there, as the
/// whole-shape triangulation leaves it: refining alone past the agreed
/// chord, as [`triangulate_face`] would, is what the neighbours cannot
/// follow.
///
/// # Errors
///
/// As [`triangulate_face`].
pub fn triangulate_face_with(
    model: &Model,
    face: &Shape,
    deflection: Deflection,
    chords: &EdgeChords,
    tol: Tolerances,
) -> OgeomResult<Triangulation> {
    triangulate_with(model, face, deflection, Some(chords), None, tol)
}

/// Make the appended face meshes traverse their shared boundaries in
/// opposite directions, flipping as few faces as possible.
///
/// A closed oriented surface walks each interior edge once each way. A file
/// whose face orientation flags disagree with each other still *pairs* every
/// edge (closed by counting, two-sided nowhere), and every integral over the
/// result is reference-dependent garbage. The face flags cannot be judged
/// from the model's wires (a synthesised seam's occurrence directions are
/// chart bookkeeping, not 3D traversal), but the meshes tell the truth:
/// shared boundary vertices are anchored to their edges' own curves, so two
/// faces meeting along an edge carry bitwise-identical positions, and the
/// direction each walks them is right there in the triangles. Faces are
/// flood-filled across those shared runs, each constrained to oppose its
/// neighbour, and each connected component keeps the polarity that flips the
/// minority; flipping means reversing the piece's windings and negating its
/// normals. Conflicts are left standing: this repairs orientation, not
/// topology.
fn orient_pieces(mesh: &mut Triangulation, pieces: &[(usize, usize)]) {
    use ogeom_core::FastMap;
    type Key = (u64, u64, u64);
    let key = |p: &Point| -> Key { (p.x.to_bits(), p.y.to_bits(), p.z.to_bits()) };

    // Directed boundary edges per piece, keyed by position: only each
    // piece's *border* edges (used once within the piece) face other pieces.
    let mut owners: FastMap<(Key, Key), Vec<(usize, bool)>> = FastMap::default();
    for (i, &(t0, t1)) in pieces.iter().enumerate() {
        let mut inside: FastMap<(u32, u32), u32> = FastMap::default();
        for t in &mesh.triangles[t0..t1] {
            for k in 0..3 {
                let (a, b) = (t[k], t[(k + 1) % 3]);
                *inside.entry((a.min(b), a.max(b))).or_default() += 1;
            }
        }
        for t in &mesh.triangles[t0..t1] {
            for k in 0..3 {
                let (a, b) = (t[k], t[(k + 1) % 3]);
                if inside.get(&(a.min(b), a.max(b))).copied() != Some(1) {
                    continue;
                }
                let (ka, kb) = (
                    key(&mesh.positions[a as usize]),
                    key(&mesh.positions[b as usize]),
                );
                // A side whose corners lie on one point (a triangle closing
                // at an apex, or a fan ruled onto a seam that shrinks to a
                // point) is walked no way at all, and says nothing of
                // which way its pieces turn.
                if ka == kb {
                    continue;
                }
                // One undirected key, direction recorded: `true` where this
                // piece walks the smaller key first.
                let (lo, hi, forward) = if ka <= kb {
                    (ka, kb, true)
                } else {
                    (kb, ka, false)
                };
                owners.entry((lo, hi)).or_default().push((i, forward));
            }
        }
    }

    // Manifold constraints: exactly two pieces sharing a run. Same recorded
    // direction means exactly one must flip. A pair of faces shares many
    // runs and slop can corrupt a few, so each pair's constraint is the
    // majority of its runs, and ties abstain rather than guess. Everything
    // is aggregated in sorted order, so the answer does not depend on hash
    // iteration.
    let mut votes: std::collections::BTreeMap<(usize, usize), (u32, u32)> =
        std::collections::BTreeMap::new();
    for list in owners.values() {
        if let [(a, fa), (b, fb)] = list[..]
            && a != b
        {
            let pair = (a.min(b), a.max(b));
            let entry = votes.entry(pair).or_insert((0, 0));
            if fa == fb {
                entry.0 += 1;
            } else {
                entry.1 += 1;
            }
        }
    }
    // Strongest constraints first: a pair vouched for by many runs outranks
    // one hanging on a sliver, so when the graph carries a contradiction
    // (an odd cycle born of slop), the weakest link is the one disbelieved.
    // Union-find with parity keeps the whole resolution order-independent.
    // A pair whose runs split evenly still says the faces touch: it joins
    // the components with the consistent-shell prior, at zero strength, so
    // a fragment cannot drift off and choose its polarity alone.
    let mut constraints: Vec<(u32, usize, usize, bool)> = votes
        .iter()
        .map(|(&(a, b), &(same, opposite))| (same.abs_diff(opposite), a, b, same > opposite))
        .collect();
    constraints.sort_by(|x, y| y.0.cmp(&x.0).then(x.1.cmp(&y.1)).then(x.2.cmp(&y.2)));

    let mut parent: Vec<usize> = (0..pieces.len()).collect();
    // Whether a piece is flipped relative to its parent.
    let mut parity: Vec<bool> = vec![false; pieces.len()];
    fn find(parent: &mut [usize], parity: &mut [bool], i: usize) -> (usize, bool) {
        if parent[i] == i {
            return (i, false);
        }
        let (root, above) = find(parent, parity, parent[i]);
        parent[i] = root;
        parity[i] ^= above;
        (root, parity[i])
    }
    for &(_, a, b, same_direction) in &constraints {
        let (ra, pa) = find(&mut parent, &mut parity, a);
        let (rb, pb) = find(&mut parent, &mut parity, b);
        // One of a same-direction pair flips: their parities must differ.
        let need = same_direction;
        if ra == rb {
            // Agreeing or contradicting, the die is cast; a contradiction
            // here lost to stronger evidence.
            continue;
        }
        parent[rb] = ra;
        parity[rb] = (pa != pb) != need;
    }

    let mut flip: Vec<bool> = vec![false; pieces.len()];
    let mut components: std::collections::BTreeMap<usize, Vec<usize>> =
        std::collections::BTreeMap::new();
    for (i, f) in flip.iter_mut().enumerate() {
        let (root, p) = find(&mut parent, &mut parity, i);
        *f = p;
        components.entry(root).or_default().push(i);
    }
    for members in components.values() {
        let flipped = members.iter().filter(|&&i| flip[i]).count();
        if flipped * 2 > members.len() {
            for &i in members {
                flip[i] = !flip[i];
            }
        }
    }

    for (i, &(t0, t1)) in pieces.iter().enumerate() {
        if !flip[i] {
            continue;
        }
        let mut flipped_vertices: Vec<u32> = Vec::new();
        for t in &mut mesh.triangles[t0..t1] {
            t.swap(1, 2);
            flipped_vertices.extend_from_slice(&t[..]);
        }
        flipped_vertices.sort_unstable();
        flipped_vertices.dedup();
        for v in flipped_vertices {
            mesh.normals[v as usize] = -mesh.normals[v as usize];
        }
    }
}

/// A face's trimming boundary, in its surface's parameter space.
///
/// The outer wire first, then any holes, each as a closed ring of `(u, v)`
/// points with no repeated closing point. A face with no wires covers its
/// surface's whole domain, and gets that rectangle as its boundary.
///
/// Public because parameter-space trimming is not only the triangulator's
/// concern: classifying a point against a face, splitting a face in a boolean,
/// and hidden-line removal all ask the same question of the same rings, and
/// each deriving them separately would be three chances to disagree.
///
/// # Errors
///
/// As [`triangulate_face`].
pub fn face_boundary(
    model: &Model,
    face: &Shape,
    deflection: Deflection,
    tol: Tolerances,
) -> OgeomResult<Vec<Vec<Point2>>> {
    drawn_boundary(model, face, deflection, false, tol)
}

/// [`face_boundary`], each ring keeping within the chord of its edges in
/// the chart too, read in space through the surface's tangent scales: a
/// point on an edge then lies within the chord of its ring however much
/// more the edge bends in the chart than in space.
///
/// The rings to place points against. [`face_boundary`]'s are drawn at the
/// points a mesh welds its faces along, and keep to their edges in space
/// only: a parabola where a cone meets a plane drafted as the cone is bends
/// little in space and much in the cone's chart, and its ring stands off it
/// there by more than the chord.
///
/// # Errors
///
/// As [`face_boundary`].
pub fn face_boundary_in_chart(
    model: &Model,
    face: &Shape,
    deflection: Deflection,
    tol: Tolerances,
) -> OgeomResult<Vec<Vec<Point2>>> {
    drawn_boundary(model, face, deflection, true, tol)
}

/// [`face_boundary_in_chart`] where every edge's polyline honours the
/// deflection, and `None` where one reached [`Deflection::max_segments`]
/// first: rings that keep within the chord are what a point can be placed
/// against at that chord.
///
/// # Errors
///
/// As [`face_boundary`].
pub fn face_boundary_in_chart_met(
    model: &Model,
    face: &Shape,
    deflection: Deflection,
    tol: Tolerances,
) -> OgeomResult<Option<Vec<Vec<Point2>>>> {
    let drawn = drawn_trimming(model, face, deflection, true, tol)?;
    Ok(drawn.met.then_some(drawn.rings))
}

/// [`face_boundary`], or with `in_chart` [`face_boundary_in_chart`].
fn drawn_boundary(
    model: &Model,
    face: &Shape,
    deflection: Deflection,
    in_chart: bool,
    tol: Tolerances,
) -> OgeomResult<Vec<Vec<Point2>>> {
    Ok(drawn_trimming(model, face, deflection, in_chart, tol)?.rings)
}

fn drawn_trimming(
    model: &Model,
    face: &Shape,
    deflection: Deflection,
    in_chart: bool,
    tol: Tolerances,
) -> OgeomResult<Trimming> {
    deflection.validate()?;
    if model.kind_of(face)? != ShapeType::Face {
        ogeom_bail!(Construction, "expected a face");
    }
    let Some(node) = model.node(face) else {
        ogeom_bail!(Dangling, "face is not in this model");
    };
    let NodeData::Face(data) = node.data() else {
        ogeom_bail!(Construction, "face node holds no face data");
    };
    let Some(surface) = model.geometry().surface(data.surface) else {
        ogeom_bail!(Dangling, "face refers to a surface not in this model");
    };

    trimming_rings(
        model,
        face,
        data.surface,
        surface,
        deflection,
        &EdgeChords::new(),
        in_chart,
        tol,
    )
}

/// The chord each edge is to be drawn to, by edge node index, where a face
/// wants its edges finer than the caller's chord: what every face sharing
/// an edge has to agree on, or their meshes disagree along it.
///
/// Keyed by the edge's node, so both faces bounding it look the same value
/// up and sample it identically. Absent means the caller's own chord.
#[allow(
    clippy::disallowed_types,
    reason = "public type; a std map keeps the API"
)]
pub type EdgeChords = std::collections::HashMap<u32, f64>;

/// How many times a face's boundary may be redrawn finer before its
/// crossing is taken to be something the chord cannot fix.
///
/// Six halvings is a chord sixty-four times tighter than asked. A sliver
/// still crossing itself there is degenerate in a way refinement does not
/// reach (two boundaries genuinely on top of one another), and drawing it
/// a seventh time only spends longer to say so.
const REFINEMENTS: usize = 6;

/// [`face_chords`], through `kept` where there is one: a face is known by
/// its mesh's key at the chords held for its edges, and its answer is kept
/// by where each edge comes in the face's walk, so the same face built
/// again in another model, its edges named afresh, finds it.
#[allow(
    clippy::too_many_arguments,
    reason = "face_chords' own arguments and the cache"
)]
fn face_chords_kept(
    model: &Model,
    face: &Shape,
    id: ogeom_topo::SurfaceId,
    surface: &SurfaceGeometry,
    deflection: Deflection,
    held: &EdgeChords,
    kept: Option<&FaceMeshCache>,
    tol: Tolerances,
) -> OgeomResult<EdgeChords> {
    let walked = || -> Option<Vec<u32>> {
        let mut edges = Vec::new();
        for wire in model.ordered_children_of(face).ok()? {
            for edge in model.ordered_children_of(&wire).ok()? {
                edges.push(edge.node().index());
            }
        }
        Some(edges)
    };
    let (Some(kept), Some(key), Some(walk)) = (
        kept,
        mesh_key(model, face, deflection, Some(held), Pass::Agreeing, tol),
        walked(),
    ) else {
        return face_chords(model, face, id, surface, deflection, held, tol);
    };
    if let Some(found) = kept.chords(&key)
        && found.len() == walk.len()
    {
        return Ok(walk
            .iter()
            .zip(found)
            .filter_map(|(&edge, chord)| chord.map(|c| (edge, c)))
            .collect());
    }
    let found = face_chords(model, face, id, surface, deflection, held, tol)?;
    // Kept only where the walk names every edge answered for.
    if found.keys().all(|e| walk.contains(e)) {
        kept.keep_chords(key, walk.iter().map(|e| found.get(e).copied()).collect());
    }
    Ok(found)
}

/// What one face needs its own edges drawn with.
///
/// A face narrower than the chord error its boundary is drawn with has a
/// boundary that crosses *itself*: on an arc sliver a few microns wide,
/// drawn at a chord whose sagitta exceeds that width, the inner polyline
/// bulges straight through the outer one, and what reaches the
/// triangulator is not a region at all but disconnected fragments with
/// holes between them, which stop the body meshing closed.
///
/// The deflection a caller asks for bounds how far the mesh may sit from
/// the surface; it is not a licence to hand the triangulator a polygon that
/// crosses itself. Refining *lowers* that error, so this never breaks the
/// caller's bound.
///
/// Only this face's own edges are named, and only when they need it.
fn face_chords(
    model: &Model,
    face: &Shape,
    id: ogeom_topo::SurfaceId,
    surface: &SurfaceGeometry,
    deflection: Deflection,
    held: &EdgeChords,
    tol: Tolerances,
) -> OgeomResult<EdgeChords> {
    let deflection = in_own_frame(model, face, deflection);
    let edges: Vec<u32> =
        ogeom_topo::explore(model, face, ogeom_topo::Filter::OfType(ShapeType::Edge))?
            .iter()
            .map(|e| e.node().index())
            .collect();
    // From the chords its edges are already held to: a face whose
    // neighbour drew one of its edges finer is asked about its rings as
    // they are drawn, one rim fine and the rest as they were, which is
    // where a narrow strip's coarse rim sags over its fine one. A chord
    // only ever comes down.
    let mut finer: EdgeChords = edges
        .iter()
        .filter_map(|e| held.get(e).map(|c| (*e, *c)))
        .collect();
    let lower = |finer: &mut EdgeChords, chord: f64| {
        for edge in &edges {
            let at = finer.entry(*edge).or_insert(chord);
            *at = at.min(chord);
        }
    };
    let mut chord = deflection.chord;
    // A face narrower than a few chords first: its edges drawn to a
    // fraction of its width, so the boundary's sag is small against it.
    let first = trimming_rings(model, face, id, surface, deflection, &finer, false, tol)?.rings;
    if let Some(want) = narrow_chord(surface, &first, deflection, tol) {
        chord = want;
        lower(&mut finer, chord);
    }
    for _ in 0..=REFINEMENTS {
        let rings = trimming_rings(model, face, id, surface, deflection, &finer, false, tol)?.rings;
        // Only the verdict is read: where the rings crossed, the mesh
        // inside them is not drawn.
        let planar = region_mesh(
            &rings,
            surface,
            deflection,
            rim_holds(model, face, deflection, tol)?,
            Upto::Crossing,
            tol,
        )?;
        if !planar.crossed && planar.triangles.len() >= planar.boundary_triangles {
            break;
        }
        chord *= 0.5;
        lower(&mut finer, chord);
    }
    Ok(finer)
}

/// How many times [`split_flat_triangles`] goes round a face: one per flat
/// triangle in the longest row of them along one line.
const FLAT_PASSES: usize = 16;

/// Take apart every triangle whose corners lie on one line in space.
///
/// A chart triangle can be sound and its lift have no area: on a cone, the
/// apex and two points on one ruling, the apex taken at another `u`. Its
/// edges pair with its neighbours', so a weld finds the mesh closed, but a
/// reader that drops triangles with no area (as a 3MF reader may) finds a
/// hole where it was. The corner on the long side splits the neighbour
/// across that side in two instead, and the flat triangle goes: the pair's
/// outline, and so every edge the face shares, is what it was.
fn split_flat_triangles(mesh: &mut Triangulation, tol: Tolerances) {
    let mut across: ogeom_core::FastMap<(u32, u32), usize> = ogeom_core::FastMap::default();
    for (i, t) in mesh.triangles.iter().enumerate() {
        for k in 0..3 {
            across.insert((t[k], t[(k + 1) % 3]), i);
        }
    }
    let at = |mesh: &Triangulation, i: u32| mesh.positions[i as usize];
    // Whether `m` lies on the segment from `p` to `q`, strictly between them.
    let between = |p: Point, q: Point, m: Point| {
        let side = q - p;
        let length = side.magnitude();
        if length <= tol.confusion() {
            return false;
        }
        let along = (m - p).dot(side) / (length * length);
        let off = (m - p).cross(side).magnitude() / length;
        off <= tol.confusion()
            && along * length > tol.confusion()
            && (1.0 - along) * length > tol.confusion()
    };
    let flat = |mesh: &Triangulation, t: [u32; 3]| {
        let [a, b, c] = t.map(|i| at(mesh, i));
        (b - a).cross(c - a).magnitude()
            <= tol.confusion() * a.distance(b).max(b.distance(c)).max(c.distance(a))
    };
    let mut gone = vec![false; mesh.triangles.len()];
    // A ruling can carry several flat triangles in a row, each the other's
    // neighbour across its long side: the one whose neighbour has area goes
    // first, and the next pass finds the one inside it a neighbour that
    // now has area too.
    for _ in 0..FLAT_PASSES {
        let mut split = false;
        for i in 0..mesh.triangles.len() {
            if gone[i] {
                continue;
            }
            let t = mesh.triangles[i];
            // The corner `m` on the side from `p` to `q`, wound `p, q, m`.
            let Some([p, q, m]) = (0..3)
                .map(|k| [t[k], t[(k + 1) % 3], t[(k + 2) % 3]])
                .find(|[p, q, m]| between(at(mesh, *p), at(mesh, *q), at(mesh, *m)))
            else {
                continue;
            };
            let Some(&n) = across.get(&(q, p)) else {
                continue;
            };
            if n == i || gone[n] {
                continue;
            }
            let other = mesh.triangles[n];
            let Some(&d) = other.iter().find(|&&v| v != p && v != q) else {
                continue;
            };
            let halves = [[q, m, d], [m, p, d]];
            if halves.iter().any(|h| flat(mesh, *h)) {
                continue;
            }
            gone[i] = true;
            mesh.triangles[n] = halves[0];
            for k in 0..3 {
                across.insert((halves[0][k], halves[0][(k + 1) % 3]), n);
            }
            across.remove(&(q, p));
            across.remove(&(p, q));
            let added = mesh.triangles.len();
            mesh.triangles.push(halves[1]);
            gone.push(false);
            for k in 0..3 {
                across.insert((halves[1][k], halves[1][(k + 1) % 3]), added);
            }
            split = true;
        }
        if !split {
            break;
        }
    }
    let mut kept = gone.iter();
    mesh.triangles
        .retain(|_| !kept.next().copied().unwrap_or(false));
}

/// How many chords wide a face must be for its edges to be drawn at the
/// caller's chord; narrower, they are drawn at the width over this.
///
/// The boundary of a face drawn at chord `c` sags up to `c` between its
/// points, and a triangle from that boundary across the face to the
/// other side stands off the surface by the sag. On a face `w` wide that
/// is `c / w` of the way to standing on end; held to a quarter, the
/// triangles lean at most fourteen degrees.
const NARROW: f64 = 4.0;

/// The chord a face narrower than [`NARROW`] chords wants its edges drawn
/// to, `None` for a face wide enough at the caller's.
///
/// Width is the rings' chart extent each way scaled by the mean tangent
/// length that way, the smaller of the two: right for a strip along a
/// chart axis, an over-estimate for one across the chart, which then
/// keeps the caller's chord.
fn narrow_chord(
    surface: &SurfaceGeometry,
    rings: &[Vec<Point2>],
    deflection: Deflection,
    tol: Tolerances,
) -> Option<f64> {
    let (lo, hi) = chart_extent(rings);
    if !(lo.x.is_finite() && hi.x.is_finite() && lo.y.is_finite() && hi.y.is_finite()) {
        return None;
    }
    let (mut du, mut dv, mut n) = (0.0_f64, 0.0_f64, 0usize);
    for i in 0..=2 {
        for j in 0..=2 {
            let u = lo.x + (hi.x - lo.x) * (0.25 + 0.25 * f64::from(i));
            let v = lo.y + (hi.y - lo.y) * (0.25 + 0.25 * f64::from(j));
            if let Ok((a, b)) = surface.d1_at(u, v, tol) {
                du += a.magnitude();
                dv += b.magnitude();
                n += 1;
            }
        }
    }
    if n == 0 {
        return None;
    }
    #[allow(clippy::cast_precision_loss)]
    let width = ((hi.x - lo.x) * du / n as f64).min((hi.y - lo.y) * dv / n as f64);
    if !width.is_finite() || width <= tol.confusion() || width >= deflection.chord * NARROW {
        return None;
    }
    Some(width / NARROW)
}

/// Whether every edge of `face` is already held to `want` or finer.
fn edges_already_at(model: &Model, face: &Shape, finer: &EdgeChords, want: f64) -> bool {
    ogeom_topo::explore(model, face, ogeom_topo::Filter::OfType(ShapeType::Edge)).is_ok_and(
        |edges| {
            edges
                .iter()
                .all(|e| finer.get(&e.node().index()).is_some_and(|&c| c <= want))
        },
    )
}

/// A face's trimming rings, walked, folded and cleaned.
struct Trimming {
    /// The outer ring first, then the holes, each closed without a repeated
    /// closing point.
    rings: Vec<Vec<Point2>>,
    /// Each ring point's position in space where its edge's own curve put
    /// it, `None` where it was made up.
    anchors: Vec<Vec<Option<Point>>>,
    /// Whether every edge's polyline honoured the deflection.
    met: bool,
    /// Each ring as its edges first walked it, chart and space together:
    /// what a face whose rings all collapse is drawn from.
    walked: Vec<Vec<(Point2, Point)>>,
}

/// One walked ring: chart points, anchors, deflection honesty, the ambiguous
/// whole-period folds taken, and the half-period ties left undecided.
type WalkedRing = (
    Vec<Point2>,
    Vec<Option<Point>>,
    bool,
    Vec<(usize, f64)>,
    Vec<usize>,
);

/// The whole periods a ring spanning `extent` along a closed direction is
/// slid by to stand in `domain`: up where it starts below, down where it
/// ends above, and none where it is already inside.
///
/// A hair past an end is fit noise, not a period: the knots take it, and
/// rounding it up to a whole turn would carry the ring a period the wrong
/// way. The same hair is taken off before counting, so a ring a whole
/// period out by rounding comes back one period, not two.
fn turns_into_domain(extent: (f64, f64), domain: (f64, f64)) -> f64 {
    let ((least, most), (lo, hi)) = (extent, domain);
    let span = hi - lo;
    let slack = span * 1e-6;
    if least < lo - slack {
        ((lo - least - slack) / span).ceil()
    } else if most > hi + slack {
        -((most - hi - slack) / span).ceil()
    } else {
        0.0
    }
}

/// A face's rings walked from its wires, each drawn by [`boundary_ring`]
/// with `in_chart`.
#[allow(
    clippy::too_many_arguments,
    reason = "the walk's settings, each read once"
)]
fn trimming_rings(
    model: &Model,
    face: &Shape,
    id: ogeom_topo::SurfaceId,
    surface: &SurfaceGeometry,
    deflection: Deflection,
    finer: &EdgeChords,
    in_chart: bool,
    tol: Tolerances,
) -> OgeomResult<Trimming> {
    let mut rings = Vec::new();
    let mut ring_anchors = Vec::new();
    let mut ring_folds: Vec<Vec<(usize, f64)>> = Vec::new();
    let mut ring_ties: Vec<Vec<usize>> = Vec::new();
    let mut ring_wires: Vec<Shape> = Vec::new();
    let mut met = true;
    let wires = model.ordered_children_of(face)?;
    let bounded = !wires.is_empty();
    let mut walked = Vec::new();
    for wire in wires {
        let (ring, anchors, ring_met, folds, ties) =
            boundary_ring(model, &wire, id, deflection, finer, in_chart, tol)?;
        met &= ring_met;
        if ring.len() >= 3 {
            let mut here = Vec::with_capacity(ring.len());
            for (uv, anchor) in ring.iter().zip(&anchors) {
                let at = match anchor {
                    Some(p) => *p,
                    None => surface.point_at(uv.x, uv.y, tol)?,
                };
                here.push((*uv, at));
            }
            walked.push(here);
            rings.push(ring);
            ring_anchors.push(anchors);
            ring_folds.push(folds);
            ring_ties.push(ties);
            ring_wires.push(wire);
        }
    }

    // A band between two *wound* rings: each boundary winds a periodic
    // direction once (a bore wall whose ends are a full rim and a
    // staircase of arcs, a torus band between two parallels, a ball's
    // belt between two latitudes), and neither ring closes on its own. In
    // the unrolled chart the face is the strip between the two chains, so
    // the pair is merged into one ring by two joining runs standing
    // exactly one period apart: their lifted points are the same 3D
    // points, and the weld closes them the way it closes a seam. Paired
    // before either is closed on its own: a rim closed alone spans the
    // whole surface, and two of those cancel where they overlap.
    {
        use ogeom_geom::Surface as _;
        let ((ua, ub), (va, vb)) = surface.domain();
        let axes = [
            (surface.is_periodic_u(), ub - ua),
            (surface.is_periodic_v(), vb - va),
        ];
        for (axis, (periodic, period)) in axes.into_iter().enumerate() {
            if !periodic || period <= 0.0 {
                continue;
            }
            let along = |p: Point2| if axis == 0 { p.x } else { p.y };
            let across = |p: Point2| if axis == 0 { p.y } else { p.x };
            let make = |a: f64, c: f64| {
                if axis == 0 {
                    Point2::new(a, c)
                } else {
                    Point2::new(c, a)
                }
            };
            let d_of = |ring: &[Point2]| -> f64 {
                ring.last().map_or(0.0, |l| along(*l) - along(ring[0]))
            };
            let open_wound: Vec<usize> = rings
                .iter()
                .enumerate()
                .filter(|(_, ring)| {
                    ring.len() >= 3
                        && (d_of(ring).abs() - period).abs() <= period * 1e-3
                        && ring
                            .last()
                            .is_some_and(|l| ring[0].distance(*l) > period * 0.5)
                })
                .map(|(i, _)| i)
                .collect();
            let [i, j] = open_wound[..] else {
                continue;
            };
            if d_of(&rings[i]).signum() == d_of(&rings[j]).signum() {
                continue;
            }
            let sign = d_of(&rings[i]).signum();
            // The column the joining runs stand on. Wherever the first rim's
            // chain happens to end is as good as anywhere for the rims
            // themselves, and no good at all when a third ring lies there:
            // a cross hole through a bore wall, whose loop the runs would
            // then cut through, two constraints refused and the face never
            // drawn whole. So the column is chosen where no other ring is,
            // in the widest gap the others leave round the period, and
            // both rims are re-cut to begin there: the same chains, the
            // same lifted points, begun a turn's fraction round.
            let column = {
                let mut spans: Vec<(f64, f64)> = Vec::new();
                for (k, ring) in rings.iter().enumerate() {
                    if k == i || k == j {
                        continue;
                    }
                    let (lo, hi) = ring
                        .iter()
                        .fold((f64::INFINITY, f64::NEG_INFINITY), |(lo, hi), p| {
                            (lo.min(along(*p)), hi.max(along(*p)))
                        });
                    if lo.is_finite() && hi - lo < period {
                        spans.push((lo.rem_euclid(period), hi - lo));
                    }
                }
                let fallback = rings[i].last().map_or(0.0, |l| along(*l));
                if spans.is_empty() {
                    fallback
                } else {
                    // The widest gap between the spans' images on one turn.
                    let mut edges: Vec<(f64, f64)> =
                        spans.iter().map(|(lo, len)| (*lo, lo + len)).collect();
                    edges.sort_by(|a, b| a.0.total_cmp(&b.0));
                    let mut best = (f64::NEG_INFINITY, fallback);
                    for k in 0..edges.len() {
                        let end = edges[k].1;
                        let next = if k + 1 < edges.len() {
                            edges[k + 1].0
                        } else {
                            edges[0].0 + period
                        };
                        if next - end > best.0 {
                            best = (next - end, f64::midpoint(end, next));
                        }
                    }
                    best.1
                }
            };
            let recut = |chain: &mut Vec<Point2>, anchors: &mut Vec<Option<Point>>| {
                // The chain runs one period from its first point, either way
                // round; it is cut where it first reaches the column in its
                // own direction, and its head carried a period on to the
                // tail, closing on the same lifted point.
                let n = chain.len();
                if n < 3 {
                    return;
                }
                let forward = d_of(chain) > 0.0;
                let start = along(chain[0]);
                let target = if forward {
                    start + (column - start).rem_euclid(period)
                } else {
                    start - (start - column).rem_euclid(period)
                };
                let reached = |k: usize| {
                    if forward {
                        along(chain[k]) >= target
                    } else {
                        along(chain[k]) <= target
                    }
                };
                let Some(k) = (1..n - 1).find(|&k| reached(k)) else {
                    return;
                };
                if k <= 1 {
                    return;
                }
                let carry = if forward { period } else { -period };
                let mut turned: Vec<Point2> = chain[k..].to_vec();
                let mut turned_anchors: Vec<Option<Point>> = anchors[k..].to_vec();
                for (p, a) in chain[1..=k].iter().zip(&anchors[1..=k]) {
                    turned.push(make(along(*p) + carry, across(*p)));
                    turned_anchors.push(*a);
                }
                *chain = turned;
                *anchors = turned_anchors;
            };
            recut(&mut rings[i], &mut ring_anchors[i]);
            recut(&mut rings[j], &mut ring_anchors[j]);
            let b = rings.remove(j);
            let b_anchors = ring_anchors.remove(j);
            ring_folds.remove(j);
            ring_ties.remove(j);
            ring_wires.remove(j);
            let a = &mut rings[i];
            let a_anchors = &mut ring_anchors[i];
            let a_last = a.last().copied().unwrap_or(a[0]);
            // Whole periods only: the second chain slides along the unrolled
            // chart until its start stands nearest the first chain's end.
            let shift = ((along(a_last) - along(b[0])) / period).round() * period;
            let b_first = make(along(b[0]) + shift, across(b[0]));
            // The run is no edge of the file's, so nothing else draws it to
            // the chord: it is drawn here, as the grid's lines are.
            let fractions = joining_run(surface, a_last, b_first, deflection, tol);
            let inner = &fractions[1..fractions.len() - 1];
            for &f in inner {
                a.push(make(
                    along(a_last) + (along(b_first) - along(a_last)) * f,
                    across(a_last) + (across(b_first) - across(a_last)) * f,
                ));
                a_anchors.push(None);
            }
            for (p, anchor) in b.iter().zip(&b_anchors) {
                a.push(make(along(*p) + shift, across(*p)));
                a_anchors.push(*anchor);
            }
            // The way back: the same run, one period over, walked the other
            // way; the two runs lift to identical points.
            for &f in inner.iter().rev() {
                a.push(make(
                    along(a_last) + (along(b_first) - along(a_last)) * f - sign * period,
                    across(a_last) + (across(b_first) - across(a_last)) * f,
                ));
                a_anchors.push(None);
            }
            // The band now stands on one stretch of the chart, a period wide
            // from where the first rim was cut; every other ring (a hole
            // through the wall) is slid by whole periods into that stretch,
            // the same lifted points, or it stands outside the band it
            // belongs in.
            let (band_lo, band_hi) = a
                .iter()
                .fold((f64::INFINITY, f64::NEG_INFINITY), |(lo, hi), p| {
                    (lo.min(along(*p)), hi.max(along(*p)))
                });
            for (k, ring) in rings.iter_mut().enumerate() {
                if k == i || ring.is_empty() {
                    continue;
                }
                #[allow(clippy::cast_precision_loss)]
                let mean = ring.iter().map(|p| along(*p)).sum::<f64>() / ring.len() as f64;
                if mean >= band_lo && mean <= band_hi {
                    continue;
                }
                let turns = ((mean - band_lo) / period).floor();
                if turns == 0.0 {
                    continue;
                }
                for p in ring.iter_mut() {
                    *p = make(along(*p) - turns * period, across(*p));
                }
            }
        }
    }

    // A single ring still winding after the pairing had no partner, and if
    // it cannot close it was mis-folded: a jump of exactly one period is the
    // same 3D point, and the greedy fold that zeroed it may have wound the
    // ring instead. Undoing one such ambiguous fold (translating the walk
    // from there on by a period) closes the ring; the join it reopens
    // still lifts to a single point.
    if geometry_winds(surface) {
        use ogeom_geom::Surface as _;
        let ((ua, ub), _) = surface.domain();
        let period = ub - ua;
        for ((ring, folds), ties) in rings.iter_mut().zip(&ring_folds).zip(&ring_ties) {
            let Some(last) = ring.last().copied() else {
                continue;
            };
            let du = last.x - ring[0].x;
            let k = (du / period).round();
            if k == 0.0
                || (du - k * period).abs() > period * 1e-3
                || ring[0].distance(last) <= period * 1e-3
            {
                continue;
            }
            // The last ambiguous fold whose applied shift matches the
            // winding is the one to undo. Where none does, the walk wound at
            // a half-period tie it left in place (a ruling into a cone's
            // apex, whose other half the ring needed): the last one is
            // where to translate from.
            let fold = folds
                .iter()
                .rev()
                .find(|(_, shift)| (shift - k * period).abs() <= period * 1e-6)
                .map(|&(fold_start, _)| fold_start)
                .filter(|&fold_start| fold_start < ring.len());
            if fold.is_none()
                && let Some(&tie) = ties.last()
                && tie < ring.len()
            {
                for p in &mut ring[tie..] {
                    p.x -= k * period;
                }
                continue;
            }
            if let Some(fold_start) = fold {
                // Translating from the last undecided tie before the fold
                // (where the walk first guessed) keeps the period jump on
                // the degenerate row, where it lifts to nothing.
                let start = ties
                    .iter()
                    .rev()
                    .find(|&&t| t < fold_start)
                    .copied()
                    .unwrap_or(fold_start);
                for p in &mut ring[start..] {
                    p.x -= k * period;
                }
            }
        }
    }

    // A ring still winding after all that closes on its own: against a
    // pole row, or its own translate a period over the other way.
    for (ring, anchors) in rings.iter_mut().zip(ring_anchors.iter_mut()) {
        close_wound_ring(ring, anchors, surface, tol);
        snap_chart_junctions(ring, anchors, surface, tol);
    }

    // Points closer than the triangulator's own resolution make it refuse
    // the constraint; a real file's chart can carry them. The consecutive
    // near-duplicates collapse, anchors staying aligned.
    for (ring, anchors) in rings.iter_mut().zip(ring_anchors.iter_mut()) {
        collapse_close_points(ring, anchors);
    }
    // A hole's ring is walked on its own and can come to rest a whole
    // period from the ring that bounds the face, its pcurves named in
    // another turn of the chart: a cross bore's loop at u = 2pi beside a
    // wall drawn from -pi to pi. Left there it lies outside the face's
    // region, which is then drawn whole and the hole drawn again as a disc
    // of its own. So each ring is slid by whole periods, the same lifted
    // points, to stand where the roomiest ring (the face's own bound) does.
    {
        use ogeom_geom::Surface as _;
        let ((ua, ub), (va, vb)) = surface.domain();
        let axes = [
            (surface.is_periodic_u(), ub - ua),
            (surface.is_periodic_v(), vb - va),
        ];
        let area = |ring: &[Point2]| -> f64 {
            let Some(&o) = ring.first() else {
                return 0.0;
            };
            ring.iter()
                .zip(ring.iter().cycle().skip(1))
                .map(|(a, b)| (*a - o).cross(*b - o))
                .sum::<f64>()
                .abs()
        };
        let bound = (0..rings.len()).max_by(|&a, &b| area(&rings[a]).total_cmp(&area(&rings[b])));
        for (axis, (periodic, period)) in axes.into_iter().enumerate() {
            let Some(bound) = bound else {
                break;
            };
            if !periodic || period <= 0.0 {
                continue;
            }
            let along = |p: &Point2| if axis == 0 { p.x } else { p.y };
            let (lo, hi) = rings[bound]
                .iter()
                .fold((f64::INFINITY, f64::NEG_INFINITY), |(lo, hi), p| {
                    (lo.min(along(p)), hi.max(along(p)))
                });
            if !(lo.is_finite() && hi.is_finite()) {
                continue;
            }
            for (k, ring) in rings.iter_mut().enumerate() {
                if k == bound || ring.is_empty() {
                    continue;
                }
                #[allow(clippy::cast_precision_loss)]
                let mean = ring.iter().map(along).sum::<f64>() / ring.len() as f64;
                if mean >= lo && mean <= hi {
                    continue;
                }
                let turns = ((mean - f64::midpoint(lo, hi)) / period).round();
                let moved = mean - turns * period;
                if turns == 0.0 || moved < lo || moved > hi {
                    continue;
                }
                for p in ring.iter_mut() {
                    if axis == 0 {
                        p.x -= turns * period;
                    } else {
                        p.y -= turns * period;
                    }
                }
            }
        }
    }
    // Folded across a join, a ring on a *closed* surface can come to rest a
    // whole period outside the domain. A periodic surface would not mind
    // (it wraps), but a surface that merely closes on itself evaluates only
    // where its knots are and refuses everywhere else, so a ring every
    // point of which lies a turn past the end does not lift at all. The
    // fold keeps the ring continuous, which is the part that matters, and
    // a rigid slide by whole periods keeps it so: the same points on the
    // surface, named inside the chart.
    {
        use ogeom_geom::Surface as _;
        let ((ua, ub), (va, vb)) = surface.domain();
        let slides = [
            (!surface.is_periodic_u() && surface.is_closed_u(tol), ua, ub),
            (!surface.is_periodic_v() && surface.is_closed_v(tol), va, vb),
        ];
        for (across, (closed, lo, hi)) in [true, false].into_iter().zip(slides) {
            if !closed || hi <= lo {
                continue;
            }
            let span = hi - lo;
            for ring in &mut rings {
                let read = |p: &Point2| if across { p.x } else { p.y };
                let (least, most) = ring
                    .iter()
                    .fold((f64::INFINITY, f64::NEG_INFINITY), |(a, b), p| {
                        (a.min(read(p)), b.max(read(p)))
                    });
                if !(least.is_finite() && most.is_finite()) || most - least > span * (1.0 + 1e-9) {
                    continue;
                }
                let turns = turns_into_domain((least, most), (lo, hi));
                if turns == 0.0 {
                    continue;
                }
                for p in ring.iter_mut() {
                    if across {
                        p.x += turns * span;
                    } else {
                        p.y += turns * span;
                    }
                }
            }
        }
    }
    if *MESH_DEBUG {
        for (i, (ring, anchors)) in rings.iter().zip(&ring_anchors).enumerate() {
            eprintln!("DBG ring {i}: {} points", ring.len());
            for (p, a) in ring.iter().zip(anchors) {
                eprintln!(
                    "DBG   uv({:.5},{:.5}) anchor {}",
                    p.x,
                    p.y,
                    a.map_or("-".to_string(), |q| format!(
                        "({:.4},{:.4},{:.4})",
                        q.x, q.y, q.z
                    ))
                );
            }
        }
    }
    // A run out along an edge and straight back along it (a ring that
    // reads `p, q, p`) is a spike into the region that bounds nothing:
    // the file's way of drawing a slit of no width at all on the face's
    // own boundary. It triangulates to two hairs and one vertex too many,
    // which is one triangle more than a boundary that encloses a region
    // has, and it is not a crossing. Off it comes, out to in, until the
    // ring reverses nowhere.
    // And two consecutive points a millionth of the ring's size apart are
    // one point: the end of the last edge and the start of the first, each
    // where its own curve put the shared vertex, a file's slop apart. Kept
    // both, the second sits on the first's next segment to the last bit,
    // and a constraint through a vertex is one the triangulation refuses.
    for (ring, anchors) in rings.iter_mut().zip(ring_anchors.iter_mut()) {
        let reach = chart_reach(ring);
        merge_near_duplicates(ring, anchors, reach);
        remove_spikes(ring, anchors, reach);
    }
    let long_enough: Vec<bool> = rings.iter().map(|r| r.len() >= 3).collect();
    let mut it = long_enough.iter();
    ring_wires.retain(|_| *it.next().unwrap_or(&true));
    rings.retain(|r| r.len() >= 3);
    ring_anchors.retain(|a| a.len() >= 3);

    // An inner ring that walks out and back again along edges through the
    // same vertices, thinner than a micron, is a slit, not a hole: a loop
    // run out along two arcs and back along two splines fitted to the same
    // arcs encloses nothing, and the triangulator can only read it as a
    // tangle, drawing holes the face does not have (see [`is_slit`]). The
    // outer ring is left alone, whatever it is, since a face that is
    // itself a slit is a different question.
    if rings.len() > 1 {
        // About each ring's first point: the chart's origin may stand far
        // off, where the sum is a difference of huge products.
        let chart_area = |ring: &[Point2]| -> f64 {
            let Some(&anchor) = ring.first() else {
                return 0.0;
            };
            let mut a = 0.0;
            for i in 0..ring.len() {
                a += (ring[i] - anchor).cross(ring[(i + 1) % ring.len()] - anchor);
            }
            a.abs()
        };
        let outer = (0..rings.len())
            .max_by(|&i, &j| chart_area(&rings[i]).total_cmp(&chart_area(&rings[j])))
            .unwrap_or(0);
        let mut keep = Vec::with_capacity(rings.len());
        for i in 0..rings.len() {
            keep.push(i == outer || !is_slit(model, &ring_wires[i], &ring_anchors[i], tol)?);
        }
        let mut it = keep.iter();
        rings.retain(|_| *it.next().unwrap_or(&true));
        let mut it = keep.iter();
        ring_anchors.retain(|_| *it.next().unwrap_or(&true));
    }

    // A face with no wires covers its surface's whole domain, so the domain
    // rectangle is the boundary. A face whose wires all collapse (a sliver
    // narrower than a millionth of its own length, its two long sides one
    // line in the chart) encloses nothing, and keeps no ring at all: taken
    // for a face without wires, it would be drawn as the whole plane.
    if rings.is_empty() && !bounded {
        let ring = domain_ring(surface, deflection, tol);
        ring_anchors.push(vec![None; ring.len()]);
        rings.push(ring);
    }
    Ok(Trimming {
        rings,
        anchors: ring_anchors,
        met,
        walked,
    })
}

/// Collapse a ring's consecutive near-duplicate points, its closing one
/// with them, anchors staying aligned. A ring of fewer than three points is
/// left as it is.
fn collapse_close_points(ring: &mut Vec<Point2>, anchors: &mut Vec<Option<Point>>) {
    if ring.len() < 3 {
        return;
    }
    let mut extent = 0.0_f64;
    for pair in ring.windows(2) {
        extent = extent.max(pair[0].distance(pair[1]));
    }
    let eps = extent.mul_add(1e-9, 1e-12);
    let mut kept_ring = Vec::with_capacity(ring.len());
    let mut kept_anchors = Vec::with_capacity(anchors.len());
    for (p, a) in ring.iter().zip(anchors.iter()) {
        if kept_ring
            .last()
            .is_some_and(|held: &Point2| held.distance(*p) <= eps)
        {
            continue;
        }
        kept_ring.push(*p);
        kept_anchors.push(*a);
    }
    if kept_ring.len() > 2
        && let (Some(first), Some(last)) = (kept_ring.first(), kept_ring.last())
        && first.distance(*last) <= eps
    {
        kept_ring.pop();
        kept_anchors.pop();
    }
    *ring = kept_ring;
    *anchors = kept_anchors;
}

/// Whether an inner ring is a slit rather than a hole: its wire walks back
/// the way it came, and it is thinner than a micron.
///
/// A slit's wire runs out along a chain of edges and home along a second
/// chain through the same vertices in the opposite order, every edge from
/// one vertex to the next answered by another from the next back to it,
/// so the two chains bound nothing between them but the rounding of two
/// fits to one curve. A hole's wire goes round: a ring of facets, or a
/// rim of edges, meets each vertex once, and is a hole however thin. A
/// wire of two edges between the same two vertices walks back by its
/// vertices whether it is a slit or a round hole cut in two halves; it is
/// a slit only where its two sides lie along each other, the ring's width
/// under a thousandth of its perimeter. Two halves of a hole bound a round
/// region, as wide as a sixth of its perimeter for a circle. Width is
/// twice the enclosed area over the perimeter, measured in space through
/// the anchors, so a chart's units do not enter into it. A ring with a
/// point not anchored is not called a slit.
///
/// # Errors
///
/// [`OgeomError::Dangling`](ogeom_core::OgeomError::Dangling) where an edge
/// of the wire is not in the model.
fn is_slit(
    model: &Model,
    wire: &Shape,
    anchors: &[Option<Point>],
    tol: Tolerances,
) -> OgeomResult<bool> {
    let measure = || -> Option<(f64, f64)> {
        let pts: Option<Vec<Point>> = anchors.iter().copied().collect();
        let pts = pts?;
        // About the ring's first point: the sum is the same about any
        // point, and about a far one it is a difference of huge products.
        let anchor = *pts.first()?;
        let mut normal = Vector::ZERO;
        let mut perimeter = 0.0;
        for i in 0..pts.len() {
            let (a, b) = (pts[i], pts[(i + 1) % pts.len()]);
            normal += (a - anchor).cross(b - anchor);
            perimeter += a.distance(b);
        }
        (perimeter > 0.0).then(|| (normal.magnitude() / perimeter, perimeter))
    };
    let Some((width, perimeter)) = measure() else {
        return Ok(false);
    };
    if width >= tol.confusion() * 1e4 || !walks_back(model, wire)? {
        return Ok(false);
    }
    Ok(model.ordered_children_of(wire)?.len() > 2 || width < perimeter * 1e-3)
}

/// Whether a wire's edges pair off, each running from one vertex to
/// another answered by a second running from that one back: the walk out
/// and the walk home of a ring that goes nowhere. An edge without
/// vertices, or one starting where it ends, has no answer.
fn walks_back(model: &Model, wire: &Shape) -> OgeomResult<bool> {
    let mut ends: Vec<(Shape, Shape)> = Vec::new();
    for edge in model.ordered_children_of(wire)? {
        let bounds = model.children_of(&edge)?;
        let (Some(first), Some(last)) = (bounds.first(), bounds.last()) else {
            return Ok(false);
        };
        if first.is_same(last) {
            return Ok(false);
        }
        ends.push(if edge.orientation() == Orientation::Reversed {
            (last.clone(), first.clone())
        } else {
            (first.clone(), last.clone())
        });
    }
    if ends.len() < 2 || !ends.len().is_multiple_of(2) {
        return Ok(false);
    }
    let mut answered = vec![false; ends.len()];
    for i in 0..ends.len() {
        if answered[i] {
            continue;
        }
        let (from, to) = &ends[i];
        let Some(j) = (i + 1..ends.len())
            .find(|&j| !answered[j] && ends[j].0.is_same(to) && ends[j].1.is_same(from))
        else {
            return Ok(false);
        };
        answered[i] = true;
        answered[j] = true;
    }
    Ok(true)
}

/// One wire's trimming ring on a face whose chart neither repeats nor
/// closes on itself, as [`face_boundary_in_chart`] draws it among the
/// face's other rings; `None` where the wire leaves no ring there.
///
/// On such a chart each ring is drawn from its own wire alone, so a caller
/// asking about a few holes of a face with hundreds draws only those. The
/// one ring drawn differently is a slit (a wire walking out and back along
/// the same vertices, thinner than a micron), which [`face_boundary`] drops
/// unless it is the face's largest ring:
/// `largest` says whether this wire's ring is that one.
///
/// # Errors
///
/// As [`face_boundary`], and
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) where
/// the face's surface repeats or closes in either chart direction.
pub fn open_chart_ring(
    model: &Model,
    face: &Shape,
    wire: &Shape,
    largest: bool,
    deflection: Deflection,
    tol: Tolerances,
) -> OgeomResult<Option<Vec<Point2>>> {
    use ogeom_geom::Surface as _;
    deflection.validate()?;
    let Some(NodeData::Face(data)) = model.node(face).map(|n| n.data()) else {
        ogeom_bail!(Construction, "expected a face");
    };
    let Some(surface) = model.geometry().surface(data.surface) else {
        ogeom_bail!(Dangling, "face refers to a surface not in this model");
    };
    if surface.is_periodic_u()
        || surface.is_periodic_v()
        || surface.is_closed_u(tol)
        || surface.is_closed_v(tol)
    {
        ogeom_bail!(
            Construction,
            "a chart that repeats or closes draws its rings together"
        );
    }
    let (mut ring, mut anchors, ..) = boundary_ring(
        model,
        wire,
        data.surface,
        deflection,
        &EdgeChords::new(),
        true,
        tol,
    )?;
    if ring.len() < 3 {
        return Ok(None);
    }
    close_wound_ring(&mut ring, &mut anchors, surface, tol);
    collapse_close_points(&mut ring, &mut anchors);
    let reach = chart_reach(&ring);
    merge_near_duplicates(&mut ring, &mut anchors, reach);
    remove_spikes(&mut ring, &mut anchors, reach);
    if ring.len() < 3 || anchors.len() < 3 || (!largest && is_slit(model, wire, &anchors, tol)?) {
        return Ok(None);
    }
    Ok(Some(ring))
}

/// A face with no area as a fan of triangles across each of its walked
/// rings, every point where its edges put it. Its normals are the
/// surface's carried through the face's placement, turned where the face
/// is reversed.
fn sliver_fan(
    walked: &[Vec<(Point2, Point)>],
    surface: &SurfaceGeometry,
    placement: &ogeom_math::Transform,
    reversed: bool,
    tol: Tolerances,
) -> Triangulation {
    let mut mesh = Triangulation::new();
    for ring in walked {
        let base = u32::try_from(mesh.positions.len()).unwrap_or(u32::MAX);
        for (uv, at) in ring {
            mesh.positions.push(*at);
            mesh.parameters.push((uv.x, uv.y));
            let normal = surface
                .normal_at(uv.x, uv.y, tol)
                .map_or(ogeom_math::Vector::Z, |n| {
                    placement.apply_vector(n.vector())
                });
            mesh.normals.push(if reversed { -normal } else { normal });
        }
        let count = u32::try_from(ring.len()).unwrap_or(0);
        for k in 1..count.saturating_sub(1) {
            mesh.triangles.push([base, base + k, base + k + 1]);
        }
    }
    mesh
}

/// How near two chart points of a ring stand to be one point: within a
/// millionth of the ring's extent along each chart direction, well under
/// any feature and well over the slop two edges leave at the vertex they
/// share. Each direction is held to its own extent: a chart's two
/// directions need not share a unit (a drum's u is an angle and its v a
/// length), and a ring a few microns wide round a drum and hundreds of
/// millimetres along it keeps its width, where a millionth of its length
/// read in radians is wider than the whole ring.
#[derive(Debug, Clone, Copy)]
struct ChartReach {
    u: f64,
    v: f64,
}

impl ChartReach {
    /// Whether `p` and `q` are one point.
    fn joins(self, p: Point2, q: Point2) -> bool {
        (p.x - q.x).abs() <= self.u && (p.y - q.y).abs() <= self.v
    }
}

/// The [`ChartReach`] of `ring`, or `None` for a ring with no extent.
fn chart_reach(ring: &[Point2]) -> Option<ChartReach> {
    let (mut lo, mut hi) = (
        Point2::new(f64::INFINITY, f64::INFINITY),
        Point2::new(f64::NEG_INFINITY, f64::NEG_INFINITY),
    );
    for p in ring {
        lo = Point2::new(lo.x.min(p.x), lo.y.min(p.y));
        hi = Point2::new(hi.x.max(p.x), hi.y.max(p.y));
    }
    let (du, dv) = (hi.x - lo.x, hi.y - lo.y);
    (du.is_finite() && dv.is_finite() && du.max(dv) > 0.0).then_some(ChartReach {
        u: du * 1e-6,
        v: dv * 1e-6,
    })
}

/// Merge consecutive ring points within `reach` of each other,
/// cyclically; the earlier point and its anchor stay.
fn merge_near_duplicates(
    ring: &mut Vec<Point2>,
    anchors: &mut Vec<Option<Point>>,
    reach: Option<ChartReach>,
) {
    let Some(reach) = reach else {
        return;
    };
    if ring.len() < 2 || anchors.len() != ring.len() {
        return;
    }
    let mut i = 0;
    while i < ring.len() && ring.len() >= 2 {
        let next = (i + 1) % ring.len();
        if reach.joins(ring[i], ring[next]) {
            // The later point goes (the last one when the ring's end
            // repeats its start), and the earlier is looked at again, in
            // case it sits next to another near-duplicate.
            if next == 0 {
                ring.remove(i);
                anchors.remove(i);
                break;
            }
            ring.remove(next);
            anchors.remove(next);
        } else {
            i += 1;
        }
    }
}

/// Strip every `p, q, p` from a ring (a point stepped out to and straight
/// back from, the two `p` within `reach` of each other), with its anchors,
/// until none is left. Cyclic: the ring's last point is its first's
/// neighbour.
fn remove_spikes(
    ring: &mut Vec<Point2>,
    anchors: &mut Vec<Option<Point>>,
    reach: Option<ChartReach>,
) {
    let Some(reach) = reach else {
        return;
    };
    loop {
        let n = ring.len();
        if n < 3 || anchors.len() != n {
            return;
        }
        let Some(tip) = (0..n).find(|&i| reach.joins(ring[(i + n - 1) % n], ring[(i + 1) % n]))
        else {
            return;
        };
        // The tip and one of its two identical neighbours go.
        let neighbour = (tip + 1) % n;
        let (first, second) = if tip < neighbour {
            (neighbour, tip)
        } else {
            (tip, neighbour)
        };
        ring.remove(first);
        anchors.remove(first);
        ring.remove(second);
        anchors.remove(second);
    }
}

/// Split a nonincident side at a boundary vertex within chart and surface
/// precision. Both occurrences carry the same coordinates and spatial anchor,
/// making the junction explicit even under exact orientation predicates.
fn snap_chart_junctions(
    ring: &mut Vec<Point2>,
    anchors: &mut Vec<Option<Point>>,
    surface: &SurfaceGeometry,
    tol: Tolerances,
) {
    let original = ring.to_vec();
    if original.len() < 3 {
        return;
    }
    // Index both axes: a seam can have thousands of samples at one u, even
    // when u is the chart's longer axis. Query the sparser band at each point.
    let indexed: Vec<_> = original.iter().map(|p| Point2::new(p.y, p.x)).collect();
    let x_bands = RingBands::over(std::slice::from_ref(&indexed));
    let y_bands = RingBands::over(std::slice::from_ref(&original));
    if x_bands.bands.is_empty() || y_bands.bands.is_empty() {
        return;
    }
    let mut junctions = Vec::new();
    for (i, p) in original.iter().copied().enumerate() {
        let nearest = junction_candidates(&x_bands, &y_bands, p, tol.parametric())
            .iter()
            .flatten()
            .map(|&k| k as usize)
            .filter(|&k| k != i && (k + 1) % original.len() != i)
            .filter_map(|k| {
                let (a, b) = (original[k], original[(k + 1) % original.len()]);
                let (dx, dy) = (b.x - a.x, b.y - a.y);
                let l2 = dx * dx + dy * dy;
                if l2 == 0.0 {
                    return None;
                }
                let along = ((p.x - a.x) * dx + (p.y - a.y) * dy) / l2;
                if !(0.0 < along && along < 1.0) {
                    return None;
                }
                let q = Point2::new(a.x + dx * along, a.y + dy * along);
                let distance = p.distance(q);
                (distance <= tol.parametric()).then_some((k, along, q, distance))
            })
            .min_by(|a, b| a.3.total_cmp(&b.3));
        if let Some((side, along, q, _)) = nearest
            && let (Ok(at), Ok(on)) = (
                surface.point_at(p.x, p.y, tol),
                surface.point_at(q.x, q.y, tol),
            )
            && at.distance(on) <= tol.confusion()
        {
            junctions.push((side, along, i));
        }
    }
    if junctions.is_empty() {
        return;
    }
    // Rebuild in side order so many junctions require one linear copy.
    // Each shared vertex bends its constraint by at most the accepted reach.
    junctions.sort_by(|a, b| a.0.cmp(&b.0).then(a.1.total_cmp(&b.1)));
    let held = std::mem::take(anchors);
    ring.clear();
    let mut junctions = junctions.into_iter().peekable();
    for (side, (&p, &anchor)) in original.iter().zip(&held).enumerate() {
        ring.push(p);
        anchors.push(anchor);
        while let Some(&(at, _, vertex)) = junctions.peek() {
            if at != side {
                break;
            }
            ring.push(original[vertex]);
            anchors.push(held[vertex]);
            junctions.next();
        }
    }
}

/// Any side within `reach` of a point overlaps its band on both axes. Either
/// axis is complete; choose the one with fewer entries to keep dense seams
/// from turning the junction search into an all-pairs walk.
fn junction_candidates<'a>(
    x: &'a RingBands,
    y: &'a RingBands,
    p: Point2,
    reach: f64,
) -> &'a [Vec<u32>] {
    let range = |bands: &RingBands, height: f64| {
        let first = RingBands::band_of(height - reach, bands.low, bands.step, bands.bands.len());
        let last = RingBands::band_of(height + reach, bands.low, bands.step, bands.bands.len());
        first..=last
    };
    let (xs, ys) = (&x.bands[range(x, p.x)], &y.bands[range(y, p.y)]);
    if xs.iter().map(Vec::len).sum::<usize>() <= ys.iter().map(Vec::len).sum::<usize>() {
        xs
    } else {
        ys
    }
}

/// Whether a surface is periodic in `u` alone: the charts on which a wound
/// ring cannot close against a translate of itself and needs a partner.
fn geometry_winds(surface: &SurfaceGeometry) -> bool {
    use ogeom_geom::Surface as _;
    surface.is_periodic_u() && !surface.is_periodic_v()
}

/// Whether a point in parameter space lies inside the region `rings` bound.
///
/// Even-odd winding: inside the outer ring and outside every hole. The rings
/// come from wires whose direction already encodes outer from inner, but
/// counting crossings does not depend on that being right, which makes it
/// robust to a wire that was built the wrong way round.
///
/// Says nothing about a point *on* a ring: the crossing count of a boundary
/// point is whichever side rounding puts it. A caller that cares has to measure
/// its distance to the boundary and decide, which is what classification does.
///
/// Decided with exact predicates. See [`inside_boundary_with`].
#[must_use]
pub fn inside_boundary(rings: &[Vec<Point2>], p: Point2) -> bool {
    inside_boundary_with::<Exact>(rings, p)
}

/// As [`inside_boundary`], with the predicate implementation named.
///
/// This is the seam `docs/DATA_MODEL.md` §9 describes, and it is here rather
/// than anywhere else because this is where the *combinatorial* decision is.
/// Whether a point is inside a ring is not a measurement that can be a little
/// wrong: it decides whether a triangle is kept or dropped, so an error near a
/// boundary is a hole in the mesh rather than a slightly misplaced one.
///
/// The naive form divides to find where an edge crosses the sampling ray, and
/// that division cancels catastrophically for a point nearly on the edge.
/// `orient2d` answers the same question with no division at all, and
/// [`Exact`] answers it correctly however close the point is.
#[must_use]
pub fn inside_boundary_with<P: Predicates>(rings: &[Vec<Point2>], p: Point2) -> bool {
    let mut inside = false;
    for ring in rings {
        if crosses_odd_times::<P>(ring, p) {
            inside = !inside;
        }
    }
    inside
}

/// The boundary of one wire, in the face's parameter space.
///
/// Each edge is drawn at its own polyline's points ([`edge_polyline`]) and
/// its pcurve evaluated at their parameters, so two faces sharing the edge
/// agree on where its points are.
/// With `in_chart`, each edge's samples are also refined until its chords
/// keep to its pcurve in the chart ([`refine_in_chart`]): a ring to place
/// points against rather than to weld a mesh along.
fn boundary_ring(
    model: &Model,
    wire: &Shape,
    surface: ogeom_topo::SurfaceId,
    deflection: Deflection,
    finer: &EdgeChords,
    in_chart: bool,
    tol: Tolerances,
) -> OgeomResult<WalkedRing> {
    let mut ring: Vec<Point2> = Vec::new();
    let mut anchors: Vec<Option<Point>> = Vec::new();
    // Ambiguous folds: where an edge was translated a whole period to
    // continue the walk from a start that already coincided modulo the
    // period; the fold was one of two defensible choices, recorded so a
    // mis-wound ring can be unwound.
    let mut folds: Vec<(usize, f64)> = Vec::new();
    // Half-period jumps whose side could not be decided when walked.
    let mut ties: Vec<usize> = Vec::new();
    let mut met = true;
    // Whether the ring's first and last points are vertices drawn one step
    // past their curves' ends.
    let mut stepped_ring = [false; 2];
    // Whether each chart direction comes back on itself: periodic, or
    // closed without repeating. Asked once here: closure on a spline is a
    // walk down a control column, too costly to repeat at every edge of
    // every face.
    let (wraps_u, wraps_v) = model
        .geometry()
        .surface(surface)
        .map_or((false, false), |g| {
            use ogeom_geom::Surface as _;
            (
                g.is_periodic_u() || g.is_closed_u(tol),
                g.is_periodic_v() || g.is_closed_v(tol),
            )
        });

    // Start the walk off a seam if the wire allows it: a seam's side is
    // chosen by continuity with the point already walked to, and continuity
    // needs something to continue from. The ring is cyclic, so rotating the
    // walk changes nothing it reports.
    let mut children = model.ordered_children_of(wire)?;
    let is_seam = |model: &Model, e: &Shape| -> bool {
        model
            .node(e)
            .and_then(|n| n.data().as_edge())
            .and_then(|d| d.pcurve_for(surface, e.location()))
            .is_some_and(|r| matches!(r, EdgeRepr::Seam { .. }))
    };
    if let Some(start) = children.iter().position(|e| !is_seam(model, e)) {
        children.rotate_left(start);
    }

    // The column each seam edge's first traversal effectively walked (after
    // folding), and whether its two sides differ in u or in v. A seam bounds
    // its face twice, and the two traversals must bracket the ring exactly
    // one period apart; the walk checks the second against this record.
    let mut seam_walked: ogeom_core::FastMap<ogeom_topo::TShapeId, Point2> =
        ogeom_core::FastMap::default();

    for edge in children {
        let Some(node) = model.node(&edge) else {
            ogeom_bail!(Dangling, "edge is not in this model");
        };
        let NodeData::Edge(data) = node.data() else {
            ogeom_bail!(Construction, "edge node holds no edge data");
        };
        // Whether this edge is a seam, and if so whether its two sides
        // differ in u (true) or in v.
        let mut seam: Option<bool> = None;
        let (pcurve_id, pcurve_range) = match data.pcurve_for(surface, edge.location()) {
            Some(EdgeRepr::PCurve { curve, range, .. }) => (*curve, *range),
            // A seam edge runs along a closed surface's join and bounds its
            // face twice: up one side of the parameter rectangle and down
            // the other. Which side this occurrence takes is decided by the
            // ring itself: the side whose oriented start continues the point
            // already walked to. Orientation flags cannot answer it (a
            // reversed face flips every occurrence while the chart columns
            // stay where they were built), but the chart can.
            Some(EdgeRepr::Seam {
                forward,
                reversed,
                range,
                ..
            }) => {
                let (f, r) = (*forward, *reversed);
                let side_start = |id: ogeom_topo::PCurveId| -> Option<Point2> {
                    model.geometry().pcurve(id)?.point_at(range.0, tol).ok()
                };
                seam = Some(match (side_start(f), side_start(r)) {
                    (Some(a), Some(b)) => (a.x - b.x).abs() >= (a.y - b.y).abs(),
                    _ => true,
                });
                let picked = if let Some(last) = ring.last().copied() {
                    let start_of = |id: ogeom_topo::PCurveId| -> Option<Point2> {
                        let pc = model.geometry().pcurve(id)?;
                        let t = if edge.orientation() == Orientation::Reversed {
                            range.1
                        } else {
                            range.0
                        };
                        pc.point_at(t, tol).ok()
                    };
                    match (start_of(f), start_of(r)) {
                        (Some(a), Some(b)) => {
                            if last.distance(a) <= last.distance(b) {
                                f
                            } else {
                                r
                            }
                        }
                        _ => f,
                    }
                } else if edge.orientation() == Orientation::Reversed {
                    r
                } else {
                    f
                };
                (picked, *range)
            }
            _ => ogeom_bail!(
                Construction,
                "edge has no pcurve on this face, so the face cannot be \
                 triangulated in its own parameter space"
            ),
        };
        let Some(pcurve) = model.geometry().pcurve(pcurve_id) else {
            ogeom_bail!(Dangling, "pcurve is not in this model");
        };

        // Sample where the *3D* curve says to, so an adjacent face lands on
        // the same points, and *anchor* the boundary vertices to that curve
        // too: two faces sharing an edge lift the same parameters through
        // different surfaces, and on an imported file those surfaces
        // disagree by the file's own slop. The edge is the shared authority,
        // so its points are the positions both faces use, and the weld is a
        // matter of identity rather than luck.
        let mut edge_anchors: Vec<Option<Point>> = Vec::new();
        // An edge drawn finer is drawn finer for *every* face that bounds
        // it, which is the whole point: the two sides must agree point for
        // point or the weld has nothing to join.
        let along = match finer.get(&edge.node().index()) {
            Some(chord) => Deflection {
                chord: *chord,
                ..deflection
            },
            None => deflection,
        };
        // The ends drawn from the edge's vertices rather than its curve.
        let mut snapped_ends: Vec<usize> = Vec::new();
        // Whether the edge's start and end are vertices drawn one step past
        // its curve's ends, in the curve's direction.
        let mut stepped_ends = [false; 2];
        let samples = match sample_parameters(model, data, along, tol)? {
            Some((parameters, edge_met)) => {
                met &= edge_met;
                let mut parameters = match (in_chart, model.geometry().surface(surface)) {
                    (true, Some(geometry)) => {
                        let (refined, chart_met) = refine_in_chart(
                            parameters,
                            |t| to_pcurve_parameter(t, data, pcurve_range),
                            pcurve,
                            geometry,
                            along,
                            tol,
                        );
                        met &= chart_met;
                        refined
                    }
                    _ => parameters,
                };
                if let Some(EdgeRepr::Curve3d { curve, .. }) = data.curve3d()
                    && let Some(geometry) = model.geometry().curve(*curve)
                    && let Ok(edge_placement) = edge.transform(model.datums())
                {
                    for t in &parameters {
                        edge_anchors.push(
                            geometry
                                .point_at(*t, tol)
                                .ok()
                                .map(|p| edge_placement.apply(p)),
                        );
                    }
                    // An edge's curve may end off its vertex by up to the
                    // vertex's tolerance, and the next edge's off it on
                    // another side: each face round the vertex would then
                    // close its own gap between two of the ends, and
                    // between them the faces leave a hole. Drawn from the
                    // vertex itself, every face round it meets at one
                    // point.
                    //
                    // Moved onto a vertex further off than the chord, an
                    // end takes the edge's last segment off its faces by
                    // that much along the whole segment: a straight edge
                    // is one segment, and a face bounded by such edges
                    // alone is drawn wholly beside its surface. The edge
                    // then takes one more point on its curve, as far in
                    // from the end as the vertex stands off it, and only
                    // the step from there to the vertex leaves the curve.
                    // The point is the edge's own, so every face along the
                    // edge draws it, and every face round the vertex still
                    // meets there. In the chart the vertex keeps the
                    // curve's end: solved onto the vertex, it can stand
                    // inside the face, where the triangulation joins it
                    // across the face and takes the face off its surface
                    // again.
                    let chord = along.chord * edge_placement.scale_factor().abs();
                    let corners: Vec<(Point, f64)> = model
                        .ordered_children_of(&edge)?
                        .iter()
                        .filter_map(|v| {
                            let data = model.node(v)?.data().as_vertex()?;
                            let at = v.transform(model.datums()).ok()?.apply(data.point);
                            Some((at, data.tolerance.get()))
                        })
                        .collect();
                    let reach = data.tolerance.get();
                    // The end last, so the start's index still holds.
                    let mut snapped = [false; 2];
                    for i in [edge_anchors.len().saturating_sub(1), 0] {
                        let Some(Some(end)) = edge_anchors.get(i).copied() else {
                            continue;
                        };
                        let Some(&(at, within)) = corners
                            .iter()
                            .min_by(|a, b| a.0.distance(end).total_cmp(&b.0.distance(end)))
                        else {
                            continue;
                        };
                        let off = at.distance(end);
                        if off > within + reach {
                            continue;
                        }
                        let next = if i == 0 { 1 } else { i - 1 };
                        let mut stepped = false;
                        if off > chord
                            && let (Some(&from), Some(&to)) =
                                (parameters.get(i), parameters.get(next))
                            && let Some(Some(inner)) = edge_anchors.get(next).copied()
                            && inner.distance(end) > 2.0 * off
                        {
                            let t = from + (to - from) * off / inner.distance(end);
                            let step = if i == 0 { 1 } else { i };
                            parameters.insert(step, t);
                            edge_anchors.insert(
                                step,
                                geometry
                                    .point_at(t, tol)
                                    .ok()
                                    .map(|p| edge_placement.apply(p)),
                            );
                            stepped = true;
                        }
                        let i = if i == 0 { 0 } else { edge_anchors.len() - 1 };
                        edge_anchors[i] = Some(at);
                        snapped[usize::from(i != 0)] = !stepped;
                        stepped_ends[usize::from(i != 0)] = stepped;
                    }
                    if snapped[0] {
                        snapped_ends.push(0);
                    }
                    if snapped[1] {
                        snapped_ends.push(edge_anchors.len() - 1);
                    }
                }
                map_to_pcurve(&parameters, data, pcurve_range)
            }
            // No 3D curve to defer to: the pcurve's own shape, measured in
            // space through the surface, so the chord tolerance means the
            // same thing it means everywhere else.
            None => {
                let Some(geometry) = model.geometry().surface(surface) else {
                    ogeom_bail!(Dangling, "face refers to a surface not in this model");
                };
                // At the chord the faces agreed for this edge, like the
                // 3D path above: a face that asked for finer gets it here
                // too.
                let (_, parameters) = crate::discretize::discretize_on_surface(
                    pcurve,
                    pcurve_range,
                    geometry,
                    along,
                    tol,
                )?;
                if in_chart {
                    let (refined, chart_met) =
                        refine_in_chart(parameters, |t| t, pcurve, geometry, along, tol);
                    met &= chart_met;
                    refined
                } else {
                    parameters
                }
            }
        };

        let mut points: Vec<Point2> = samples
            .iter()
            .map(|u| pcurve.point_at(*u, tol))
            .collect::<OgeomResult<_>>()?;
        // A point drawn from the vertex stands where the vertex does in the
        // chart too, or the boundary runs from the curve's end to the
        // vertex and back, a fold the triangulation covers twice. Solved
        // onto the vertex from the pcurve's own end.
        if !snapped_ends.is_empty()
            && let Some(geometry) = model.geometry().surface(surface)
            && let Ok(placement) = wire.transform(model.datums())
            && let Ok(local) = placement.inverse()
        {
            use ogeom_geom::Surface as _;
            for &i in &snapped_ends {
                let (Some(Some(target)), Some(start)) = (edge_anchors.get(i), points.get(i)) else {
                    continue;
                };
                let target = local.apply(*target);
                let mut uv = *start;
                let mut solved = false;
                for _ in 0..8 {
                    let Ok(jet) = geometry.jet_at(uv.x, uv.y, tol) else {
                        break;
                    };
                    let r = jet.point - target;
                    if r.magnitude() <= tol.confusion() {
                        solved = true;
                        break;
                    }
                    let (a, b, c) = (jet.du.dot(jet.du), jet.du.dot(jet.dv), jet.dv.dot(jet.dv));
                    let det = a.mul_add(c, -(b * b));
                    if det.abs() <= f64::MIN_POSITIVE {
                        break;
                    }
                    let (gu, gv) = (jet.du.dot(r), jet.dv.dot(r));
                    uv = Point2::new(
                        uv.x - (c * gu - b * gv) / det,
                        uv.y - (a * gv - b * gu) / det,
                    );
                }
                // Short of the vertex, the point that came nearest it.
                if !solved
                    && let Ok(p) = geometry.point_at(uv.x, uv.y, tol)
                    && let Ok(q) = geometry.point_at(start.x, start.y, tol)
                {
                    solved = p.distance(target) < q.distance(target);
                }
                if solved && let Some(slot) = points.get_mut(i) {
                    *slot = uv;
                }
            }
        }
        if edge.orientation() == Orientation::Reversed {
            points.reverse();
            edge_anchors.reverse();
            stepped_ends.reverse();
        }
        if edge_anchors.len() != points.len() {
            edge_anchors = vec![None; points.len()];
        }
        // The ends of an edge belong to its *vertices* (the one authority
        // every face and every neighbouring edge shares), but only within the
        // tolerance the vertex itself records. An imported curve ends within
        // the vertex's widened tolerance of it, and lifting the ends through
        // the curve alone would leave each corner split as many ways as there
        // are curves meeting there. A vertex that sits *beyond* its stated
        // tolerance from the curve is not describing the curve's end at all,
        // and the curve stays the authority.
        if !points.is_empty()
            && let Ok(vs) = model.children_of(&edge)
            && vs.len() >= 2
            && let Ok(edge_placement) = edge.transform(model.datums())
        {
            let point_of = |v: &Shape| -> Option<(Point, f64)> {
                let data = model.node(v)?.data().as_vertex()?;
                Some((edge_placement.apply(data.point), data.tolerance.get()))
            };
            let (from, to) = if edge.orientation() == Orientation::Reversed {
                (&vs[vs.len() - 1], &vs[0])
            } else {
                (&vs[0], &vs[vs.len() - 1])
            };
            if let Some((p, within)) = point_of(from)
                && let Some(a) = edge_anchors.first_mut()
                && a.is_none_or(|end| end.distance(p) <= within + tol.confusion())
            {
                *a = Some(p);
            }
            if let Some((p, within)) = point_of(to)
                && let Some(a) = edge_anchors.last_mut()
                && a.is_none_or(|end| end.distance(p) <= within + tol.confusion())
            {
                *a = Some(p);
            }
        }
        // Fold onto the branch that continues the ring. Two faces sharing an
        // edge share its pcurve, and on a periodic surface the pcurve sits in
        // *one* face's window: a cylinder split into two halves has a ruling
        // at u = 0 that the other half needs at u = 2pi. The chart cannot
        // store both. Continuity with the ring being walked recovers the
        // right branch, exactly as the seam sides are chosen.
        if let Some(last) = ring.last().copied()
            && let Some(first) = points.first().copied()
            && let Some(geometry) = model.geometry().surface(surface)
        {
            use ogeom_geom::Surface as _;
            let ((ua, ub), (va, vb)) = geometry.domain();
            // Nearest whole period, ties broken toward *not moving*: a jump
            // of exactly half a period is what a boundary crossing a
            // degenerate row looks like (two rulings into an apex stand
            // half a turn apart, and the connecting run along the apex row
            // lifts to nothing), and folding it would drag the edge a full
            // period from the column its own projection put it on.
            let whole_periods = |gap: f64, span: f64| -> f64 {
                let r = gap / span;
                if (r.fract().abs() - 0.5).abs() <= 1e-9 {
                    r.trunc() * span
                } else {
                    r.round() * span
                }
            };
            let mut shift = Point2::new(0.0, 0.0);
            // A repeated seam edge folds like any other, but never onto its
            // own first traversal. The two traversals bound the face up one
            // side of the chart and down the other, one period apart, and at
            // a degenerate row continuity cannot say so: a cone walked to
            // its apex reaches a corner that maps to the whole row, both
            // sides continue it equally, and folding by nearness closes the
            // ring over nothing. The record decides instead: land exactly a
            // period from the first walk, on the side the ring occupies.
            let prior = seam.and_then(|_| seam_walked.get(&edge.node())).copied();
            if wraps_u && (ub - ua) > 0.0 {
                let span = ub - ua;
                let gap = last.x - first.x;
                shift.x = whole_periods(gap, span);
                let mut bracketed = false;
                if seam == Some(true)
                    && let Some(prior) = prior
                    && (first.x + shift.x - prior.x).abs() < span * 0.5
                {
                    let (lo, hi) = ring
                        .iter()
                        .fold((f64::INFINITY, f64::NEG_INFINITY), |(lo, hi), p| {
                            (lo.min(p.x), hi.max(p.x))
                        });
                    // Bracketing is for the face that wraps the period: a
                    // band whose rims run the whole way round, its two seam
                    // columns a period apart. A *slit* uses one edge twice
                    // without wrapping: the ring stays in a fraction of the
                    // chart, both traversals stand on one column, and
                    // forcing them apart winds the ring, invites a pole row
                    // it never touches, and meshes the complement of the
                    // face. The ring's own reach
                    // says which face this is.
                    if hi - lo >= span * 0.5 {
                        let side = if f64::midpoint(lo, hi) >= prior.x {
                            1.0
                        } else {
                            -1.0
                        };
                        shift.x = prior.x + side * span - first.x;
                        bracketed = true;
                    }
                }
                if !bracketed {
                    if shift.x != 0.0 && (gap - shift.x).abs() <= span * 1e-6 {
                        // The start already stood a whole period from the walk
                        // (the same 3D point), so this fold is a choice, not a
                        // repair; recorded so a mis-wound ring can be unwound.
                        folds.push((ring.len(), shift.x));
                    } else if ((gap / span).fract().abs() - 0.5).abs() <= 1e-9 {
                        // A half-period tie: either side of the degenerate row
                        // was defensible, and if the ring comes out wound the
                        // unwinding starts here rather than at the later fold.
                        ties.push(ring.len());
                    }
                }
            }
            if wraps_v && (vb - va) > 0.0 {
                let span = vb - va;
                shift.y = whole_periods(last.y - first.y, span);
                if seam == Some(false)
                    && let Some(prior) = prior
                    && (first.y + shift.y - prior.y).abs() < span * 0.5
                {
                    let (lo, hi) = ring
                        .iter()
                        .fold((f64::INFINITY, f64::NEG_INFINITY), |(lo, hi), p| {
                            (lo.min(p.y), hi.max(p.y))
                        });
                    let side = if f64::midpoint(lo, hi) >= prior.y {
                        1.0
                    } else {
                        -1.0
                    };
                    shift.y = prior.y + side * span - first.y;
                }
            }
            if shift.x != 0.0 || shift.y != 0.0 {
                for p in &mut points {
                    p.x += shift.x;
                    p.y += shift.y;
                }
            }
        }
        if seam.is_some()
            && let Some(first) = points.first().copied()
        {
            seam_walked.entry(edge.node()).or_insert(first);
        }
        // The previous edge already contributed the shared vertex, but only
        // where the chart agrees it is shared. Two rulings meeting at an
        // apex share the *vertex* while standing apart in the chart, and the
        // run between them along the degenerate row is boundary the ring
        // needs: dropping its start would cut the corner straight through
        // the face's interior.
        //
        // Two things have to hold for a gap to be such a run, and neither
        // alone is enough. It has to be wide: two ends that disagree by the
        // file's slop stand a micron over a radius apart, and a thousandth
        // of the period is three orders above that. And what lies between
        // has to be degenerate: a run along an apex row lifts to one point
        // the whole way, so its chart midpoint lands on the shared vertex,
        // where a wide gap on a live row lifts to somewhere the width of the
        // gap away. Width alone keeps slop on fitted splines. The lift alone
        // keeps every gap too small for its midpoint to land anywhere else.
        // And no fraction of the period alone is right at all: rulings into
        // an apex can stand exactly a quarter turn apart, and a threshold at
        // that fraction reads them as not apart and loses the triangle at
        // the apex.
        let keep_gap = if let (Some(last), Some(first)) = (ring.last(), points.first()) {
            model.geometry().surface(surface).is_some_and(|geometry| {
                use ogeom_geom::Surface as _;
                let ((ua, ub), (va, vb)) = geometry.domain();
                let wide = (wraps_u && (last.x - first.x).abs() > (ub - ua) * 1e-3)
                    || (wraps_v && (last.y - first.y).abs() > (vb - va) * 1e-3);
                if !wide {
                    return false;
                }
                let mid = Point2::new(
                    f64::midpoint(last.x, first.x),
                    f64::midpoint(last.y, first.y),
                );
                let reach = tol.confusion() * 1e4;
                match (
                    geometry.point_at(last.x, last.y, tol),
                    geometry.point_at(mid.x, mid.y, tol),
                    geometry.point_at(first.x, first.y, tol),
                ) {
                    (Ok(a), Ok(m), Ok(b)) => a.distance(m) <= reach && m.distance(b) <= reach,
                    _ => false,
                }
            })
        } else {
            false
        };
        if !ring.is_empty() && !points.is_empty() && !keep_gap {
            points.remove(0);
            edge_anchors.remove(0);
        }
        if ring.is_empty() {
            stepped_ring[0] = stepped_ends[0];
        }
        stepped_ring[1] = stepped_ends[1];
        ring.extend(points);
        anchors.extend(edge_anchors);
    }

    // A closed ring repeats its first point at the end; the triangulator wants
    // it named once.
    if ring.len() > 2
        && let (Some(first), Some(last)) = (ring.first().copied(), ring.last().copied())
    {
        // Equal in the chart, or the same vertex in space: the last edge's
        // curve ends where the first edge's begins to within the file's
        // slop: up to ten microns in a real assembly, recorded on the
        // vertex as its widened tolerance. Kept as two points, the ring
        // closes with a fold back over its own first segment, a crossing
        // a fraction of a micron deep that the triangulation refuses as a
        // constraint and the face is then drawn six times finer for.
        // The same vertex in space is not enough on its own: a ring that
        // winds a periodic chart ends a whole period from where it began
        // and lifts to the same point, and that closing is a seam, not
        // slop. Close in the chart too (within a thousandth of the ring's
        // own extent), or the ring is left to the winding rule below.
        let extent = ring.iter().fold(
            (
                Point2::new(f64::INFINITY, f64::INFINITY),
                Point2::new(f64::NEG_INFINITY, f64::NEG_INFINITY),
            ),
            |(lo, hi), p| {
                (
                    Point2::new(lo.x.min(p.x), lo.y.min(p.y)),
                    Point2::new(hi.x.max(p.x), hi.y.max(p.y)),
                )
            },
        );
        let extent = (extent.1.x - extent.0.x).max(extent.1.y - extent.0.y);
        let near_in_chart = (first.x - last.x).hypot(first.y - last.y) <= extent * 1e-3;
        let same_vertex = near_in_chart
            && match (anchors.first(), anchors.last()) {
                (Some(Some(a)), Some(Some(b))) => a.distance(*b) <= tol.confusion() * 1e5,
                _ => false,
            };
        // A vertex drawn one step past both the last edge's curve and the
        // first's stands at each curve's end in the chart: two points for
        // one vertex, closed as the walk closes between any two edges. Ends
        // a whole period apart lift to one point and stay to be wound.
        let stepped = stepped_ring == [true, true]
            && match (anchors.first(), anchors.last()) {
                (Some(Some(a)), Some(Some(b))) => a.distance(*b) <= tol.confusion(),
                _ => false,
            }
            && model.geometry().surface(surface).is_some_and(|geometry| {
                use ogeom_geom::Surface as _;
                match (
                    geometry.point_at(first.x, first.y, tol),
                    geometry.point_at(last.x, last.y, tol),
                ) {
                    (Ok(a), Ok(b)) => a.distance(b) > tol.confusion(),
                    _ => false,
                }
            });
        if first.is_equal(last, tol) || same_vertex || stepped {
            ring.pop();
            anchors.pop();
        }
    }
    Ok((ring, anchors, met, folds, ties))
}

/// Close a ring that winds a periodic direction of its chart and found no
/// partner to pair with.
///
/// A ring that winds one periodic direction of a doubly-periodic surface
/// (a diagonal loop on a torus) ends a whole period from where it
/// began, and the face is the band between the chain and its own
/// translate one period over in the *other* periodic direction, joined
/// at the ends by columns that lift to one 3D circle. The translate's
/// anchors are the same 3D points, and the joining columns' two copies
/// lift identically, so the weld closes them exactly as it closes a seam.
/// On a cone or a sphere the ring closes against the row where the
/// surface collapses to a point (the apex or the pole), which every `u`
/// reaches: that closure costs no area, because the row has none.
///
/// A ring that has a partner (the other rim of a band) is not closed
/// here but paired with it, or each rim would close the whole surface on
/// its own and the two would cancel where they overlap.
fn close_wound_ring(
    ring: &mut Vec<Point2>,
    anchors: &mut Vec<Option<Point>>,
    geometry: &SurfaceGeometry,
    tol: Tolerances,
) {
    let (Some(first), Some(last)) = (ring.first().copied(), ring.last().copied()) else {
        return;
    };
    use ogeom_geom::Surface as _;
    let ((ua, ub), (va, vb)) = geometry.domain();
    let du = last.x - first.x;
    let dv = last.y - first.y;
    let winds_u = geometry.is_periodic_u()
        && (du.abs() - (ub - ua)).abs() <= (ub - ua) * 1e-3
        && dv.abs() <= (vb - va).max(1.0) * 1e-3;
    let winds_v = geometry.is_periodic_u()
        && geometry.is_periodic_v()
        && (dv.abs() - (vb - va)).abs() <= (vb - va) * 1e-3
        && du.abs() <= (ub - ua).max(1.0) * 1e-3;
    // Where does a u-winding ring close against? On a doubly
    // periodic surface, its own translate one v-period over. On a
    // cone or sphere, the row where the surface collapses to a point
    // (the apex or the pole), which every u reaches: the closure
    // costs no area error because the row has none.
    let degenerate_row = |v: f64| -> bool {
        let (Ok(p), Ok(q), Ok(r)) = (
            geometry.point_at(ua, v, tol),
            geometry.point_at(f64::midpoint(ua, ub), v, tol),
            geometry.point_at(ub, v, tol),
        ) else {
            return false;
        };
        p.distance(q) <= tol.confusion() * 10.0 && p.distance(r) <= tol.confusion() * 10.0
    };
    let target_v = if winds_u && !geometry.is_periodic_v() {
        // The nearer degenerate row, if either end has one.
        let mid_v = f64::midpoint(first.y, last.y);
        if degenerate_row(va) && (mid_v - va).abs() <= (mid_v - vb).abs() {
            Some(va)
        } else if degenerate_row(vb) {
            Some(vb)
        } else if degenerate_row(va) {
            Some(va)
        } else {
            None
        }
    } else {
        None
    };
    let column_steps = 8;
    if let Some(v_apex) = target_v {
        // Down the seam column to the apex row, across it, and back
        // up: the row has no length in space, so the closure adds no
        // area and its lifted points weld to the one apex.
        let row_steps = ring.len().max(8);
        for k in 1..=column_steps {
            let f = f64::from(k) / f64::from(column_steps);
            ring.push(Point2::new(last.x, last.y + (v_apex - last.y) * f));
            anchors.push(None);
        }
        for k in 1..row_steps {
            #[allow(clippy::cast_precision_loss)]
            let f = k as f64 / row_steps as f64;
            ring.push(Point2::new(last.x + (first.x - last.x) * f, v_apex));
            anchors.push(None);
        }
        for k in 0..column_steps {
            let f = f64::from(column_steps - k) / f64::from(column_steps);
            ring.push(Point2::new(first.x, first.y + (v_apex - first.y) * f));
            anchors.push(None);
        }
    } else if (winds_u && geometry.is_periodic_v()) || winds_v {
        let shift = if winds_u {
            Point2::new(0.0, -(vb - va))
        } else {
            Point2::new(-(ub - ua), 0.0)
        };
        let chain: Vec<Point2> = ring.clone();
        let chain_anchors = anchors.clone();
        // Down from the chain's end to its translate's end.
        for k in 1..=column_steps {
            let f = f64::from(k) / f64::from(column_steps);
            ring.push(Point2::new(last.x + shift.x * f, last.y + shift.y * f));
            anchors.push(None);
        }
        // The translate, walked back.
        for (p, a) in chain.iter().rev().zip(chain_anchors.iter().rev()).skip(1) {
            ring.push(Point2::new(p.x + shift.x, p.y + shift.y));
            anchors.push(*a);
        }
        // Up from the translate's start back to the chain's start,
        // stopping one step short of closing.
        for k in 1..column_steps {
            let f = f64::from(column_steps - k) / f64::from(column_steps);
            ring.push(Point2::new(first.x + shift.x * f, first.y + shift.y * f));
            anchors.push(None);
        }
    }
}

/// Parameters at which to sample an edge, taken from its 3D curve.
fn sample_parameters(
    model: &Model,
    data: &ogeom_topo::EdgeData,
    deflection: Deflection,
    tol: Tolerances,
) -> OgeomResult<Option<(Vec<f64>, bool)>> {
    Ok(edge_polyline(model, data, deflection, tol)?
        .map(|line| (line.parameters, line.deflection_met)))
}

/// An edge's polyline in its curve's own frame: its 3D curve discretized
/// at `deflection`, refined where a chart chord of it on a face it bounds
/// sags past the repair's reach ([`refine_in_every_chart`]). Every face
/// bounding the edge draws its boundary through these points, and the
/// edge's own polyline is these points; `None` for an edge with no 3D
/// curve.
pub(crate) fn edge_polyline(
    model: &Model,
    data: &ogeom_topo::EdgeData,
    deflection: Deflection,
    tol: Tolerances,
) -> OgeomResult<Option<crate::discretize::Polyline>> {
    let Some(EdgeRepr::Curve3d { curve, range, .. }) = data.curve3d() else {
        return Ok(None);
    };
    let Some(geometry) = model.geometry().curve(*curve) else {
        ogeom_bail!(Dangling, "curve is not in this model");
    };
    let line = discretize(geometry, *range, deflection, tol)?;
    let drawn = line.parameters.len();
    let (parameters, met) =
        refine_in_every_chart(model, data, line.parameters.clone(), deflection, tol);
    if parameters.len() == drawn {
        return Ok(Some(crate::discretize::Polyline {
            deflection_met: line.deflection_met && met,
            ..line
        }));
    }
    let points = parameters
        .iter()
        .map(|&t| geometry.point_at(t, tol))
        .collect::<OgeomResult<Vec<Point>>>()?;
    Ok(Some(crate::discretize::Polyline {
        points,
        parameters,
        deflection_met: line.deflection_met && met,
    }))
}

/// Map parameters on the 3D curve onto the pcurve's own range.
///
/// The two are parameterized over their own intervals; `same_parameter` means
/// they agree *proportionally*, which is what this converts.
fn map_to_pcurve(
    parameters: &[f64],
    data: &ogeom_topo::EdgeData,
    pcurve_range: (f64, f64),
) -> Vec<f64> {
    let Some(EdgeRepr::Curve3d { range, .. }) = data.curve3d() else {
        return parameters.to_vec();
    };
    let (ca, cb) = *range;
    let (pa, pb) = pcurve_range;
    if (cb - ca).abs() <= f64::MIN_POSITIVE {
        return parameters.to_vec();
    }
    parameters
        .iter()
        .map(|u| pa + (pb - pa) * (u - ca) / (cb - ca))
        .collect()
}

/// [`map_to_pcurve`] for one parameter.
fn to_pcurve_parameter(u: f64, data: &ogeom_topo::EdgeData, pcurve_range: (f64, f64)) -> f64 {
    map_to_pcurve(&[u], data, pcurve_range)[0]
}

/// An edge's sample `parameters` refined until its polyline in the chart
/// keeps within the chord of its pcurve, read in space through the
/// surface's tangent scales, as a classifier's band reads a chart step:
/// `(du |S_u|, dv |S_v|)`. Samples placed by the chord in space alone keep
/// the edge's chords near the edge in space, but where the edge bends more
/// in the chart than in space (a parabola on a cone) the chart's straight
/// chords stand off its pcurve by many chords.
///
/// Each span is read at its quarters and its middle, the greatest of the
/// three distances from its chart chord, and split at its middle where that
/// exceeds the chord. `to_pcurve` maps a sample parameter onto the pcurve.
/// A span whose surface has no tangents where read is left as drawn; the
/// flag is false where `max_segments` stopped the refinement short.
fn refine_in_chart(
    parameters: Vec<f64>,
    to_pcurve: impl Fn(f64) -> f64,
    pcurve: &ogeom_geom::PlanarCurve,
    surface: &SurfaceGeometry,
    deflection: Deflection,
    tol: Tolerances,
) -> (Vec<f64>, bool) {
    let at = |t: f64| pcurve.point_at(to_pcurve(t), tol).ok();
    let stands_off = |a: f64, b: f64| -> bool {
        let (Some(pa), Some(pb)) = (at(a), at(b)) else {
            return false;
        };
        [0.25, 0.5, 0.75].into_iter().any(|share| {
            let t = (b - a).mul_add(share, a);
            let Some(q) = at(t) else {
                return false;
            };
            let Ok((du, dv)) = surface.d1_at(q.x, q.y, tol) else {
                return false;
            };
            let (su, sv) = (du.magnitude(), dv.magnitude());
            let scaled = |p: Point2| Point2::new((p.x - q.x) * su, (p.y - q.y) * sv);
            let (a, b) = (scaled(pa), scaled(pb));
            let d = b - a;
            let squared = d.dot(d);
            let w = if squared <= f64::MIN_POSITIVE {
                0.0
            } else {
                (-(a.x * d.x + a.y * d.y) / squared).clamp(0.0, 1.0)
            };
            let foot = a + d * w;
            foot.x.hypot(foot.y) > deflection.chord
        })
    };
    let mut parameters = parameters;
    let ceiling = deflection.max_segments.max(parameters.len());
    let mut i = 0;
    while i + 1 < parameters.len() {
        let (a, b) = (parameters[i], parameters[i + 1]);
        let mid = f64::midpoint(a, b);
        if mid > a && mid < b && stands_off(a, b) {
            if parameters.len() > ceiling {
                return (parameters, false);
            }
            parameters.insert(i + 1, mid);
        } else {
            i += 1;
        }
    }
    (parameters, true)
}

/// An edge's sample `parameters` refined until, on every curved surface
/// it has a pcurve on, no span's chart chord lifts farther off the line
/// through its lifted ends than the repair lets a chart edge sag
/// ([`REPAIR_SAG`] chords). A face is drawn from its rings' chart chords,
/// and the repair splits every triangle along one that sags so, but never
/// a span the face shares with its neighbour: where an edge bends in the
/// chart and not in space (a curve where a plane cuts a torus near the
/// tube's rim), the triangles along such a span draw the repair's points
/// towards it, off the face by the sag. Each span is read at its middle and split there, the
/// pcurves taken in the order the edge holds them and the round repeated
/// until none splits a span: the samples depend on the edge alone, so
/// every face bounding it draws the same points. A plane's chart is its
/// space and is not asked. The flag is false where `max_segments` stopped
/// a refinement short.
fn refine_in_every_chart(
    model: &Model,
    data: &ogeom_topo::EdgeData,
    parameters: Vec<f64>,
    deflection: Deflection,
    tol: Tolerances,
) -> (Vec<f64>, bool) {
    use ogeom_geom::Surface as _;
    let charts: Vec<(&ogeom_geom::PlanarCurve, &SurfaceGeometry, (f64, f64))> = data
        .representations
        .iter()
        .filter_map(|repr| {
            let (curve, surface, range) = match repr {
                EdgeRepr::PCurve {
                    curve,
                    surface,
                    range,
                    ..
                } => (*curve, *surface, *range),
                EdgeRepr::Seam {
                    forward,
                    surface,
                    range,
                    ..
                } => (*forward, *surface, *range),
                _ => return None,
            };
            let surface = model.geometry().surface(surface)?;
            if matches!(surface, SurfaceGeometry::Plane(_)) {
                return None;
            }
            Some((model.geometry().pcurve(curve)?, surface, range))
        })
        .collect();
    let mut parameters = parameters;
    let ceiling = deflection.max_segments.max(parameters.len());
    loop {
        let before = parameters.len();
        for &(pcurve, surface, range) in &charts {
            let at = |t: f64| {
                pcurve
                    .point_at(to_pcurve_parameter(t, data, range), tol)
                    .ok()
            };
            let lift = |p: Point2| surface.point_at(p.x, p.y, tol).ok();
            let bows = |a: f64, b: f64| -> bool {
                let (Some(pa), Some(pb)) = (at(a), at(b)) else {
                    return false;
                };
                let (Some(start), Some(end)) = (lift(pa), lift(pb)) else {
                    return false;
                };
                lift(pa + (pb - pa) * 0.5)
                    .is_some_and(|m| chord_sag(start, end, m, tol) > deflection.chord * REPAIR_SAG)
            };
            let mut i = 0;
            while i + 1 < parameters.len() {
                let (a, b) = (parameters[i], parameters[i + 1]);
                let mid = f64::midpoint(a, b);
                if mid > a && mid < b && bows(a, b) {
                    if parameters.len() > ceiling {
                        return (parameters, false);
                    }
                    parameters.insert(i + 1, mid);
                } else {
                    i += 1;
                }
            }
        }
        if parameters.len() == before {
            return (parameters, true);
        }
    }
}

/// The edge of a surface's domain, as a boundary ring.
///
/// Refined, not just the four corners. A triangulation only ever connects the
/// points it is given, so a boundary named by its corners alone forces long
/// triangles reaching right across the domain to find one: on a sphere, a
/// sliver from the equator to the pole. The interior refinement cannot fix that;
/// the missing points are on the boundary.
fn domain_ring(surface: &SurfaceGeometry, deflection: Deflection, tol: Tolerances) -> Vec<Point2> {
    let ((ua, ub), (va, vb)) = surface.domain();
    if ![ua, ub, va, vb].iter().all(|x| x.is_finite()) {
        return Vec::new();
    }

    let along_u = |v: f64| refine_line(surface, (ua, ub), |u| (u, v), deflection, tol).0;
    let along_v = |u: f64| refine_line(surface, (va, vb), |v| (u, v), deflection, tol).0;

    // Counter-clockwise around the rectangle. Each side drops its final point,
    // which the next side contributes: a repeated vertex would be a
    // zero-length boundary edge, and a constraint of zero length is not one.
    let (bottom, top) = (along_u(va), along_u(vb));
    let (left, right) = (along_v(ua), along_v(ub));
    let mut ring = Vec::new();
    ring.extend(
        bottom[..bottom.len() - 1]
            .iter()
            .map(|u| Point2::new(*u, va)),
    );
    ring.extend(right[..right.len() - 1].iter().map(|v| Point2::new(ub, *v)));
    ring.extend(top[1..].iter().rev().map(|u| Point2::new(*u, vb)));
    ring.extend(left[1..].iter().rev().map(|v| Point2::new(ua, *v)));
    ring
}

/// A triangulation still in parameter space, before it is lifted onto the
/// surface.
struct PlanarMesh {
    /// The `(u, v)` of each vertex.
    parameters: Vec<(f64, f64)>,
    /// Triangles as indices into `parameters`.
    triangles: Vec<[u32; 3]>,
    /// Minimum triangle count, accounting for repeated junction vertices.
    boundary_triangles: usize,
    /// The boundary alone did not triangulate to the count a boundary that
    /// encloses a region gives: it crosses itself somewhere.
    crossed: bool,
}

/// Triangulate a region in parameter space, given its boundary rings.
///
/// The first ring is the outer boundary; the rest are holes.
fn triangulate_region(
    rings: &[Vec<Point2>],
    surface: &SurfaceGeometry,
    deflection: Deflection,
    rows_kept: Rows,
    tol: Tolerances,
) -> OgeomResult<PlanarMesh> {
    region_mesh(rings, surface, deflection, rows_kept, Upto::Whole, tol)
}

/// How far [`region_mesh`] goes.
#[derive(Clone, Copy, PartialEq, Eq)]
enum Upto {
    /// The whole mesh.
    Whole,
    /// Only as far as whether the rings crossed: where they did, the mesh
    /// comes back empty, its boundary count and verdict as they are.
    Crossing,
}

/// [`triangulate_region`], stopping where `upto` says.
fn region_mesh(
    rings: &[Vec<Point2>],
    surface: &SurfaceGeometry,
    deflection: Deflection,
    rows_kept: Rows,
    upto: Upto,
    tol: Tolerances,
) -> OgeomResult<PlanarMesh> {
    // A chart mangled enough (trims fitted through millimetres of boundary
    // error) can drive the triangulation library past its own asserts, and
    // a panic in a dependency is a crash in every consumer. The rings are
    // screened for what provably breaks it, and whatever still slips
    // through is caught at this boundary and spoken as the refusal it is:
    // the kernel's contract is refusal by name, never a crash on bad input.
    for ring in rings {
        let mut extent = 0.0_f64;
        for p in ring {
            if !p.x.is_finite() || !p.y.is_finite() {
                ogeom_bail!(
                    NotDone,
                    "a boundary ring carries a non-finite chart coordinate; \
                     the face's trim does not describe a region"
                );
            }
            extent = extent.max(p.x.abs()).max(p.y.abs());
        }
        if extent > 1e12 {
            ogeom_bail!(
                NotDone,
                "a boundary ring reaches {extent:.1e} in the chart; a trim \
                 that far out describes no face"
            );
        }
    }
    match std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
        triangulate_region_inner(rings, surface, deflection, rows_kept, upto, tol)
    })) {
        Ok(result) => result,
        Err(_) => ogeom_bail!(
            NotDone,
            "the face's boundary broke the triangulation; the chart is too \
             degenerate to mesh"
        ),
    }
}

/// The chart stretched per axis to the surface's own metric, so that the
/// triangulation (Delaunay in the chart) sees distances as space does.
///
/// A fitted strip a fifth of a millimetre wide and a centimetre long may
/// carry `u` over a fiftieth of a unit and `v` over one: in the chart the
/// long way is the short way, sixty times over, and Delaunay, which
/// connects nearest neighbours *in the chart*, joins points along the
/// strip across several columns rather than to the row beside them. The
/// triangles it makes are slivers in space that lift folded (a flat
/// triangle spanning a bend the surface takes in between, its normal
/// pointing where none of its vertices' do), and the face shades as a
/// quilt of creases. Scaled by the mean tangent length each way, the
/// chart is the surface to first order, and Delaunay in it is Delaunay
/// on the surface.
#[derive(Debug, Clone, Copy)]
struct ChartScale {
    su: f64,
    sv: f64,
}

impl ChartScale {
    /// The mean tangent length each way over the rings' extent, the longer
    /// normalized to one; `(1, 1)` where the surface will not say.
    fn of(surface: &SurfaceGeometry, rings: &[Vec<Point2>], tol: Tolerances) -> Self {
        let (lo, hi) = chart_extent(rings);
        if !(lo.x.is_finite() && hi.x.is_finite() && lo.y.is_finite() && hi.y.is_finite()) {
            return Self { su: 1.0, sv: 1.0 };
        }
        let (mut du, mut dv, mut n) = (0.0_f64, 0.0_f64, 0usize);
        for i in 0..=2 {
            for j in 0..=2 {
                let u = lo.x + (hi.x - lo.x) * (0.25 + 0.25 * f64::from(i));
                let v = lo.y + (hi.y - lo.y) * (0.25 + 0.25 * f64::from(j));
                if let Ok((a, b)) = surface.d1_at(u, v, tol) {
                    du += a.magnitude();
                    dv += b.magnitude();
                    n += 1;
                }
            }
        }
        if n == 0 || !(du > 0.0 && dv > 0.0) || !du.is_finite() || !dv.is_finite() {
            return Self { su: 1.0, sv: 1.0 };
        }
        let longer = du.max(dv);
        Self {
            su: du / longer,
            sv: dv / longer,
        }
    }

    /// A chart point into the scaled chart.
    fn to(self, u: f64, v: f64) -> SpadePoint<f64> {
        SpadePoint::new(u * self.su, v * self.sv)
    }

    /// A scaled-chart point back into the chart.
    fn from(self, x: f64, y: f64) -> (f64, f64) {
        (x / self.su, y / self.sv)
    }
}

/// The rings' bounding box in the chart.
fn chart_extent(rings: &[Vec<Point2>]) -> (Point2, Point2) {
    rings.iter().flatten().fold(
        (
            Point2::new(f64::INFINITY, f64::INFINITY),
            Point2::new(f64::NEG_INFINITY, f64::NEG_INFINITY),
        ),
        |(lo, hi), p| {
            (
                Point2::new(lo.x.min(p.x), lo.y.min(p.y)),
                Point2::new(hi.x.max(p.x), hi.y.max(p.y)),
            )
        },
    )
}

fn triangulate_region_inner(
    rings: &[Vec<Point2>],
    surface: &SurfaceGeometry,
    deflection: Deflection,
    rows_kept: Rows,
    upto: Upto,
    tol: Tolerances,
) -> OgeomResult<PlanarMesh> {
    let mut cdt: ConstrainedDelaunayTriangulation<SpadePoint<f64>> =
        ConstrainedDelaunayTriangulation::new();
    let sub = Instant::now();
    let mut refused_total = 0usize;
    let mut junctions = 0usize;

    // Everything the triangulation sees is in the scaled chart; the rings'
    // own chart points are kept by their scaled bits so a boundary vertex
    // comes back with exactly the parameters its anchor was keyed by.
    let scale = ChartScale::of(surface, rings, tol);
    let scaled: Vec<Vec<Point2>> = rings
        .iter()
        .map(|ring| {
            ring.iter()
                .map(|p| {
                    let q = scale.to(p.x, p.y);
                    Point2::new(q.x, q.y)
                })
                .collect()
        })
        .collect();
    let mut exact: ogeom_core::FastMap<(u64, u64), (f64, f64)> = ogeom_core::FastMap::default();

    // The boundary edges are constraints, so the triangulation respects the
    // trimming rather than spanning across a hole.
    for (ring, chart) in scaled.iter().zip(rings) {
        let mut ring_handles = Vec::with_capacity(ring.len());
        for (p, uv) in ring.iter().zip(chart) {
            // A chart coordinate can come out subnormal-tiny (the sine of
            // a fold angle, the residue of an exact cancellation), and the
            // triangulation refuses what is, for every purpose, zero. The
            // point is known by where the triangulation holds it, so its
            // exact parameters, and the anchor they key, are found again.
            let held = mitigate_underflow(SpadePoint::new(p.x, p.y));
            exact.insert((held.x.to_bits(), held.y.to_bits()), (uv.x, uv.y));
            let handle = cdt
                .insert(held)
                .map_err(|e| ogeom_core::ogeom_err!(NotDone, "boundary insertion failed: {e}"))?;
            ring_handles.push(handle);
        }
        let distinct: ogeom_core::FastSet<_> = ring_handles.iter().copied().collect();
        junctions += ring_handles.len() - distinct.len();
        let mut refused = 0usize;
        let mut same = 0usize;
        for i in 0..ring_handles.len() {
            let (a, b) = (ring_handles[i], ring_handles[(i + 1) % ring_handles.len()]);
            if a == b {
                same += 1;
            } else if cdt.can_add_constraint(a, b) {
                cdt.add_constraint(a, b);
            } else {
                refused += 1;
                if *MESH_DEBUG_REFINE {
                    let n = ring.len();
                    let partners: Vec<usize> = (0..n)
                        .filter(|&j| j != i && j != (i + 1) % n && j != (i + n - 1) % n)
                        .filter(|&j| {
                            segments_cross(ring[i], ring[(i + 1) % n], ring[j], ring[(j + 1) % n])
                        })
                        .collect();
                    let len = ring[i].distance(ring[(i + 1) % n]);
                    eprintln!(
                        "REFUSED segment {i} of {n} (chart length {len:.3e}) crosses {partners:?}"
                    );
                    for j in [(i + n - 1) % n, i, (i + 1) % n, (i + 2) % n, (i + 3) % n] {
                        eprintln!("   ring[{j}] = ({:.12}, {:.12})", ring[j].x, ring[j].y);
                    }
                }
            }
        }
        if *MESH_DEBUG_REFINE && (refused > 0 || same > 0) {
            eprintln!(
                "CONSTRAINTS ring of {}: {refused} refused, {same} zero-length",
                ring.len()
            );
        }
        refused_total += refused;
    }

    // Before interior points go in, the boundary gives an exact triangle
    // count. A simple polygon with holes gives `b + 2w - 4`, where `b` is
    // its distinct vertices and `w` its rings. A repeated vertex within a
    // ring joins otherwise separate lobes and reduces that count by one.
    let bands = RingBands::over(&scaled);
    // A refused constraint crosses one already present, regardless of count.
    let boundary_triangles = (cdt.num_vertices() + 2 * rings.len()).saturating_sub(4 + junctions);
    let crossed = refused_total > 0
        || inside_by_parity(&cdt).is_none_or(|inside| inside != boundary_triangles);
    if *MESH_DEBUG_REFINE && crossed {
        let points: usize = rings.iter().map(Vec::len).sum();
        eprintln!(
            "PARITY inside {:?} vertices {} points {points} rings {} inner faces {}",
            inside_by_parity(&cdt),
            cdt.num_vertices(),
            rings.len(),
            cdt.num_inner_faces()
        );
    }
    if crossed && upto == Upto::Crossing {
        return Ok(PlanarMesh {
            parameters: Vec::new(),
            triangles: Vec::new(),
            boundary_triangles,
            crossed,
        });
    }

    // Interior points where the surface bends away from the flat triangle. A
    // planar face needs none, which is why this is driven by measured
    // deflection rather than by a fixed grid.
    let boundary_ms = sub.elapsed().as_secs_f64() * 1e3;
    let sub = Instant::now();
    add_interior_points(&mut cdt, rings, surface, deflection, scale, rows_kept, tol)?;
    let interior_ms = sub.elapsed().as_secs_f64() * 1e3;
    let interior_points = cdt.num_vertices();
    let sub = Instant::now();
    let mut rounds_run = 0usize;

    // The scale a degenerate chart triangle is measured against: the
    // region's own span, the longer way. Its *position* is not its size:
    // a face on a cylinder whose axis point sits half a metre away has
    // `v` near −500 000 and a span of twenty, and a scale taken from where
    // the ring sits rather than how far it reaches would call every honest
    // cell a hair.
    let (lo, hi) = chart_extent(&scaled);
    let extent = (hi.x - lo.x).max(hi.y - lo.y);
    let degenerate_area = extent.max(1.0).powi(2) * 1e-12;

    // The grid rows guarantee the deflection along their own lines, but a
    // hole in the face punches a gap through a row, and where the surface is
    // flat in one direction (a cylinder along its axis) there may be no
    // other row for the mesher to reach. The band around the hole then fans
    // from the rim to the far side of the gap, in triangles that sag through
    // the solid by far more than the deflection while every one of their
    // vertices sits exactly on the surface: a mesh whose faces all have the
    // right area and the wrong volume.
    //
    // The repair measures the truth: any kept triangle whose midpoints sag
    // beyond the chord gets its centre inserted, and the loop runs until the
    // mesh is honest or the cap says the surface is being unreasonable.
    // The rings do not change while the mesh is refined, so the containment
    // test they answer is indexed once and reused by every round below and by
    // the output pass.
    if !matches!(surface.kind(), ogeom_geom::SurfaceKind::Plane) {
        // Whether a chart edge sags past the repair's threshold depends on
        // its two ends alone: measured once, however many triangles share
        // it and however many rounds it survives.
        let mut sags: ogeom_core::FastMap<([u64; 2], [u64; 2]), bool> =
            ogeom_core::FastMap::default();
        // And a vertex's point on the surface is asked once, however many
        // chart edges end at it.
        let mut lifted: ogeom_core::FastMap<[u64; 2], Option<Point>> =
            ogeom_core::FastMap::default();
        for _ in 0..REFINEMENT_ROUNDS {
            rounds_run += 1;
            let before = cdt.num_vertices();
            let mut worst: Vec<SpadePoint<f64>> = Vec::new();
            for triangle in cdt.inner_faces() {
                let vertices = triangle.vertices();
                let centre = triangle.center();
                let at = Point2::new(centre.x, centre.y);
                if !bands.holds(at) {
                    continue;
                }
                let corners: [(f64, f64); 3] = [
                    (vertices[0].position().x, vertices[0].position().y),
                    (vertices[1].position().x, vertices[1].position().y),
                    (vertices[2].position().x, vertices[2].position().y),
                ];
                // A chart-degenerate hair is dropped from the mesh, not
                // refined: its 3D chord can sag enormously, and feeding its
                // centre back in only breeds more hairs along the same line.
                let area = ((corners[1].0 - corners[0].0) * (corners[2].1 - corners[0].1)
                    - (corners[1].1 - corners[0].1) * (corners[2].0 - corners[0].0))
                    .abs()
                    / 2.0;
                if area < degenerate_area {
                    continue;
                }
                // The grid already bounds sag along rows and columns, and a
                // grid triangle's diagonal spanning one cell each way may
                // legitimately sag up to the sum (twice the chord). The
                // threshold sits clear above that band so the repair fires
                // only on the fan triangles it exists for, which sag through
                // a hole's gap by tens of chords, and an honest grid
                // (including a perfectly symmetric one, whose mesh must stay
                // symmetric) is left untouched.
                let sagged = (0..3).any(|i| {
                    let (a, b) = (corners[i], corners[(i + 1) % 3]);
                    let bits = |p: (f64, f64)| [p.0.to_bits(), p.1.to_bits()];
                    let key = if bits(a) <= bits(b) {
                        (bits(a), bits(b))
                    } else {
                        (bits(b), bits(a))
                    };
                    *sags.entry(key).or_insert_with(|| {
                        let (from, to) = (scale.from(a.0, a.1), scale.from(b.0, b.1));
                        let mut lift = |p: (f64, f64), uv: (f64, f64)| {
                            *lifted
                                .entry(bits(p))
                                .or_insert_with(|| GridSample::point_at(surface, uv, tol))
                        };
                        let (start, end) = (lift(a, from), lift(b, to));
                        let mid = (f64::midpoint(from.0, to.0), f64::midpoint(from.1, to.1));
                        let sag = match (start, end, GridSample::point_at(surface, mid, tol)) {
                            (Some(start), Some(end), Some(middle)) => {
                                chord_sag(start, end, middle, tol)
                            }
                            // Off the surface's domain; nothing to refine towards.
                            _ => 0.0,
                        };
                        sag > deflection.chord * REPAIR_SAG
                    })
                });
                // A surface curved both ways twists between its rows and
                // columns, and a triangle whose edges each keep to the chord
                // can still stand off it in its middle. A ruled face is not
                // safe either: a triangle joining points on rulings far
                // apart (at a cone's apex, where every ruling meets) cuts
                // inside its curve. So on every curved face a triangle's
                // middle and the middles of its sides, lifted, are held to
                // the chord against its plane.
                let middle_off = || {
                    let lift =
                        |p: (f64, f64)| GridSample::point_at(surface, scale.from(p.0, p.1), tol);
                    let [a, b, c] = corners.map(lift);
                    let (Some(a), Some(b), Some(c)) = (a, b, c) else {
                        return false;
                    };
                    let Ok(n) = (b - a).cross(c - a).normalized(tol) else {
                        return false;
                    };
                    let half = |i: usize, j: usize| {
                        (
                            f64::midpoint(corners[i].0, corners[j].0),
                            f64::midpoint(corners[i].1, corners[j].1),
                        )
                    };
                    [(centre.x, centre.y), half(0, 1), half(1, 2), half(2, 0)]
                        .into_iter()
                        .filter_map(lift)
                        .any(|m| (m - a).dot(n).abs() > deflection.chord)
                };
                if sagged || middle_off() {
                    worst.push(SpadePoint::new(centre.x, centre.y));
                }
            }
            if worst.is_empty() {
                break;
            }
            let mut inserted = 0usize;
            for point in worst {
                // A centre that lands on a vertex already there is not a new
                // point. A sliver whose apex sits on its own base (three
                // grid points on a diagonal, the middle one a rounding off
                // the line) has its centre at that apex to the last bits,
                // and inserting it breeds a hair a few ulps wide, whose
                // centre is the same point again: round after round, a
                // stack of hairs the degenerate filter then drops, and a
                // hole in the face where they were.
                if lands_on_the_mesh(&cdt, point, extent * 1e-9, (1.0, 1.0)) {
                    continue;
                }
                cdt.insert(mitigate_underflow(point)).map_err(|e| {
                    ogeom_core::ogeom_err!(NotDone, "refinement insertion failed: {e}")
                })?;
                inserted += 1;
            }
            if inserted == 0 {
                // Everything that sagged was a hair on a vertex; another
                // round would find the same hairs.
                break;
            }
            if *MESH_DEBUG_REFINE {
                eprintln!(
                    "ROUND {rounds_run}: +{} vertices",
                    cdt.num_vertices() - before
                );
            }
        }
    }

    let refine_ms = sub.elapsed().as_secs_f64() * 1e3;
    if *MESH_DEBUG_REFINE && boundary_ms + interior_ms + refine_ms > 50.0 {
        eprintln!(
            "SUB boundary {boundary_ms:.0}ms interior {interior_ms:.0}ms ({interior_points} verts) refine {refine_ms:.0}ms ({rounds_run} rounds, {} verts)",
            cdt.num_vertices()
        );
    }
    let mut parameters = Vec::new();
    let mut index_of = ogeom_core::FastMap::default();
    for (i, vertex) in cdt.vertices().enumerate() {
        let p = vertex.position();
        index_of.insert(vertex.fix(), i);
        parameters.push(
            exact
                .get(&(p.x.to_bits(), p.y.to_bits()))
                .copied()
                .unwrap_or_else(|| scale.from(p.x, p.y)),
        );
    }

    let mut triangles = Vec::new();
    let (mut dbg_total, mut dbg_outside, mut dbg_degenerate) = (0usize, 0usize, 0usize);
    for triangle in cdt.inner_faces() {
        dbg_total += 1;
        let vertices = triangle.vertices();
        let centre = triangle.center();
        // A constrained Delaunay covers the convex hull of its input, so
        // triangles outside the trimmed region (across a concavity, or inside
        // a hole) have to be discarded. Winding tells them apart.
        if !bands.holds(Point2::new(centre.x, centre.y)) {
            dbg_outside += 1;
            continue;
        }
        // A boundary run whose points differ by last-bit noise (a chart row
        // whose corner came off a different pcurve than its interior) lets
        // the triangulation weave a hair of a triangle along it: chart area
        // measured in ulps, a centroid *on* the boundary that even-odd
        // counting places wherever rounding falls, and a lifted sliver that
        // spans the run in one spurious stroke. It bounds nothing; drop it.
        let area = {
            let [a, b, c] = [
                vertices[0].position(),
                vertices[1].position(),
                vertices[2].position(),
            ];
            ((b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x)).abs() / 2.0
        };
        if area < degenerate_area {
            dbg_degenerate += 1;
            continue;
        }
        #[allow(clippy::cast_possible_truncation)]
        let indices = [
            index_of[&vertices[0].fix()] as u32,
            index_of[&vertices[1].fix()] as u32,
            index_of[&vertices[2].fix()] as u32,
        ];
        triangles.push(indices);
    }

    if *MESH_DEBUG_REFINE {
        eprintln!(
            "FILTER {dbg_total} triangles: {dbg_outside} outside, {dbg_degenerate} degenerate (area < {degenerate_area:.3e}), {} kept",
            triangles.len()
        );
    }
    // No triangles is not refused here: a boundary drawn coarsely enough
    // to cross itself can enclose nothing at all (an annulus narrower
    // than the sag of its rims' polygons, two arcs a hair apart), and the
    // caller's answer to a crossing is to draw the edges finer and ask
    // again. Empty counts as short; only a face still empty after that is
    // refused, by [`triangulate_with`].
    Ok(PlanarMesh {
        parameters,
        triangles,
        boundary_triangles,
        crossed,
    })
}

/// Add interior points wherever the surface deviates from flat by more than the
/// deflection allows.
///
/// The two parameter directions are refined *independently*, and that is not an
/// optimization; it is the difference between converging and not. A uniform
/// grid on a sphere puts as many meridians through the pole as through the
/// equator, so the triangles there become arbitrarily thin slivers in space.
/// The summed area of such a mesh does not approach the sphere's; it grows
/// without bound as the grid tightens. (Schwarz's lantern is the standard
/// example: an inscribed polyhedron whose area diverges under refinement.)
///
/// Refining each direction by its own measured sag fixes it at the source. Near
/// a pole the circle of latitude has almost no radius, so a chord right across
/// it sags by almost nothing and the direction stops subdividing after one or
/// two steps, while the meridian direction, whose curvature does not change,
/// keeps refining. The mesh degenerates into a fan, which is the right shape.
fn add_interior_points(
    cdt: &mut ConstrainedDelaunayTriangulation<SpadePoint<f64>>,
    rings: &[Vec<Point2>],
    surface: &SurfaceGeometry,
    deflection: Deflection,
    scale: ChartScale,
    rows_kept: Rows,
    tol: Tolerances,
) -> OgeomResult<()> {
    // A plane is flat everywhere; sampling it would add points that buy nothing.
    if matches!(surface.kind(), ogeom_geom::SurfaceKind::Plane) {
        return Ok(());
    }
    // Up to half a million grid points each asking the whole boundary: the
    // banded index answers the same exact test from one band's edges.
    let region = RingBands::over(rings);

    let bound = rings
        .iter()
        .flatten()
        .fold(ogeom_math::Aabb::EMPTY, |acc, p| {
            acc.with_point(Point::new(p.x, p.y, 0.0))
        });
    let (Some(low), Some(high)) = (bound.low(), bound.high()) else {
        return Ok(());
    };
    // Within this of a boundary vertex or segment is on it: a hair's width
    // at the scaled chart's scale, the same reach the repair pass keeps.
    let reach = ((high.x - low.x) * scale.su)
        .max((high.y - low.y) * scale.sv)
        .max(1.0)
        * 1e-9;

    // The v resolution has to hold everywhere the region reaches, so its sag is
    // the worst over a spread of u probes rather than the sag along one line.
    // A surface of revolution is the same at every u and a lofted one is not.
    #[allow(clippy::cast_precision_loss)]
    let probes: Vec<f64> = (0..=U_PROBES)
        .map(|i| low.x + (high.x - low.x) * i as f64 / U_PROBES as f64)
        .collect();
    let (rows, row_samples) = refine_lines(
        surface,
        (low.y, high.y),
        &probes,
        |u, v| (u, v),
        deflection,
        tol,
    );

    // Sag alone leaves a cylinder one row: it is straight along its axis, so
    // nothing along `v` ever sags. But the triangulation is Delaunay in the
    // chart, and a bore four hundred millimetres long with one row in the
    // middle hands it two-hundred-millimetre spans from each rim to that
    // row. Delaunay bridges those however it likes, and the repair below
    // fires only at three chords: a triangle a quarter turn wide on a two
    // millimetre bore sags less than that, so it stays, and the bore draws
    // as a square between its holes. Cells are held to a bounded aspect
    // instead: rows close enough, measured in space through the surface,
    // that no triangle between two rows can reach across more than a few
    // columns. Rows are added, never removed, and spread evenly, so a grid
    // that was symmetric stays symmetric. The same the other way round.
    // How far apart two chart points are in space, where the surface says.
    let span = |p: (f64, f64), q: (f64, f64)| -> f64 {
        use ogeom_geom::Surface as _;
        match (
            surface.point_at(p.0, p.1, tol),
            surface.point_at(q.0, q.1, tol),
        ) {
            (Ok(a), Ok(b)) => a.distance(b),
            _ => 0.0,
        }
    };
    // The columns along one row: the u resolution measured at that row.
    let columns_at = |v: f64| refine_line(surface, (low.x, high.x), |u| (u, v), deflection, tol);
    // A cone's circumference runs from nothing at its apex to its widest
    // at the other end, and Delaunay in a chart scaled by the mean bridges
    // the wide end's columns three at a time between rows held to six
    // cells: its rows are held to a cell and a half.
    let aspect = if matches!(surface.kind(), ogeom_geom::SurfaceKind::Cone) {
        CONE_CELL_ASPECT
    } else {
        CELL_ASPECT
    };
    // A surface straight along `v` (a drum, a cone, an extrusion) has its
    // columns held as constraints below instead: no triangle crosses one,
    // so none bridges columns, and rows are wanted only where it sags.
    let ruled = matches!(
        surface.kind(),
        ogeom_geom::SurfaceKind::Cylinder
            | ogeom_geom::SurfaceKind::Cone
            | ogeom_geom::SurfaceKind::Extrusion
    );
    // A drum's or cone's rows are its two ends where its rims allow it
    // (see `rim_holds`), its rulings drawn between them as constraints
    // below: a row inside adds a point per column and asks nothing of the
    // chord. Not where one end closes on an apex: its
    // columns all meet there, and Delaunay fans across the face from the
    // apex's few chart points.
    let apex = || {
        let um = f64::midpoint(low.x, high.x);
        let across = |v: f64| span((low.x, v), (um, v)) + span((um, v), (high.x, v));
        let (bottom, top) = (across(low.y), across(high.y));
        bottom.min(top) <= bottom.max(top) * 1e-3
    };
    let round = matches!(
        surface.kind(),
        ogeom_geom::SurfaceKind::Cylinder | ogeom_geom::SurfaceKind::Cone
    );
    let rows = if round && rows_kept == Rows::Ends && rows.len() > 2 && !apex() {
        vec![rows[0], rows[rows.len() - 1]]
    } else if ruled {
        rows
    } else {
        spread_to_aspect(
            rows,
            aspect,
            low.x,
            high.x,
            |v| columns_at(v).0.len(),
            |i| {
                row_samples[i]
                    .iter()
                    .zip(&row_samples[i + 1])
                    .map(|(a, b)| a.distance(b))
                    .fold(0.0_f64, f64::max)
            },
            |a, b, v| span((a, v), (b, v)),
        )
    };

    if *MESH_DEBUG_REFINE {
        let v = f64::midpoint(low.y, high.y);
        let (columns, _) = columns_at(v);
        eprintln!(
            "GRID u [{:.3},{:.3}] v [{:.3},{:.3}]: {} rows after aspect, {} columns at the middle row, chord {}",
            low.x,
            high.x,
            low.y,
            high.y,
            rows.len(),
            columns.len(),
            deflection.chord
        );
    }
    // Each row and each column is held to the chord along its own line,
    // and a cell's middle stands off by both sags together: on a surface
    // curved both ways, a cell whose sides each sag the chord stands off
    // twice it across its diagonal. Where any cell's diagonal sags past
    // the chord, every cell gets its centre, which halves the diagonals a
    // triangle can take, and the grid's lines between neighbouring points
    // are held as constraints, so no triangle crosses a cell's side and
    // each cell is drawn alone: four triangles round its centre. A grid
    // line inside the chart then holds the triangles as the seam does, and
    // a symmetric grid stays symmetric. A surface straight one way never
    // asks.
    let lifted = |u: f64, v: f64| {
        use ogeom_geom::Surface as _;
        surface.point_at(u, v, tol).ok()
    };
    let cells: Vec<(f64, f64, f64, f64)> = rows
        .windows(2)
        .flat_map(|pair| {
            let (v0, v1) = (pair[0], pair[1]);
            let (columns, _) = columns_at(f64::midpoint(v0, v1));
            columns
                .windows(2)
                .map(|cell| (cell[0], cell[1], v0, v1))
                .collect::<Vec<_>>()
        })
        .collect();
    // A surface straight one way sags only across its rulings, which the
    // columns already hold.
    let sags = !ruled
        && cells.len() >= MIN_CENTRED_CELLS
        && cells.iter().any(|&(u0, u1, v0, v1)| {
            let (um, vm) = (f64::midpoint(u0, u1), f64::midpoint(v0, v1));
            match (
                lifted(u0, v0),
                lifted(u1, v1),
                lifted(u1, v0),
                lifted(u0, v1),
                lifted(um, vm),
            ) {
                (Some(a), Some(b), Some(c), Some(d), Some(m)) => {
                    chord_sag(a, b, m, tol).max(chord_sag(c, d, m, tol)) > deflection.chord
                }
                _ => false,
            }
        });
    let round_ruled = matches!(
        surface.kind(),
        ogeom_geom::SurfaceKind::Cylinder | ogeom_geom::SurfaceKind::Cone
    );
    let mut placed: Vec<Vec<(usize, f64, FixedVertexHandle)>> = Vec::new();
    let mut row_of: Vec<usize> = Vec::new();
    for (row, &v) in rows
        .iter()
        .enumerate()
        .take(rows.len().saturating_sub(1))
        .skip(1)
    {
        let mut here: Vec<(usize, f64, FixedVertexHandle)> = Vec::new();
        // The row gap either side of this row, the finer of the two: the
        // chart scale a keep-out band is measured against along `v`.
        let dv = (v - rows[row - 1]).abs().min((rows[row + 1] - v).abs());
        // Each row gets its own u resolution, measured at that row.
        let (columns, column_samples) = columns_at(v);
        // The same the other way round: a surface straight along `u` gets
        // two columns from sag, and a row a hundred millimetres wide would
        // bridge across the rows as badly as the bore bridged its columns.
        let columns = spread_to_aspect(
            columns,
            CELL_ASPECT,
            low.y,
            high.y,
            |_| rows.len(),
            |i| column_samples[i].distance(&column_samples[i + 1]),
            |a, b, u| span((u, a), (u, b)),
        );
        for (column, &u) in columns
            .iter()
            .enumerate()
            .take(columns.len().saturating_sub(1))
            .skip(1)
        {
            if round_ruled && rows_kept == Rows::Sparse && column % 4 != 0 {
                continue;
            }
            // Interior points only: the boundary is already constrained, and a
            // point landing just off a constraint would split it.
            if !region.holds(Point2::new(u, v)) {
                continue;
            }
            // Inside, and not *on* the boundary: a grid point can fall
            // exactly on a ring segment that runs diagonally across the
            // chart (the midpoint of two grid corners the ring happens to
            // join), and even-odd counting calls it inside. Inserted, it
            // splits that constraint on this face alone, and the face
            // across the edge is drawn to the unsplit segment.
            let point = scale.to(u, v);
            if lands_on_the_mesh(cdt, point, reach, (1.0, 1.0)) {
                continue;
            }
            // Nor *near* it, measured in cells. A grid point a sliver's
            // width from a boundary chord makes a triangle with that
            // chord's two ends that is thin in the chart and, lifted, is
            // not thin at all: the chord cuts across the curvature by its
            // sag and the point sits on the surface, so the triangle
            // stands off the surface as a fin whose normal is tangent to
            // it and whose sign is whichever way the sliver leaned. A face
            // shades with a crease along every such chord. The point is
            // left out and the boundary's own row of triangles reaches
            // to the next grid line instead.
            let du = (u - columns[column - 1])
                .abs()
                .min((columns[column + 1] - u).abs());
            if du > 0.0
                && dv > 0.0
                && lands_on_the_mesh(
                    cdt,
                    point,
                    KEEP_OUT,
                    (1.0 / (du * scale.su), 1.0 / (dv * scale.sv)),
                )
            {
                continue;
            }
            let handle = cdt
                .insert(mitigate_underflow(point))
                .map_err(|e| ogeom_core::ogeom_err!(NotDone, "interior insertion failed: {e}"))?;
            if sags {
                here.push((column, u, handle));
            }
        }
        if sags {
            // Along the row, between neighbouring columns.
            for pair in here.windows(2) {
                if pair[1].0 == pair[0].0 + 1 && cdt.can_add_constraint(pair[0].2, pair[1].2) {
                    cdt.add_constraint(pair[0].2, pair[1].2);
                }
            }
            // Across to the row before, at the same column.
            if let (Some(before), Some(&last)) = (placed.last(), row_of.last())
                && last + 1 == row
            {
                for &(_, u, handle) in &here {
                    if let Some(&(_, _, above)) =
                        before.iter().find(|p| p.1.to_bits() == u.to_bits())
                        && cdt.can_add_constraint(above, handle)
                    {
                        cdt.add_constraint(above, handle);
                    }
                }
            }
            // To the boundary where the grid line meets it at a vertex: a
            // row's ends on the chart's end columns, the first and last
            // rows' points on the end rows.
            let edge = |cdt: &ConstrainedDelaunayTriangulation<SpadePoint<f64>>, u: f64, v: f64| {
                cdt.locate_vertex(mitigate_underflow(scale.to(u, v)))
                    .map(|h| h.fix())
            };
            let mut ties: Vec<(FixedVertexHandle, FixedVertexHandle)> = Vec::new();
            if let (Some(&(first, _, start)), Some(&(last, _, end))) = (here.first(), here.last()) {
                if first == 1
                    && let Some(b) = edge(cdt, columns[0], v)
                {
                    ties.push((b, start));
                }
                if last + 2 == columns.len()
                    && let Some(b) = edge(cdt, columns[columns.len() - 1], v)
                {
                    ties.push((end, b));
                }
            }
            for &(_, u, handle) in &here {
                if row == 1
                    && let Some(b) = edge(cdt, u, rows[0])
                {
                    ties.push((b, handle));
                }
                if row + 2 == rows.len()
                    && let Some(b) = edge(cdt, u, rows[rows.len() - 1])
                {
                    ties.push((handle, b));
                }
            }
            for (a, b) in ties {
                if a != b && cdt.can_add_constraint(a, b) {
                    cdt.add_constraint(a, b);
                }
            }
            placed.push(here);
            row_of.push(row);
        }
    }
    if sags {
        for (u0, u1, v0, v1) in cells {
            let (um, vm) = (f64::midpoint(u0, u1), f64::midpoint(v0, v1));
            if !region.holds(Point2::new(um, vm)) {
                continue;
            }
            // A cell whose side collapses to a point (a pole, an apex) is a
            // fan's sliver already, and a centre there lands wherever the
            // collapse left the columns.
            let collapsed = |v: f64| match (lifted(u0, v), lifted(u1, v)) {
                (Some(a), Some(b)) => a.distance(b) <= tol.confusion(),
                _ => true,
            };
            if collapsed(v0) || collapsed(v1) {
                continue;
            }
            let point = scale.to(um, vm);
            let (du, dv) = ((u1 - u0) * 0.5, (v1 - v0) * 0.5);
            if lands_on_the_mesh(cdt, point, reach, (1.0, 1.0))
                || (du > 0.0
                    && dv > 0.0
                    && lands_on_the_mesh(
                        cdt,
                        point,
                        KEEP_OUT,
                        (1.0 / (du * scale.su), 1.0 / (dv * scale.sv)),
                    ))
            {
                continue;
            }
            cdt.insert(mitigate_underflow(point))
                .map_err(|e| ogeom_core::ogeom_err!(NotDone, "interior insertion failed: {e}"))?;
        }
    }
    if ruled {
        // Each column, where it runs inside the face, a constraint from a
        // keep-out's width inside the boundary to as far inside the other
        // end: the triangles between two columns stay between them.
        let (bottom, _) = columns_at(low.y);
        let (top, _) = columns_at(high.y);
        let columns = if bottom.len() >= top.len() {
            bottom
        } else {
            top
        };
        // The points already in, by chart position across: a column runs
        // through the grid's points on it, a rounding off each, and is
        // held through each of them rather than past it, which would leave
        // the point on the column's far side a hair from its line.
        let mut across: Vec<(f64, f64, FixedVertexHandle)> = cdt
            .vertices()
            .map(|v| (v.position().x, v.position().y, v.fix()))
            .collect();
        across.sort_by(|a, b| a.0.total_cmp(&b.0));
        for (k, &u) in columns
            .iter()
            .enumerate()
            .take(columns.len().saturating_sub(1))
            .skip(1)
        {
            let du = (u - columns[k - 1]).abs().min((columns[k + 1] - u).abs());
            let margin = KEEP_OUT * du * scale.su / scale.sv;
            // Where the column crosses the boundary, and the boundary's
            // vertex there when the column runs through one: a column drawn
            // a rounding off a vertex still ends on it.
            let snap = du * 1e-6;
            let mut crossings: Vec<(f64, Option<Point2>)> = Vec::new();
            for ring in rings {
                for (i, &p) in ring.iter().enumerate() {
                    let q = ring[(i + 1) % ring.len()];
                    if (p.x < u) != (q.x < u) {
                        let y = p.y + (u - p.x) / (q.x - p.x) * (q.y - p.y);
                        let on = [p, q].into_iter().find(|w| (w.x - u).abs() <= snap);
                        crossings.push((on.map_or(y, |w| w.y), on));
                    }
                }
            }
            crossings.sort_by(|a, b| a.0.total_cmp(&b.0));
            for pair in crossings.windows(2) {
                let (low_end, high_end) = (pair[0], pair[1]);
                if !region.holds(Point2::new(u, f64::midpoint(low_end.0, high_end.0))) {
                    continue;
                }
                // A stretch short against the columns' spacing is one Delaunay
                // draws between neighbouring columns anyway.
                if (high_end.0 - low_end.0) * scale.sv <= CELL_ASPECT * du * scale.su {
                    continue;
                }
                let pair = [low_end.0, high_end.0];
                let mut ends = Vec::with_capacity(2);
                // From the boundary's own vertex where the column meets one,
                // else from a keep-out's width inside it.
                for ((at, on), inward) in [(low_end, margin), (high_end, -margin)] {
                    let vertex = on.unwrap_or(Point2::new(u, at));
                    if let Some(h) =
                        cdt.locate_vertex(mitigate_underflow(scale.to(vertex.x, vertex.y)))
                    {
                        ends.push(h.fix());
                        continue;
                    }
                    let v = at + inward;
                    if (v - f64::midpoint(pair[0], pair[1])).abs()
                        >= (pair[1] - pair[0]).abs() * 0.5
                    {
                        break;
                    }
                    // A boundary whose pcurves reach past the surface's
                    // domain leaves no point there to end on.
                    if surface.point_at(u, v, tol).is_err() {
                        break;
                    }
                    let point = mitigate_underflow(scale.to(u, v));
                    let handle = match cdt.locate_vertex(point) {
                        Some(h) => h.fix(),
                        None if !lands_on_the_mesh(cdt, point, reach, (1.0, 1.0)) => {
                            cdt.insert(point).map_err(|e| {
                                ogeom_core::ogeom_err!(NotDone, "column insertion failed: {e}")
                            })?
                        }
                        None => break,
                    };
                    ends.push(handle);
                }
                let [a, b] = ends[..] else {
                    continue;
                };
                if a == b {
                    continue;
                }
                let (ya, yb) = (cdt.vertex(a).position().y, cdt.vertex(b).position().y);
                let (x, near) = (u * scale.su, 1e-3 * du * scale.su);
                let first = across.partition_point(|p| p.0 < x - near);
                let mut through: Vec<(f64, FixedVertexHandle)> = across[first..]
                    .iter()
                    .take_while(|p| p.0 <= x + near)
                    .filter(|p| p.1 > ya.min(yb) && p.1 < ya.max(yb))
                    .map(|p| (p.1, p.2))
                    .collect();
                through.sort_by(|p, q| p.0.total_cmp(&q.0));
                if ya > yb {
                    through.reverse();
                }
                let chain: Vec<FixedVertexHandle> = core::iter::once(a)
                    .chain(through.into_iter().map(|p| p.1))
                    .chain(core::iter::once(b))
                    .collect();
                for link in chain.windows(2) {
                    if link[0] != link[1] && cdt.can_add_constraint(link[0], link[1]) {
                        cdt.add_constraint(link[0], link[1]);
                    }
                }
            }
        }
    }
    Ok(())
}

/// Whether segments `a..b` and `c..d` cross properly: at a point interior
/// to both, neither touching the other's end.
pub(crate) fn segments_cross(a: Point2, b: Point2, c: Point2, d: Point2) -> bool {
    let orient =
        |p: Point2, q: Point2, r: Point2| (q.x - p.x) * (r.y - p.y) - (q.y - p.y) * (r.x - p.x);
    let (o1, o2) = (orient(a, b, c), orient(a, b, d));
    let (o3, o4) = (orient(c, d, a), orient(c, d, b));
    o1 != 0.0
        && o2 != 0.0
        && o3 != 0.0
        && o4 != 0.0
        && (o1 > 0.0) != (o2 > 0.0)
        && (o3 > 0.0) != (o4 > 0.0)
}

/// How many of the triangulation's faces lie inside its constraints, told
/// by parity rather than by geometry.
///
/// Walking from a face on the convex hull, which is outside, every
/// constraint edge crossed flips inside for outside. Asked of the
/// triangulation of a boundary and nothing else, this is exact where the
/// even-odd test of a triangle's centre is not: a sliver face triangulates
/// to hairs whose centres sit on the boundary to the last bit, and which
/// side rounding puts them is a coin toss. `None` when the walk reaches a
/// face both ways with different answers: the constraints do not enclose
/// consistently, which is a crossing by another name.
fn inside_by_parity(cdt: &ConstrainedDelaunayTriangulation<SpadePoint<f64>>) -> Option<usize> {
    use ogeom_core::FastMap;
    let mut parity: FastMap<spade::handles::FixedFaceHandle<spade::handles::InnerTag>, bool> =
        FastMap::with_capacity_and_hasher(cdt.num_inner_faces(), Default::default());
    let mut queue = Vec::new();
    for hull in cdt.convex_hull() {
        // The hull edge's far side is the outer face; its near side is a
        // face of the triangulation, outside unless the hull edge itself
        // is a boundary.
        let Some(face) = hull.rev().face().as_inner() else {
            continue;
        };
        let inside = hull.is_constraint_edge();
        match parity.get(&face.fix()) {
            Some(&known) if known != inside => return None,
            Some(_) => {}
            None => {
                parity.insert(face.fix(), inside);
                queue.push((face.fix(), inside));
            }
        }
    }
    while let Some((face, inside)) = queue.pop() {
        for edge in cdt.face(face).adjacent_edges() {
            let Some(next) = edge.rev().face().as_inner() else {
                continue;
            };
            let next_inside = inside != edge.is_constraint_edge();
            match parity.get(&next.fix()) {
                Some(&known) if known != next_inside => return None,
                Some(_) => {}
                None => {
                    parity.insert(next.fix(), next_inside);
                    queue.push((next.fix(), next_inside));
                }
            }
        }
    }
    Some(parity.values().filter(|&&inside| inside).count())
}

/// Whether `point` sits within `reach` of a vertex the triangulation has,
/// or of one of its constraint edges.
///
/// Asked of whatever the point lands on (a vertex, an edge's two ends, a
/// face's three corners and whichever of its sides are constraints),
/// which is where anything that close must be. A point on a vertex is not
/// a new point; a point on a constraint would split it, and a boundary
/// split on one face only is a crack against the face across it.
fn lands_on_the_mesh(
    cdt: &ConstrainedDelaunayTriangulation<SpadePoint<f64>>,
    point: SpadePoint<f64>,
    reach: f64,
    scale: (f64, f64),
) -> bool {
    use spade::PositionInTriangulation as At;
    // Distances in a chart scaled per axis: `scale` is one over the local
    // grid step each way, so `reach` reads in cells, whatever the chart's
    // own units: one face's `u` runs over a fiftieth of its `v`.
    let scaled = |p: SpadePoint<f64>| ((p.x - point.x) * scale.0, (p.y - point.y) * scale.1);
    let near = |v: SpadePoint<f64>| {
        let (x, y) = scaled(v);
        x.hypot(y) <= reach
    };
    let along = |a: SpadePoint<f64>, b: SpadePoint<f64>| {
        // Distance to the segment `a..b`, the point at the origin.
        let (ax, ay) = scaled(a);
        let (bx, by) = scaled(b);
        let (dx, dy) = (bx - ax, by - ay);
        let len2 = dx * dx + dy * dy;
        let t = if len2 > 0.0 {
            ((-ax * dx - ay * dy) / len2).clamp(0.0, 1.0)
        } else {
            0.0
        };
        (ax + t * dx).hypot(ay + t * dy) <= reach
    };
    match cdt.locate(point) {
        At::OnVertex(_) => true,
        At::OnEdge(edge) => {
            let edge = cdt.directed_edge(edge);
            edge.is_constraint_edge() || edge.vertices().iter().any(|v| near(v.position()))
        }
        At::OnFace(face) => {
            let face = cdt.face(face);
            face.vertices().iter().any(|v| near(v.position()))
                || face.adjacent_edges().iter().any(|e| {
                    e.is_constraint_edge() && along(e.from().position(), e.to().position())
                })
        }
        At::OutsideOfConvexHull(_) | At::NoTriangulation => false,
    }
}

/// How close to the boundary, in grid cells, an interior point may sit.
///
/// Closer than this the triangle between the point and a boundary chord
/// is a sliver in the chart and a fin in space; at this and beyond the
/// boundary's own row of triangles is at least this tall against the
/// chord, and lifts as a facet on the surface rather than off it.
const KEEP_OUT: f64 = 0.35;

/// How many column widths a grid cell may be tall before rows are added.
///
/// Delaunay in the chart connects nearest neighbours in the chart; held to
/// this aspect, a cell's nearest neighbours across a row gap are the same
/// columns, not columns several away, and the triangles between rows stay
/// as narrow as the columns are.
const CELL_ASPECT: f64 = 6.0;

/// [`CELL_ASPECT`] for a cone's rows: see where it is used.
const CONE_CELL_ASPECT: f64 = 1.5;

/// Grid lines in one direction spread so that no cell is longer, in
/// space, than `aspect` times its width the other way: extra
/// lines added evenly between the ones sag chose.
///
/// Written for rows against columns and used both ways round. `lines`
/// are the parameters sag chose in this direction; `lo..hi` is the chart
/// range the other way; `crossings_at(t)` counts the lines the other way
/// at parameter `t` of this one; `length(i)` is the extent in space
/// between lines `i` and `i + 1` of this direction; `width(a, b, t)` is
/// the extent in space between two parameters of the other direction,
/// along this one's line `t`.
fn spread_to_aspect(
    lines: Vec<f64>,
    aspect: f64,
    lo: f64,
    hi: f64,
    crossings_at: impl Fn(f64) -> usize,
    length: impl Fn(usize) -> f64,
    width: impl Fn(f64, f64, f64) -> f64,
) -> Vec<f64> {
    if lines.len() < 2 {
        return lines;
    }
    let mid = f64::midpoint(lines[0], lines[lines.len() - 1]);
    let crossings = crossings_at(mid);
    if crossings < 4 {
        // Flat the other way as well: a plane in all but name, and nothing
        // to hold an aspect against.
        return lines;
    }
    // One cell's width, not the whole range's: across a closed direction
    // the whole range comes back to its own start and measures nothing.
    // Read at the middle line, and never narrower than the mean of the
    // widths a quarter, half and three quarters of the way along: the
    // triangulation is Delaunay in a chart scaled by the mean tangent
    // lengths, so the mean is what it bridges against. The middle line
    // alone can be a parallel shrinking to a point, such as the inner
    // equator of a torus whose hole nearly closes, and would spread every
    // line to that parallel's width.
    #[allow(clippy::cast_precision_loss)]
    let step = (hi - lo) / (crossings - 1) as f64;
    let (first, last) = (lines[0], lines[lines.len() - 1]);
    let quarter = (last - first) * 0.25;
    let widths = [first + quarter, mid, last - quarter].map(|t| width(lo, lo + step, t));
    let cell = widths[1].max(widths.iter().sum::<f64>() / 3.0);
    if cell.partial_cmp(&0.0) != Some(std::cmp::Ordering::Greater) {
        return lines;
    }
    let mut out = Vec::with_capacity(lines.len());
    for (i, pair) in lines.windows(2).enumerate() {
        let (a, b) = (pair[0], pair[1]);
        out.push(a);
        let tall = length(i);
        let pieces = (tall / (aspect * cell)).ceil();
        if pieces.is_finite() && pieces > 1.0 {
            #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
            let n = (pieces as usize).min(MAX_DIRECTION_STEPS);
            #[allow(clippy::cast_precision_loss)]
            for i in 1..n {
                out.push(a + (b - a) * i as f64 / n as f64);
            }
        }
    }
    out.push(lines[lines.len() - 1]);
    out
}

/// How many places across the domain the v resolution is measured at.
const U_PROBES: usize = 8;

/// How many rounds of sag-driven refinement a region may take.
///
/// Each round halves the worst sag roughly; a surface not honest after this
/// many is degenerate, and the cap makes that a coarse mesh rather than an
/// exhausted allocator.
const REFINEMENT_ROUNDS: usize = 12;

/// How many chords a chart edge's lift may stand off the line through its
/// ends before the repair splits the triangles along it.
const REPAIR_SAG: f64 = 3.0;

/// The fewest grid cells a face has for its cells to be given centres.
///
/// A face of fewer (a fillet's sliver, a small corner ball) keeps its grid's
/// triangles as drawn, and its middle may stand off by up to twice the
/// chord: such a face is drawn mostly by the rims it shares, and its
/// neighbours, and a mesh read back into faces, are drawn to agree with the
/// grid's own triangles there.
const MIN_CENTRED_CELLS: usize = 128;

/// The fewest subdivisions one parameter direction may take before its cap
/// stops it: the caller's [`Deflection::max_segments`] where that is more.
///
/// A surface that has not converged by then has a singularity, not a
/// resolution problem, and the ceiling is what makes that a coarse mesh
/// rather than an exhausted allocator.
const MAX_DIRECTION_STEPS: usize = 512;

/// How many subdivisions one parameter direction may take at `deflection`.
fn direction_steps(deflection: Deflection) -> usize {
    deflection.max_segments.max(MAX_DIRECTION_STEPS)
}

/// Subdivide `[lo, hi]` until no sub-interval errs further than `chord`,
/// or there are `cap` values.
///
/// The same adaptive bisection [`discretize`] uses on a curve, applied to a
/// line through parameter space. Each value is sampled once, by `sample`;
/// an interval is judged by `error` from its two ends' samples and what
/// `middle` gives at its midpoint, and a split keeps that middle, finished
/// by `promote`.
///
/// An interval's verdict depends on its own two ends alone, so the order
/// splits are made in changes nothing until the cap: the values come out
/// the same however the work is ordered. The splits go in rounds, every
/// interval that errs split at once, so a direction the cap stops is
/// refined evenly along its whole length; a round the cap cuts short
/// splits the intervals erring most.
fn refine_direction<S, M>(
    (lo, hi): (f64, f64),
    chord: f64,
    cap: usize,
    sample: impl Fn(f64) -> S,
    middle: impl Fn(f64) -> M,
    error: impl Fn(&S, &S, &M) -> f64,
    promote: impl Fn(f64, M) -> S,
) -> (Vec<f64>, Vec<S>) {
    let centre = f64::midpoint(lo, hi);
    let mut values = vec![lo, centre, hi];
    let mut samples: Vec<S> = vec![sample(lo), sample(centre), sample(hi)];
    // Whether the interval starting at each value is known to be within
    // the chord: its ends never move, so it stays known.
    let mut settled = vec![false; values.len() - 1];
    loop {
        // Each interval still open, judged once: (index, error, middle).
        let mut erring: Vec<(usize, f64, f64, M)> = Vec::new();
        for i in 0..settled.len() {
            if settled[i] {
                continue;
            }
            let mid = f64::midpoint(values[i], values[i + 1]);
            let at_mid = middle(mid);
            let err = error(&samples[i], &samples[i + 1], &at_mid);
            // A split that does not divide the interval means the
            // parameters have reached the resolution of f64.
            if err <= chord || mid <= values[i] || mid >= values[i + 1] {
                settled[i] = true;
                continue;
            }
            erring.push((i, err, mid, at_mid));
        }
        if erring.is_empty() {
            break;
        }
        let room = cap.saturating_sub(values.len());
        if room == 0 {
            break;
        }
        if erring.len() > room {
            erring.sort_by(|x, y| y.1.total_cmp(&x.1).then(x.0.cmp(&y.0)));
            erring.truncate(room);
            erring.sort_by_key(|e| e.0);
        }
        let mut grown_values = Vec::with_capacity(values.len() + erring.len());
        let mut grown_samples = Vec::with_capacity(values.len() + erring.len());
        let mut grown_settled = Vec::with_capacity(settled.len() + erring.len());
        let mut splits = erring.into_iter().peekable();
        for (i, (value, sampled)) in values.into_iter().zip(samples).enumerate() {
            grown_values.push(value);
            grown_samples.push(sampled);
            if i == settled.len() {
                continue;
            }
            match splits.next_if(|e| e.0 == i) {
                Some((_, _, mid, at_mid)) => {
                    grown_settled.extend([false, false]);
                    grown_values.push(mid);
                    grown_samples.push(promote(mid, at_mid));
                }
                None => grown_settled.push(settled[i]),
            }
        }
        values = grown_values;
        samples = grown_samples;
        settled = grown_settled;
    }
    (values, samples)
}

/// What a grid line knows of the surface at one of its parameters: the
/// point and the normal there, where the surface gives them.
#[derive(Clone, Copy)]
struct GridSample {
    point: Option<Point>,
    normal: Option<Direction>,
}

impl GridSample {
    fn at(surface: &SurfaceGeometry, uv: (f64, f64), tol: Tolerances) -> Self {
        Self {
            point: Self::point_at(surface, uv, tol),
            normal: Self::normal_at(surface, uv, tol),
        }
    }

    fn point_at(surface: &SurfaceGeometry, (u, v): (f64, f64), tol: Tolerances) -> Option<Point> {
        surface.point_at(u, v, tol).ok()
    }

    fn normal_at(
        surface: &SurfaceGeometry,
        (u, v): (f64, f64),
        tol: Tolerances,
    ) -> Option<Direction> {
        surface.normal_at(u, v, tol).ok()
    }

    /// How far apart two samples are in space; nothing where either has no
    /// point.
    fn distance(&self, other: &Self) -> f64 {
        match (self.point, other.point) {
            (Some(a), Some(b)) => a.distance(b),
            _ => 0.0,
        }
    }
}

/// A grid line through `range` of one chart direction, refined until no
/// cell edge along it is further from honest than the chord (see
/// [`cell_error`]), with the line's samples at every value. `at(t)` is the
/// chart point at `t` along the line.
fn refine_line(
    surface: &SurfaceGeometry,
    range: (f64, f64),
    at: impl Fn(f64) -> (f64, f64),
    deflection: Deflection,
    tol: Tolerances,
) -> (Vec<f64>, Vec<GridSample>) {
    refine_direction(
        range,
        deflection.chord,
        direction_steps(deflection),
        |t| GridSample::at(surface, at(t), tol),
        |t| GridSample::point_at(surface, at(t), tol),
        |a, b, m| cell_error(a, b, *m, deflection, tol),
        |t, point| GridSample {
            point,
            normal: GridSample::normal_at(surface, at(t), tol),
        },
    )
}

/// Fractions along the straight chart run from `from` to `to`, from 0 to
/// 1: eight even steps, each halved until no step along the run is
/// further from honest than the chord (see [`cell_error`]), or the run
/// holds as many steps as one grid direction may.
fn joining_run(
    surface: &SurfaceGeometry,
    from: Point2,
    to: Point2,
    deflection: Deflection,
    tol: Tolerances,
) -> Vec<f64> {
    const EVEN_STEPS: u32 = 8;
    let at = |f: f64| (from.x + (to.x - from.x) * f, from.y + (to.y - from.y) * f);
    let mut fractions: Vec<f64> = (0..=EVEN_STEPS)
        .map(|k| f64::from(k) / f64::from(EVEN_STEPS))
        .collect();
    let mut samples: Vec<GridSample> = fractions
        .iter()
        .map(|&f| GridSample::at(surface, at(f), tol))
        .collect();
    let cap = direction_steps(deflection);
    loop {
        let mut grown = Vec::with_capacity(fractions.len() * 2);
        let mut grown_samples = Vec::with_capacity(fractions.len() * 2);
        let mut split = false;
        for i in 0..fractions.len() {
            grown.push(fractions[i]);
            grown_samples.push(samples[i]);
            let Some(&next) = fractions.get(i + 1) else {
                continue;
            };
            let mid = f64::midpoint(fractions[i], next);
            if mid <= fractions[i] || mid >= next {
                continue;
            }
            let middle = GridSample::at(surface, at(mid), tol);
            if cell_error(&samples[i], &samples[i + 1], middle.point, deflection, tol)
                > deflection.chord
            {
                grown.push(mid);
                grown_samples.push(middle);
                split = true;
            }
        }
        if !split || grown.len() > cap {
            break;
        }
        fractions = grown;
        samples = grown_samples;
    }
    fractions
}

/// As [`refine_line`] over several lines at once, sharing their values:
/// `lines` are the lines by their parameter in the other direction,
/// `at(line, t)` the chart point at `t` along one of them, and an
/// interval's error the worst across them.
fn refine_lines(
    surface: &SurfaceGeometry,
    range: (f64, f64),
    lines: &[f64],
    at: impl Fn(f64, f64) -> (f64, f64),
    deflection: Deflection,
    tol: Tolerances,
) -> (Vec<f64>, Vec<Vec<GridSample>>) {
    refine_direction(
        range,
        deflection.chord,
        direction_steps(deflection),
        |t| {
            lines
                .iter()
                .map(|&line| GridSample::at(surface, at(line, t), tol))
                .collect::<Vec<_>>()
        },
        |t| {
            lines
                .iter()
                .map(|&line| GridSample::point_at(surface, at(line, t), tol))
                .collect::<Vec<_>>()
        },
        |a, b, m| {
            a.iter()
                .zip(b)
                .zip(m)
                .map(|((a, b), m)| cell_error(a, b, *m, deflection, tol))
                .fold(0.0_f64, f64::max)
        },
        |t, points| {
            lines
                .iter()
                .zip(points)
                .map(|(&line, point)| GridSample {
                    point,
                    normal: GridSample::normal_at(surface, at(line, t), tol),
                })
                .collect()
        },
    )
}

/// How far a grid cell's edge is from honest, as a sag: the chord sag
/// itself, or the normal's turn across it scaled so that a turn of the
/// angular deflection weighs the same as a sag of the chord, whichever
/// is worse. From the edge's two ends and the surface's point at its
/// middle; the sag is nothing where any of the three is missing.
///
/// The chord alone is what the boundary's edges are *not* drawn to: a
/// curve is discretized to both deflections, so a bore's rims come out
/// round at the angular limit while columns held to the chord alone come
/// out a polygon of far fewer sides, and the bore changes shape a chord's
/// length in from each rim. The interior is held to the same two limits
/// the boundary is.
fn cell_error(
    from: &GridSample,
    to: &GridSample,
    middle: Option<Point>,
    deflection: Deflection,
    tol: Tolerances,
) -> f64 {
    let sag = match (from.point, to.point, middle) {
        (Some(a), Some(b), Some(m)) => chord_sag(a, b, m, tol),
        // Off the surface's domain; nothing to refine towards.
        _ => 0.0,
    };
    let turn = match (from.normal, to.normal) {
        (Some(a), Some(b)) => a.angle(b),
        // A pole or an apex has no normal to compare; the sag still governs.
        _ => 0.0,
    };
    sag.max(turn / deflection.angular * deflection.chord)
}

/// How far `m` stands off the chord from `a` to `b`.
///
/// Measured in space, which is the only place the number means anything: the
/// same step in `u` covers a metre at a sphere's equator and a millimetre near
/// its pole.
fn chord_sag(a: Point, b: Point, m: Point, tol: Tolerances) -> f64 {
    ogeom_math::Axis::through(a, b, tol).map_or_else(|_| a.distance(m), |axis| axis.distance_to(m))
}

/// Whether a point is inside the region the rings bound.
#[cfg(test)]
fn inside_region(rings: &[Vec<Point2>], p: Point2) -> bool {
    inside_boundary_with::<Exact>(rings, p)
}

/// Ray-crossing count for one ring, decided by orientation rather than by
/// arithmetic.
///
/// A horizontal ray in `+u`. The half-open `y` comparison counts a vertex lying
/// exactly on the ray once rather than twice or not at all; which side of the
/// edge the point falls on is then an `orient2d` sign.
///
/// Deliberately *not* "solve for where the edge crosses the ray, then compare".
/// That form divides by the edge's `y` extent, which is near zero for a nearly
/// horizontal edge, and subtracts two nearly equal numbers to compare, so for a
/// point close to the boundary it can answer either way. Here there is no
/// division and the comparison is a determinant's sign, which an exact predicate
/// gets right at any separation.
fn crosses_odd_times<P: Predicates>(ring: &[Point2], p: Point2) -> bool {
    let mut inside = false;
    let n = ring.len();
    let at = |q: Point2| [q.x, q.y];
    for i in 0..n {
        let (a, b) = (ring[i], ring[(i + 1) % n]);
        if (a.y > p.y) == (b.y > p.y) {
            continue;
        }
        // The edge crosses the ray's line. Whether it crosses the ray *itself*
        // (to the right of `p`) is which side of the directed edge `p` is on,
        // read the right way round for the edge's direction in `y`.
        let side = P::orient2d(at(a), at(b), [p.x, p.y]);
        let rightwards = if b.y > a.y {
            side == ogeom_core::Sign::Positive
        } else {
            side == ogeom_core::Sign::Negative
        };
        if rightwards {
            inside = !inside;
        }
    }
    inside
}

/// Whether to report, per shape, how many faces were drawn again finer.
///
/// Read once: `env::var` takes a process-wide lock and allocates its answer,
/// and this is asked per face.
static MESH_DEBUG_REFINE: std::sync::LazyLock<bool> =
    std::sync::LazyLock::new(|| std::env::var("OGEOM_MESH_DEBUG_REFINE").is_ok());

/// Whether the mesh debug dump is on, read once for the same reason.
static MESH_DEBUG: std::sync::LazyLock<bool> =
    std::sync::LazyLock::new(|| std::env::var("OGEOM_MESH_DEBUG").is_ok());

/// The ring edges that can cross a horizontal ray, bucketed by height.
///
/// [`crosses_odd_times`] walks every edge of every ring for each point it is
/// asked about. That is fine for a handful of queries and ruinous for the
/// refinement loop, which asks once per triangle per round while the rings
/// themselves never change: a face with 544 boundary points and 41 000
/// triangles pays a quarter of a billion edge visits, nearly all on edges
/// nowhere near the point.
///
/// An edge can only straddle a ray at height `y` if `y` lies within the edge's
/// own `y` span, so bucketing edges by that span and querying one bucket tests
/// a conservative superset of the edges that could contribute. **The answer is
/// therefore identical**: the same straddle test and the same exact predicate
/// decide each candidate; the index only declines to visit edges that could
/// not have counted.
///
/// Parity is taken over all rings at once, which is what
/// [`inside_boundary_with`] computes as an exclusive-or of per-ring parities:
/// the two agree because the parity of the total crossing count is the sum of
/// the rings' parities.
struct RingBands {
    /// Every ring's edges, flattened.
    edges: Vec<(Point2, Point2)>,
    /// Edge indices per band, low `y` first.
    bands: Vec<Vec<u32>>,
    low: f64,
    high: f64,
    /// Band height. Zero when every point shares one `y`, which leaves a
    /// single band holding everything.
    step: f64,
}

impl RingBands {
    /// Index the rings. Cheap enough to build per face and paid back by the
    /// first few hundred queries.
    fn over(rings: &[Vec<Point2>]) -> Self {
        let mut edges = Vec::new();
        for ring in rings {
            for i in 0..ring.len() {
                edges.push((ring[i], ring[(i + 1) % ring.len()]));
            }
        }
        let (mut low, mut high) = (f64::INFINITY, f64::NEG_INFINITY);
        for (a, b) in &edges {
            low = low.min(a.y).min(b.y);
            high = high.max(a.y).max(b.y);
        }
        if edges.is_empty() || !low.is_finite() || !high.is_finite() {
            return Self {
                edges,
                bands: Vec::new(),
                low: 0.0,
                high: 0.0,
                step: 0.0,
            };
        }
        // About four edges to a band: enough to keep the per-query walk short
        // without spreading a long edge across a table of mostly empty bands.
        let count = edges.len().div_ceil(4).clamp(1, 4096);
        #[allow(
            clippy::cast_precision_loss,
            reason = "a band count, far below the integers f64 represents exactly"
        )]
        let step = (high - low) / count as f64;
        let mut bands: Vec<Vec<u32>> = vec![Vec::new(); count];
        for (i, (a, b)) in edges.iter().enumerate() {
            let (lo, hi) = (a.y.min(b.y), a.y.max(b.y));
            let first = Self::band_of(lo, low, step, count);
            let last = Self::band_of(hi, low, step, count);
            for band in &mut bands[first..=last] {
                #[allow(
                    clippy::cast_possible_truncation,
                    reason = "an edge index, bounded by the ring lengths"
                )]
                band.push(i as u32);
            }
        }
        Self {
            edges,
            bands,
            low,
            high,
            step,
        }
    }

    /// Which band a height falls in, clamped to the table.
    fn band_of(y: f64, low: f64, step: f64, count: usize) -> usize {
        if step <= 0.0 {
            return 0;
        }
        let raw = (y - low) / step;
        if raw <= 0.0 {
            return 0;
        }
        #[allow(
            clippy::cast_possible_truncation,
            clippy::cast_sign_loss,
            reason = "clamped to the band count on the next line"
        )]
        let index = raw as usize;
        index.min(count - 1)
    }

    /// Whether the point lies inside the region the rings bound.
    fn holds(&self, p: Point2) -> bool {
        if self.bands.is_empty() || p.y < self.low || p.y > self.high {
            return false;
        }
        let band = Self::band_of(p.y, self.low, self.step, self.bands.len());
        let at = |q: Point2| [q.x, q.y];
        let mut inside = false;
        for &i in &self.bands[band] {
            let (a, b) = self.edges[i as usize];
            if (a.y > p.y) == (b.y > p.y) {
                continue;
            }
            let side = Exact::orient2d(at(a), at(b), [p.x, p.y]);
            let rightwards = if b.y > a.y {
                side == ogeom_core::Sign::Positive
            } else {
                side == ogeom_core::Sign::Negative
            };
            if rightwards {
                inside = !inside;
            }
        }
        inside
    }
}

/// The unit normal of a triangle, or `None` if it is degenerate.
#[must_use]
pub fn triangle_normal(a: Point, b: Point, c: Point, tol: Tolerances) -> Option<Direction> {
    Direction::from_cross(b - a, c - a, tol).ok()
}

/// Discretize an edge into a polyline in space, for display or coarse queries.
///
/// # Errors
///
/// As [`discretize`].
pub fn polyline_of_edge(
    model: &Model,
    edge: &Shape,
    deflection: Deflection,
    tol: Tolerances,
) -> OgeomResult<Vec<Point>> {
    if model.kind_of(edge)? != ShapeType::Edge {
        ogeom_bail!(Construction, "expected an edge");
    }
    let Some(node) = model.node(edge) else {
        ogeom_bail!(Dangling, "edge is not in this model");
    };
    let NodeData::Edge(data) = node.data() else {
        ogeom_bail!(Construction, "edge node holds no edge data");
    };
    let Some(line) = edge_polyline(model, data, deflection, tol)? else {
        return Ok(Vec::new());
    };
    let placement = edge.transform(model.datums())?;
    let mut points: Vec<Point> = line.points.iter().map(|p| placement.apply(*p)).collect();
    if edge.orientation() == Orientation::Reversed {
        points.reverse();
    }
    Ok(points)
}

#[cfg(test)]
#[allow(clippy::unwrap_used)]
mod tests {
    use super::*;

    /// Two pieces wound alike about +z share one side, walked once each
    /// way, and each closes a triangle at the same two points, a side of
    /// which runs from a point to itself. Those sides are walked no way, so
    /// they cast no vote, and neither piece is flipped.
    #[test]
    fn a_side_from_a_point_to_itself_turns_no_piece() {
        let (far, farther) = (Point::new(5.0, 5.0, 5.0), Point::new(6.0, 5.0, 5.0));
        let mut mesh = Triangulation::new();
        mesh.positions = vec![
            Point::new(0.0, 0.0, 0.0),
            Point::new(1.0, 0.0, 0.0),
            Point::new(0.0, 1.0, 0.0),
            Point::new(1.0, 1.0, 0.0),
            far,
            far,
            farther,
            farther,
            far,
            far,
            farther,
            farther,
        ];
        mesh.normals = vec![Vector::new(0.0, 0.0, 1.0); mesh.positions.len()];
        let first = [[0, 1, 2], [0, 4, 5], [0, 6, 7]];
        let second = [[1, 3, 2], [3, 8, 9], [3, 10, 11]];
        mesh.triangles = first.iter().chain(&second).copied().collect();
        let before = mesh.triangles.clone();
        orient_pieces(&mut mesh, &[(0, 3), (3, 6)]);
        assert_eq!(mesh.triangles, before);
    }

    /// A ring wrapping a closed chart once, read a whole turn low: the
    /// extent a patch's full band reads before it is slid, on a chart of
    /// [-3.14159265358979, 3.14159265358979] as a file states it. The slide
    /// is one turn up, though the turn's length over the chart's span
    /// rounds to a hair over one; a hair past an end is no turn, and a ring
    /// inside stays.
    #[test]
    #[allow(
        clippy::excessive_precision,
        clippy::approx_constant,
        reason = "the values as a file states them, a hair short of pi"
    )]
    fn a_ring_a_whole_turn_out_slides_back_one_turn() {
        let domain = (-3.141_592_653_589_79, 3.141_592_653_589_79);
        let span = domain.1 - domain.0;
        let low = (-9.424_777_960_769_370_47, -3.141_592_653_589_744_27);
        assert!((domain.0 - low.0) / span > 1.0);
        assert_eq!(turns_into_domain(low, domain), 1.0);
        let high = (low.0 + 2.0 * span, low.1 + 2.0 * span);
        assert_eq!(turns_into_domain(high, domain), -1.0);
        assert_eq!(turns_into_domain((domain.0 - 1e-9, domain.1), domain), 0.0);
        assert_eq!(turns_into_domain((-1.0, 2.0), domain), 0.0);
        // Two turns out is two turns back.
        let two = (low.0 - span, low.1 - span);
        assert_eq!(turns_into_domain(two, domain), 2.0);
    }

    /// A plane strip 0.03 wide between two concentric arcs of a 15 mm
    /// circle, its inner rim already held to a fine chord by the face
    /// across it: the outer rim, drawn at the coarse chord, sags across the
    /// inner one, so the face asks for its outer rim drawn finer. Asked
    /// from nothing, both rims drawn coarse at the same angles, it asks
    /// for nothing: the narrow strip is whole until one rim is refined.
    #[test]
    fn a_strip_with_one_rim_held_fine_asks_for_its_other_rim_finer() {
        use ogeom_geom::{CircleCurve, Curve, PlaneSurface};
        use ogeom_math::{Circle, Frame, Plane};
        let mut model = Model::new();
        let (inner, outer, angle) = (15.0, 15.03, core::f64::consts::FRAC_PI_3);
        let at = |r: f64, a: f64| Point::new(r * a.cos(), r * a.sin(), 0.0);
        let corners = [
            at(outer, 0.0),
            at(outer, angle),
            at(inner, angle),
            at(inner, 0.0),
        ];
        let vertices: Vec<Shape> = corners
            .iter()
            .map(|p| ogeom_algo::make_vertex(&mut model, *p).shape)
            .collect();
        let arc =
            |r: f64| -> Curve { CircleCurve::new(Circle::new(Frame::WORLD, r, T).unwrap()).into() };
        let line = |a: Point, b: Point| -> Curve {
            ogeom_geom::LineCurve::segment(a, b, T).unwrap().into()
        };
        let mut edge = |curve: Curve, range: (f64, f64), from: usize, to: usize| {
            ogeom_algo::make_edge_between(
                &mut model,
                curve,
                range,
                &vertices[from],
                &vertices[to],
                T,
            )
            .unwrap()
            .shape
        };
        let outer_edge = edge(arc(outer), (0.0, angle), 0, 1);
        let top = edge(
            line(corners[1], corners[2]),
            (0.0, corners[1].distance(corners[2])),
            1,
            2,
        );
        let inner_edge = edge(arc(inner), (0.0, angle), 3, 2);
        let bottom = edge(
            line(corners[3], corners[0]),
            (0.0, corners[3].distance(corners[0])),
            3,
            0,
        );
        let plane: SurfaceGeometry = PlaneSurface::new(Plane::new(Frame::WORLD)).into();
        let ring = vec![outer_edge.clone(), top, inner_edge.reversed(), bottom];
        let face = ogeom_algo::make_face_with_pcurves(&mut model, plane, &[ring], T)
            .unwrap()
            .shape;
        let data = model.node(&face).unwrap().data().as_face().unwrap().clone();
        let surface = model.geometry().surface(data.surface).unwrap().clone();
        let coarse = Deflection::default();
        let alone = face_chords(
            &model,
            &face,
            data.surface,
            &surface,
            coarse,
            &EdgeChords::new(),
            T,
        )
        .unwrap();
        assert!(!alone.contains_key(&outer_edge.node().index()), "{alone:?}");
        let held: EdgeChords = [(inner_edge.node().index(), coarse.chord / 128.0)].into();
        let asked = face_chords(&model, &face, data.surface, &surface, coarse, &held, T).unwrap();
        let outer_chord = asked.get(&outer_edge.node().index()).copied();
        assert!(outer_chord.is_some_and(|c| c < coarse.chord), "{asked:?}");
        assert!(asked[&inner_edge.node().index()] <= coarse.chord / 128.0);
    }

    #[test]
    fn a_wild_chart_refuses_instead_of_crashing() {
        let plane: SurfaceGeometry = ogeom_geom::PlaneSurface::over(
            ogeom_math::Plane::through(
                ogeom_math::Point::new(0.0, 0.0, 0.0),
                ogeom_math::Direction::Z,
            ),
            (-10.0, 10.0),
            (-10.0, 10.0),
        )
        .unwrap()
        .into();
        let tol = Tolerances::millimetres();
        // Coordinates wilder than any face: refused by the screen.
        let wild = vec![vec![
            Point2::new(0.0, 0.0),
            Point2::new(1e15, 0.0),
            Point2::new(0.0, 1e15),
        ]];
        let Err(err) = triangulate_region(&wild, &plane, Deflection::default(), Rows::Ends, tol)
        else {
            panic!("a wild chart must refuse");
        };
        assert!(err.to_string().contains("describes no face"), "{err}");
        // A non-finite coordinate: refused by name, not fed to the library.
        let nan = vec![vec![
            Point2::new(0.0, 0.0),
            Point2::new(f64::NAN, 1.0),
            Point2::new(1.0, 1.0),
        ]];
        let Err(err) = triangulate_region(&nan, &plane, Deflection::default(), Rows::Ends, tol)
        else {
            panic!("a NaN chart must refuse");
        };
        assert!(err.to_string().contains("non-finite"), "{err}");
    }

    use approx::assert_relative_eq;

    #[test]
    fn chart_tangency_retains_both_lobes_and_all_boundary_sides() {
        let tol = Tolerances::millimetres();
        let frame = ogeom_math::Frame::new(
            Point::ORIGIN,
            ogeom_math::Direction::Z,
            ogeom_math::Direction::X,
            tol,
        )
        .unwrap();
        let plane: SurfaceGeometry =
            ogeom_geom::PlaneSurface::new(ogeom_math::Plane::new(frame)).into();
        for angle in [0.0_f64, 0.1, 0.7] {
            let (sin, cos) = angle.sin_cos();
            let rotate = |p: Point2| Point2::new(cos * p.x - sin * p.y, sin * p.x + cos * p.y);
            for offset in [-1e-12, 0.0, 1e-12] {
                let mut ring: Vec<_> = [
                    Point2::new(0.0, 0.0),
                    Point2::new(0.0, 10.0),
                    Point2::new(1.0, 10.0),
                    Point2::new(offset, 3.0),
                    Point2::new(1.0, 0.0),
                ]
                .into_iter()
                .map(rotate)
                .collect();
                let original = ring.clone();
                let mut anchors = vec![None; ring.len()];
                snap_chart_junctions(&mut ring, &mut anchors, &plane, tol);
                assert_eq!(
                    ring[1], original[3],
                    "the target side shares the touching vertex"
                );
                let boundary = ring.clone();
                let mesh = triangulate_region(
                    &[ring],
                    &plane,
                    Deflection::with_chord(0.1).unwrap(),
                    Rows::Ends,
                    tol,
                )
                .unwrap();
                assert!(!mesh.crossed, "angle {angle}, offset {offset}");
                assert_eq!(mesh.triangles.len(), 2);
                for (&a, &b) in boundary.iter().zip(boundary.iter().cycle().skip(1)) {
                    assert!(
                        mesh.triangles.iter().any(|t| (0..3).any(|i| {
                            let x = mesh.parameters[t[i] as usize];
                            let y = mesh.parameters[t[(i + 1) % 3] as usize];
                            (x == (a.x, a.y) && y == (b.x, b.y))
                                || (x == (b.x, b.y) && y == (a.x, a.y))
                        })),
                        "missing boundary side {a:?} -- {b:?}"
                    );
                }
                // The two lobes have total area 5. The tolerance covers the
                // sub-parametric boundary bend and floating-point rotation.
                let area: f64 = mesh
                    .triangles
                    .iter()
                    .map(|t| {
                        let [a, b, c] = t.map(|i| {
                            let (u, v) = mesh.parameters[i as usize];
                            Point2::new(u, v)
                        });
                        (b - a).cross(c - a).abs() * 0.5
                    })
                    .sum();
                assert!((area - 5.0).abs() < 1e-9, "{area}");
                let mut model = Model::new();
                let points: Vec<_> = original.iter().map(|p| Point::new(p.x, p.y, 0.0)).collect();
                let wire = ogeom_algo::make_polygon(&mut model, &points, true, tol)
                    .unwrap()
                    .shape;
                let face = ogeom_algo::make_face(&mut model, plane.clone(), &[wire], tol)
                    .unwrap()
                    .shape;
                ogeom_algo::build::trimmed_where_bare(&mut model, &face, tol).unwrap();
                let (mesh, verdict) = triangulate_reporting(
                    &model,
                    &face,
                    Deflection::with_chord(0.1).unwrap(),
                    Some(&EdgeChords::new()),
                    tol,
                )
                .unwrap();
                assert!(verdict != Verdict::Short, "angle {angle}, offset {offset}");
                let area: f64 = mesh
                    .triangles
                    .iter()
                    .map(|t| {
                        let [a, b, c] = t.map(|i| mesh.positions[i as usize]);
                        (b - a).cross(c - a).magnitude() * 0.5
                    })
                    .sum();
                assert!((area - 5.0).abs() < 1e-9, "face area {area}");
            }
        }
        let mut separate = vec![
            Point2::new(0.0, 0.0),
            Point2::new(0.0, 10.0),
            Point2::new(1.0, 10.0),
            Point2::new(2.0 * tol.parametric(), 5.0),
            Point2::new(1.0, 0.0),
        ];
        let before = separate.clone();
        snap_chart_junctions(&mut separate, &mut vec![None; before.len()], &plane, tol);
        assert_eq!(separate, before, "a resolvable gap must stay separate");
    }
    #[test]
    fn chart_junction_candidates_stay_bounded_on_dense_seams() {
        let samples = 1024;
        for transpose in [false, true] {
            let mut ring = Vec::new();
            for side in 0..4 {
                for i in 0..samples {
                    let t = i as f64 / samples as f64;
                    let p = match side {
                        0 => Point2::new(10.0 * t, 0.0),
                        1 => Point2::new(10.0, t),
                        2 => Point2::new(10.0 * (1.0 - t), 1.0),
                        _ => Point2::new(0.0, 1.0 - t),
                    };
                    ring.push(if transpose { Point2::new(p.y, p.x) } else { p });
                }
            }
            let swapped: Vec<_> = ring.iter().map(|p| Point2::new(p.y, p.x)).collect();
            let x = RingBands::over(&[swapped]);
            let y = RingBands::over(std::slice::from_ref(&ring));
            let mut visited = 0;
            for (i, &p) in ring.iter().enumerate() {
                let candidates: Vec<_> = junction_candidates(&x, &y, p, 1e-9)
                    .iter()
                    .flatten()
                    .copied()
                    .collect();
                // Both incident sides must be found, including band boundaries.
                let current = u32::try_from(i).unwrap();
                let previous = u32::try_from((i + ring.len() - 1) % ring.len()).unwrap();
                assert!(candidates.contains(&current));
                assert!(candidates.contains(&previous));
                visited += candidates.len();
            }
            // Count work rather than wall time: one axis would visit over a
            // million entries on either the horizontal or vertical seam.
            assert!(visited < 32 * ring.len(), "visited {visited} sides");
        }
    }

    use ogeom_algo::make_box;
    use ogeom_math::Frame;
    use ogeom_topo::explore_unique;

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

    /// A cap no profile below reaches but the one written to.
    const UNCAPPED: usize = 1 << 16;

    /// `refine_direction` by rescan from zero after every split: the
    /// reference the rounds are held to.
    fn refine_by_rescan<F: Fn(f64, f64) -> f64>(lo: f64, hi: f64, chord: f64, sag: F) -> Vec<f64> {
        let mut values = vec![lo, f64::midpoint(lo, hi), hi];
        while values.len() < UNCAPPED {
            let Some(i) = (0..values.len() - 1).find(|&i| sag(values[i], values[i + 1]) > chord)
            else {
                break;
            };
            let mid = f64::midpoint(values[i], values[i + 1]);
            if mid <= values[i] || mid >= values[i + 1] {
                break;
            }
            values.insert(i + 1, mid);
        }
        values
    }

    #[test]
    fn rounds_split_where_rescanning_did() {
        // Rounds are only sound because splitting an interval cannot change
        // whether another one sags. Held to the rescan's output exactly,
        // over sag profiles that bite in different places: flat, steep at one
        // end, periodic; and one savage enough to reach the step cap, which
        // must leave it refined evenly.
        /// A named sag profile to hold both forms to.
        type Profile = (&'static str, Box<dyn Fn(f64, f64) -> f64>);
        let cases: Vec<Profile> = vec![
            ("flat", Box::new(|_a: f64, _b: f64| 0.0)),
            ("width", Box::new(|a: f64, b: f64| (b - a).abs())),
            (
                "steep at the low end",
                Box::new(|a: f64, b: f64| (b - a).abs() / a.abs().max(1e-3)),
            ),
            (
                "periodic",
                Box::new(|a: f64, b: f64| (b - a).abs() * (a * 12.0).sin().abs()),
            ),
            (
                "beyond the cap",
                Box::new(|a: f64, b: f64| (b - a).abs() * 1e6),
            ),
        ];
        for (name, sag) in cases {
            for chord in [1.0, 0.1, 0.01, 1e-3] {
                let capped = name == "beyond the cap";
                let cap = if capped {
                    MAX_DIRECTION_STEPS
                } else {
                    UNCAPPED
                };
                let (walked, samples) = refine_direction(
                    (0.0, 1.0),
                    chord,
                    cap,
                    |t| t,
                    |t| t,
                    |a, b, _| sag(*a, *b),
                    |t, _| t,
                );
                assert_eq!(walked, samples, "{name}: each value keeps its own sample");
                if capped {
                    // Stopped by the cap, and refined evenly all along.
                    assert!(walked.len() <= cap, "{name}: {} values", walked.len());
                    let gaps: Vec<f64> = walked.windows(2).map(|w| w[1] - w[0]).collect();
                    let (least, most) = gaps
                        .iter()
                        .fold((f64::INFINITY, 0.0_f64), |(l, m), g| (l.min(*g), m.max(*g)));
                    assert!(most <= least * 2.0, "{name}: gaps from {least} to {most}");
                    continue;
                }
                let rescanned = refine_by_rescan(0.0, 1.0, chord, &sag);
                assert_eq!(
                    walked, rescanned,
                    "{name} at chord {chord}: the rounds split somewhere the rescan did not"
                );
            }
        }
    }

    #[test]
    fn banded_containment_answers_what_the_full_scan_answers() {
        // The index may only decline to visit edges that could not have
        // counted. Held to the unindexed predicate over a ring with a hole,
        // on a grid that straddles both boundaries and the vertices themselves.
        let outer: Vec<Point2> = vec![
            Point2::new(0.0, 0.0),
            Point2::new(4.0, 0.0),
            Point2::new(4.0, 3.0),
            Point2::new(2.0, 1.5),
            Point2::new(0.0, 3.0),
        ];
        let hole: Vec<Point2> = vec![
            Point2::new(1.0, 0.5),
            Point2::new(1.0, 1.0),
            Point2::new(1.5, 1.0),
            Point2::new(1.5, 0.5),
        ];
        let rings = vec![outer, hole];
        let bands = RingBands::over(&rings);
        for i in 0..=80 {
            for j in 0..=60 {
                #[allow(clippy::cast_precision_loss)]
                let p = Point2::new(f64::from(i) * 0.05 - 0.1, f64::from(j) * 0.05 - 0.1);
                assert_eq!(
                    bands.holds(p),
                    inside_region(&rings, p),
                    "the index disagreed with the full scan at {p:?}"
                );
            }
        }
    }

    fn fine() -> Deflection {
        Deflection {
            chord: 1e-3,
            angular: 0.05,
            ..Deflection::default()
        }
    }

    /// A vertex at a cone's apex or a sphere's pole carries the normal the
    /// surface tends to there along its own column, not nothing.
    #[test]
    fn apex_and_pole_vertices_carry_the_limit_normal() {
        let mut model = Model::new();
        let cone = ogeom_algo::make_cone(&mut model, Frame::WORLD, 2.0, 0.0, 3.0, T).unwrap();
        let sphere = ogeom_algo::make_sphere(&mut model, Frame::WORLD, 1.5, T).unwrap();
        for (shape, what) in [(&cone.shape, "cone"), (&sphere.shape, "sphere")] {
            for face in explore_unique(&model, shape, ShapeType::Face).unwrap() {
                let mesh = triangulate_face(&model, &face, fine(), T).unwrap();
                for (i, n) in mesh.normals.iter().enumerate() {
                    assert!(
                        (n.magnitude() - 1.0).abs() < 1e-9,
                        "{what} vertex {i} at {:?} has normal {n:?}",
                        mesh.positions[i]
                    );
                }
            }
        }
        // At the apex the limit normal along a ruling is that ruling's
        // normal: on a cone of half-angle atan(2/3) it leans out by that
        // much from the axis, the same as every other normal in its column.
        let faces = explore_unique(&model, &cone.shape, ShapeType::Face).unwrap();
        let lean = (2.0_f64 / 3.0).atan();
        for face in &faces {
            let mesh = triangulate_face(&model, face, fine(), T).unwrap();
            for (i, p) in mesh.positions.iter().enumerate() {
                if p.distance(Point::new(0.0, 0.0, 3.0)) < 1e-9 {
                    let n = mesh.normals[i];
                    let from_axis = n.z.abs().acos();
                    assert!(
                        ((std::f64::consts::FRAC_PI_2 - from_axis) - lean).abs() < 1e-6,
                        "apex normal {n:?} leans {from_axis} from the axis"
                    );
                }
            }
        }
    }

    #[test]
    fn a_box_face_triangulates_into_two_triangles() {
        // A planar square needs no interior points at all, which is what makes
        // deflection-driven refinement worth having: a fixed grid would add
        // dozens that buy nothing.
        let mut model = Model::new();
        let built = make_box(&mut model, Frame::WORLD, (2.0, 3.0, 4.0), T).unwrap();
        let faces = explore_unique(&model, &built.shape, ShapeType::Face).unwrap();

        for face in &faces {
            let mesh = triangulate_face(&model, face, fine(), T).unwrap();
            assert_eq!(mesh.triangle_count(), 2, "a rectangle is two triangles");
            assert_eq!(mesh.vertex_count(), 4);
            assert!(mesh.deflection_met);
        }
    }

    #[test]
    fn a_boxs_mesh_is_closed_and_reports_the_right_volume() {
        // The end-to-end check: triangulate, weld, and ask the mesh what it
        // encloses. A volume that comes out negative would mean the faces are
        // wound inward; one that is wrong in magnitude would mean the
        // triangulation is not covering the boundary.
        let mut model = Model::new();
        let size = (2.0, 3.0, 4.0);
        let built = make_box(&mut model, Frame::WORLD, size, T).unwrap();
        let mesh = triangulate(&model, &built.shape, fine(), T).unwrap();

        assert_eq!(mesh.triangle_count(), 12, "six faces, two triangles each");
        assert_eq!(mesh.vertex_count(), 8, "welding merged the shared corners");
        assert!(
            mesh.is_closed(),
            "every triangle edge should be shared by two"
        );
        assert_relative_eq!(mesh.volume(), size.0 * size.1 * size.2, epsilon = 1e-9);
        assert_relative_eq!(
            mesh.area(),
            2.0 * (size.0 * size.1 + size.1 * size.2 + size.2 * size.0),
            epsilon = 1e-9
        );
    }

    /// The segments a full circle of radius `r` takes at chord `c`: the
    /// sag of a segment spanning `θ` is `r (1 − cos(θ/2))`, so the circle
    /// wants `π / acos(1 − c/r)` of them, and the bisection reaches it from
    /// its floor of two by halving every segment alike.
    fn circle_segments(r: f64, c: f64) -> usize {
        let wanted = std::f64::consts::PI / (1.0 - c / r).acos();
        let mut n = Deflection::default().min_segments;
        #[allow(clippy::cast_precision_loss, reason = "a segment count")]
        while (n as f64) < wanted {
            n *= 2;
        }
        n
    }

    /// No triangle is longer for want of a length cap: a box is two
    /// triangles a face whether it is a 20 mm cube or a bar a thousand
    /// millimetres long and a tenth across, and every one of them has area
    /// and faces out with its face.
    #[test]
    fn a_box_is_two_triangles_a_face_however_long_its_faces() {
        let deflection = Deflection::curvature_only(0.01, 0.5).unwrap();
        for size in [(20.0, 20.0, 20.0), (1000.0, 0.1, 0.1)] {
            let mut model = Model::new();
            let built = make_box(&mut model, Frame::WORLD, size, T).unwrap();
            let mesh = triangulate(&model, &built.shape, deflection, T).unwrap();
            assert_eq!(mesh.triangle_count(), 12, "{size:?}");
            assert!(mesh.is_closed(), "{size:?}");
            // Two triangles on each face of the box, exact to rounding.
            assert_relative_eq!(
                mesh.volume(),
                size.0 * size.1 * size.2,
                max_relative = 1e-12
            );
            for face in explore_unique(&model, &built.shape, ShapeType::Face).unwrap() {
                assert_faces_out(&triangulate_face(&model, &face, deflection, T).unwrap());
            }
        }
    }

    /// Every triangle has area and turns the way the surface's normals
    /// at its corners do.
    fn assert_faces_out(mesh: &Triangulation) {
        for t in &mesh.triangles {
            let [a, b, c] = t.map(|i| mesh.positions[i as usize]);
            let along = t
                .iter()
                .map(|&i| mesh.normals[i as usize])
                .fold(Vector::ZERO, |sum, n| sum + n);
            let normal = (b - a).cross(c - a);
            assert!(normal.magnitude() > 0.0, "a triangle with no area");
            assert!(normal.dot(along) > 0.0, "a triangle turned over");
        }
    }

    /// A cylinder of radius 10 at a chord of a hundredth: each rim takes
    /// the segments the chord asks, the side two triangles a segment and
    /// each cap the triangulation of its polygon, at a height of 20 and of
    /// 2000 alike.
    #[test]
    fn a_cylinder_takes_the_segments_its_chord_asks_at_any_height() {
        let (r, c) = (10.0, 0.01);
        let n = circle_segments(r, c);
        // π / acos(0.999) = 70.24: seventy-one segments would do, and the
        // halving from two reaches 128.
        assert_eq!(n, 128);
        let deflection = Deflection::curvature_only(c, 0.5).unwrap();
        for height in [20.0, 2000.0] {
            let mut model = Model::new();
            let built = ogeom_algo::make_cylinder(&mut model, Frame::WORLD, r, height, T).unwrap();
            let mesh = triangulate(&model, &built.shape, deflection, T).unwrap();
            assert!(mesh.is_closed(), "height {height}");
            assert_eq!(
                mesh.triangle_count(),
                2 * n + 2 * (n - 2),
                "height {height}"
            );
            let mut counts: Vec<usize> = explore_unique(&model, &built.shape, ShapeType::Face)
                .unwrap()
                .iter()
                .map(|face| {
                    let drawn = triangulate_face(&model, face, deflection, T).unwrap();
                    assert_faces_out(&drawn);
                    drawn.triangle_count()
                })
                .collect();
            counts.sort_unstable();
            assert_eq!(counts, [n - 2, n - 2, 2 * n], "height {height}");
        }
    }

    /// A tube a twentieth of a millimetre thick: each cap is an annulus far
    /// longer round than across, drawn as its two polygons' triangulation,
    /// one triangle per boundary point, none of them folded. Its side is
    /// ruled between the rims, so the mesh encloses exactly the prism on
    /// the two polygons.
    #[test]
    fn a_thin_annulus_is_its_outline_s_triangulation() {
        use ogeom_geom::Curve3d as _;
        let (outer, inner, height, c) = (10.05, 10.0, 5.0, 0.01);
        let n = circle_segments(inner, c);
        assert_eq!(n, circle_segments(outer, c));
        let mut model = Model::new();
        let ring = |model: &mut Model, frame: Frame, radius: f64| {
            let circle = ogeom_math::Circle::new(frame, radius, T).unwrap();
            let curve: ogeom_geom::Curve = ogeom_geom::CircleCurve::new(circle).into();
            let domain = curve.domain();
            ogeom_algo::make_edge(model, curve, domain, T)
                .unwrap()
                .shape
        };
        let down =
            Frame::from_axes(Point::ORIGIN, Direction::X, -Direction::Y, -Direction::Z, T).unwrap();
        let rim = ring(&mut model, Frame::WORLD, outer);
        let hole = ring(&mut model, down, inner);
        let plane = ogeom_math::Plane::new(Frame::WORLD);
        let surface: SurfaceGeometry =
            ogeom_geom::PlaneSurface::over(plane, (-20.0, 20.0), (-20.0, 20.0))
                .unwrap()
                .into();
        let cap =
            ogeom_algo::make_face_with_pcurves(&mut model, surface, &[vec![rim], vec![hole]], T)
                .unwrap()
                .shape;
        let tube =
            ogeom_algo::make_prism(&mut model, &cap, Vector::new(0.0, 0.0, height), T).unwrap();

        let deflection = Deflection::curvature_only(c, 0.5).unwrap();
        let mesh = triangulate(&model, &tube.shape, deflection, T).unwrap();
        assert!(mesh.is_closed());
        #[allow(clippy::cast_precision_loss, reason = "a segment count")]
        let polygons = n as f64 / 2.0 * (std::f64::consts::TAU / n as f64).sin();
        // The prism on the two inscribed polygons, exact to rounding.
        assert_relative_eq!(
            mesh.volume(),
            polygons * (outer * outer - inner * inner) * height,
            max_relative = 1e-9
        );
        for face in explore_unique(&model, &tube.shape, ShapeType::Face).unwrap() {
            let drawn = triangulate_face(&model, &face, deflection, T).unwrap();
            assert_faces_out(&drawn);
            assert_eq!(drawn.triangle_count(), 2 * n);
        }
    }

    #[test]
    fn welding_is_what_closes_the_mesh() {
        // Without it each face brings its own copy of every boundary vertex, so
        // no triangle edge is shared and the surface is a pile of loose squares.
        let mut model = Model::new();
        let built = make_box(&mut model, Frame::WORLD, (1.0, 1.0, 1.0), T).unwrap();

        let mut loose = Triangulation::new();
        for face in ogeom_topo::explore(
            &model,
            &built.shape,
            ogeom_topo::Filter::OfType(ShapeType::Face),
        )
        .unwrap()
        {
            loose.append(&triangulate_face(&model, &face, fine(), T).unwrap());
        }
        assert_eq!(loose.vertex_count(), 24, "four corners per face, unmerged");
        assert!(!loose.is_closed());

        let welded = loose.welded(T);
        assert_eq!(welded.vertex_count(), 8);
        assert!(welded.is_closed());
    }

    #[test]
    fn a_reversed_face_presents_the_other_side() {
        // A renderer or a volume computation that ignored orientation would
        // have the solid inside out, and nothing about the positions says so.
        let mut model = Model::new();
        let built = make_box(&mut model, Frame::WORLD, (1.0, 1.0, 1.0), T).unwrap();
        let face = explore_unique(&model, &built.shape, ShapeType::Face).unwrap()[0].clone();

        let forward = triangulate_face(&model, &face, fine(), T).unwrap();
        let backward = triangulate_face(&model, &face.reversed(), fine(), T).unwrap();

        assert_eq!(forward.triangle_count(), backward.triangle_count());
        for (a, b) in forward.normals.iter().zip(&backward.normals) {
            assert!(a.is_equal(-*b, T), "normals did not flip: {a:?} vs {b:?}");
        }
        // And the winding flipped with them, so the two agree.
        let winding = |m: &Triangulation, i: usize| {
            let [a, b, c] = m.triangles[i].map(|k| m.positions[k as usize]);
            (b - a).cross(c - a)
        };
        assert!(winding(&forward, 0).dot(winding(&backward, 0)) < 0.0);
    }

    /// Triangles whose winding points against the sum of their vertex
    /// normals.
    fn against(mesh: &Triangulation) -> usize {
        mesh.triangles
            .iter()
            .filter(|t| {
                let [a, b, c] = t.map(|k| mesh.positions[k as usize]);
                let facing = (b - a).cross(c - a);
                let normal = t
                    .iter()
                    .fold(Vector::ZERO, |sum, &k| sum + mesh.normals[k as usize]);
                facing.dot(normal) < 0.0
            })
            .collect::<Vec<_>>()
            .len()
    }

    #[test]
    fn a_mirrored_face_is_shaded_the_way_its_triangles_face() {
        // A reflection turns the chart's natural normal over, so the
        // winding has to turn to keep facing the way the image faces; the
        // normals are the original's carried through the mirror and do not.
        use ogeom_geom::{CircleCurve, Curve};
        use ogeom_math::{Circle, Transform};
        let mut model = Model::new();
        let frame = Frame::new(Point::ORIGIN, Direction::Z, Direction::X, T).unwrap();
        let circle = Circle::new(frame, 10.0, T).unwrap();
        let arc = ogeom_algo::make_edge(
            &mut model,
            Curve::Circle(CircleCurve::new(circle)),
            (0.0, std::f64::consts::PI),
            T,
        )
        .unwrap()
        .shape;
        let wall = ogeom_algo::make_prism(&mut model, &arc, Vector::new(0.0, 0.0, 5.0), T)
            .unwrap()
            .shape;
        let face = explore_unique(&model, &wall, ShapeType::Face)
            .unwrap()
            .remove(0);
        let copy = ogeom_algo::copied(&mut model, &face).unwrap().shape;
        let mirror = Transform::plane_mirror(Point::ORIGIN, Direction::Y);
        let image = ogeom_algo::transformed(&mut model, &copy, mirror)
            .unwrap()
            .shape;
        for f in [&face, &image, &image.reversed()] {
            let mesh = triangulate_face(&model, f, Deflection::default(), T).unwrap();
            assert!(mesh.triangle_count() > 0);
            assert_eq!(against(&mesh), 0, "{} triangles", mesh.triangle_count());
        }
        // The image's normal at (10, 0, 0) is the original's there, which
        // the mirror leaves where it was.
        let original = triangulate_face(&model, &face, Deflection::default(), T).unwrap();
        let mirrored = triangulate_face(&model, &image, Deflection::default(), T).unwrap();
        let at = |mesh: &Triangulation| {
            let i = mesh
                .positions
                .iter()
                .position(|p| p.distance(Point::new(10.0, 0.0, 0.0)) < 1e-9)
                .unwrap();
            mesh.normals[i]
        };
        assert!(at(&mirrored).is_equal(at(&original), T));
        // The whole-shape path agrees.
        let whole = triangulate(&model, &image, Deflection::default(), T).unwrap();
        assert_eq!(against(&whole), 0);
    }

    #[test]
    fn every_triangle_vertex_lies_on_the_surface_it_came_from() {
        let mut model = Model::new();
        let built = make_box(&mut model, Frame::WORLD, (2.0, 1.0, 3.0), T).unwrap();

        for face in explore_unique(&model, &built.shape, ShapeType::Face).unwrap() {
            let mesh = triangulate_face(&model, &face, fine(), T).unwrap();
            let data = model.node(&face).unwrap().data().as_face().unwrap().clone();
            let surface = model.geometry().surface(data.surface).unwrap();
            for (position, (u, v)) in mesh.positions.iter().zip(&mesh.parameters) {
                let exact = surface.point_at(*u, *v, T).unwrap();
                assert!(
                    position.is_equal(exact, T),
                    "{position:?} is not on its surface"
                );
            }
        }
    }

    #[test]
    fn a_finer_deflection_never_gives_fewer_triangles() {
        let mut model = Model::new();
        let built = make_box(&mut model, Frame::WORLD, (1.0, 1.0, 1.0), T).unwrap();
        let coarse = triangulate(&model, &built.shape, Deflection::default(), T).unwrap();
        let detailed = triangulate(&model, &built.shape, fine(), T).unwrap();
        assert!(detailed.triangle_count() >= coarse.triangle_count());
        // Both enclose the same volume, since the faces are flat.
        assert_relative_eq!(coarse.volume(), 1.0, epsilon = 1e-9);
        assert_relative_eq!(detailed.volume(), 1.0, epsilon = 1e-9);
    }

    #[test]
    fn a_sphere_gets_interior_points_and_converges_on_its_true_area() {
        // The whole point of measuring deflection in space rather than in
        // parameter space. A sphere's parameter rectangle is uniform; the
        // surface it maps to is not, and a fixed grid would be dense at the
        // poles and coarse at the equator.
        use ogeom_algo::make_natural_face;
        use ogeom_geom::SphereSurface;
        use ogeom_math::Sphere;

        let radius = 10.0;
        let exact = 4.0 * std::f64::consts::PI * radius * radius;
        let mut previous = 0.0;

        for chord in [1.0_f64, 0.25, 0.05] {
            let mut model = Model::new();
            let surface = SphereSurface::new(Sphere::new(Frame::WORLD, radius, T).unwrap());
            let face = make_natural_face(&mut model, surface.into()).unwrap().shape;
            let deflection = Deflection {
                chord,
                ..Deflection::default()
            };
            let mesh = triangulate_face(&model, &face, deflection, T).unwrap();

            assert!(
                mesh.triangle_count() > 2,
                "a curved face needs interior points, got {} triangles",
                mesh.triangle_count()
            );
            // Every triangle chord-cuts the sphere, so the area comes in under
            // the truth and climbs as the tolerance tightens.
            let area = mesh.area();
            assert!(area < exact, "a chord-cut area cannot exceed the surface's");
            assert!(
                area > previous,
                "tightening the chord from the previous step lost area: \
                 {area} after {previous}"
            );
            previous = area;
        }
        assert!(
            previous > exact * 0.99,
            "at a chord of 0.05 on a radius of 10 the area should be within a \
             percent, got {previous} against {exact}"
        );
    }

    #[test]
    fn every_sphere_vertex_sits_at_the_right_radius() {
        // Lifting through the surface is what makes the mesh curved at all; a
        // vertex left in parameter space, or lifted with the wrong parameters,
        // would land nowhere near the sphere.
        use ogeom_algo::make_natural_face;
        use ogeom_geom::SphereSurface;
        use ogeom_math::Sphere;

        let mut model = Model::new();
        let surface = SphereSurface::new(Sphere::new(Frame::WORLD, 3.0, T).unwrap());
        let face = make_natural_face(&mut model, surface.into()).unwrap().shape;
        let mesh = triangulate_face(&model, &face, Deflection::default(), T).unwrap();

        for p in &mesh.positions {
            assert_relative_eq!(p.to_vector().magnitude(), 3.0, epsilon = 1e-9);
        }
    }

    #[test]
    fn a_curved_domains_boundary_is_refined_not_just_its_corners() {
        // The corners alone would leave the triangulation nothing to connect to
        // along a side, so it reaches right across the domain for one, and a
        // sliver from a sphere's equator to its pole makes the summed area
        // diverge under refinement rather than converge.
        use ogeom_geom::{PlaneSurface, SphereSurface};
        use ogeom_math::{Plane, Sphere};

        let sphere: SurfaceGeometry =
            SphereSurface::new(Sphere::new(Frame::WORLD, 10.0, T).unwrap()).into();
        let coarse = domain_ring(&sphere, Deflection::default(), T);
        let fine_ring = domain_ring(&sphere, fine(), T);
        assert!(coarse.len() > 4, "a sphere's domain edge is curved");
        assert!(
            fine_ring.len() > coarse.len(),
            "a tighter chord should place more boundary points"
        );

        // No side may repeat a corner: a zero-length constraint is not one.
        for w in fine_ring.windows(2) {
            assert!(!w[0].is_equal(w[1], T), "the ring repeats a point");
        }

        // A plane is flat, so its domain needs only what closes the rectangle.
        let plane: SurfaceGeometry = PlaneSurface::new(Plane::new(Frame::WORLD)).into();
        assert!(domain_ring(&plane, fine(), T).len() >= 4);
    }

    #[test]
    fn winding_detects_points_inside_and_outside_a_ring() {
        let square = vec![
            Point2::new(0.0, 0.0),
            Point2::new(1.0, 0.0),
            Point2::new(1.0, 1.0),
            Point2::new(0.0, 1.0),
        ];
        assert!(inside_region(
            std::slice::from_ref(&square),
            Point2::new(0.5, 0.5)
        ));
        assert!(!inside_region(
            std::slice::from_ref(&square),
            Point2::new(1.5, 0.5)
        ));
        assert!(!inside_region(
            std::slice::from_ref(&square),
            Point2::new(0.5, -0.5)
        ));

        // With a hole, the middle is outside again.
        let hole = vec![
            Point2::new(0.4, 0.4),
            Point2::new(0.6, 0.4),
            Point2::new(0.6, 0.6),
            Point2::new(0.4, 0.6),
        ];
        let with_hole = vec![square, hole];
        assert!(!inside_region(&with_hole, Point2::new(0.5, 0.5)));
        assert!(inside_region(&with_hole, Point2::new(0.2, 0.2)));
    }

    #[test]
    fn a_polyline_of_an_edge_runs_in_the_edges_direction() {
        let mut model = Model::new();
        let built = make_box(&mut model, Frame::WORLD, (1.0, 1.0, 1.0), T).unwrap();
        let edge = explore_unique(&model, &built.shape, ShapeType::Edge).unwrap()[0].clone();

        let forward = polyline_of_edge(&model, &edge, fine(), T).unwrap();
        let backward = polyline_of_edge(&model, &edge.reversed(), fine(), T).unwrap();
        assert!(forward.len() >= 2);
        assert!(forward[0].is_equal(backward[backward.len() - 1], T));
        assert!(forward[forward.len() - 1].is_equal(backward[0], T));
    }

    #[test]
    fn triangulating_something_that_is_not_a_face_is_refused() {
        let mut model = Model::new();
        let built = make_box(&mut model, Frame::WORLD, (1.0, 1.0, 1.0), T).unwrap();
        assert!(triangulate_face(&model, &built.shape, fine(), T).is_err());

        let vertex = explore_unique(&model, &built.shape, ShapeType::Vertex).unwrap()[0].clone();
        assert!(triangulate_face(&model, &vertex, fine(), T).is_err());
        assert!(polyline_of_edge(&model, &vertex, fine(), T).is_err());
    }

    #[test]
    fn a_triangles_normal_is_none_when_it_is_degenerate() {
        let a = Point::ORIGIN;
        let b = Point::new(1.0, 0.0, 0.0);
        assert!(triangle_normal(a, b, Point::new(0.0, 1.0, 0.0), T).is_some());
        assert!(triangle_normal(a, b, Point::new(2.0, 0.0, 0.0), T).is_none());
        assert!(triangle_normal(a, a, a, T).is_none());
    }

    /// A ten millimetre square plate at `z = 0` moved by `at`, with the
    /// inner wire `hole` builds into its model given the same move, and its
    /// mesh's Euler characteristic and area.
    fn plate_with(at: Vector, hole: impl Fn(&mut Model, Vector) -> Shape) -> (i64, f64) {
        use ogeom_algo::{make_face, make_polygon};
        let mut model = Model::new();
        let p = |x: f64, y: f64| Point::new(x, y, 0.0) + at;
        let outer = make_polygon(
            &mut model,
            &[p(0.0, 0.0), p(10.0, 0.0), p(10.0, 10.0), p(0.0, 10.0)],
            true,
            T,
        )
        .unwrap()
        .shape;
        let inner = hole(&mut model, at);
        let plane: SurfaceGeometry = ogeom_geom::PlaneSurface::new(ogeom_math::Plane::through(
            Point::ORIGIN + at,
            Direction::Z,
        ))
        .into();
        let face = make_face(&mut model, plane, &[outer, inner], T)
            .unwrap()
            .shape;
        ogeom_algo::build::trimmed_where_bare(&mut model, &face, T).unwrap();
        let mesh = triangulate_face(&model, &face, Deflection::default(), T).unwrap();
        let mut edges = ogeom_core::FastSet::default();
        let mut area = 0.0;
        for t in &mesh.triangles {
            for i in 0..3 {
                let (a, b) = (t[i], t[(i + 1) % 3]);
                edges.insert((a.min(b), a.max(b)));
            }
            let [a, b, c] = t.map(|i| mesh.positions[i as usize]);
            area += 0.5 * (b - a).cross(c - a).magnitude();
        }
        #[allow(clippy::cast_possible_wrap)]
        let euler = mesh.positions.len() as i64 - edges.len() as i64 + mesh.triangles.len() as i64;
        (euler, area)
    }

    #[test]
    fn a_hole_thinner_than_a_micron_stays_a_hole() {
        // A slot two millimetres long and six tenths of a micron across: a
        // ring of four distinct vertices, which goes round rather than out
        // and back, so it is a hole however thin, and the plate meshes
        // with one hole in it.
        let (euler, area) = plate_with(Vector::ZERO, |model, _| {
            let p = |x: f64, y: f64| Point::new(x, y, 0.0);
            ogeom_algo::make_polygon(
                model,
                &[p(4.0, 5.0), p(4.0, 5.0006), p(6.0, 5.0006), p(6.0, 5.0)],
                true,
                T,
            )
            .unwrap()
            .shape
        });
        assert_eq!(euler, 0, "one hole: V - E + F is 1 - 1");
        assert_relative_eq!(area, 100.0 - 2.0 * 0.0006, epsilon = 1e-7);
    }

    #[test]
    fn a_ring_walked_out_and_back_thinner_than_a_micron_is_a_slit() {
        // Out from one vertex to a second and a third along two lines, and
        // home through the same three along two arcs bulging six tenths of
        // a micron off them: every edge is answered by one walking it back,
        // and the ring is a few tenths of a micron wide on average. It is a
        // slit, and the plate meshes whole, up to 2e7 from the origin too.
        for offset in [0.0, 1.0e6, 1.0e7, 2.0e7] {
            let (euler, area) = slit_plate(Vector::new(offset, -0.7 * offset, 0.3 * offset));
            assert_eq!(euler, 1, "no hole at {offset:e}: V - E + F is 1");
            assert_relative_eq!(area, 100.0, epsilon = 1e-6);
        }
    }

    #[test]
    fn a_small_round_hole_in_two_halves_stays_a_hole_far_from_the_origin() {
        // A hole five microns in radius and one four tenths of a micron in
        // radius, each circle cut into two halves between the same two
        // vertices: the wire walks back by its vertices, and the smaller
        // is under a micron wide. Both halves bound a round region whose
        // width is a sixth of its perimeter, not two sides lying along
        // each other, so the plate meshes with one hole in it up to 2e7
        // from the origin.
        for radius in [0.005, 0.0004] {
            for offset in [0.0, 1.0e6, 1.0e7, 2.0e7] {
                let at = Vector::new(offset, -0.7 * offset, 0.3 * offset);
                let (euler, _) = plate_with(at, |model, by| {
                    use ogeom_algo::{make_edge_between, make_wire};
                    let frame = ogeom_math::Frame::new(
                        Point::new(5.0, 5.0, 0.0) + by,
                        Direction::Z,
                        Direction::X,
                        T,
                    )
                    .unwrap();
                    let circle = ogeom_math::Circle::new(frame, radius, T).unwrap();
                    let ends = [
                        frame.origin() + frame.x().vector() * radius,
                        frame.origin() - frame.x().vector() * radius,
                    ];
                    let v: Vec<Shape> = ends
                        .iter()
                        .map(|q| model.add_vertex(ogeom_topo::VertexData::new(*q)))
                        .collect();
                    let half = |model: &mut Model, range: (f64, f64), from: &Shape, to: &Shape| {
                        let curve: ogeom_geom::Curve = ogeom_geom::CircleCurve::new(circle).into();
                        make_edge_between(model, curve, range, from, to, T)
                            .unwrap()
                            .shape
                    };
                    let pi = core::f64::consts::PI;
                    let edges = [
                        half(model, (0.0, pi), &v[0], &v[1]),
                        half(model, (pi, 2.0 * pi), &v[1], &v[0]),
                    ];
                    make_wire(model, &edges, T).unwrap().shape.reversed()
                });
                assert_eq!(
                    euler, 0,
                    "radius {radius}, one hole at {offset:e}: V - E + F is 1 - 1"
                );
            }
        }
    }

    #[test]
    fn two_edges_lying_along_each_other_thinner_than_a_micron_are_a_slit() {
        // Out from one vertex to a second along a line two millimetres
        // long, home along an arc bulging six tenths of a micron off it:
        // two edges between the same two vertices, like a round hole cut
        // in two halves, but the ring is a ten-thousandth as wide as it is
        // long. It is a slit, and the plate meshes whole.
        let (euler, area) = plate_with(Vector::ZERO, |model, _| {
            use ogeom_algo::{Spacing, interpolate, make_edge_between, make_wire};
            use ogeom_geom::Curve3d as _;
            let at = [Point::new(4.0, 5.0, 0.0), Point::new(6.0, 5.0, 0.0)];
            let v: Vec<Shape> = at
                .iter()
                .map(|q| model.add_vertex(ogeom_topo::VertexData::new(*q)))
                .collect();
            let line: ogeom_geom::Curve = ogeom_geom::LineCurve::segment(at[0], at[1], T)
                .unwrap()
                .into();
            let out = make_edge_between(model, line, (0.0, 2.0), &v[0], &v[1], T)
                .unwrap()
                .shape;
            let mid = at[0].midpoint(at[1]) + Vector::new(0.0, 0.0006, 0.0);
            let arc: ogeom_geom::Curve =
                interpolate(&[at[1], mid, at[0]], 2, Spacing::Centripetal, T)
                    .unwrap()
                    .into();
            let range = arc.domain();
            let back = make_edge_between(model, arc, range, &v[1], &v[0], T)
                .unwrap()
                .shape;
            make_wire(model, &[out, back], T).unwrap().shape
        });
        assert_eq!(euler, 1, "no hole: V - E + F is 1");
        // The plate's own area, to the rounding of a hundred square
        // millimetres summed over its triangles.
        assert_relative_eq!(area, 100.0, epsilon = 1e-6);
    }

    /// The plate of the test above with its slit, moved by `at`.
    fn slit_plate(at: Vector) -> (i64, f64) {
        plate_with(at, |model, by| {
            use ogeom_algo::{Spacing, interpolate, make_edge_between, make_wire};
            use ogeom_geom::Curve3d as _;
            let p = |x: f64, y: f64| Point::new(x, y, 0.0) + by;
            let at = [p(3.0, 5.0), p(5.0, 5.0), p(7.0, 5.0)];
            let v: Vec<Shape> = at
                .iter()
                .map(|q| model.add_vertex(ogeom_topo::VertexData::new(*q)))
                .collect();
            let mut edges = Vec::new();
            for i in 0..2 {
                let line: ogeom_geom::Curve = ogeom_geom::LineCurve::segment(at[i], at[i + 1], T)
                    .unwrap()
                    .into();
                let range = (0.0, at[i].distance(at[i + 1]));
                edges.push(
                    make_edge_between(model, line, range, &v[i], &v[i + 1], T)
                        .unwrap()
                        .shape,
                );
            }
            for i in (0..2).rev() {
                let mid = at[i].midpoint(at[i + 1]) + Vector::new(0.0, 0.0006, 0.0);
                let arc: ogeom_geom::Curve =
                    interpolate(&[at[i + 1], mid, at[i]], 2, Spacing::Centripetal, T)
                        .unwrap()
                        .into();
                let range = arc.domain();
                edges.push(
                    make_edge_between(model, arc, range, &v[i + 1], &v[i], T)
                        .unwrap()
                        .shape,
                );
            }
            make_wire(model, &edges, T).unwrap().shape
        })
    }
}

#[cfg(test)]
#[allow(clippy::unwrap_used)]
mod predicate_tests {
    use super::*;
    use ogeom_core::Fast;

    /// A triangle with one very long, very nearly diagonal edge.
    ///
    /// The configuration where a naive crossing test goes wrong. Subtracting
    /// the edge's ends from a query point close to the line cancels almost
    /// every significant digit, and what is left is rounding rather than
    /// geometry.
    fn sliver() -> Vec<Vec<Point2>> {
        vec![vec![
            Point2::new(0.5, 0.5),
            Point2::new(1000.0, 1000.0),
            Point2::new(1000.0, 0.5),
        ]]
    }

    #[test]
    fn the_two_implementations_are_a_real_choice_and_not_a_decoration() {
        // The point of the seam. In a band of near-degenerate queries the two
        // answer differently, and if they never did, routing through the trait
        // would be ceremony rather than a design.
        //
        // The query points straddle the long edge at a spacing far below what
        // the subtraction can resolve, which is exactly the case a mesh hits
        // when a face's boundary passes close to a triangulation vertex.
        let rings = sliver();
        let step = 2.0_f64.powi(-48);
        let mut disagreements = 0;
        for i in 0..256i32 {
            for j in 0..256i32 {
                let p = Point2::new(
                    250.0 + f64::from(i - 128) * step,
                    250.0 + f64::from(j - 128) * step,
                );
                if inside_boundary_with::<Exact>(&rings, p)
                    != inside_boundary_with::<Fast>(&rings, p)
                {
                    disagreements += 1;
                }
            }
        }
        assert!(
            disagreements > 0,
            "the exact and fast predicates never disagreed, so the seam is not \
             carrying anything"
        );
    }

    #[test]
    fn the_exact_predicate_is_the_one_that_is_right_near_the_edge() {
        // Disagreeing is not enough; the exact one has to be the *correct*
        // one, and that needs points whose side is known without asking either
        // implementation.
        //
        // Stepped in units of the last place rather than by a small distance.
        // A distance below the spacing of `f64` at 250 rounds away entirely,
        // leaving two points that are both exactly *on* the diagonal, where
        // either answer is defensible and the test would be asserting nothing.
        let rings = sliver();
        let base = 250.0_f64;
        let mut off = base;
        for k in 1..64 {
            off = off.next_up();
            assert_ne!(off, base, "step {k} did not move the point at all");
            // Inside this triangle is below the diagonal `y = x`: larger x.
            assert!(
                inside_boundary_with::<Exact>(&rings, Point2::new(off, base)),
                "a point {k} ulps below the diagonal was reported outside"
            );
            assert!(
                !inside_boundary_with::<Exact>(&rings, Point2::new(base, off)),
                "a point {k} ulps above the diagonal was reported inside"
            );
        }
    }

    #[test]
    fn the_exact_answer_is_the_one_the_geometry_supports() {
        // Points placed by construction, so the right answer is known without
        // asking either implementation.
        let square = vec![vec![
            Point2::new(0.0, 0.0),
            Point2::new(1.0, 0.0),
            Point2::new(1.0, 1.0),
            Point2::new(0.0, 1.0),
        ]];
        assert!(inside_boundary_with::<Exact>(
            &square,
            Point2::new(0.5, 0.5)
        ));
        assert!(!inside_boundary_with::<Exact>(
            &square,
            Point2::new(1.5, 0.5)
        ));
        assert!(!inside_boundary_with::<Exact>(
            &square,
            Point2::new(-0.5, 0.5)
        ));
        assert!(!inside_boundary_with::<Exact>(
            &square,
            Point2::new(0.5, 1.5)
        ));

        // A vertex exactly on the sampling ray is counted once, not twice or
        // not at all, which is what the half-open comparison is for.
        let diamond = vec![vec![
            Point2::new(0.0, 0.0),
            Point2::new(1.0, 1.0),
            Point2::new(2.0, 0.0),
            Point2::new(1.0, -1.0),
        ]];
        assert!(inside_boundary_with::<Exact>(
            &diamond,
            Point2::new(1.0, 0.0)
        ));
        assert!(!inside_boundary_with::<Exact>(
            &diamond,
            Point2::new(3.0, 0.0)
        ));
        assert!(!inside_boundary_with::<Exact>(
            &diamond,
            Point2::new(-1.0, 0.0)
        ));
    }

    /// A cylinder's grid has as many rows as its length needs, not as
    /// many as its sag asks for.
    #[test]
    fn rows_are_spread_until_no_cell_is_taller_than_its_aspect() {
        // Sixteen columns a millimetre wide over a region a hundred long:
        // sag left one interior row, the aspect wants cells six tall.
        let sagged = vec![0.0, 50.0, 100.0];
        let spread = spread_to_aspect(
            sagged.clone(),
            CELL_ASPECT,
            0.0,
            16.0,
            |_| 17,
            |i| (sagged[i + 1] - sagged[i]).abs(),
            |a, b, _| (b - a).abs(),
        );
        assert_eq!(
            spread.len(),
            2 * 9 + 1,
            "each 50 mm half in nine 6 mm pieces: {spread:?}"
        );
        assert_eq!(spread[0], 0.0);
        assert_eq!(spread[9], 50.0, "the rows sag chose stay where they were");
        assert_eq!(spread[18], 100.0);
        assert!(spread.windows(2).all(|w| w[1] > w[0]));

        // Three columns is flat along `u` as well; nothing to hold to.
        let flat = spread_to_aspect(
            sagged.clone(),
            CELL_ASPECT,
            0.0,
            16.0,
            |_| 3,
            |i| (sagged[i + 1] - sagged[i]).abs(),
            |a, b, _| (b - a).abs(),
        );
        assert_eq!(flat, sagged);

        // Cells already shorter than the aspect are left alone.
        let fine = vec![0.0, 4.0, 8.0];
        let same = spread_to_aspect(
            fine.clone(),
            CELL_ASPECT,
            0.0,
            16.0,
            |_| 17,
            |i| (fine[i + 1] - fine[i]).abs(),
            |a, b, _| (b - a).abs(),
        );
        assert_eq!(same, fine);
    }

    /// A ring's slop-duplicates and spikes come off before it is
    /// triangulated.
    #[test]
    fn a_ring_is_cleaned_of_duplicates_and_spikes() {
        let p = |x: f64, y: f64| Point2::new(x, y);
        // A square whose closing point repeats its first a hair off, and
        // whose right side steps out to a point and straight back.
        let mut ring = vec![
            p(0.0, 0.0),
            p(10.0, 0.0),
            p(10.0, 5.0),
            p(10.2, 5.0),
            p(10.0, 5.0 + 1e-9),
            p(10.0, 10.0),
            p(0.0, 10.0),
            p(1e-9, 1e-9),
        ];
        let mut anchors = vec![None; ring.len()];
        let reach = chart_reach(&ring);
        merge_near_duplicates(&mut ring, &mut anchors, reach);
        assert_eq!(ring.len(), 7, "the closing duplicate is gone: {ring:?}");
        remove_spikes(&mut ring, &mut anchors, reach);
        assert_eq!(
            ring,
            vec![
                p(0.0, 0.0),
                p(10.0, 0.0),
                p(10.0, 5.0),
                p(10.0, 10.0),
                p(0.0, 10.0)
            ],
            "the spike is gone and the ring is the square with a point on one side"
        );
        assert_eq!(anchors.len(), ring.len());
    }

    /// A ring a ten-thousandth of a radian across a drum's chart and 365
    /// along it (a strip of a radius-32 drum a few microns wide) keeps its
    /// width: the points across it are a millionth of the ring's extent
    /// across apart many times over, though not a millionth of its length.
    /// Its closing point repeated a hair off along both directions still
    /// goes.
    #[test]
    fn a_ring_narrow_in_one_chart_direction_keeps_its_width() {
        let p = |x: f64, y: f64| Point2::new(x, y);
        let w = 1.77e-4;
        let strip = vec![
            p(0.0, 0.0),
            p(w / 2.0, 0.0),
            p(w, 0.0),
            p(w, 365.0),
            p(w / 2.0, 365.0),
            p(0.0, 365.0),
        ];
        let mut ring = strip.clone();
        ring.push(p(w * 1e-9, 1e-9));
        let mut anchors = vec![None; ring.len()];
        let reach = chart_reach(&ring);
        merge_near_duplicates(&mut ring, &mut anchors, reach);
        remove_spikes(&mut ring, &mut anchors, reach);
        assert_eq!(ring, strip, "the strip keeps every corner");
        assert_eq!(anchors.len(), ring.len());
    }

    #[test]
    fn a_hole_is_outside_however_either_ring_is_wound() {
        // Even-odd does not depend on the wires being wound consistently, which
        // is what makes it survive imported geometry.
        let outer = vec![
            Point2::new(0.0, 0.0),
            Point2::new(4.0, 0.0),
            Point2::new(4.0, 4.0),
            Point2::new(0.0, 4.0),
        ];
        let hole: Vec<Point2> = vec![
            Point2::new(1.0, 1.0),
            Point2::new(3.0, 1.0),
            Point2::new(3.0, 3.0),
            Point2::new(1.0, 3.0),
        ];
        let backwards: Vec<Point2> = hole.iter().rev().copied().collect();
        for inner in [hole, backwards] {
            let rings = vec![outer.clone(), inner];
            assert!(!inside_boundary_with::<Exact>(
                &rings,
                Point2::new(2.0, 2.0)
            ));
            assert!(inside_boundary_with::<Exact>(&rings, Point2::new(0.5, 0.5)));
        }
    }
}