geop_core_topology/edit/

splice_edge_into_face.rs

use geop_core_math::{
    geop_error::{GeopError, GeopResult, WithContext},
    scalars::Scalar,
    solvers::bisect::bisect,
    vector::Vector2,
};

use crate::{
    Coedge, CoedgeGeometry, CoedgeId, Curve2, EdgeId, Face, FaceId, Model, Sense, VertexId,
    boundary::{BoundaryIndex, BoundaryType},
    contains::face::{PointClassification, face_contains},
    uv_loop::loop_signed_area,
};

impl Model {
    /// Splice an **already-existing** `edge_id` into `face_id`'s boundary
    /// structure, as a forward/reversed coedge pair whose pcurves it fits onto
    /// the face's own surface.
    ///
    /// The edge already exists because two faces are meant to *share* it
    /// (the whole point of imprinting an intersection curve), so this
    /// rewires loops around the edge it's handed rather than making one.
    ///
    /// What it does depends on how much of the edge's own topology the face
    /// already knows about, which is the only thing that determines what a
    /// correct result even looks like:
    ///
    /// | edge's start/end vertex on this face's boundary | result |
    /// |---|---|
    /// | neither | a new self-contained ring floating inside the face — a **hole** |
    /// | exactly one | a spur (out along the edge and back) inserted into that vertex's own loop |
    /// | both, same hole | that hole is **divided into two holes** |
    /// | both, two different holes | those holes are **merged into one** |
    /// | both, a hole and the outer loop | the hole is **absorbed** into the outer loop |
    /// | both, the outer loop | the face is **split into two faces** |
    ///
    /// Splitting the face, and dividing a hole (whose material side becomes
    /// a face of its own), are the cases that create a face; its id is
    /// returned, `None` otherwise. The face split is the only one that has to
    /// reclassify the old face's holes — each now lies inside exactly one of
    /// the two halves (see [`split_face`]). A spur changes no connectivity at
    /// all, being "wire" topology rather than a real trim boundary.
    pub fn splice_edge_into_face(
        &mut self,
        edge_id: EdgeId,
        face_id: FaceId,
        max_nodes: usize,
        min_subdivision_size: Scalar,
    ) -> GeopResult<Option<FaceId>> {
        let model = self;
        let ctx = |e: GeopError| {
            e.with_context(format!(
                "Model::splice_edge_into_face(edge={edge_id}, face={face_id})"
            ))
        };

        let edge = model.get_edge(edge_id).with_context(&ctx)?;
        let (start_vertex, end_vertex) = (edge.start_vertex, edge.end_vertex);
        let curve = edge.curve.clone();

        // The coedge each of the edge's endpoints would attach *after*: one
        // already arriving at that vertex. Where the face already has one,
        // its pcurve end is the authoritative `(u, v)` there — the new pcurve
        // is pinned to it so the loop stays exactly continuous.
        let surface = model.get_face(face_id).with_context(&ctx)?.surface.clone();
        let fit = |model: &Model, at_start: Option<CoedgeId>, at_end: Option<CoedgeId>| {
            surface
                .fit_pcurve(
                    &curve,
                    pin(model, &surface, at_start, start_vertex)?,
                    pin(model, &surface, at_end, end_vertex)?,
                    max_nodes,
                    min_subdivision_size,
                )
                .with_context(&ctx)
        };
        let mut arriving_start =
            coedges_ending_at(model, face_id, start_vertex).with_context(&ctx)?;
        let mut arriving_end = coedges_ending_at(model, face_id, end_vertex).with_context(&ctx)?;
        let first = |arriving: &[CoedgeId]| arriving.first().copied();
        // A pole is a whole row of the surface at one vertex: no one corner
        // of the loop is where an edge arrives at it, so the first fit pins
        // neither end there, and the end's own `(u, v)` says where along the
        // row it is (see [`split_pole_row`]).
        let pole = |model: &Model, vertex: VertexId, arriving: &[CoedgeId]| {
            arriving.iter().any(|&c| {
                model
                    .get_coedge(c)
                    .is_ok_and(|c| c.geometry == CoedgeGeometry::Vertex(vertex))
            })
        };
        let (pole_start, pole_end) = (
            pole(model, start_vertex, &arriving_start),
            pole(model, end_vertex, &arriving_end),
        );
        let mut pcurve_fwd = fit(
            model,
            first(&arriving_start).filter(|_| !pole_start),
            first(&arriving_end).filter(|_| !pole_end),
        )?;
        if pole_start || pole_end {
            let (t0, t1) = pcurve_fwd.domain();
            for (is_pole, vertex, t, arriving) in [
                (pole_start, start_vertex, t0, &mut arriving_start),
                (pole_end, end_vertex, t1, &mut arriving_end),
            ] {
                if !is_pole {
                    continue;
                }
                let uv = pcurve_fwd.evaluate(t).with_context(&ctx)?;
                if let Some(corner) =
                    split_pole_row(model, vertex, arriving, uv).with_context(&ctx)?
                {
                    *arriving = vec![corner];
                }
            }
            pcurve_fwd = fit(model, first(&arriving_start), first(&arriving_end))?;
        }
        // A vertex the loop passes more than once — the root of a spur, or a
        // point where the face touches itself — has as many corners there,
        // and the edge belongs in the one it leaves into. Its pcurve says
        // which, so it is chosen after a first fit and refitted if that
        // moved the pin.
        let at_start = choose_corner(
            model,
            &arriving_start,
            Leg {
                pcurve: &pcurve_fwd,
                from_end: false,
            },
        )
        .with_context(&ctx)?;
        let at_end = choose_corner(
            model,
            &arriving_end,
            Leg {
                pcurve: &pcurve_fwd,
                from_end: true,
            },
        )
        .with_context(&ctx)?;
        if at_start != first(&arriving_start) || at_end != first(&arriving_end) {
            pcurve_fwd = fit(model, at_start, at_end)?;
        }
        let pcurve_rev = pcurve_fwd.reverse();

        // A spur dangles into the face: its free end has to lie strictly
        // inside it. Were it outside, the face's loop would run out of the
        // face and back — every check of the loop's own structure passes, and
        // the face is left straddling whatever the edge was meant to divide
        // (`narrow_groove` surfaced that far downstream, in classification).
        // On the boundary, it would be a vertex of the loop already.
        if let (Some(_), None) | (None, Some(_)) = (at_start, at_end) {
            let (t0, t1) = pcurve_fwd.domain();
            let (free, t) = if at_start.is_some() {
                (end_vertex, t1)
            } else {
                (start_vertex, t0)
            };
            let uv = pcurve_fwd.evaluate(t).with_context(&ctx)?;
            let class = face_contains(
                model,
                face_id,
                uv[0],
                uv[1],
                max_nodes,
                min_subdivision_size,
                HOLE_CLASSIFY_SEED,
            )
            .with_context(&ctx)?;
            if class != PointClassification::Inside {
                let point = model.get_vertex(free).with_context(&ctx)?.point;
                return Err(ctx(GeopError::new(format!(
                    "a spur's free end must lie inside the face, but vertex {free} at {point:?} (uv={uv:?}) is {class:?}; the face is bounded by (corners, outer loop first) {:?}",
                    model.face_corners(face_id).with_context(&ctx)?
                ))));
            }
        }

        let fwd = model.insert_coedge(Coedge {
            geometry: CoedgeGeometry::Edge(edge_id),
            sense: Sense::Forward,
            pcurve: pcurve_fwd,
            next: CoedgeId(0),
            prev: CoedgeId(0),
            face: face_id,
        });
        let rev = model.insert_coedge(Coedge {
            geometry: CoedgeGeometry::Edge(edge_id),
            sense: Sense::Reversed,
            pcurve: pcurve_rev,
            next: CoedgeId(0),
            prev: CoedgeId(0),
            face: face_id,
        });

        let mut new_face = None;
        match (at_start, at_end) {
            // Neither endpoint is on this face yet: a self-contained
            // two-coedge ring, floating inside the face — that is a hole.
            (None, None) => {
                link(model, fwd, rev).with_context(&ctx)?;
                link(model, rev, fwd).with_context(&ctx)?;
                model
                    .get_face_mut(face_id)
                    .with_context(&ctx)?
                    .holes
                    .push(BoundaryType::Loop(fwd));
            }
            // Only one endpoint is on the face: splice a spur into that loop —
            // out along the edge and straight back, so the loop still closes.
            (Some(a), None) => {
                splice_spur(model, a, fwd, rev).with_context(&ctx)?;
            }
            (None, Some(a)) => {
                splice_spur(model, a, rev, fwd).with_context(&ctx)?;
            }
            // Both endpoints are already on the face, so the new edge joins
            // two points of its boundary structure. *Which* boundaries they
            // sit on decides what that means topologically, and the four
            // cases are genuinely different operations.
            (Some(a), Some(b)) => {
                let loop_a = model
                    .find_boundary_containing(face_id, a)
                    .with_context(&ctx)?;
                let loop_b = model
                    .find_boundary_containing(face_id, b)
                    .with_context(&ctx)?;

                // `fwd` runs start->end, so it leaves the coedge arriving at
                // `start_vertex` and lands on whatever leaves `end_vertex`;
                // `rev` closes the complementary path the other way round.
                let a_next = model.get_coedge(a).with_context(&ctx)?.next;
                let b_next = model.get_coedge(b).with_context(&ctx)?.next;

                link(model, a, fwd).with_context(&ctx)?;
                link(model, fwd, b_next).with_context(&ctx)?;
                link(model, b, rev).with_context(&ctx)?;
                link(model, rev, a_next).with_context(&ctx)?;

                match (loop_a, loop_b) {
                    // Two points of the *outer* loop: the edge cuts the face
                    // itself in two. Both rings bound material, so neither can
                    // be a hole of the other — this is the only case that
                    // creates a face.
                    (BoundaryIndex::Outer, BoundaryIndex::Outer) => {
                        // Both halves must enclose material. A half of zero
                        // area means the new edge runs along the stretch of
                        // boundary between its own endpoints — it duplicates
                        // a path the face already has, so "splitting" there
                        // carves off a sliver rather than two faces. Reported
                        // rather than built: a zero-area face is structurally
                        // perfect and every other check accepts it, so it
                        // would surface much later as a face nothing can be
                        // classified against.
                        for ring in [fwd, rev] {
                            let area = loop_signed_area(model, ring).with_context(&ctx)?;
                            if !area.abs().definitely_greater(Scalar::ZERO) {
                                let widest = widest_pcurve(model, ring).with_context(&ctx)?;
                                return Err(ctx(GeopError::new(format!(
                                    "{DEGENERATE_SPLIT}: splicing edge {edge_id} into face {face_id} would split its outer loop into a ring of signed area {area:?} — the edge runs along the boundary it is being spliced into, so one side encloses nothing (the ring's widest pcurve: {widest})"
                                ))));
                            }
                        }
                        new_face = Some(
                            split_face(model, face_id, fwd, rev, max_nodes, min_subdivision_size)
                                .with_context(&ctx)?,
                        );
                    }
                    // Two points of the *same* hole. The edge runs through
                    // material, so together with one of the two arcs it just
                    // cut the hole into, it encloses a patch of material that
                    // is now bounded on its own — a new face — while the other
                    // arc still bounds a void and stays a hole.
                    //
                    // Which is which is exactly the winding: by this kernel's
                    // convention an outer loop runs counter-clockwise in
                    // `(u, v)` and a hole runs clockwise, so the ring with
                    // positive signed area bounds material and the one with
                    // negative area bounds a void. Nothing weaker will do —
                    // both rings pass through the same vertices and contain
                    // the same edge, so no purely topological test can tell
                    // them apart.
                    (BoundaryIndex::Hole(i), BoundaryIndex::Hole(j)) if i == j => {
                        let area = loop_signed_area(model, fwd).with_context(&ctx)?;
                        let (material, void) = if area.definitely_greater(Scalar::ZERO) {
                            (fwd, rev)
                        } else if area.definitely_less(Scalar::ZERO) {
                            (rev, fwd)
                        } else {
                            return Err(ctx(GeopError::new(format!(
                                "splice_edge_into_face: splitting hole {i} of face {face_id} with edge {edge_id} produced a ring of signed area {area:?}, which is not decidably clockwise or counter-clockwise — the split would be degenerate"
                            ))));
                        };
                        model
                            .get_face_mut(face_id)
                            .with_context(&ctx)?
                            .set_boundary(BoundaryIndex::Hole(i), BoundaryType::Loop(void));
                        new_face = Some(
                            new_face_from_ring(
                                model,
                                face_id,
                                material,
                                max_nodes,
                                min_subdivision_size,
                            )
                            .with_context(&ctx)?,
                        );
                    }
                    // Two *different* holes: the bridge merges them into one
                    // hole.
                    (BoundaryIndex::Hole(i), BoundaryIndex::Hole(j)) => {
                        let face = model.get_face_mut(face_id).with_context(&ctx)?;
                        face.set_boundary(BoundaryIndex::Hole(i), BoundaryType::Loop(fwd));
                        face.holes.remove(j);
                    }
                    // A hole bridged to the outer loop: the hole stops being a
                    // separate boundary and its coedges become part of the one
                    // ring that now bounds the face.
                    (BoundaryIndex::Outer, BoundaryIndex::Hole(j))
                    | (BoundaryIndex::Hole(j), BoundaryIndex::Outer) => {
                        let face = model.get_face_mut(face_id).with_context(&ctx)?;
                        face.outer = BoundaryType::Loop(fwd);
                        face.holes.remove(j);
                    }
                }
            }
        }

        Ok(new_face)
    }
}

/// Where an edge arrives at a pole — `vertex`, at which a face's surface
/// collapses a whole row — if that is in the middle of the row's coedge: the
/// coedge among `arriving` that runs along the row (see
/// [`CoedgeGeometry::Vertex`]) is cut at `uv`, where the edge's pcurve ends,
/// into two, and the first is returned, which now ends there. The edge is
/// spliced after it, between the two halves of the row.
///
/// `None` where no such coedge holds `uv` strictly inside: the edge arrives
/// at one of the corners of the row, which the usual choice of corner
/// finds. A pole has a corner for every end of its row and every meridian
/// that arrives at it, each at a `(u, v)` of its own, though they are one
/// vertex — so which of them an edge arrives at is told by its `(u, v)`,
/// not by the way it leaves the vertex.
///
/// Which `u` the edge arrives at, along a row that is all one point, is the
/// edge's own: taken from where its pcurve was fitted to end — sharp, since
/// the row is cut there, and the edge is then pinned to that cut.
fn split_pole_row(
    model: &mut Model,
    vertex: VertexId,
    arriving: &[CoedgeId],
    uv: Vector2,
) -> GeopResult<Option<CoedgeId>> {
    for &row in arriving {
        let coedge = model.get_coedge(row)?;
        if coedge.geometry != CoedgeGeometry::Vertex(vertex) {
            continue;
        }
        let (t0, t1) = coedge.pcurve.domain();
        let (a, b) = (coedge.pcurve.evaluate(t0)?, coedge.pcurve.evaluate(t1)?);
        let (face, next) = (coedge.face, coedge.next);
        if !uv[1].could_be_equal(a[1]) || !uv[1].could_be_equal(b[1]) {
            continue;
        }
        let u = uv[0].midpoint().sharpen();
        let (lo, hi) = (a[0].min(b[0]), a[0].max(b[0]));
        if !u.definitely_greater(lo) || !u.definitely_less(hi) {
            continue;
        }
        // The row runs straight along `u`, at a speed of its own.
        let t = t0
            .add(t1.sub(t0).mul(u.sub(a[0]).div(b[0].sub(a[0]))?))
            .sharpen();
        let (before, after) = coedge.pcurve.split(t)?;
        model.get_coedge_mut(row)?.pcurve = before;
        let rest = model.insert_coedge(Coedge {
            geometry: CoedgeGeometry::Vertex(vertex),
            sense: Sense::Forward,
            pcurve: after,
            next,
            prev: row,
            face,
        });
        link(model, row, rest)?;
        link(model, rest, next)?;
        return Ok(Some(row));
    }
    Ok(None)
}

/// The `(u, v)` a coedge's pcurve ends at, if there is such a coedge.
fn pcurve_end_uv(model: &Model, coedge: Option<CoedgeId>) -> GeopResult<Option<Vector2>> {
    let Some(coedge) = coedge else {
        return Ok(None);
    };
    let pcurve = &model.get_coedge(coedge)?.pcurve;
    let (_, t1) = pcurve.domain();
    Ok(Some(pcurve.evaluate(t1)?))
}

/// The `(u, v)` to pin a new pcurve's end to at `vertex`, where `coedge`
/// arrives on `surface`: that coedge's end — what keeps the loop continuous
/// — intersected with `vertex`'s own projection onto the surface. Both are
/// enclosures of the one `(u, v)` the vertex is at, so their intersection is
/// too, and it is no wider than the narrower: an end widened by its own
/// pcurve's history (a pin of its own, the drift it encloses) does not pass
/// that width on — a pcurve is widened to hold its pins whole, all along
/// (see `NurbCurve2D::interpolate_enclosing`), so a wide pin is a wide
/// pcurve. A projection that fails, or misses the end altogether, leaves
/// the end as it is.
fn pin(
    model: &Model,
    surface: &geop_core_geometry::nurb_surface::NurbSurface3D,
    coedge: Option<CoedgeId>,
    vertex: VertexId,
) -> GeopResult<Option<Vector2>> {
    let Some(end) = pcurve_end_uv(model, coedge)? else {
        return Ok(None);
    };
    let point = model.get_vertex(vertex)?.point;
    let Ok((u, v)) = surface.project(point, end[0].midpoint(), end[1].midpoint(), PIN_ITERATIONS)
    else {
        return Ok(Some(end));
    };
    let overlaps = end[0].could_be_equal(u) && end[1].could_be_equal(v);
    Ok(Some(if overlaps {
        Vector2::from_array([end[0].intersect(u), end[1].intersect(v)])
    } else {
        end
    }))
}

/// Newton iterations refining a pin (see [`pin`]): from a seed already on
/// the foot point's doorstep, a handful converge.
const PIN_ITERATIONS: usize = 8;

/// Every coedge of `face_id` that *arrives* at `vertex`: one per time its
/// boundary passes through it. A new edge leaving `vertex` has to be spliced
/// in after one of them — see [`choose_corner`].
fn coedges_ending_at(
    model: &Model,
    face_id: FaceId,
    vertex: VertexId,
) -> GeopResult<Vec<CoedgeId>> {
    let mut arriving = Vec::new();
    for coedge_id in model.iterate_face_coedges(face_id) {
        if model.coedge_end_vertex_id(coedge_id)? == vertex {
            arriving.push(coedge_id);
        }
    }
    Ok(arriving)
}

/// Of the coedges `arriving` at one vertex, the one whose corner a new edge
/// leaving the vertex along `leaving` runs into.
///
/// Each pass of the boundary through the vertex makes a corner: the face's
/// material fills the angle swept counter-clockwise from the way the boundary
/// leaves to the way it arrived from — material lies to the left of every
/// loop, outer and hole alike. One pass leaves no choice. Several do — a
/// spur's root, where the boundary passes once on each side of it — and
/// splicing into the wrong corner puts the edge on the wrong side of the
/// spur: the face split it makes then carries the spur into the half it does
/// not lie in. That happened on `three_inscribed_bores_remesh_imprints_every_arc`:
/// an arc spliced as a spur from a point of a bore's rim, and the rim arc
/// spliced at that same point a moment later, landed the spur in the outer
/// piece of the wall.
///
/// The curves are ordered around the vertex by where they first leave a
/// small circle around it, not by their tangents: see [`Leg`].
fn choose_corner(
    model: &Model,
    arriving: &[CoedgeId],
    leaving: Leg<'_>,
) -> GeopResult<Option<CoedgeId>> {
    if arriving.len() <= 1 {
        return Ok(arriving.first().copied());
    }
    let mut corners = Vec::new();
    for &coedge_id in arriving {
        let coedge = model.get_coedge(coedge_id)?;
        let goes_to = Leg {
            pcurve: &model.get_coedge(coedge.next)?.pcurve,
            from_end: false,
        };
        let came_from = Leg {
            pcurve: &coedge.pcurve,
            from_end: true,
        };
        corners.push((coedge_id, goes_to, came_from));
    }

    // One radius for every leg, small enough that no leg meets another
    // vertex inside it (see `Leg`).
    let mut radius = leaving.reach()?;
    for (_, goes_to, came_from) in &corners {
        for leg in [goes_to, came_from] {
            let reach = leg.reach()?;
            if reach.definitely_less(radius) {
                radius = reach;
            }
        }
    }
    let radius = radius.div(Scalar::TWO)?.sharpen();
    if !radius.definitely_greater(Scalar::ZERO) {
        return Err(GeopError::new(format!(
            "splice_edge_into_face: a curve at the vertex returns to it at once, so no circle around the vertex orders the corners after coedges {arriving:?}"
        )));
    }

    let d = leaving.exit(radius)?;
    let mut chosen = Vec::new();
    for (coedge_id, goes_to, came_from) in &corners {
        let (from, to) = (goes_to.exit(radius)?, came_from.exit(radius)?);
        match within_corner(from, to, d) {
            Some(true) => chosen.push(*coedge_id),
            Some(false) => {}
            None => {
                return Err(GeopError::new(format!(
                    "splice_edge_into_face: cannot tell whether an edge leaving towards {d:?} runs into the corner after coedge {coedge_id}, from {from:?} counter-clockwise to {to:?} (at radius {radius:?})"
                )));
            }
        }
    }
    match chosen[..] {
        [coedge_id] => Ok(Some(coedge_id)),
        _ => Err(GeopError::new(format!(
            "splice_edge_into_face: an edge leaving towards {d:?} runs into {} of the corners after coedges {arriving:?}, where it must run into exactly one",
            chosen.len()
        ))),
    }
}

/// A pcurve run away from the vertex at one of its ends.
///
/// Several curves leaving one vertex are ordered around it by where each
/// first leaves a circle of some radius around the vertex. Within a radius
/// that no curve meets another vertex inside, their pieces in the disk are
/// paths from its centre to its rim that do not cross — curves of a face's
/// boundary, and an edge being spliced into it, meet only at vertices — and
/// such paths leave the disk in the same cyclic order they leave its centre.
/// So the order is exact at any such radius, and a larger one only separates
/// the curves further.
///
/// Tangents would answer the same question only where the curves leave in
/// different directions. Intersection curves of tangent configurations often
/// leave a vertex tangent to each other, and then only curvature separates
/// them — which, differentiated twice off a pcurve that encloses a traced
/// curve, comes out too wide to decide anything.
struct Leg<'a> {
    pcurve: &'a Curve2,
    /// Whether the vertex is at the pcurve's end (and it is run backwards).
    from_end: bool,
}

impl Leg<'_> {
    /// The parameters of the vertex end and of the pcurve's midpoint.
    fn ends(&self) -> GeopResult<(Scalar, Scalar)> {
        let (t0, t1) = self.pcurve.domain();
        let mid = t0.add(t1).div(Scalar::TWO)?;
        Ok((if self.from_end { t1 } else { t0 }, mid))
    }

    /// A radius this leg reaches, before meeting any vertex other than the
    /// one it leaves: its distance to its midpoint, or to its far end where
    /// that is nearer. A curve returning to the vertex it leaves has no far
    /// end to meet.
    fn reach(&self) -> GeopResult<Scalar> {
        let (t0, t1) = self.pcurve.domain();
        let (t_vertex, t_mid) = self.ends()?;
        let t_far = if self.from_end { t0 } else { t1 };
        let origin = self.pcurve.evaluate(t_vertex)?;
        let mid = self.pcurve.evaluate(t_mid)?.sub(&origin).norm();
        let far = self.pcurve.evaluate(t_far)?.sub(&origin).norm();
        Ok(
            if far.definitely_greater(Scalar::ZERO) && far.definitely_less(mid) {
                far
            } else {
                mid
            },
        )
    }

    /// Where the leg leaves the circle of `radius` around its vertex,
    /// relative to the vertex. Located by bisection between the vertex and
    /// the midpoint, which lies outside a circle of at most half its reach;
    /// the parameters bisected are the search's own choice, so they are sharp.
    ///
    /// Bisection finds *a* crossing of the circle, which is the first exit
    /// the ordering argument needs as long as the leg does not leave the disk
    /// and come back before its midpoint — true of the short, smooth trim
    /// curves at half their reach, and an assumption for anything else.
    fn exit(&self, radius: Scalar) -> GeopResult<Vector2> {
        let (t_vertex, t_mid) = self.ends()?;
        let origin = self.pcurve.evaluate(t_vertex)?;
        let radius_sq = radius.mul(radius);
        let mut point = self.pcurve.evaluate(t_mid)?.sub(&origin);
        bisect(
            t_vertex,
            t_mid,
            EXIT_BISECTIONS,
            false,
            |t| -> GeopResult<_> {
                point = self.pcurve.evaluate(t)?.sub(&origin);
                // Inside the circle is the vertex's side.
                Ok(point
                    .norm_sq()
                    .sub(radius_sq)
                    .definite_sign()
                    .map(|out| !out))
            },
        )?;
        Ok(point)
    }
}

/// Bisection steps for [`Leg::exit`]. Bounds effort only: each step halves
/// the stretch of curve the exit point is known to lie on, and the point
/// returned is on the curve either way.
const EXIT_BISECTIONS: usize = 60;

/// Whether `d` lies strictly inside the angle swept counter-clockwise from
/// `from` to `to`; `None` where the intervals cannot decide it. `from` and
/// `to` pointing the same way sweep the full turn — the tip of a spur, whose
/// corner is everything around it.
fn within_corner(from: Vector2, to: Vector2, d: Vector2) -> Option<bool> {
    let sign = |x: Scalar| {
        if x.definitely_greater(Scalar::ZERO) {
            Some(true)
        } else if x.definitely_less(Scalar::ZERO) {
            Some(false)
        } else {
            None
        }
    };
    let after_from = sign(from.prod_cross(&d));
    let before_to = sign(d.prod_cross(&to));
    let turn = from.prod_cross(&to);
    if turn.definitely_greater(Scalar::ZERO) {
        // Less than half a turn: inside means past `from` and short of `to`.
        match (after_from, before_to) {
            (Some(false), _) | (_, Some(false)) => Some(false),
            (Some(true), Some(true)) => Some(true),
            _ => None,
        }
    } else if turn.definitely_less(Scalar::ZERO) {
        // More than half a turn: everything but the complementary angle.
        match (after_from, before_to) {
            (Some(true), _) | (_, Some(true)) => Some(true),
            (Some(false), Some(false)) => Some(false),
            _ => None,
        }
    } else if from.prod_dot(&to).definitely_greater(Scalar::ZERO) {
        // The full turn, unless `d` leaves along it.
        after_from.map(|_| true)
    } else {
        None
    }
}

/// Make `to` follow `from` in loop order, keeping `prev` consistent.
fn link(model: &mut Model, from: CoedgeId, to: CoedgeId) -> GeopResult<()> {
    model.get_coedge_mut(from)?.next = to;
    model.get_coedge_mut(to)?.prev = from;
    Ok(())
}

/// Insert `out`/`back` (an out-and-return pair on the same edge) into the
/// loop right after `after`, leaving the rest of the loop untouched.
fn splice_spur(
    model: &mut Model,
    after: CoedgeId,
    out: CoedgeId,
    back: CoedgeId,
) -> GeopResult<()> {
    let after_next = model.get_coedge(after)?.next;
    link(model, after, out)?;
    link(model, out, back)?;
    link(model, back, after_next)?;
    Ok(())
}

/// Marks the error raised when a splice would divide a face into a piece of
/// no area. Callers that are *imprinting* an edge they already know follows
/// this face's boundary can treat it as "already represented" rather than a
/// failure — see `booleans::remesh`'s tracing. Matched on the message because
/// `GeopError` carries no code; that is fragile enough to be worth the
/// constant rather than a literal at both ends.
pub const DEGENERATE_SPLIT: &str = "degenerate split";

/// Cut `face_id` in two along the newly inserted edge: `keep` and `moved`
/// anchor the two rings the old outer loop just split into.
///
/// `keep`'s ring stays with `face_id`; `moved`'s ring becomes the outer loop
/// of a brand-new face, and the holes are re-sorted between the two
/// ([`new_face_from_ring`]). Both rings bound material — that is what
/// distinguishes this from every other case in `splice_edge_into_face` — so
/// neither can become a hole of the other.
fn split_face(
    model: &mut Model,
    face_id: FaceId,
    keep: CoedgeId,
    moved: CoedgeId,
    max_nodes: usize,
    min_subdivision_size: Scalar,
) -> GeopResult<FaceId> {
    model
        .get_face_mut(face_id)?
        .set_boundary(BoundaryIndex::Outer, BoundaryType::Loop(keep));
    new_face_from_ring(model, face_id, moved, max_nodes, min_subdivision_size)
}

/// Split `ring` off `face_id` as the outer loop of a brand-new face on the
/// same surface and in the same shell.
///
/// `face_id`'s own holes are then re-sorted, since each now lies inside
/// exactly one of the two faces and the split has no idea which. Each is
/// classified by taking a point on it and asking `face_contains` — a hole the
/// new face contains moves there, everything else stays. A hole with no
/// usable point (a bare `Vertex` boundary, or a loop whose `(u, v)` cannot be
/// read) stays put rather than being dropped: leaving it on the wrong face is
/// a recoverable error, losing it silently fills in a hole that should exist.
fn new_face_from_ring(
    model: &mut Model,
    face_id: FaceId,
    ring: CoedgeId,
    max_nodes: usize,
    min_subdivision_size: Scalar,
) -> GeopResult<FaceId> {
    let old = model.get_face(face_id)?.clone();
    let new_face_id = model.insert_face(Face {
        surface: old.surface.clone(),
        outer: BoundaryType::Loop(ring),
        holes: Vec::new(),
        shell: old.shell,
    });
    model.get_shell_mut(old.shell)?.faces.push(new_face_id);

    // Every coedge of the moved ring now belongs to the new face. Walked with
    // a hard cap: a ring that never returns to its anchor would otherwise spin
    // here forever, and a corrupted `next` chain is exactly the kind of thing
    // a restructuring like this can introduce.
    for c in ring_coedges(model, ring)? {
        model.get_coedge_mut(c)?.face = new_face_id;
    }

    let mut stay = Vec::new();
    let mut move_over = Vec::new();
    for hole in old.holes {
        let Some(uv) = boundary_uv(model, hole)? else {
            stay.push(hole);
            continue;
        };
        let inside_new = matches!(
            face_contains(
                model,
                new_face_id,
                uv[0],
                uv[1],
                max_nodes,
                min_subdivision_size,
                HOLE_CLASSIFY_SEED,
            )?,
            PointClassification::Inside
        );
        if inside_new {
            move_over.push(hole);
        } else {
            stay.push(hole);
        }
    }
    for hole in &move_over {
        if let BoundaryType::Loop(anchor) = hole {
            for c in ring_coedges(model, *anchor)? {
                model.get_coedge_mut(c)?.face = new_face_id;
            }
        }
    }
    model.get_face_mut(face_id)?.holes = stay;
    model.get_face_mut(new_face_id)?.holes = move_over;
    Ok(new_face_id)
}

/// Which coedge of the ring anchored at `anchor` has the widest pcurve, at
/// `WIDTH_SAMPLES` points along each, and how wide: what an area too
/// uncertain to be told from zero comes from — a ring that really encloses
/// nothing, or a pcurve carrying too much uncertainty to say. Only for the
/// error message, never for a decision.
fn widest_pcurve(model: &Model, anchor: CoedgeId) -> GeopResult<String> {
    const WIDTH_SAMPLES: usize = 8;
    let width =
        |uv: &geop_core_math::vector::Vector2| uv[0].width().to_f64().max(uv[1].width().to_f64());
    let mut widest = (0.0, None);
    let mut each = Vec::new();
    for coedge_id in ring_coedges(model, anchor)? {
        let coedge = model.get_coedge(coedge_id)?;
        let (t0, t1) = coedge.pcurve.domain();
        let mut samples = Vec::with_capacity(WIDTH_SAMPLES);
        for i in 0..WIDTH_SAMPLES {
            let frac = Scalar::from_ratio(i as i64, (WIDTH_SAMPLES - 1) as i64)?;
            let uv = coedge.pcurve.evaluate(t0.add(t1.sub(t0).mul(frac)))?;
            if width(&uv) > widest.0 {
                widest = (width(&uv), Some((coedge_id, coedge.geometry, uv)));
            }
            samples.push(width(&uv));
        }
        let mid = samples[samples.len() / 2];
        each.push(format!(
            "{:?}: {:.1e} / {mid:.1e} / {:.1e}",
            coedge.geometry,
            samples[0],
            samples[samples.len() - 1]
        ));
    }
    let each = each.join(", ");
    Ok(match widest.1 {
        Some((coedge, geometry, uv)) => {
            // How fast the surface moves there: a slow parametrization
            // turns a narrow point into a wide `(u, v)`.
            let face = model.get_coedge(coedge)?.face;
            let surface_width = model
                .get_face(face)?
                .surface
                .control_points
                .iter()
                .flat_map(|cp| (0..4).map(move |k| cp[k].width().to_f64()))
                .fold(0.0, f64::max);
            let speed = model
                .get_face(face)?
                .surface
                .derivatives(uv[0].midpoint(), uv[1].midpoint())
                .map(|(du, dv)| format!("|dS/du| = {:?}, |dS/dv| = {:?}", du.norm(), dv.norm()))
                .unwrap_or_else(|e| format!("no derivatives: {}", e.root_message()));
            // Whether its ends land on the surface where its vertices are: a
            // pcurve pinned to a `(u, v)` off its vertex is widened to span
            // the gap all along.
            let coedge_data = model.get_coedge(coedge)?;
            let surface = &model.get_face(face)?.surface;
            let (t0, t1) = coedge_data.pcurve.domain();
            let mut ends = Vec::new();
            for (t, vertex) in [
                (t0, model.coedge_start_vertex(coedge)?),
                (t1, model.coedge_end_vertex(coedge)?),
            ] {
                let uv = coedge_data.pcurve.evaluate(t)?;
                let width = (0..3)
                    .map(|k| vertex.point[k].width().to_f64())
                    .fold(0.0, f64::max);
                ends.push(match surface.evaluate(uv[0].midpoint(), uv[1].midpoint()) {
                    Ok(p) => format!(
                        "{:.1e} (vertex {width:.1e} wide)",
                        p.sub(&vertex.point).norm().to_f64()
                    ),
                    Err(_) => {
                        format!("(outside the surface's domain at {uv:?}; vertex {width:.1e} wide)")
                    }
                });
            }
            let ends = ends.join(" and ");
            // And whether it inherits that from its edge's 3-D curve.
            let curve_width = match geometry {
                CoedgeGeometry::Edge(edge) => model
                    .get_edge(edge)?
                    .curve
                    .control_points
                    .iter()
                    .flat_map(|cp| (0..4).map(move |k| cp[k].width().to_f64()))
                    .fold(0.0, f64::max),
                CoedgeGeometry::Vertex(_) => 0.0,
            };
            format!(
                "coedge {coedge} on {geometry:?}, {:e} wide at uv={uv:?}, where {speed} and the surface's control points are up to {surface_width:e} wide; its ends land {ends} away from their vertices; its 3-D curve's control points up to {curve_width:e} wide. Every coedge's width at its start / middle / end: {each}",
                widest.0
            )
        }
        None => "every sample sharp".into(),
    })
}

/// Every coedge of the ring anchored at `anchor`, erroring rather than
/// spinning if the `next` chain never returns to it.
fn ring_coedges(model: &Model, anchor: CoedgeId) -> GeopResult<Vec<CoedgeId>> {
    let cap = model.coedges.len() + 1;
    let ring: Vec<CoedgeId> = model.iterate_loop_coedges(anchor).take(cap).collect();
    if ring.len() >= cap {
        return Err(GeopError::new(format!(
            "splice_edge_into_face: the ring anchored at coedge {anchor} never returns to its anchor"
        )));
    }
    Ok(ring)
}

/// Fixed seed for the ray casting behind hole classification — `face_contains`
/// retries until it finds a ray grazing nothing, so the answer is
/// seed-independent and a constant keeps a split reproducible run to run.
const HOLE_CLASSIFY_SEED: u64 = 0x1234_5678_9ABC_DEF0;

/// A `(u, v)` on `boundary`, for classifying which side of a split it falls
/// on. `None` for a bare `Vertex` boundary, which has no pcurve to read one
/// from.
fn boundary_uv(
    model: &Model,
    boundary: BoundaryType,
) -> GeopResult<Option<geop_core_math::vector::Vector2>> {
    let BoundaryType::Loop(anchor) = boundary else {
        return Ok(None);
    };
    let pcurve = &model.get_coedge(anchor)?.pcurve;
    let (t0, _) = pcurve.domain();
    Ok(Some(pcurve.evaluate(t0)?))
}