ogeom_heal/

divide.rs

//! Dividing a shape into smaller pieces: faces cut along lines of their own
//! parameters, edges cut at parameters of their own.
//!
//! Every piece keeps the curve or surface it was cut from, and so its
//! parameters: nothing is refitted, and a pcurve that held on the whole
//! holds on each piece. A cut face's new boundary is the surface's own
//! iso-curve, exact for every surface whose iso-curves have a closed form
//! (planes, drums, cones, balls, tori, extrusions, revolutions and
//! splines).
//!
//! The divisions on top pick where to cut: at the knots where a spline is
//! less smooth than asked ([`divide_by_continuity`]), wherever a face turns
//! further than an angle ([`divide_by_angle`]), until no face is larger than
//! an area ([`divide_by_area`]), and at every knot, each piece's geometry
//! then restated as the single Bézier span it covers ([`to_bezier`]).

use ogeom_core::{FastMap, FastSet};

use ogeom_algo::{Built, History, edge_vertices, surface_properties};
use ogeom_core::{OgeomResult, Tolerances, ogeom_bail};
use ogeom_geom::{
    BSpline2d, CircleCurve, Continuity, Curve, Curve2d as _, Curve3d as _, LineCurve, PlanarCurve,
    Surface as _, SurfaceGeometry, Transformable as _,
};
use ogeom_math::{Axis, Circle, Direction, Frame, KnotVector, Point, Point2, Transform};
use ogeom_mesh::Deflection;
use ogeom_topo::{
    EdgeData, EdgeRepr, Filter, Location, Model, NodeData, Orientation, Shape, ShapeType,
    SurfaceId, TShapeId, VertexData, explore, explore_unique,
};

use crate::Reshape;

/// A line of a face's parameters: `u = at` or `v = at`.
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum IsoLine {
    /// The line `u = at`, running along `v`.
    U(f64),
    /// The line `v = at`, running along `u`.
    V(f64),
}

impl IsoLine {
    const fn at(self) -> f64 {
        match self {
            Self::U(c) | Self::V(c) => c,
        }
    }

    /// The parameter the line holds fixed, read off a chart point.
    const fn fixed(self, p: Point2) -> f64 {
        match self {
            Self::U(_) => p.x,
            Self::V(_) => p.y,
        }
    }

    /// The parameter the line runs along.
    const fn free(self, p: Point2) -> f64 {
        match self {
            Self::U(_) => p.y,
            Self::V(_) => p.x,
        }
    }

    const fn point(self, free: f64) -> Point2 {
        match self {
            Self::U(c) => Point2::new(c, free),
            Self::V(c) => Point2::new(free, c),
        }
    }
}

/// Cut `face` of `shape` in two or more along `line`, and rebuild `shape`
/// around the pieces. The edges the line crosses are cut where it crosses
/// them, in every face that holds them.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if
/// the line does not cross the face, runs along one of its edges, or the
/// surface has no closed-form iso-curve.
pub fn divide_face(
    model: &mut Model,
    shape: &Shape,
    face: &Shape,
    line: IsoLine,
    tol: Tolerances,
) -> OgeomResult<Built> {
    let Some(reshape) = cut_face(model, face, line, tol)? else {
        ogeom_bail!(Construction, "the line {line:?} does not cross the face");
    };
    reshape.apply(model, shape)
}

/// Cut every edge and face of `shape` where its spline is less smooth than
/// `at_least`, so each piece is at least that smooth throughout. A knot is
/// judged by its multiplicity, so G1 and G2 ask what C1 and C2 ask.
///
/// Edges are cut at their curves' knots; faces along their surfaces' knot
/// lines, for a spline, an extrusion or revolution of a spline, and a
/// trimmed or offset surface over one.
///
/// # Errors
///
/// As [`divide_face`], for a cut that cannot be made.
pub fn divide_by_continuity(
    model: &mut Model,
    shape: &Shape,
    at_least: Continuity,
    tol: Tolerances,
) -> OgeomResult<Built> {
    let order = match at_least {
        Continuity::C0 => return Ok(Built::new(shape.clone(), History::identity())),
        Continuity::G1 | Continuity::C1 => 1,
        Continuity::G2 | Continuity::C2 => 2,
        Continuity::CInfinity => usize::MAX,
    };
    let start = placed_baked(model, shape, tol)?;
    let edges = divide_edges(model, &start.shape, order, tol)?;
    let faces = divide_faces(model, &edges.shape, tol, |_, _, surface, _| {
        Ok(knot_lines(surface, order))
    })?;
    Ok(Built::new(
        faces.shape,
        start.history.then(&edges.history).then(&faces.history),
    ))
}

/// Cut every face of `shape` whose angular parameters (a drum's, cone's,
/// ball's, torus's or revolution's turn) sweep more than `max_angle`, into
/// equal pieces sweeping no more than it. A closed face comes apart at its
/// seam into open ones.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if
/// `max_angle` is not a positive angle, and as [`divide_face`].
pub fn divide_by_angle(
    model: &mut Model,
    shape: &Shape,
    max_angle: f64,
    tol: Tolerances,
) -> OgeomResult<Built> {
    if !max_angle.is_finite() || max_angle <= tol.angular() {
        ogeom_bail!(Construction, "{max_angle} is not an angle to divide by");
    }
    let start = placed_baked(model, shape, tol)?;
    let faces = divide_faces(model, &start.shape, tol, |_, _, surface, bounds| {
        let (u_turns, v_turns) = angular(surface);
        let mut lines = Vec::new();
        for (turns, (lo, hi), iso) in [
            (u_turns, bounds.0, IsoLine::U as fn(f64) -> IsoLine),
            (v_turns, bounds.1, IsoLine::V),
        ] {
            let span = hi - lo;
            if turns && span > max_angle * (1.0 + 1e-9) {
                let pieces = (span / max_angle - 1e-6).ceil();
                lines.push(iso(lo + span / pieces));
            }
        }
        Ok(lines)
    })?;
    Ok(Built::new(faces.shape, start.history.then(&faces.history)))
}

/// Cut every face of `shape` larger than `max_area` in half, across its
/// longer extent, until none is.
///
/// # Errors
///
/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if
/// `max_area` is not a positive area, and as [`divide_face`].
pub fn divide_by_area(
    model: &mut Model,
    shape: &Shape,
    max_area: f64,
    tol: Tolerances,
) -> OgeomResult<Built> {
    if !max_area.is_finite() || max_area <= tol.confusion() * tol.confusion() {
        ogeom_bail!(Construction, "{max_area} is not an area to divide by");
    }
    let start = placed_baked(model, shape, tol)?;
    let faces = divide_faces(model, &start.shape, tol, |model, face, surface, bounds| {
        let area = surface_properties(model, face, Deflection::default(), tol)?.mass;
        if area <= max_area {
            return Ok(Vec::new());
        }
        let ((u0, u1), (v0, v1)) = bounds;
        let (um, vm) = (0.5 * (u0 + u1), 0.5 * (v0 + v1));
        // How long the face runs each way, along the lines through its
        // middle.
        let length = |from: (f64, f64), to: (f64, f64)| -> OgeomResult<f64> {
            let mut total = 0.0;
            let mut last = surface.point_at(from.0, from.1, tol)?;
            for i in 1..=16 {
                let s = f64::from(i) / 16.0;
                let p = surface.point_at(
                    from.0 + (to.0 - from.0) * s,
                    from.1 + (to.1 - from.1) * s,
                    tol,
                )?;
                total += p.distance(last);
                last = p;
            }
            Ok(total)
        };
        let across_u = length((u0, vm), (u1, vm))?;
        let across_v = length((um, v0), (um, v1))?;
        Ok(if across_u >= across_v {
            vec![IsoLine::U(um), IsoLine::V(vm)]
        } else {
            vec![IsoLine::V(vm), IsoLine::U(um)]
        })
    })?;
    Ok(Built::new(faces.shape, start.history.then(&faces.history)))
}

/// `shape` with every curve and surface a single Bézier span: converted to
/// splines, cut at every knot, and each piece's geometry restated as the
/// span it covers, parameters kept.
///
/// What `shape` shares with other shapes is copied before it is
/// worked on ([`ogeom_algo::on_own_nodes`]), which leaves them as they are.
///
/// # Errors
///
/// As [`ogeom_algo::to_nurbs`] and [`divide_by_continuity`].
pub fn to_bezier(model: &mut Model, shape: &Shape, tol: Tolerances) -> OgeomResult<Built> {
    ogeom_algo::on_own_nodes(model, shape, &[], |model, shape, _| {
        to_bezier_own(model, shape, tol)
    })
}

/// [`to_bezier`] on a shape every node below which is its own.
fn to_bezier_own(model: &mut Model, shape: &Shape, tol: Tolerances) -> OgeomResult<Built> {
    let nurbs = ogeom_algo::to_nurbs(model, shape, tol)?;
    let divided = divide_by_continuity(model, &nurbs.shape, Continuity::CInfinity, tol)?;

    // Each edge's curve restated as its span.
    let mut reshape = Reshape::new();
    for edge in explore_unique(model, &divided.shape, ShapeType::Edge)? {
        let data = edge_data(model, &edge)?;
        let Some(EdgeRepr::Curve3d { curve, range, .. }) = data.curve3d().cloned() else {
            continue;
        };
        let Some(Curve::BSpline(spline)) = model.geometry().curve(curve).cloned() else {
            continue;
        };
        if spline.knots().distinct().len() <= 2 {
            continue;
        }
        let span = model
            .geometry_mut()
            .add_curve(Curve::BSpline(spline.segment(range, tol)?));
        let mut fresh = data.clone();
        for repr in &mut fresh.representations {
            if let EdgeRepr::Curve3d { curve, .. } = repr {
                *curve = span;
            }
        }
        let bounds = model.children_of(&forward(&edge))?;
        let rebuilt = model.add_edge(fresh, &bounds)?;
        reshape.replace(&forward(&edge), rebuilt);
    }
    let spans = if reshape.is_empty() {
        Built::new(divided.shape.clone(), History::identity())
    } else {
        reshape.apply(model, &divided.shape)?
    };

    // Each face's surface restated as its patch, its edges' trims kept
    // under the patch's id: the patch keeps the parameters they are in.
    let mut reshape = Reshape::new();
    for face in explore_unique(model, &spans.shape, ShapeType::Face)? {
        let data = face_data(model, &face)?;
        let Some(SurfaceGeometry::BSpline(spline)) =
            model.geometry().surface(data.surface).cloned()
        else {
            continue;
        };
        if spline.u_knots().distinct().len() <= 2 && spline.v_knots().distinct().len() <= 2 {
            continue;
        }
        let rings = rings(model, &forward(&face), data.surface, tol)?;
        let ((u0, u1), (v0, v1)) = chart_bounds(&rings, tol);
        let ((du0, du1), (dv0, dv1)) = spline.domain();
        // A cut's trims stand on its knot line up to their fit; read a
        // hair off it, the patch would keep a sliver of the next span.
        let snap = |x: f64, knots: &[(f64, usize)]| {
            knots
                .iter()
                .map(|&(k, _)| k)
                .find(|k| (k - x).abs() <= tol.parametric())
                .unwrap_or(x)
        };
        let (us, vs) = (spline.u_knots().distinct(), spline.v_knots().distinct());
        let (u0, u1) = (snap(u0, &us), snap(u1, &us));
        let (v0, v1) = (snap(v0, &vs), snap(v1, &vs));
        let patch = spline.segment((u0.max(du0), u1.min(du1)), (v0.max(dv0), v1.min(dv1)), tol)?;
        let id = model
            .geometry_mut()
            .add_surface(SurfaceGeometry::BSpline(patch));
        // A seam of the whole surface that a cut left bounding this patch
        // once is no seam of the patch, which does not close across it: it
        // keeps only the side the face's ring walks, the one inside the
        // ring's bounds.
        let walked = explore(model, &forward(&face), Filter::OfType(ShapeType::Edge))?;
        let outside = |p: Point2| {
            (u0 - p.x)
                .max(p.x - u1)
                .max(0.0)
                .hypot((v0 - p.y).max(p.y - v1).max(0.0))
        };
        for edge in explore_unique(model, &face, ShapeType::Edge)? {
            let repr = edge_data(model, &edge)?
                .pcurve_for(data.surface, edge.location())
                .cloned();
            let Some(mut repr) = repr else { continue };
            let once = walked.iter().filter(|e| e.node() == edge.node()).count() == 1;
            if let (
                true,
                EdgeRepr::Seam {
                    forward: ahead,
                    reversed: back,
                    range,
                    location,
                    ..
                },
            ) = (once, &repr)
            {
                let middle = |id| -> OgeomResult<Point2> {
                    let Some(pcurve) = model.geometry().pcurve(id) else {
                        ogeom_bail!(Dangling, "pcurve is not in this model");
                    };
                    pcurve.point_at(0.5 * (range.0 + range.1), tol)
                };
                let curve = if outside(middle(*back)?) < outside(middle(*ahead)?) {
                    *back
                } else {
                    *ahead
                };
                repr = EdgeRepr::PCurve {
                    curve,
                    range: *range,
                    surface: data.surface,
                    location: location.clone(),
                };
            }
            match &mut repr {
                EdgeRepr::PCurve { surface, .. } | EdgeRepr::Seam { surface, .. } => *surface = id,
                _ => continue,
            }
            let Some(node) = model.node_mut(&edge) else {
                ogeom_bail!(Dangling, "edge is not in this model");
            };
            if let NodeData::Edge(e) = node.data_mut() {
                let agreed = e.same_parameter();
                e.add(repr);
                e.assert_same_parameter(agreed);
            }
        }
        let mut fresh = data.clone();
        fresh.surface = id;
        fresh.triangulation = None;
        let wires = model.children_of(&forward(&face))?;
        let rebuilt = model.add_face(fresh, &wires)?;
        reshape.replace(&forward(&face), rebuilt);
    }
    let patches = if reshape.is_empty() {
        Built::new(spans.shape.clone(), History::identity())
    } else {
        reshape.apply(model, &spans.shape)?
    };
    Ok(Built::new(
        patches.shape,
        nurbs
            .history
            .then(&divided.history)
            .then(&spans.history)
            .then(&patches.history),
    ))
}

// --- drivers ---------------------------------------------------------------

pub(crate) type Bounds = ((f64, f64), (f64, f64));

/// Cut faces until `lines` asks for no cut any face takes: each round asks
/// every face not yet settled for the lines it wants, in order, and makes
/// the first that crosses it.
fn divide_faces(
    model: &mut Model,
    shape: &Shape,
    tol: Tolerances,
    mut lines: impl FnMut(&Model, &Shape, &SurfaceGeometry, Bounds) -> OgeomResult<Vec<IsoLine>>,
) -> OgeomResult<Built> {
    let mut current = Built::new(shape.clone(), History::identity());
    let mut settled: FastSet<TShapeId> = FastSet::default();
    'rounds: for _ in 0..100_000 {
        for face in explore_unique(model, &current.shape, ShapeType::Face)? {
            if settled.contains(&face.node()) {
                continue;
            }
            let data = face_data(model, &face)?;
            if data.natural_restriction {
                // Bounded first, by its chart's own sides: then it cuts as
                // any face does.
                match bounded_natural(model, &face, tol)? {
                    Some(bounded) => {
                        let mut reshape = Reshape::new();
                        reshape.replace(&forward(&face), bounded);
                        let next = reshape.apply(model, &current.shape)?;
                        current = Built::new(next.shape, current.history.then(&next.history));
                        continue 'rounds;
                    }
                    None => {
                        settled.insert(face.node());
                        continue;
                    }
                }
            }
            let Some(surface) = model.geometry().surface(data.surface).cloned() else {
                ogeom_bail!(Dangling, "surface is not in this model");
            };
            let rings = rings(model, &forward(&face), data.surface, tol)?;
            let bounds = chart_bounds(&rings, tol);
            for line in lines(&*model, &face, &surface, bounds)? {
                let (lo, hi) = match line {
                    IsoLine::U(_) => bounds.0,
                    IsoLine::V(_) => bounds.1,
                };
                let margin = (hi - lo).abs() * 1e-7 + tol.parametric();
                if line.at() <= lo + margin || line.at() >= hi - margin {
                    continue;
                }
                if let Some(reshape) = cut_face(model, &face, line, tol)? {
                    let next = reshape.apply(model, &current.shape)?;
                    current = Built::new(next.shape, current.history.then(&next.history));
                    continue 'rounds;
                }
            }
            settled.insert(face.node());
        }
        return Ok(current);
    }
    ogeom_bail!(Construction, "the division did not settle");
}

/// A face covering its surface's whole chart, rebuilt with that chart's
/// sides as its boundary: a periodic direction's two sides one seam edge,
/// a side the surface pinches to a point (a ball's pole) a degenerate edge.
/// `None` where the chart is unbounded or a side has no closed-form curve.
pub(crate) fn bounded_natural(
    model: &mut Model,
    face: &Shape,
    tol: Tolerances,
) -> OgeomResult<Option<Shape>> {
    let face_fwd = forward(face);
    let data = face_data(model, &face_fwd)?;
    let Some(surface) = model.geometry().surface(data.surface).cloned() else {
        ogeom_bail!(Dangling, "surface is not in this model");
    };
    let ((u0, u1), (v0, v1)) = surface.domain();
    if ![u0, u1, v0, v1]
        .iter()
        .all(|x| x.is_finite() && x.abs() < 1e6)
    {
        return Ok(None);
    }
    let (pu, pv) = (surface.is_periodic_u(), surface.is_periodic_v());
    let corner = |u: f64, v: f64| surface.point_at(u, v, tol);
    // Each side as (line, from, to) in the free parameter, walked so the
    // ring runs counter-clockwise round the chart.
    let sides = [
        (IsoLine::V(v0), u0, u1),
        (IsoLine::U(u1), v0, v1),
        (IsoLine::V(v1), u1, u0),
        (IsoLine::U(u0), v1, v0),
    ];
    // Vertices by chart corner, one where the corners meet in space.
    let corners = [(u0, v0), (u1, v0), (u1, v1), (u0, v1)];
    let mut vertices: Vec<Shape> = Vec::with_capacity(4);
    for (i, (u, v)) in corners.iter().enumerate() {
        let p = corner(*u, *v)?;
        let found = (0..i).find(|j| {
            let (a, b) = corners[*j];
            corner(a, b).is_ok_and(|q| q.distance(p) <= tol.confusion())
        });
        vertices.push(match found {
            Some(j) => vertices[j].clone(),
            None => model.add_vertex(VertexData::new(p)),
        });
    }
    let chart_line = |line: IsoLine, a: f64, b: f64, t: (f64, f64)| -> OgeomResult<PlanarCurve> {
        let knots = KnotVector::new(vec![t.0, t.0, t.1, t.1], 1)?;
        let (p, q) = (line.point(a), line.point(b));
        let (p, q) = if t.0 <= t.1 { (p, q) } else { (q, p) };
        Ok(PlanarCurve::BSpline(BSpline2d::new(
            knots,
            vec![p, q],
            tol,
        )?))
    };
    // One edge per side, a periodic pair shared; each side's occurrence.
    let mut ring: Vec<Shape> = Vec::with_capacity(4);
    let mut seam_u: Option<Shape> = None;
    let mut seam_v: Option<Shape> = None;
    for (index, (line, a, b)) in sides.into_iter().enumerate() {
        let (from, to) = (vertices[index].clone(), vertices[(index + 1) % 4].clone());
        let periodic_pair = match line {
            IsoLine::U(_) => pu,
            IsoLine::V(_) => pv,
        };
        // The second side of a periodic pair is the first's seam, walked
        // back.
        if periodic_pair {
            let held = match line {
                IsoLine::U(_) => &seam_u,
                IsoLine::V(_) => &seam_v,
            };
            if let Some(seam) = held {
                ring.push(seam.reversed());
                continue;
            }
        }
        let (lo, hi) = (a.min(b), a.max(b));
        let mid = line.point(0.5 * (lo + hi));
        let pinched = corner(line.point(lo).x, line.point(lo).y)?.distance(corner(mid.x, mid.y)?)
            <= tol.confusion()
            && corner(line.point(hi).x, line.point(hi).y)?.distance(corner(mid.x, mid.y)?)
                <= tol.confusion();
        let edge = if pinched {
            let mut edge_data = EdgeData::new();
            edge_data.degenerate = true;
            let trim = chart_line(line, lo, hi, (lo, hi))?;
            let id = model.geometry_mut().add_pcurve(trim);
            edge_data.add(EdgeRepr::PCurve {
                curve: id,
                surface: data.surface,
                location: Location::identity(),
                range: (lo, hi),
            });
            let pole = model.add_edge(edge_data, &[from.clone(), from.clone()])?;
            // Its trim runs up the chart; the side may walk it down.
            if a <= b { pole } else { pole.reversed() }
        } else {
            let Some((curve, scale, offset)) =
                iso_curve(&surface, line, &[0.5 * (lo + hi)], (lo, hi), tol)?
            else {
                return Ok(None);
            };
            let (t_lo, t_hi) = (scale * lo + offset, scale * hi + offset);
            let (range, rising) = if t_lo <= t_hi {
                ((t_lo, t_hi), true)
            } else {
                ((t_hi, t_lo), false)
            };
            let curve_id = model.geometry_mut().add_curve(curve);
            let mut edge_data = EdgeData::on_curve(curve_id, Location::identity(), range);
            // The edge runs with its curve; the side walks from `a` to `b`.
            let (start_v, end_v) = {
                let low_end = if a <= b { &from } else { &to };
                let high_end = if a <= b { &to } else { &from };
                if rising {
                    (low_end.clone(), high_end.clone())
                } else {
                    (high_end.clone(), low_end.clone())
                }
            };
            let pcurve_at = |line: IsoLine| -> OgeomResult<PlanarCurve> {
                let (p, q) = if rising {
                    (line.point(lo), line.point(hi))
                } else {
                    (line.point(hi), line.point(lo))
                };
                let knots = KnotVector::new(vec![range.0, range.0, range.1, range.1], 1)?;
                Ok(PlanarCurve::BSpline(BSpline2d::new(
                    knots,
                    vec![p, q],
                    tol,
                )?))
            };
            if periodic_pair {
                // Its two chart sides: this one, and the one a period on.
                let other = match line {
                    IsoLine::U(c) => IsoLine::U(if (c - u1).abs() < (c - u0).abs() {
                        u0
                    } else {
                        u1
                    }),
                    IsoLine::V(c) => IsoLine::V(if (c - v1).abs() < (c - v0).abs() {
                        v0
                    } else {
                        v1
                    }),
                };
                let here = model.geometry_mut().add_pcurve(pcurve_at(line)?);
                let there = model.geometry_mut().add_pcurve(pcurve_at(other)?);
                edge_data.add(EdgeRepr::Seam {
                    forward: here,
                    reversed: there,
                    surface: data.surface,
                    location: Location::identity(),
                    range,
                });
            } else {
                let id = model.geometry_mut().add_pcurve(pcurve_at(line)?);
                edge_data.add(EdgeRepr::PCurve {
                    curve: id,
                    surface: data.surface,
                    location: Location::identity(),
                    range,
                });
            }
            edge_data.assert_same_parameter(true);
            let edge = model.add_edge(edge_data, &[start_v, end_v])?;
            let walked_with_curve = (a <= b) == rising;
            let occurrence = if walked_with_curve {
                edge.clone()
            } else {
                edge.reversed()
            };
            if periodic_pair {
                match line {
                    IsoLine::U(_) => seam_u = Some(occurrence.clone()),
                    IsoLine::V(_) => seam_v = Some(occurrence.clone()),
                }
            }
            ring.push(occurrence);
            continue;
        };
        ring.push(edge);
    }
    let wire = model.add_wire(&ring)?;
    let mut fresh = data.clone();
    fresh.natural_restriction = false;
    fresh.triangulation = None;
    let bounded = model.add_face(fresh, &[wire])?;
    Ok(Some(bounded))
}

/// Cut every edge whose spline is less than `order` times differentiable at
/// an interior knot, there.
fn divide_edges(
    model: &mut Model,
    shape: &Shape,
    order: usize,
    tol: Tolerances,
) -> OgeomResult<Built> {
    let mut reshape = Reshape::new();
    for edge in explore_unique(model, shape, ShapeType::Edge)? {
        let edge = forward(&edge);
        let data = edge_data(model, &edge)?;
        let Some(EdgeRepr::Curve3d { curve, range, .. }) = data.curve3d().cloned() else {
            continue;
        };
        let Some(geometry) = model.geometry().curve(curve) else {
            ogeom_bail!(Dangling, "curve is not in this model");
        };
        let margin = (range.1 - range.0) * 1e-9 + tol.parametric();
        let at: Vec<f64> = curve_knots(geometry)
            .map(|knots| weak_knots(&knots, order))
            .unwrap_or_default()
            .into_iter()
            .filter(|t| *t > range.0 + margin && *t < range.1 - margin)
            .collect();
        if at.is_empty() {
            continue;
        }
        let fractions: Vec<f64> = at
            .iter()
            .map(|t| (t - range.0) / (range.1 - range.0))
            .collect();
        let (pieces, _) = split_edge(model, &edge, &fractions, None, tol)?;
        reshape.split(&edge, pieces);
    }
    if reshape.is_empty() {
        return Ok(Built::new(shape.clone(), History::identity()));
    }
    reshape.apply(model, shape)
}

/// The shape, its placements baked in where any edge or face is placed:
/// cuts make new nodes in the chart of the node they cut, and a node placed
/// twice has two.
pub(crate) fn placed_baked(
    model: &mut Model,
    shape: &Shape,
    tol: Tolerances,
) -> OgeomResult<Built> {
    let mut placed = false;
    for kind in [ShapeType::Edge, ShapeType::Face] {
        for s in explore_unique(model, shape, kind)? {
            placed |= !s.location().is_identity();
        }
    }
    for face in explore_unique(model, shape, ShapeType::Face)? {
        placed |= !face_data(model, &face)?.location.is_identity();
    }
    if placed {
        ogeom_algo::baked_shape(model, shape, tol)
    } else {
        Ok(Built::new(shape.clone(), History::identity()))
    }
}

// --- where to cut -----------------------------------------------------------

fn curve_knots(curve: &Curve) -> Option<KnotVector> {
    match curve {
        Curve::BSpline(b) => Some(b.knots().clone()),
        Curve::Trimmed(t) => curve_knots(t.basis()),
        _ => None,
    }
}

/// Interior knots where the spline is less than `order` times
/// differentiable.
fn weak_knots(knots: &KnotVector, order: usize) -> Vec<f64> {
    let (a, b) = knots.domain();
    let p = knots.degree();
    knots
        .distinct()
        .into_iter()
        .filter(|(value, m)| *value > a && *value < b && p.saturating_sub(*m) < order)
        .map(|(value, _)| value)
        .collect()
}

/// The knot lines of a surface where it is less than `order` times
/// differentiable.
fn knot_lines(surface: &SurfaceGeometry, order: usize) -> Vec<IsoLine> {
    match surface {
        SurfaceGeometry::BSpline(b) => weak_knots(b.u_knots(), order)
            .into_iter()
            .map(IsoLine::U)
            .chain(weak_knots(b.v_knots(), order).into_iter().map(IsoLine::V))
            .collect(),
        SurfaceGeometry::Extrusion(e) => curve_knots(e.curve())
            .map(|k| weak_knots(&k, order).into_iter().map(IsoLine::U).collect())
            .unwrap_or_default(),
        SurfaceGeometry::Revolution(r) => curve_knots(r.curve())
            .map(|k| weak_knots(&k, order).into_iter().map(IsoLine::V).collect())
            .unwrap_or_default(),
        SurfaceGeometry::Trimmed(t) => knot_lines(t.basis(), order),
        // An offset is one order rougher than its basis at a knot.
        SurfaceGeometry::Offset(o) => knot_lines(o.basis(), order.saturating_add(1)),
        _ => Vec::new(),
    }
}

/// Whether each parameter of a surface is an angle.
const fn angular(surface: &SurfaceGeometry) -> (bool, bool) {
    match surface {
        SurfaceGeometry::Cylinder(_)
        | SurfaceGeometry::Cone(_)
        | SurfaceGeometry::Revolution(_) => (true, false),
        SurfaceGeometry::Sphere(_) | SurfaceGeometry::Torus(_) => (true, true),
        _ => (false, false),
    }
}

// --- one cut ------------------------------------------------------------------

/// One traversal of an edge around a face, in the chart.
#[derive(Debug, Clone)]
pub(crate) struct Occurrence {
    pub(crate) edge: Shape,
    pcurve: PlanarCurve,
    range: (f64, f64),
}

impl Occurrence {
    pub(crate) fn reversed(&self) -> bool {
        self.edge.orientation() == Orientation::Reversed
    }

    /// The chart point a fraction `s` of the way along the traversal.
    pub(crate) fn at(&self, s: f64, tol: Tolerances) -> OgeomResult<Point2> {
        let s = if self.reversed() { 1.0 - s } else { s };
        self.pcurve
            .point_at(self.range.0 + (self.range.1 - self.range.0) * s, tol)
    }

    pub(crate) fn polyline(&self, samples: usize, tol: Tolerances) -> OgeomResult<Vec<Point2>> {
        (0..=samples)
            .map(|i| {
                #[allow(clippy::cast_precision_loss)]
                let s = i as f64 / samples as f64;
                self.at(s, tol)
            })
            .collect()
    }
}

/// An edge's trim on a face: one pcurve, or a seam's two, over a range.
type Trims = (PlanarCurve, Option<PlanarCurve>, (f64, f64));

/// A ring's edges, its chart outline, and the holes it holds.
type Outer = (Vec<Shape>, Vec<Point2>, Vec<Vec<Shape>>);

/// A face's rings as traversals, each seam traversal on the side of the
/// chart its ring continues on.
pub(crate) fn rings(
    model: &Model,
    face: &Shape,
    surface: SurfaceId,
    tol: Tolerances,
) -> OgeomResult<Vec<Vec<Occurrence>>> {
    let mut out = Vec::new();
    for wire in model.ordered_children_of(face)? {
        let edges = model.ordered_children_of(&wire)?;
        let mut choices: Vec<Option<Trims>> = Vec::new();
        for edge in &edges {
            let data = edge_data(model, edge)?;
            let pcurve =
                |id| -> OgeomResult<PlanarCurve> {
                    model.geometry().pcurve(id).cloned().ok_or_else(|| {
                        ogeom_core::ogeom_err!(Dangling, "pcurve is not in this model")
                    })
                };
            choices.push(match data.pcurve_for(surface, edge.location()) {
                Some(EdgeRepr::PCurve { curve, range, .. }) => {
                    Some((pcurve(*curve)?, None, *range))
                }
                Some(EdgeRepr::Seam {
                    forward,
                    reversed,
                    range,
                    ..
                }) => Some((pcurve(*forward)?, Some(pcurve(*reversed)?), *range)),
                _ => None,
            });
        }
        // Walked from an edge off the seam where there is one. A ring of
        // seams alone (a torus's) has no such edge, and where its seams are
        // cut into pieces every corner is one vertex, so the side its first
        // occurrence takes is tried both ways: the right one closes the
        // ring in the chart.
        let first_plain = choices
            .iter()
            .position(|c| c.as_ref().is_some_and(|c| c.1.is_none()));
        let start = first_plain.unwrap_or(0);
        let walk = |flip: bool| -> OgeomResult<(Vec<Occurrence>, f64)> {
            let n = edges.len();
            let mut ring: Vec<Option<Occurrence>> = vec![None; n];
            let mut last: Option<Point2> = None;
            let mut gap = 0.0;
            for k in 0..n {
                let i = (start + k) % n;
                let Some((a, b, range)) = choices[i].clone() else {
                    ogeom_bail!(Construction, "an edge has no trim on this face");
                };
                let mut occurrence = Occurrence {
                    edge: edges[i].clone(),
                    pcurve: a,
                    range,
                };
                if let Some(b) = b {
                    let other = Occurrence {
                        pcurve: b,
                        ..occurrence.clone()
                    };
                    let take_other = match last {
                        Some(last) => {
                            other.at(0.0, tol)?.distance(last)
                                < occurrence.at(0.0, tol)?.distance(last)
                        }
                        None => occurrence.reversed() != flip,
                    };
                    if take_other {
                        occurrence = other;
                    }
                }
                if let Some(last) = last {
                    gap += occurrence.at(0.0, tol)?.distance(last);
                }
                last = Some(occurrence.at(1.0, tol)?);
                ring[i] = Some(occurrence);
            }
            let ring: Vec<Occurrence> = ring.into_iter().flatten().collect();
            if let (Some(first), Some(last)) = (ring.first(), last) {
                gap += first.at(0.0, tol)?.distance(last);
            }
            Ok((ring, gap))
        };
        let (mut ring, gap) = walk(false)?;
        if first_plain.is_none() && gap > tol.parametric() {
            let (flipped, flipped_gap) = walk(true)?;
            if flipped_gap < gap {
                ring = flipped;
            }
        }
        out.push(ring);
    }
    Ok(out)
}

pub(crate) fn chart_bounds(rings: &[Vec<Occurrence>], tol: Tolerances) -> Bounds {
    let mut b = (
        (f64::INFINITY, f64::NEG_INFINITY),
        (f64::INFINITY, f64::NEG_INFINITY),
    );
    for o in rings.iter().flatten() {
        for i in 0..=32 {
            if let Ok(p) = o.at(f64::from(i) / 32.0, tol) {
                b.0 = (b.0.0.min(p.x), b.0.1.max(p.x));
                b.1 = (b.1.0.min(p.y), b.1.1.max(p.y));
            }
        }
    }
    b
}

/// Where a line meets a face's boundary.
#[derive(Debug, Clone, Copy)]
enum Meeting {
    /// At an existing vertex.
    Vertex(TShapeId),
    /// Inside an edge, a fraction of the way along it (forward).
    Inside(TShapeId, u64),
}

/// The substitutions that cut `face` along `line`, or `None` where the
/// line does not cross it.
#[allow(clippy::too_many_lines)]
fn cut_face(
    model: &mut Model,
    face: &Shape,
    line: IsoLine,
    tol: Tolerances,
) -> OgeomResult<Option<Reshape>> {
    let face_fwd = forward(face);
    let data = face_data(model, &face_fwd)?;
    let Some(surface) = model.geometry().surface(data.surface).cloned() else {
        ogeom_bail!(Dangling, "surface is not in this model");
    };
    let rings = rings(model, &face_fwd, data.surface, tol)?;
    let bounds = chart_bounds(&rings, tol);
    let span = (bounds.0.1 - bounds.0.0).max(bounds.1.1 - bounds.1.0);
    let eps = span * 1e-9 + tol.parametric();
    let c = line.at();

    // Where the line crosses each edge, as fractions of the edge's forward
    // range; and every point where it meets the boundary, along the line.
    let mut cuts: FastMap<TShapeId, Vec<f64>> = FastMap::default();
    let mut meetings: Vec<(f64, Meeting)> = Vec::new();
    for occurrence in rings.iter().flatten() {
        let f = |s: f64| -> OgeomResult<f64> { Ok(line.fixed(occurrence.at(s, tol)?) - c) };
        const N: usize = 64;
        #[allow(clippy::cast_precision_loss)]
        let values: Vec<f64> = (0..=N)
            .map(|i| f(i as f64 / N as f64))
            .collect::<OgeomResult<_>>()?;
        if values.iter().all(|v| v.abs() <= eps) {
            ogeom_bail!(Construction, "an edge of the face runs along {line:?}");
        }
        let end = occurrence.at(1.0, tol)?;
        if (line.fixed(end) - c).abs() <= eps {
            let Some((_, v)) = edge_vertices(model, &occurrence.edge)? else {
                continue;
            };
            meetings.push((line.free(end), Meeting::Vertex(v.node())));
        }
        for i in 0..N {
            let (a, b) = (values[i], values[i + 1]);
            if a.abs() <= eps || b.abs() <= eps || a.signum() == b.signum() {
                // Through a sample exactly: a crossing there only where the
                // sign changes across it.
                if i + 1 < N
                    && b.abs() <= eps
                    && a.abs() > eps
                    && values[i + 2].abs() > eps
                    && a.signum() != values[i + 2].signum()
                {
                    #[allow(clippy::cast_precision_loss)]
                    let s = (i + 1) as f64 / N as f64;
                    record(
                        model,
                        &surface,
                        occurrence,
                        s,
                        line,
                        &mut cuts,
                        &mut meetings,
                        tol,
                    )?;
                }
                continue;
            }
            #[allow(clippy::cast_precision_loss)]
            let (mut lo, mut hi) = (i as f64 / N as f64, (i + 1) as f64 / N as f64);
            let mut f_lo = a;
            for _ in 0..80 {
                let mid = 0.5 * (lo + hi);
                let f_mid = f(mid)?;
                if f_mid.signum() == f_lo.signum() {
                    (lo, f_lo) = (mid, f_mid);
                } else {
                    hi = mid;
                }
            }
            record(
                model,
                &surface,
                occurrence,
                0.5 * (lo + hi),
                line,
                &mut cuts,
                &mut meetings,
                tol,
            )?;
        }
    }
    meetings.sort_by(|a, b| a.0.total_cmp(&b.0));
    meetings.dedup_by(|b, a| {
        (a.0 - b.0).abs() <= eps
            && match (a.1, b.1) {
                (Meeting::Vertex(x), Meeting::Vertex(y)) => x == y,
                (Meeting::Inside(x, s), Meeting::Inside(y, t)) => x == y && s == t,
                _ => false,
            }
    });

    // The stretches of the line inside the face.
    let outlines: Vec<Vec<Point2>> = rings
        .iter()
        .map(|ring| -> OgeomResult<Vec<Point2>> {
            let mut out = Vec::new();
            for o in ring {
                out.extend(o.polyline(32, tol)?);
            }
            Ok(out)
        })
        .collect::<OgeomResult<_>>()?;
    let mut stretches: Vec<(usize, usize)> = Vec::new();
    for i in 0..meetings.len().saturating_sub(1) {
        let (w0, w1) = (meetings[i].0, meetings[i + 1].0);
        if w1 - w0 <= eps * 10.0 {
            continue;
        }
        if inside(&outlines, line.point(0.5 * (w0 + w1))) {
            stretches.push((i, i + 1));
        }
    }
    if stretches.is_empty() {
        return Ok(None);
    }
    let outer_sign = outlines
        .iter()
        .map(|o| signed_area(o))
        .max_by(|a, b| a.abs().total_cmp(&b.abs()))
        .unwrap_or(1.0)
        .signum();

    // Cut the edges, and every meeting inside one becomes its new vertex.
    let mut reshape = Reshape::new();
    let mut pieces_of: FastMap<TShapeId, Vec<(Shape, f64, f64)>> = FastMap::default();
    let mut vertex_at: FastMap<(TShapeId, u64), Shape> = FastMap::default();
    let mut vertex_of_node: FastMap<TShapeId, Shape> = FastMap::default();
    for v in explore_unique(model, &face_fwd, ShapeType::Vertex)? {
        vertex_of_node.insert(v.node(), v);
    }
    for occurrence in rings.iter().flatten() {
        let edge = forward(&occurrence.edge);
        if pieces_of.contains_key(&edge.node()) {
            continue;
        }
        let Some(fractions) = cuts.get(&edge.node()) else {
            continue;
        };
        let mut fractions = fractions.clone();
        fractions.sort_by(f64::total_cmp);
        fractions.dedup_by(|b, a| (*a - *b).abs() <= 1e-12);
        let (pieces, vertices) = split_edge(model, &edge, &fractions, Some(data.surface), tol)?;
        for (s, v) in fractions.iter().zip(&vertices) {
            vertex_at.insert((edge.node(), s.to_bits()), v.clone());
        }
        let mut bounds = vec![0.0];
        bounds.extend(&fractions);
        bounds.push(1.0);
        pieces_of.insert(
            edge.node(),
            pieces
                .iter()
                .zip(bounds.windows(2))
                .map(|(p, w)| (p.clone(), w[0], w[1]))
                .collect(),
        );
        reshape.split(&edge, pieces);
    }
    let meeting_vertex = |m: Meeting| -> OgeomResult<Shape> {
        match m {
            Meeting::Vertex(node) => vertex_of_node
                .get(&node)
                .cloned()
                .ok_or_else(|| ogeom_core::ogeom_err!(Construction, "a meeting vertex is lost")),
            Meeting::Inside(edge, s) => vertex_at
                .get(&(edge, s))
                .cloned()
                .ok_or_else(|| ogeom_core::ogeom_err!(Construction, "a cut vertex is lost")),
        }
    };

    // The boundary's pieces, each on its side of the line.
    let mut sides: [Vec<Item>; 2] = [Vec::new(), Vec::new()];
    for occurrence in rings.iter().flatten() {
        let edge = forward(&occurrence.edge);
        let parts: Vec<(Shape, f64, f64)> = match pieces_of.get(&edge.node()) {
            Some(parts) => parts.clone(),
            None => vec![(edge.clone(), 0.0, 1.0)],
        };
        let mut ordered: Vec<(Shape, f64, f64)> = if occurrence.reversed() {
            parts
                .into_iter()
                .rev()
                .map(|(p, a, b)| (p.reversed(), 1.0 - b, 1.0 - a))
                .collect()
        } else {
            parts
        };
        for (piece, a, b) in ordered.drain(..) {
            let sub = |s: f64| occurrence.at(a + (b - a) * s, tol);
            let polyline: Vec<Point2> = (0..=16)
                .map(|i| sub(f64::from(i) / 16.0))
                .collect::<OgeomResult<_>>()?;
            let side = usize::from(line.fixed(polyline[8]) > c);
            let Some((from, to)) = edge_vertices(model, &piece)? else {
                ogeom_bail!(Construction, "a boundary piece has no vertices");
            };
            sides[side].push(Item {
                edge: piece,
                from: from.node(),
                to: to.node(),
                polyline,
            });
        }
    }

    // The new edges along the line, each walked one way by each side.
    let probes: Vec<f64> = stretches
        .iter()
        .map(|&(i, j)| 0.5 * (meetings[i].0 + meetings[j].0))
        .collect();
    let span = (
        meetings.first().map_or(0.0, |m| m.0),
        meetings.last().map_or(0.0, |m| m.0),
    );
    let Some((curve, scale, offset)) = iso_curve(&surface, line, &probes, span, tol)? else {
        ogeom_bail!(
            Construction,
            "the surface has no closed-form iso-curve to cut along"
        );
    };
    for (i, j) in stretches {
        let (w0, w1) = (meetings[i].0, meetings[j].0);
        let (a, b) = (
            meeting_vertex(meetings[i].1)?,
            meeting_vertex(meetings[j].1)?,
        );
        let (t0, t1) = (scale * w0 + offset, scale * w1 + offset);
        let rising = t1 > t0;
        let (range, (start, end), (p0, p1)) = if rising {
            ((t0, t1), (&a, &b), (line.point(w0), line.point(w1)))
        } else {
            ((t1, t0), (&b, &a), (line.point(w1), line.point(w0)))
        };
        let curve_id = model.geometry_mut().add_curve(curve.clone());
        let mut edge_data = EdgeData::on_curve(curve_id, Location::identity(), range);
        let trim = BSpline2d::new(
            KnotVector::new(vec![range.0, range.0, range.1, range.1], 1)?,
            vec![p0, p1],
            tol,
        )?;
        let trim_id = model.geometry_mut().add_pcurve(PlanarCurve::BSpline(trim));
        edge_data.add(EdgeRepr::PCurve {
            curve: trim_id,
            surface: data.surface,
            location: Location::identity(),
            range,
        });
        edge_data.assert_same_parameter(true);
        // The vertices were placed by the edges they cut; the new edge's
        // ends reach them within the gap.
        let reach = [(start, range.0), (end, range.1)]
            .iter()
            .map(|(v, t)| -> OgeomResult<f64> {
                Ok(vertex_point(model, v)?.distance(curve.point_at(*t, tol)?))
            })
            .collect::<OgeomResult<Vec<f64>>>()?
            .into_iter()
            .fold(0.0_f64, f64::max);
        let edge = model.add_edge(edge_data, &[start.clone(), end.clone()])?;
        model.widen(&edge, ogeom_core::Tolerance::new(reach + tol.confusion())?)?;
        let polyline: Vec<Point2> = (0..=16)
            .map(|k| {
                let w = w0 + (w1 - w0) * f64::from(k) / 16.0;
                line.point(w)
            })
            .collect();
        // Walked up the line (increasing free parameter), the material on
        // the left is the low side of a `u` line and the high side of a `v`
        // line, for a boundary wound counter-clockwise.
        let left = match line {
            IsoLine::U(_) => 0,
            IsoLine::V(_) => 1,
        };
        let up_side = if outer_sign >= 0.0 { left } else { 1 - left };
        let up = if rising {
            edge.clone()
        } else {
            edge.reversed()
        };
        let down = up.reversed();
        sides[up_side].push(Item {
            edge: up,
            from: a.node(),
            to: b.node(),
            polyline: polyline.clone(),
        });
        sides[1 - up_side].push(Item {
            edge: down,
            from: b.node(),
            to: a.node(),
            polyline: polyline.into_iter().rev().collect(),
        });
    }

    // Each side's pieces chained into rings, and the rings into faces.
    let mut faces: Vec<Shape> = Vec::new();
    for items in sides {
        let loops = chain(items, eps.max(tol.parametric() * 1e3))?;
        let mut outers: Vec<Outer> = Vec::new();
        let mut holes: Vec<(Vec<Shape>, Vec<Point2>)> = Vec::new();
        for (edges, outline) in loops {
            if signed_area(&outline) * outer_sign > 0.0 {
                outers.push((edges, outline, Vec::new()));
            } else {
                holes.push((edges, outline));
            }
        }
        for (edges, outline) in holes {
            let probe = outline[outline.len() / 2];
            let Some(host) = outers
                .iter_mut()
                .filter(|o| inside(std::slice::from_ref(&o.1), probe))
                .min_by(|a, b| signed_area(&a.1).abs().total_cmp(&signed_area(&b.1).abs()))
            else {
                ogeom_bail!(Construction, "a hole of the cut face lies in no piece");
            };
            host.2.push(edges);
        }
        for (edges, _, hole_edges) in outers {
            let mut wires = vec![model.add_wire(&edges)?];
            for h in hole_edges {
                wires.push(model.add_wire(&h)?);
            }
            let mut fresh = data.clone();
            fresh.triangulation = None;
            fresh.natural_restriction = false;
            faces.push(model.add_face(fresh, &wires)?);
        }
    }
    if faces.len() < 2 {
        return Ok(None);
    }
    reshape.split(&face_fwd, faces);
    Ok(Some(reshape))
}

/// Record a crossing a fraction `s` along `occurrence`. A crossing that
/// lands within an end vertex's radius of it, lifted through the face's
/// surface, meets that vertex: the line passes through it as the face
/// states it, and cutting the edge there would leave a sliver shorter
/// than the vertex is wide. A degenerate edge, every point of which is its
/// vertex, is cut where the line crosses its row.
#[allow(clippy::too_many_arguments)]
fn record(
    model: &Model,
    surface: &SurfaceGeometry,
    occurrence: &Occurrence,
    s: f64,
    line: IsoLine,
    cuts: &mut FastMap<TShapeId, Vec<f64>>,
    meetings: &mut Vec<(f64, Meeting)>,
    tol: Tolerances,
) -> OgeomResult<()> {
    let at = occurrence.at(s, tol)?;
    let degenerate = model
        .node(&occurrence.edge)
        .and_then(|n| n.data().as_edge())
        .is_none_or(|d| d.degenerate);
    if !degenerate && let Some((first, last)) = edge_vertices(model, &occurrence.edge)? {
        let lifted = surface.point_at(at.x, at.y, tol)?;
        for (end, vertex) in [(0.0, first), (1.0, last)] {
            let held = model
                .node(&vertex)
                .and_then(|n| n.data().as_vertex())
                .map_or(0.0, |d| d.tolerance.get());
            if vertex_point(model, &vertex)?.distance(lifted) <= held + tol.confusion() {
                let corner = occurrence.at(end, tol)?;
                meetings.push((line.free(corner), Meeting::Vertex(vertex.node())));
                return Ok(());
            }
        }
    }
    // As a fraction of the forward edge, rounded so two traversals of one
    // seam name the same point.
    let forward_s = if occurrence.reversed() { 1.0 - s } else { s };
    let forward_s = (forward_s * 1e12).round() / 1e12;
    let node = occurrence.edge.node();
    cuts.entry(node).or_default().push(forward_s);
    meetings.push((line.free(at), Meeting::Inside(node, forward_s.to_bits())));
    Ok(())
}

/// One directed piece of a ring being assembled.
#[derive(Debug, Clone)]
struct Item {
    edge: Shape,
    from: TShapeId,
    to: TShapeId,
    polyline: Vec<Point2>,
}

/// Chain directed pieces end to start into closed rings: at a vertex, the
/// piece whose chart start is nearest where the ring stands.
fn chain(mut items: Vec<Item>, snap: f64) -> OgeomResult<Vec<(Vec<Shape>, Vec<Point2>)>> {
    let mut loops = Vec::new();
    while let Some(first) = items.pop() {
        let start_node = first.from;
        let start_at = first.polyline[0];
        let mut here = *first.polyline.last().unwrap_or(&start_at);
        let mut at_node = first.to;
        let mut edges = vec![first.edge];
        let mut outline = first.polyline;
        loop {
            if at_node == start_node && here.distance(start_at) <= snap {
                break;
            }
            let next = items
                .iter()
                .enumerate()
                .filter(|(_, it)| it.from == at_node)
                .min_by(|a, b| {
                    a.1.polyline[0]
                        .distance(here)
                        .total_cmp(&b.1.polyline[0].distance(here))
                })
                .map(|(i, _)| i);
            let Some(i) = next else {
                ogeom_bail!(Construction, "a piece of the cut face does not close");
            };
            let it = items.swap_remove(i);
            here = *it.polyline.last().unwrap_or(&here);
            at_node = it.to;
            edges.push(it.edge);
            outline.extend(it.polyline);
        }
        loops.push((edges, outline));
    }
    Ok(loops)
}

/// Split `edge` (forward) at fractions of its range, in order: the pieces
/// in order, and the new vertex at each fraction.
///
/// The fractions are measured along one description: the edge's pcurve on
/// `along` where one is named (a cut found in a face's chart), its curve
/// otherwise. Every other description is cut where it passes the same
/// point, which is the same fraction of its own range only where the edge
/// states its descriptions run at one parameter over one range: a pcurve
/// paced its own way (an arc's rational spline against its patch's chart)
/// reaches the same fraction somewhere else along the edge.
pub(crate) fn split_edge(
    model: &mut Model,
    edge: &Shape,
    fractions: &[f64],
    along: Option<SurfaceId>,
    tol: Tolerances,
) -> OgeomResult<(Vec<Shape>, Vec<Shape>)> {
    let data = edge_data(model, edge)?;
    let Some((start, end)) = edge_vertices(model, edge)? else {
        ogeom_bail!(Construction, "an edge with no vertices cannot be cut");
    };
    let measured =
        along.and_then(
            |surface| match data.pcurve_for(surface, &Location::identity()) {
                Some(
                    EdgeRepr::PCurve { curve, range, .. }
                    | EdgeRepr::Seam {
                        forward: curve,
                        range,
                        ..
                    },
                ) => Some((surface, *curve, *range)),
                _ => None,
            },
        );
    let curve = match data.curve3d() {
        Some(EdgeRepr::Curve3d { curve, range, .. }) if !data.degenerate => {
            let Some(geometry) = model.geometry().curve(*curve) else {
                ogeom_bail!(Dangling, "curve is not in this model");
            };
            Some((geometry.clone(), *range))
        }
        _ => None,
    };
    // Whether a description runs at the measured one's parameters over
    // its range, as the edge states, so the same fractions cut it.
    let same = |range: (f64, f64)| {
        let reference = match (measured, &curve) {
            (Some((_, _, r)), _) => Some(r),
            (None, Some((_, r))) => Some(*r),
            (None, None) => None,
        };
        data.same_parameter() && reference == Some(range)
    };
    let points: Vec<Option<Point>> = fractions
        .iter()
        .map(|&s| -> OgeomResult<Option<Point>> {
            Ok(match (measured, &curve) {
                (Some((surface, pcurve, range)), _) => {
                    let (Some(surface), Some(pcurve)) = (
                        model.geometry().surface(surface),
                        model.geometry().pcurve(pcurve),
                    ) else {
                        ogeom_bail!(Dangling, "an edge names geometry not in this model");
                    };
                    let uv = pcurve.point_at(range.0 + (range.1 - range.0) * s, tol)?;
                    Some(surface.point_at(uv.x, uv.y, tol)?)
                }
                (None, Some((c, range))) => {
                    Some(c.point_at(range.0 + (range.1 - range.0) * s, tol)?)
                }
                (None, None) => None,
            })
        })
        .collect::<OgeomResult<_>>()?;
    // Each description's own fractions, in the order the edge holds them.
    let mut own: Vec<Vec<f64>> = Vec::with_capacity(data.representations.len());
    for repr in &data.representations {
        let cut = match (repr, &curve) {
            (EdgeRepr::Curve3d { range, .. }, Some((c, _)))
                if measured.is_some() && !same(*range) =>
            {
                passing(fractions, &points, |s| {
                    c.point_at(range.0 + (range.1 - range.0) * s, tol).ok()
                })?
            }
            (
                EdgeRepr::PCurve {
                    curve: id,
                    range,
                    surface,
                    location,
                }
                | EdgeRepr::Seam {
                    forward: id,
                    range,
                    surface,
                    location,
                    ..
                },
                _,
            ) => {
                let measuring = measured.is_some_and(|(on, _, _)| on == *surface);
                if measuring || same(*range) || !location.is_identity() {
                    fractions.to_vec()
                } else {
                    let (Some(surface), Some(pcurve)) = (
                        model.geometry().surface(*surface),
                        model.geometry().pcurve(*id),
                    ) else {
                        ogeom_bail!(Dangling, "an edge names geometry not in this model");
                    };
                    passing(fractions, &points, |s| {
                        let uv = pcurve
                            .point_at(range.0 + (range.1 - range.0) * s, tol)
                            .ok()?;
                        surface.point_at(uv.x, uv.y, tol).ok()
                    })?
                }
            }
            _ => fractions.to_vec(),
        };
        own.push(cut);
    }
    // The new vertices stand on the curve, where it is cut.
    let curve_cut = data
        .representations
        .iter()
        .position(|r| matches!(r, EdgeRepr::Curve3d { .. }))
        .map_or_else(|| fractions.to_vec(), |k| own[k].clone());
    let mut vertices = Vec::with_capacity(fractions.len());
    for &s in &curve_cut {
        let vertex = match &curve {
            Some((c, range)) => {
                let p = c.point_at(range.0 + (range.1 - range.0) * s, tol)?;
                model.add_vertex(VertexData::with_tolerance(p, data.tolerance.get())?)
            }
            // A degenerate edge is one point, however it is cut.
            _ => start.clone(),
        };
        vertices.push(vertex);
    }
    let mut ends = vec![start];
    ends.extend(vertices.iter().cloned());
    ends.push(end);
    let mut pieces = Vec::with_capacity(fractions.len() + 1);
    for k in 0..=fractions.len() {
        let mut piece = data.clone();
        for (repr, cut) in piece.representations.iter_mut().zip(&own) {
            let (a, b) = (
                if k == 0 { 0.0 } else { cut[k - 1] },
                cut.get(k).copied().unwrap_or(1.0),
            );
            match repr {
                EdgeRepr::Curve3d { range, .. }
                | EdgeRepr::PCurve { range, .. }
                | EdgeRepr::Seam { range, .. } => {
                    *range = (
                        range.0 + (range.1 - range.0) * a,
                        range.0 + (range.1 - range.0) * b,
                    );
                }
                _ => {}
            }
        }
        piece.representations.retain(|r| {
            matches!(
                r,
                EdgeRepr::Curve3d { .. } | EdgeRepr::PCurve { .. } | EdgeRepr::Seam { .. }
            )
        });
        pieces.push(model.add_edge(piece, &[ends[k].clone(), ends[k + 1].clone()])?);
    }
    Ok((pieces, vertices))
}

/// The fractions at which a description, read through `at`, passes each
/// of `points`: for each, the nearest of evenly spaced samples, refined
/// between its neighbours. The fraction given stands where there is no
/// point (a degenerate edge) or the description cannot be read.
fn passing(
    fractions: &[f64],
    points: &[Option<Point>],
    at: impl Fn(f64) -> Option<Point>,
) -> OgeomResult<Vec<f64>> {
    const SAMPLES: u32 = 64;
    let mut out = Vec::with_capacity(fractions.len());
    for (&s, p) in fractions.iter().zip(points) {
        let Some(p) = p else {
            out.push(s);
            continue;
        };
        let off = |t: f64| at(t).map_or(f64::INFINITY, |q| q.distance(*p));
        let nearest = (0..=SAMPLES)
            .map(|i| (i, off(f64::from(i) / f64::from(SAMPLES))))
            .min_by(|a, b| a.1.total_cmp(&b.1));
        let Some((i, _)) = nearest.filter(|(_, d)| d.is_finite()) else {
            out.push(s);
            continue;
        };
        let (lo, hi) = (
            f64::from(i.saturating_sub(1)) / f64::from(SAMPLES),
            f64::from((i + 1).min(SAMPLES)) / f64::from(SAMPLES),
        );
        out.push(ogeom_math::solve::minimize(off, lo, hi, ogeom_math::Criteria::default())?.value);
    }
    Ok(out)
}

/// The surface's iso-curve along `line`, and the map from the line's free
/// parameter `w` to the curve's own: `t = scale * w + offset`.
fn iso_curve(
    surface: &SurfaceGeometry,
    line: IsoLine,
    probes: &[f64],
    span: (f64, f64),
    tol: Tolerances,
) -> OgeomResult<Option<(Curve, f64, f64)>> {
    let c = line.at();
    let circle = |centre: Point, z: Direction, x: Point, radius: f64| -> OgeomResult<Curve> {
        let x = Direction::new(x - centre, tol)?;
        let frame = Frame::new(centre, z, x, tol)?;
        Ok(CircleCurve::new(Circle::new(frame, radius, tol)?).into())
    };
    let straight = |at: Point, along: ogeom_math::Vector| -> OgeomResult<(Curve, f64)> {
        let length = along.magnitude();
        let direction = Direction::new(along, tol)?;
        Ok((
            LineCurve::new(Axis {
                location: at,
                direction,
            })
            .into(),
            length,
        ))
    };
    let found: (Curve, f64, f64) = match (surface, line) {
        (SurfaceGeometry::BSpline(b), IsoLine::U(_)) => {
            (Curve::BSpline(b.iso_u_curve(c, tol)?), 1.0, 0.0)
        }
        (SurfaceGeometry::BSpline(b), IsoLine::V(_)) => {
            (Curve::BSpline(b.iso_v_curve(c, tol)?), 1.0, 0.0)
        }
        (
            SurfaceGeometry::Plane(_) | SurfaceGeometry::Cylinder(_) | SurfaceGeometry::Cone(_),
            IsoLine::U(_),
        )
        | (SurfaceGeometry::Plane(_) | SurfaceGeometry::Extrusion(_), _) => {
            // Straight in the free parameter: through the point at the
            // span's start, along the rate the surface moves with w, both
            // read inside the span (a short drum's window ends before
            // w = 1, and the surface refuses to be read past it).
            let (lo, hi) = if span.0.is_finite() && span.1.is_finite() && span.1 > span.0 {
                span
            } else {
                (0.0, 1.0)
            };
            let at = surface.point_at(line.point(lo).x, line.point(lo).y, tol)?;
            let next = surface.point_at(line.point(hi).x, line.point(hi).y, tol)?;
            match (surface, line) {
                (SurfaceGeometry::Extrusion(e), IsoLine::V(_)) => {
                    let shift = Transform::translation(e.direction().vector() * c);
                    (e.curve().transformed(&shift, tol)?, 1.0, 0.0)
                }
                _ => {
                    let (curve, length) = straight(at, next - at)?;
                    let speed = length / (hi - lo);
                    (curve, speed, -speed * lo)
                }
            }
        }
        (
            SurfaceGeometry::Cylinder(_)
            | SurfaceGeometry::Cone(_)
            | SurfaceGeometry::Sphere(_)
            | SurfaceGeometry::Torus(_)
            | SurfaceGeometry::Revolution(_),
            IsoLine::V(_),
        ) => {
            // A parallel: the circle the surface's `u = 0` point turns on.
            let axis = match surface {
                SurfaceGeometry::Cylinder(s) => frame_axis(s.cylinder().frame()),
                SurfaceGeometry::Cone(s) => frame_axis(s.cone().frame()),
                SurfaceGeometry::Sphere(s) => frame_axis(s.sphere().frame()),
                SurfaceGeometry::Torus(s) => frame_axis(s.torus().frame()),
                SurfaceGeometry::Revolution(s) => s.axis(),
                _ => return Ok(None),
            };
            let start = surface.point_at(0.0, c, tol)?;
            let centre = axis.project(start);
            let radius = start.distance(centre);
            if radius <= tol.confusion() {
                return Ok(None);
            }
            (circle(centre, axis.direction, start, radius)?, 1.0, 0.0)
        }
        (SurfaceGeometry::Sphere(_) | SurfaceGeometry::Torus(_), IsoLine::U(_)) => {
            // A meridian or tube circle, its angle the surface's `v`.
            let at0 = surface.point_at(c, 0.0, tol)?;
            let quarter = surface.point_at(c, core::f64::consts::FRAC_PI_2, tol)?;
            let opposite = surface.point_at(c, -core::f64::consts::FRAC_PI_2, tol)?;
            let centre = Point::from_vector((quarter.to_vector() + opposite.to_vector()) * 0.5);
            let radius = at0.distance(centre);
            let z = Direction::new((at0 - centre).cross(quarter - centre), tol)?;
            (circle(centre, z, at0, radius)?, 1.0, 0.0)
        }
        (SurfaceGeometry::Revolution(r), IsoLine::U(_)) => {
            let turn = Transform::rotation(r.axis(), c);
            (r.curve().transformed(&turn, tol)?, 1.0, 0.0)
        }
        (SurfaceGeometry::Trimmed(t), _) => return iso_curve(t.basis(), line, probes, span, tol),
        // No closed form (an offset of a spline): the iso-curve fitted to
        // the surface along the line at the line's own parameters,
        // same-parameter with the straight chart line it runs along, and
        // held to the confusion distance between its samples as well as at
        // them.
        _ => {
            const N: usize = 128;
            let (lo, hi) = span;
            if hi <= lo || !hi.is_finite() || !lo.is_finite() {
                return Ok(None);
            }
            #[allow(clippy::cast_precision_loss)]
            let params: Vec<f64> = (0..=N)
                .map(|i| lo + (hi - lo) * i as f64 / N as f64)
                .collect();
            let fitted = ogeom_geom::fit::fit_curve_sampled(
                |w| {
                    let on = line.point(w);
                    surface.point_at(on.x, on.y, tol)
                },
                &params,
                false,
                3,
                tol.confusion(),
                tol,
            )?;
            if !fitted.met {
                return Ok(None);
            }
            (Curve::BSpline(fitted.curve), 1.0, 0.0)
        }
    };
    // The closed forms above assume the surfaces' own conventions. Hold
    // them to it.
    let (curve, scale, offset) = &found;
    for &w in probes {
        let on = line.point(w);
        let expected = surface.point_at(on.x, on.y, tol)?;
        let got = curve.point_at(scale * w + offset, tol)?;
        if expected.distance(got) > tol.confusion() {
            return Ok(None);
        }
    }
    Ok(Some(found))
}

fn frame_axis(frame: Frame) -> Axis {
    Axis {
        location: frame.origin(),
        direction: frame.z(),
    }
}

// --- helpers -------------------------------------------------------------------

pub(crate) fn forward(shape: &Shape) -> Shape {
    if shape.orientation() == Orientation::Reversed {
        shape.reversed()
    } else {
        shape.clone()
    }
}

pub(crate) fn edge_data(model: &Model, edge: &Shape) -> OgeomResult<EdgeData> {
    match model.node(edge).map(ogeom_topo::TShape::data) {
        Some(NodeData::Edge(data)) => Ok((**data).clone()),
        Some(_) => ogeom_bail!(Construction, "expected an edge"),
        None => ogeom_bail!(Dangling, "edge is not in this model"),
    }
}

pub(crate) fn face_data(model: &Model, face: &Shape) -> OgeomResult<ogeom_topo::FaceData> {
    match model.node(face).map(ogeom_topo::TShape::data) {
        Some(NodeData::Face(data)) => Ok((**data).clone()),
        Some(_) => ogeom_bail!(Construction, "expected a face"),
        None => ogeom_bail!(Dangling, "face is not in this model"),
    }
}

pub(crate) fn vertex_point(model: &Model, vertex: &Shape) -> OgeomResult<Point> {
    let Some(data) = model.node(vertex).and_then(|n| n.data().as_vertex()) else {
        ogeom_bail!(Construction, "a vertex holds no point");
    };
    Ok(data.point)
}

/// The area a ring encloses in its chart, positive counter-clockwise.
///
/// Summed about the ring's first point: the same sum as about the chart's
/// origin, which may stand far off, where it is a difference of huge
/// products.
pub(crate) fn signed_area(ring: &[Point2]) -> f64 {
    let n = ring.len();
    let Some(&anchor) = ring.first() else {
        return 0.0;
    };
    (0..n)
        .map(|i| (ring[i] - anchor).cross(ring[(i + 1) % n] - anchor))
        .sum::<f64>()
        * 0.5
}

/// Even-odd containment against every ring.
pub(crate) fn inside(rings: &[Vec<Point2>], p: Point2) -> bool {
    let mut odd = false;
    for ring in rings {
        let n = ring.len();
        for i in 0..n {
            let (a, b) = (ring[i], ring[(i + 1) % n]);
            if (a.y > p.y) != (b.y > p.y) {
                let x = a.x + (p.y - a.y) / (b.y - a.y) * (b.x - a.x);
                if x > p.x {
                    odd = !odd;
                }
            }
        }
    }
    odd
}

#[cfg(test)]
#[allow(clippy::unwrap_used, clippy::expect_used)]
mod tests {
    use super::*;
    use ogeom_geom::{BSplineSurface, OffsetSurface};
    use ogeom_math::ControlGrid;

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

    /// A strip one long and a hundredth wide encloses its own area wherever
    /// it stands in its chart, up to 2e7 from the chart's origin.
    #[test]
    fn a_ring_far_out_in_its_chart_encloses_its_own_area() {
        for o in [0.0, 1e6, 1e7, 2e7] {
            let (x, y) = (o, -0.7 * o);
            let ring = [(0.0, 0.0), (1.0, 0.0), (1.0, 0.01), (0.0, 0.01)]
                .map(|(dx, dy)| Point2::new(x + dx, y + dy));
            let area = signed_area(&ring);
            assert!((area - 0.01).abs() < 1e-9, "o {o:e}: area {area}");
        }
    }

    /// An iso-curve with no closed form keeps to the surface between the
    /// points it was fitted at. The surface is the offset of a patch
    /// carrying a ridge about as wide as a hundred and twenty-eighth of the
    /// line, so a spline through the line's even samples swings off the
    /// ridge between them; the curve stands within the confusion distance
    /// of the surface everywhere along the line.
    #[test]
    fn an_iso_curve_without_a_closed_form_keeps_to_the_surface_between_its_samples() {
        const ROWS: u32 = 401;
        let mut control = Vec::new();
        for i in 0..ROWS {
            let x = f64::from(i) / f64::from(ROWS - 1);
            let ridge = 0.05 * (-((x - 0.503) / 0.006).powi(2)).exp();
            control.extend([Point::new(x, 0.0, ridge), Point::new(x, 1.0, ridge)]);
        }
        let mut u_knots = vec![0.0; 4];
        u_knots.extend((1..ROWS - 3).map(|k| f64::from(k) / f64::from(ROWS - 3)));
        u_knots.extend([1.0; 4]);
        let patch = BSplineSurface::new(
            KnotVector::new(u_knots, 3).unwrap(),
            KnotVector::new(vec![0.0, 0.0, 1.0, 1.0], 1).unwrap(),
            &ControlGrid::new(control, ROWS as usize, 2).unwrap(),
            T,
        )
        .unwrap();
        let surface =
            SurfaceGeometry::Offset(Box::new(OffsetSurface::new(patch.into(), 5e-4).unwrap()));
        let line = IsoLine::V(0.5);
        let (curve, scale, offset) = iso_curve(&surface, line, &[0.5], (0.0, 1.0), T)
            .unwrap()
            .expect("the ridge is followed within the confusion distance");
        let mut worst = 0.0_f64;
        for k in 0..=20_000 {
            let w = f64::from(k) / 20_000.0;
            let on = line.point(w);
            let expected = surface.point_at(on.x, on.y, T).unwrap();
            let got = curve.point_at(scale * w + offset, T).unwrap();
            worst = worst.max(expected.distance(got));
        }
        assert!(
            worst <= 2.0 * T.confusion(),
            "the iso-curve leaves the surface by {worst:.3e}"
        );
    }
}