brep_kernel/healing/

sew.rs

//! Standalone heal/sew — Golovanov's "sewing": assemble a valid oriented
//! shell out of faces whose boundary representations were built independently
//! (imported shells, detached face groups, healing after hand edits).
//!
//! The boolean assembler and the offset pipeline each carry private sewing
//! passes specialised to their own invariants (`sew_coincident_one_use_edges`
//! requires already-opposed traversals; the offset welds know their rims).
//! This module is the GENERAL entry: it pairs coincident one-use boundary
//! edges by geometry alone — open chains by matched endpoints + locus
//! agreement, closed rims by mutual locus agreement — WITHOUT any orientation
//! precondition, because a face soup's components may be arbitrarily flipped.
//! Orientation is repaired after pairing: coedge-direction coherence is
//! propagated across the now-shared edges and the whole solid is flipped
//! outward by signed volume.
//!
//! Sewing is BEST-EFFORT and honest: pairs that cannot be joined within
//! tolerance stay open and are counted in the report; nothing is force-welded.

use crate::mass_properties::solid_signed_volume;
use crate::offset_shell::{flip_all_faces, orient_open_solid_faces};
use crate::topology::{BrepSolid, CoedgeRecord, EdgeRecord, FaceRecord, ShellRecord, VertexRecord};
use crate::{build_pcurve_on_surface, project_point_to_curve, NurbsCurve, Vec3};
use rustc_hash::{FxHashMap as HashMap, FxHashSet as HashSet};
use serde::{Deserialize, Serialize};

#[derive(Clone, Debug, Serialize, Deserialize)]
pub struct SewReport {
    pub edges_sewn: usize,
    pub shells_merged: usize,
    pub open_edges_before: usize,
    pub open_edges_after: usize,
    pub oriented_outward: bool,
    pub issues: Vec<String>,
}

fn one_use_edge_ids(solid: &BrepSolid) -> HashSet<u64> {
    let mut counts = HashMap::<u64, usize>::default();
    for coedge in solid
        .shells
        .iter()
        .flat_map(|shell| &shell.faces)
        .flat_map(|face| &face.loops)
        .flat_map(|loop_record| &loop_record.coedges)
    {
        *counts.entry(coedge.edge_id).or_default() += 1;
    }
    counts
        .into_iter()
        .filter(|(_, count)| *count == 1)
        .map(|(id, _)| id)
        .collect()
}

fn curve_point(edge: &EdgeRecord, fraction: f64) -> Result<Vec3, String> {
    edge.curve
        .evaluate(edge.t0 + (edge.t1 - edge.t0) * fraction)
}

/// Worst distance from sampled points of `piece` to the locus of `carrier`.
fn locus_deviation(
    piece: &EdgeRecord,
    carrier: &EdgeRecord,
    samples: usize,
) -> Result<f64, String> {
    let mut worst = 0.0f64;
    for index in 0..=samples {
        let point = curve_point(piece, index as f64 / samples as f64)?;
        worst = worst.max(project_point_to_curve(&carrier.curve, point)?.distance);
    }
    Ok(worst)
}

/// Pointwise agreement at matched fractions — the two edges are not merely
/// the same locus but the SAME parametrization (the duplicated-edge case),
/// letting a rebind keep the existing pcurve exactly.
fn same_parametrization(
    first: &EdgeRecord,
    second: &EdgeRecord,
    tolerance: f64,
) -> Result<bool, String> {
    for index in 0..=8 {
        let fraction = index as f64 / 8.0;
        if curve_point(first, fraction)?
            .sub(curve_point(second, fraction)?)
            .length()
            > tolerance
        {
            return Ok(false);
        }
    }
    Ok(true)
}

/// The 3D walk this coedge's loop takes along its edge, sampled from its
/// pcurve through its face surface. Fractions 0 and 0.35: an interior second
/// sample avoids the antipodal-projection ambiguity a closed rim has at 0.5.
fn walk_points(face: &FaceRecord, coedge: &CoedgeRecord) -> Result<[Vec3; 2], String> {
    let [start, end] = coedge.pcurve.domain()?;
    let sample = |fraction: f64| -> Result<Vec3, String> {
        let uv = coedge.pcurve.evaluate(start + (end - start) * fraction)?;
        face.surface.evaluate(uv.x, uv.y)
    };
    Ok([sample(0.0)?, sample(0.35)?])
}

/// One planned coedge rebind, resolved before any mutation so a failed plan
/// (an unprojectable pcurve, say) leaves the pair untouched and open.
struct RebindPatch {
    shell: usize,
    face: usize,
    loop_index: usize,
    coedge: usize,
    forward: bool,
    pcurve: Option<NurbsCurve>,
}

/// Plan the rebind of every coedge referencing `remove` onto `keep`. The
/// coedge's `forward` is derived from the 3D walk its existing pcurve
/// produces (geometry is authoritative — components may be flipped), the
/// pcurve is kept when the parametrizations agree and refit otherwise.
fn plan_rebind(
    solid: &BrepSolid,
    keep: &EdgeRecord,
    remove: &EdgeRecord,
    tolerance: f64,
) -> Result<Vec<RebindPatch>, String> {
    let identical = same_parametrization(keep, remove, tolerance)?;
    let closed = keep.start_vertex_id == keep.end_vertex_id;
    let period = keep.t1 - keep.t0;
    let mut patches = Vec::new();
    for (shell_index, shell) in solid.shells.iter().enumerate() {
        for (face_index, face) in shell.faces.iter().enumerate() {
            for (loop_index, loop_record) in face.loops.iter().enumerate() {
                for (coedge_index, coedge) in loop_record.coedges.iter().enumerate() {
                    if coedge.edge_id != remove.id {
                        continue;
                    }
                    let [walk_start, walk_next] = walk_points(face, coedge)?;
                    let t_start = project_point_to_curve(&keep.curve, walk_start)?.u;
                    let t_next = project_point_to_curve(&keep.curve, walk_next)?.u;
                    let forward = if closed {
                        let mut delta = t_next - t_start;
                        while delta > period * 0.5 {
                            delta -= period;
                        }
                        while delta < -period * 0.5 {
                            delta += period;
                        }
                        delta > 0.0
                    } else {
                        t_next > t_start
                    };
                    let pcurve = if identical {
                        None
                    } else {
                        let traversal = if forward {
                            keep.curve.clone()
                        } else {
                            keep.curve.reversed()?
                        };
                        Some(build_pcurve_on_surface(&face.surface, &traversal)?)
                    };
                    patches.push(RebindPatch {
                        shell: shell_index,
                        face: face_index,
                        loop_index,
                        coedge: coedge_index,
                        forward,
                        pcurve,
                    });
                }
            }
        }
    }
    Ok(patches)
}

fn apply_rebind(
    solid: &mut BrepSolid,
    keep: &EdgeRecord,
    remove: &EdgeRecord,
    patches: Vec<RebindPatch>,
    tolerance: f64,
) {
    for patch in patches {
        let coedge = &mut solid.shells[patch.shell].faces[patch.face].loops[patch.loop_index]
            .coedges[patch.coedge];
        coedge.edge_id = keep.id;
        coedge.forward = patch.forward;
        if let Some(pcurve) = patch.pcurve {
            coedge.pcurve = pcurve;
        }
    }
    solid.edges.retain(|edge| edge.id != remove.id);
    // Weld the removed edge's endpoint vertices into the kept edge's, so
    // OTHER edges of the removed component (a loop mixes sewn rim edges with
    // still-duplicate interior edges) chain through the shared vertices.
    // Only coincident endpoints weld — a rotated closed rim keeps its own
    // seam vertex.
    let point_of = |vertex_id: u64| {
        solid
            .vertices
            .iter()
            .find(|vertex| vertex.id == vertex_id)
            .map(|vertex| vertex.point)
    };
    let mut welds = Vec::<(u64, u64)>::new();
    for from in [remove.start_vertex_id, remove.end_vertex_id] {
        let Some(from_point) = point_of(from) else {
            continue;
        };
        let target = [keep.start_vertex_id, keep.end_vertex_id]
            .into_iter()
            .filter_map(|candidate| {
                point_of(candidate).map(|point| (candidate, point.sub(from_point).length()))
            })
            .min_by(|first, second| first.1.total_cmp(&second.1));
        if let Some((to, distance)) = target {
            if distance <= tolerance && from != to {
                welds.push((from, to));
            }
        }
    }
    for (from, to) in welds {
        for edge in &mut solid.edges {
            if edge.start_vertex_id == from {
                edge.start_vertex_id = to;
            }
            if edge.end_vertex_id == from {
                edge.end_vertex_id = to;
            }
        }
    }
}

/// Decide which of an accepted pair is KEPT and which is REMOVED.
///
/// For an open chain, and for a closed rim whose two copies share a seam
/// point, either direction works and the scan's own order stands — `first` is
/// kept, exactly as before.
///
/// A seam-ROTATED closed rim is different, and the difference is not cosmetic.
/// The removed rim takes its seam vertex out of the solid with it, so every
/// OTHER edge that chained through that vertex is left chaining through
/// nothing. `apply_rebind` deliberately does not drag it onto the kept rim's
/// seam — a rotated closed rim keeps its own seam vertex, because moving it
/// would pull a real edge endpoint an arbitrary arc round the circle, off its
/// own curve — so the tear is not repaired, it is inherited. Sewn that way a
/// cylinder comes back with `edges_sewn: 1, open_edges_after: 0` over a wall
/// loop torn at the seam: a report that reads like success.
///
/// The direction that works is the one whose REMOVED rim's seam vertex bounds
/// no other edge. When neither qualifies — both seams anchor real geometry —
/// no rebind keeps both loops closed, and the pair is left OPEN and counted
/// rather than sewn into a solid whose own report contradicts its issues.
fn orient_pair<'a>(
    solid: &BrepSolid,
    first: &'a EdgeRecord,
    second: &'a EdgeRecord,
    points: &HashMap<u64, Vec3>,
    tolerance: f64,
) -> Option<(&'a EdgeRecord, &'a EdgeRecord)> {
    // The caller has already rejected mismatched closedness, so testing
    // `first` tests both.
    if first.start_vertex_id != first.end_vertex_id {
        return Some((first, second));
    }
    let (Some(&first_seam), Some(&second_seam)) = (
        points.get(&first.start_vertex_id),
        points.get(&second.start_vertex_id),
    ) else {
        return Some((first, second));
    };
    if first_seam.sub(second_seam).length() <= tolerance {
        // Seams coincide: `apply_rebind`'s weld joins them and every loop
        // through either one still chains. Order is free.
        return Some((first, second));
    }
    let seam_is_exclusive = |rim: &EdgeRecord| {
        !solid.edges.iter().any(|edge| {
            edge.id != rim.id
                && (edge.start_vertex_id == rim.start_vertex_id
                    || edge.end_vertex_id == rim.start_vertex_id)
        })
    };
    if seam_is_exclusive(second) {
        Some((first, second))
    } else if seam_is_exclusive(first) {
        Some((second, first))
    } else {
        None
    }
}

/// Merge shells that now share an edge into single shells (union-find).
fn merge_connected_shells(solid: &mut BrepSolid) -> usize {
    let mut owner_of_edge = HashMap::<u64, usize>::default();
    let mut parent = (0..solid.shells.len()).collect::<Vec<_>>();
    fn root(parent: &mut [usize], index: usize) -> usize {
        if parent[index] != index {
            parent[index] = root(parent, parent[index]);
        }
        parent[index]
    }
    for (shell_index, shell) in solid.shells.iter().enumerate() {
        for coedge in shell
            .faces
            .iter()
            .flat_map(|face| &face.loops)
            .flat_map(|loop_record| &loop_record.coedges)
        {
            if let Some(&other) = owner_of_edge.get(&coedge.edge_id) {
                let first = root(&mut parent, other);
                let second = root(&mut parent, shell_index);
                if first != second {
                    parent[second] = first;
                }
            } else {
                owner_of_edge.insert(coedge.edge_id, shell_index);
            }
        }
    }
    let before = solid.shells.len();
    let original = std::mem::take(&mut solid.shells);
    let mut merged = Vec::<ShellRecord>::new();
    let mut group_of = HashMap::<usize, usize>::default();
    for (index, shell) in original.into_iter().enumerate() {
        let group = root(&mut parent, index);
        if let Some(&target) = group_of.get(&group) {
            merged[target].faces.extend(shell.faces);
        } else {
            group_of.insert(group, merged.len());
            merged.push(shell);
        }
    }
    solid.shells = merged;
    before - solid.shells.len()
}

/// Which `second` candidates the pair search offers for a given `first`.
///
/// The two strategies feed the SAME loop body in the SAME ascending order, so
/// "the prefilter is correct" reduces to one claim: the offered list is a
/// SUPERSET of the candidates that could pass. Nothing else in the search can
/// differ, which is what makes the parity test (`sew_prefilter_parity`) a real
/// proof rather than a spot check — first-match-wins makes iteration order
/// part of the RESULT, not just the cost.
#[derive(Clone, Copy, PartialEq, Eq, Debug)]
enum PairSearch {
    /// Production: neighbourhood query for open chains, closed-rim shortlist
    /// for closed rims.
    Hashed,
    /// The reference brute-force scan the prefilter must reproduce exactly.
    #[cfg_attr(not(test), allow(dead_code))]
    Exhaustive,
}

/// The integer cell of `point`, or `None` when the coordinates cannot be
/// bucketed at all (non-finite, or so large the index overflows the i64 key).
/// A `None` anywhere disables the prefilter for the whole pass rather than
/// silently dropping a candidate — exact parity with the exhaustive scan is
/// the contract, and one lost candidate would break it.
fn grid_cell(point: Vec3, cell: f64) -> Option<[i64; 3]> {
    let mut key = [0i64; 3];
    for (slot, value) in key.iter_mut().zip([point.x, point.y, point.z]) {
        if !value.is_finite() {
            return None;
        }
        let scaled = (value / cell).floor();
        // 9e15 keeps the cast well inside i64 and inside f64's exact-integer
        // range, so `as i64` is a faithful conversion, not a saturating one.
        if !scaled.is_finite() || scaled.abs() > 9.0e15 {
            return None;
        }
        *slot = scaled as i64;
    }
    Some(key)
}

/// Uniform spatial hash over the ENDPOINTS of the open one-use candidates —
/// the prefilter that turns the pair search from a full rescan per accepted
/// pair into a neighbourhood query.
///
/// Cells are `2·tolerance` wide, not `tolerance`. A partner endpoint is at
/// most `tolerance = cell/2` away, i.e. at most HALF a cell, so it lands in
/// this cell or an immediate neighbour with a full half-cell of slack — slack
/// that absorbs the rounding of the `point/cell` division. At `cell =
/// tolerance` two points exactly `tolerance` apart can straddle two cell
/// boundaries and the 3×3×3 block would still hold them, but with zero margin
/// for the division's last bit; the factor of two costs one extra ring of
/// mostly-empty buckets and buys certainty.
struct EndpointGrid {
    cell: f64,
    buckets: HashMap<[i64; 3], Vec<usize>>,
}

impl EndpointGrid {
    /// Bucket every open candidate under BOTH of its endpoints.
    ///
    /// Closed rims are deliberately absent. They are paired by mutual locus
    /// agreement alone — `sew_solid`'s closed branch never looks at a vertex —
    /// and a seam-ROTATED duplicate rim shares its partner's locus while
    /// sharing no point position at all, so any point-keyed bucket (midpoint
    /// included) would drop exactly the pair the closed lane exists to catch.
    /// Closed rims therefore keep a linear shortlist; see `offer_seconds`.
    ///
    /// Candidates whose vertex records are missing are absent too: the
    /// exhaustive scan `continue`s such a pair outright, so dropping it here
    /// changes nothing.
    fn build(endpoints: &[Option<[Vec3; 2]>], closed: &[bool], tolerance: f64) -> Option<Self> {
        let cell = 2.0 * tolerance;
        if !(cell.is_finite() && cell > 0.0) {
            return None;
        }
        let mut buckets = HashMap::<[i64; 3], Vec<usize>>::default();
        for (position, ends) in endpoints.iter().enumerate() {
            if closed[position] {
                continue;
            }
            let Some(ends) = ends else {
                continue;
            };
            for point in ends {
                buckets
                    .entry(grid_cell(*point, cell)?)
                    .or_default()
                    .push(position);
            }
        }
        Some(EndpointGrid { cell, buckets })
    }

    /// Every open candidate with an endpoint in the 3×3×3 block around
    /// `point`, ascending and deduplicated. A superset of "has an endpoint
    /// within tolerance of `point`", which is all the pair search needs.
    /// `false` means the point itself is unbucketable and the caller must not
    /// trust the (empty) answer.
    fn near(&self, point: Vec3, out: &mut Vec<usize>) -> bool {
        out.clear();
        let Some(centre) = grid_cell(point, self.cell) else {
            return false;
        };
        for dx in -1..=1 {
            for dy in -1..=1 {
                for dz in -1..=1 {
                    if let Some(bucket) =
                        self.buckets
                            .get(&[centre[0] + dx, centre[1] + dy, centre[2] + dz])
                    {
                        out.extend_from_slice(bucket);
                    }
                }
            }
        }
        out.sort_unstable();
        out.dedup();
        true
    }
}

/// Fill `offered` with the `second` positions to test against
/// `candidates[first_position]`, ASCENDING — the order the exhaustive scan
/// used, which first-match-wins makes load-bearing.
///
/// Three lanes, each a provable superset of what can pass:
/// - **No grid** (prefilter disabled, or `PairSearch::Exhaustive`): every
///   later position, i.e. the original scan verbatim.
/// - **Closed `first`**: the later CLOSED positions. An open `second` is
///   rejected by the `first_closed != second_closed` test before any
///   geometry is touched, so shortlisting them away changes nothing — but
///   nothing finer is available, because closed-rim pairing reads the locus
///   and not any point (see `EndpointGrid::build`). Closed rims stay
///   quadratic in their own (small) count; that is honest, not hidden.
/// - **Open `first`**: one neighbourhood query around `first`'s START point.
///   An open pair is accepted only when `first`'s start matches `second`'s
///   start (direct) or `second`'s end (reversed), so EVERY pair that can pass
///   has an endpoint within tolerance of that single point, and one query is
///   a complete superset. A `first` whose vertex records are missing offers
///   nothing, matching the exhaustive scan's `continue`.
#[allow(clippy::too_many_arguments)]
fn offer_seconds(
    grid: Option<&EndpointGrid>,
    first_position: usize,
    first_closed: bool,
    endpoints: &[Option<[Vec3; 2]>],
    closed_positions: &[usize],
    candidate_count: usize,
    neighbours: &mut Vec<usize>,
    offered: &mut Vec<usize>,
) {
    offered.clear();
    let Some(grid) = grid else {
        offered.extend(first_position + 1..candidate_count);
        return;
    };
    if first_closed {
        offered.extend(
            closed_positions
                .iter()
                .copied()
                .filter(|&position| position > first_position),
        );
        return;
    }
    let Some([start, _]) = endpoints[first_position] else {
        return;
    };
    if grid.near(start, neighbours) {
        offered.extend(
            neighbours
                .iter()
                .copied()
                .filter(|&position| position > first_position),
        );
    }
}

/// Best-effort sew of a solid's open boundary edges.
///
/// Pairs coincident one-use edges (open chains by matched endpoints + locus
/// agreement, closed rims by mutual locus agreement) and rebinds each pair to
/// one shared edge, merging the shells they join. Orientation carries NO
/// precondition: after pairing, coedge-direction coherence is propagated
/// across the shared edges and the result is flipped outward by signed
/// volume. Unsewable gaps stay open and are reported, never force-welded.
pub fn sew_solid(solid: &BrepSolid, tolerance: f64) -> Result<(BrepSolid, SewReport), String> {
    let (sewn, report, _examined) = sew_solid_with_search(solid, tolerance, PairSearch::Hashed)?;
    Ok((sewn, report))
}

/// [`sew_solid`] with the pair-search strategy exposed, plus the number of
/// (first, second) pairs the search actually examined.
///
/// The count is what the scaling smoke asserts on: it is a deterministic
/// function of the input, where wall time is not, so a band set from it is a
/// real bound rather than a machine-speed lottery.
fn sew_solid_with_search(
    solid: &BrepSolid,
    tolerance: f64,
    search: PairSearch,
) -> Result<(BrepSolid, SewReport, usize), String> {
    if !(tolerance.is_finite() && tolerance > 0.0) {
        return Err("sew_solid: tolerance must be positive".into());
    }
    let mut result = solid.clone();
    let open_edges_before = one_use_edge_ids(&result).len();
    let mut edges_sewn = 0usize;
    let mut pairs_examined = 0usize;
    // Pairs whose rebind plan failed (unprojectable pcurve) — left open
    // rather than retried forever.
    //
    // Keyed UNORDERED, because `orient_pair` may swap which of the two is kept
    // and which is removed. Keyed by (keep, remove) the block would be looked
    // up under (first, second), miss on a swapped pair, and the scan would
    // re-offer it, re-swap it, and fail to plan it again on every outer
    // iteration — nothing having been mutated in between, that is a hang, not
    // a retry. A pair that cannot be rebound is blocked whichever way round.
    let mut blocked = HashSet::<(u64, u64)>::default();
    let unordered = |first: u64, second: u64| (first.min(second), first.max(second));
    let mut offered = Vec::<usize>::new();
    let mut neighbours = Vec::<usize>::new();
    loop {
        let one_use = one_use_edge_ids(&result);
        // Candidate POSITIONS into `result.edges`, in edge order — the same
        // sequence the exhaustive scan walked, without cloning a NURBS curve
        // per candidate per accepted pair (that clone was the other O(k·n)
        // term hiding behind the O(k·n²) comparison count).
        let candidates = result
            .edges
            .iter()
            .enumerate()
            .filter(|(_, edge)| !edge.degenerate && one_use.contains(&edge.id))
            .map(|(index, _)| index)
            .collect::<Vec<_>>();
        let points = result
            .vertices
            .iter()
            .map(|vertex| (vertex.id, vertex.point))
            .collect::<HashMap<_, _>>();
        let closed = candidates
            .iter()
            .map(|&index| {
                let edge = &result.edges[index];
                edge.start_vertex_id == edge.end_vertex_id
            })
            .collect::<Vec<_>>();
        // Both endpoint positions, or `None` when a vertex record is missing —
        // exactly the tuple the exhaustive scan destructures, hoisted out of
        // the inner loop so each candidate's two lookups happen once.
        let endpoints = candidates
            .iter()
            .map(|&index| {
                let edge = &result.edges[index];
                match (
                    points.get(&edge.start_vertex_id),
                    points.get(&edge.end_vertex_id),
                ) {
                    (Some(&start), Some(&end)) => Some([start, end]),
                    _ => None,
                }
            })
            .collect::<Vec<_>>();
        let closed_positions = (0..candidates.len())
            .filter(|&position| closed[position])
            .collect::<Vec<_>>();
        let grid = match search {
            PairSearch::Hashed => EndpointGrid::build(&endpoints, &closed, tolerance),
            PairSearch::Exhaustive => None,
        };

        let mut chosen = None;
        'pairs: for first_position in 0..candidates.len() {
            let first = &result.edges[candidates[first_position]];
            let first_closed = closed[first_position];
            offer_seconds(
                grid.as_ref(),
                first_position,
                first_closed,
                &endpoints,
                &closed_positions,
                candidates.len(),
                &mut neighbours,
                &mut offered,
            );
            for second_position in offered.iter().copied() {
                let second = &result.edges[candidates[second_position]];
                pairs_examined += 1;
                if blocked.contains(&unordered(first.id, second.id)) {
                    continue;
                }
                if first_closed != closed[second_position] {
                    continue;
                }
                if !first_closed {
                    let (Some([fs, fe]), Some([ss, se])) =
                        (endpoints[first_position], endpoints[second_position])
                    else {
                        continue;
                    };
                    let direct =
                        fs.sub(ss).length() <= tolerance && fe.sub(se).length() <= tolerance;
                    let reversed =
                        fs.sub(se).length() <= tolerance && fe.sub(ss).length() <= tolerance;
                    if !direct && !reversed {
                        continue;
                    }
                }
                if locus_deviation(first, second, 8)? > tolerance
                    || locus_deviation(second, first, 8)? > tolerance
                {
                    continue;
                }
                let Some((keep, remove)) = orient_pair(&result, first, second, &points, tolerance)
                else {
                    // A seam-rotated rim pair with no safe direction: leave it
                    // open, do not sew a torn loop shut.
                    continue;
                };
                chosen = Some((keep.clone(), remove.clone()));
                break 'pairs;
            }
        }
        let Some((keep, remove)) = chosen else {
            break;
        };
        match plan_rebind(&result, &keep, &remove, tolerance) {
            Ok(patches) => {
                apply_rebind(&mut result, &keep, &remove, patches, tolerance);
                edges_sewn += 1;
            }
            Err(_) => {
                blocked.insert(unordered(keep.id, remove.id));
            }
        }
    }
    // Drop vertices only referenced by removed duplicate edges.
    let used_vertices = result
        .edges
        .iter()
        .flat_map(|edge| [edge.start_vertex_id, edge.end_vertex_id])
        .collect::<HashSet<_>>();
    result
        .vertices
        .retain(|vertex| used_vertices.contains(&vertex.id));

    let shells_merged = merge_connected_shells(&mut result);
    let mut oriented_outward = false;
    if edges_sewn > 0 {
        // Components can arrive arbitrarily flipped; make the coedge
        // directions coherent across the shared edges, then restore outward
        // normals by signed volume.
        orient_open_solid_faces(&mut result)?;
        if let Ok(volume) = solid_signed_volume(&result) {
            if volume < 0.0 {
                flip_all_faces(&mut result)?;
            }
            oriented_outward = volume.abs() > tolerance * tolerance * tolerance;
        }
    }
    // Re-derive genus from the Euler characteristic (V − E + F − H = 2 − 2g,
    // H counting inner loops, matching validate()) so validation sees the
    // merged topology, not the input's bookkeeping. Only when sewing actually
    // changed the topology — a no-op must not touch a valid genus.
    if edges_sewn > 0 || shells_merged > 0 {
        let vertex_count = result.vertices.len() as i64;
        let edge_count = result.edges.iter().filter(|edge| !edge.degenerate).count() as i64;
        let face_count = result
            .shells
            .iter()
            .map(|shell| shell.faces.len())
            .sum::<usize>() as i64;
        let ring_count = result
            .shells
            .iter()
            .flat_map(|shell| &shell.faces)
            .map(|face| face.loops.len().saturating_sub(1))
            .sum::<usize>() as i64;
        let numerator = 2 - (vertex_count - edge_count + face_count - ring_count);
        if numerator >= 0 && numerator % 2 == 0 {
            result.genus = numerator / 2;
        }
    }
    let open_edges_after = one_use_edge_ids(&result).len();
    let issues = result
        .validate()
        .into_iter()
        .map(|issue| issue.message)
        .collect();
    Ok((
        result,
        SewReport {
            edges_sewn,
            shells_merged,
            open_edges_before,
            open_edges_after,
            oriented_outward,
            issues,
        },
        pairs_examined,
    ))
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::topology::VertexRecord;
    use crate::{make_arc, make_box_brep, make_cylinder_brep, solid_mass_properties, Vec4};

    fn max_id(solid: &BrepSolid) -> u64 {
        let vertex_max = solid.vertices.iter().map(|vertex| vertex.id).max();
        let edge_max = solid.edges.iter().map(|edge| edge.id).max();
        let face_max = solid
            .shells
            .iter()
            .flat_map(|shell| &shell.faces)
            .flat_map(|face| {
                std::iter::once(face.id).chain(face.loops.iter().flat_map(|loop_record| {
                    std::iter::once(loop_record.id)
                        .chain(loop_record.coedges.iter().map(|coedge| coedge.id))
                }))
            })
            .max();
        [vertex_max, edge_max, face_max, Some(solid.id)]
            .into_iter()
            .flatten()
            .max()
            .unwrap_or(0)
    }

    /// Move the given faces into a NEW shell whose edges and vertices are
    /// fresh duplicates — the un-sewn state an importer or a face-group
    /// detachment produces.
    fn detach_faces(solid: &BrepSolid, face_ids: &[u64]) -> BrepSolid {
        let mut result = solid.clone();
        let mut next = max_id(&result) + 1;
        let mut moved = Vec::new();
        for shell in &mut result.shells {
            let mut kept = Vec::new();
            for face in shell.faces.drain(..) {
                if face_ids.contains(&face.id) {
                    moved.push(face);
                } else {
                    kept.push(face);
                }
            }
            shell.faces = kept;
        }
        assert_eq!(moved.len(), face_ids.len(), "all faces found");
        let mut edge_map = HashMap::<u64, u64>::default();
        let mut vertex_map = HashMap::<u64, u64>::default();
        for face in &mut moved {
            for coedge in face
                .loops
                .iter_mut()
                .flat_map(|loop_record| &mut loop_record.coedges)
            {
                if let Some(&mapped) = edge_map.get(&coedge.edge_id) {
                    coedge.edge_id = mapped;
                    continue;
                }
                let mut duplicate = result
                    .edges
                    .iter()
                    .find(|edge| edge.id == coedge.edge_id)
                    .expect("edge exists")
                    .clone();
                for vertex_id in [&mut duplicate.start_vertex_id, &mut duplicate.end_vertex_id] {
                    if let Some(&mapped) = vertex_map.get(vertex_id) {
                        *vertex_id = mapped;
                        continue;
                    }
                    let point = result
                        .vertices
                        .iter()
                        .find(|vertex| vertex.id == *vertex_id)
                        .expect("vertex exists")
                        .point;
                    result.vertices.push(VertexRecord { id: next, point });
                    vertex_map.insert(*vertex_id, next);
                    *vertex_id = next;
                    next += 1;
                }
                duplicate.id = next;
                next += 1;
                edge_map.insert(coedge.edge_id, duplicate.id);
                coedge.edge_id = duplicate.id;
                result.edges.push(duplicate);
            }
        }
        result.shells.push(ShellRecord {
            id: next,
            faces: moved,
        });
        // Originals whose every use moved (edges interior to the moved group)
        // are now orphaned — purge them as a real detachment would.
        let used_edges = result
            .shells
            .iter()
            .flat_map(|shell| &shell.faces)
            .flat_map(|face| &face.loops)
            .flat_map(|loop_record| &loop_record.coedges)
            .map(|coedge| coedge.edge_id)
            .collect::<HashSet<_>>();
        result.edges.retain(|edge| used_edges.contains(&edge.id));
        let used_vertices = result
            .edges
            .iter()
            .flat_map(|edge| [edge.start_vertex_id, edge.end_vertex_id])
            .collect::<HashSet<_>>();
        result
            .vertices
            .retain(|vertex| used_vertices.contains(&vertex.id));
        result
    }

    /// Translate every face surface, edge curve, and vertex belonging to the
    /// given shell — its geometry is fully duplicated by `detach_faces`, so
    /// nothing shared moves.
    fn translate_shell(solid: &mut BrepSolid, shell_index: usize, delta: Vec3) {
        let mut edge_ids = HashSet::default();
        for face in &mut solid.shells[shell_index].faces {
            for row in &mut face.surface.control_points {
                for control in row.iter_mut() {
                    let point = control.point().unwrap().add(delta);
                    *control = Vec4::from_point(point, control.w);
                }
            }
            for coedge in face
                .loops
                .iter_mut()
                .flat_map(|loop_record| &mut loop_record.coedges)
            {
                edge_ids.insert(coedge.edge_id);
            }
        }
        let mut vertex_ids = HashSet::default();
        for edge in &mut solid.edges {
            if !edge_ids.contains(&edge.id) {
                continue;
            }
            for control in &mut edge.curve.control_points {
                let point = control.point().unwrap().add(delta);
                *control = Vec4::from_point(point, control.w);
            }
            vertex_ids.insert(edge.start_vertex_id);
            vertex_ids.insert(edge.end_vertex_id);
        }
        for vertex in &mut solid.vertices {
            if vertex_ids.contains(&vertex.id) {
                vertex.point = vertex.point.add(delta);
            }
        }
    }

    fn one_use_count(solid: &BrepSolid) -> usize {
        one_use_edge_ids(solid).len()
    }

    /// Faces whose sampled centroid x exceeds the split coordinate.
    fn faces_beyond_x(solid: &BrepSolid, x: f64) -> Vec<u64> {
        solid
            .shells
            .iter()
            .flat_map(|shell| &shell.faces)
            .filter(|face| {
                face.surface
                    .evaluate(0.5, 0.5)
                    .map(|point| point.x > x)
                    .unwrap_or(false)
            })
            .map(|face| face.id)
            .collect()
    }

    #[test]
    fn sews_detached_box_faces_into_watertight_solid() {
        let solid = make_box_brep(Vec3::default(), 4.0, 3.0, 2.0).unwrap();
        let moved = faces_beyond_x(&solid, 3.9); // the +x face
        assert_eq!(moved.len(), 1);
        let detached = detach_faces(&solid, &moved);
        assert_eq!(detached.shells.len(), 2);
        assert_eq!(one_use_count(&detached), 8, "4 rim edges duplicated");

        let (sewn, report) = sew_solid(&detached, 1e-6).unwrap();
        assert_eq!(report.edges_sewn, 4, "{report:?}");
        assert_eq!(report.shells_merged, 1);
        assert_eq!(report.open_edges_after, 0);
        assert!(report.issues.is_empty(), "{:?}", report.issues);
        assert_eq!(sewn.shells.len(), 1);
        assert_eq!(one_use_count(&sewn), 0);
        let volume = solid_mass_properties(&sewn).unwrap().volume;
        assert!((volume - 24.0).abs() < 1e-9, "{volume}");
    }

    #[test]
    fn sews_a_flipped_component_and_orients_outward() {
        let solid = make_box_brep(Vec3::default(), 4.0, 3.0, 2.0).unwrap();
        let moved = faces_beyond_x(&solid, 3.9);
        let mut detached = detach_faces(&solid, &moved);
        // The detached component arrives with the opposite orientation
        // convention — sewing must not require any pre-agreement.
        crate::offset_shell::flip_shell_faces(&mut detached.shells[1]).unwrap();

        let (sewn, report) = sew_solid(&detached, 1e-6).unwrap();
        assert_eq!(report.edges_sewn, 4, "{report:?}");
        assert!(report.issues.is_empty(), "{:?}", report.issues);
        assert_eq!(one_use_count(&sewn), 0);
        let volume = solid_signed_volume(&sewn).unwrap();
        assert!((volume - 24.0).abs() < 1e-9, "outward positive: {volume}");
    }

    #[test]
    fn sews_a_detached_cap_along_its_closed_rim() {
        let axis = Vec3::new(0.0, 0.0, 1.0);
        let solid = make_cylinder_brep(Vec3::default(), axis, 1.5, 3.0).unwrap();
        // Top cap: the face whose centroid sits highest along the axis.
        let top = solid.shells[0]
            .faces
            .iter()
            .max_by(|first, second| {
                let height = |face: &FaceRecord| {
                    face.surface
                        .evaluate(0.5, 0.5)
                        .map(|point| point.z)
                        .unwrap_or(f64::NEG_INFINITY)
                };
                height(first).total_cmp(&height(second))
            })
            .unwrap()
            .id;
        let detached = detach_faces(&solid, &[top]);
        assert_eq!(one_use_count(&detached), 2, "one closed rim duplicated");

        let (sewn, report) = sew_solid(&detached, 1e-6).unwrap();
        assert_eq!(report.edges_sewn, 1, "{report:?}");
        assert_eq!(report.open_edges_after, 0);
        assert!(report.issues.is_empty(), "{:?}", report.issues);
        assert_eq!(sewn.shells.len(), 1);
        let volume = solid_mass_properties(&sewn).unwrap().volume;
        let expected = std::f64::consts::PI * 1.5 * 1.5 * 3.0;
        assert!((volume - expected).abs() < 1e-6, "{volume} vs {expected}");
    }

    #[test]
    fn sews_a_multi_face_group_with_shared_corner_vertices() {
        // Three faces meeting at a corner: their duplicated rims share
        // duplicated vertices, so several sews must chain consistently.
        let solid = make_box_brep(Vec3::default(), 4.0, 3.0, 2.0).unwrap();
        let moved = solid.shells[0]
            .faces
            .iter()
            .filter(|face| {
                face.surface
                    .evaluate(0.5, 0.5)
                    .map(|point| point.x > 3.9 || point.y > 2.9 || point.z > 1.9)
                    .unwrap_or(false)
            })
            .map(|face| face.id)
            .collect::<Vec<_>>();
        assert_eq!(moved.len(), 3);
        let detached = detach_faces(&solid, &moved);

        let (sewn, report) = sew_solid(&detached, 1e-6).unwrap();
        assert_eq!(report.open_edges_after, 0, "{report:?}");
        assert!(report.issues.is_empty(), "{:?}", report.issues);
        assert_eq!(sewn.shells.len(), 1);
        let volume = solid_mass_properties(&sewn).unwrap().volume;
        assert!((volume - 24.0).abs() < 1e-9, "{volume}");
    }

    #[test]
    fn sews_reparametrized_duplicates_via_pcurve_refit() {
        // The duplicated rim edges arrive with REVERSED curves (same locus,
        // different parametrization): the fast keep-the-pcurve path cannot
        // apply, so the rebind must refit pcurves and derive `forward` from
        // the 3D walk alone.
        let solid = make_box_brep(Vec3::default(), 4.0, 3.0, 2.0).unwrap();
        let moved = faces_beyond_x(&solid, 3.9);
        let mut detached = detach_faces(&solid, &moved);
        let duplicated = detached.shells[1]
            .faces
            .iter()
            .flat_map(|face| &face.loops)
            .flat_map(|loop_record| &loop_record.coedges)
            .map(|coedge| coedge.edge_id)
            .collect::<HashSet<_>>();
        for edge in &mut detached.edges {
            if duplicated.contains(&edge.id) {
                edge.curve = edge.curve.reversed().unwrap();
                [edge.t0, edge.t1] = edge.curve.domain().unwrap();
                std::mem::swap(&mut edge.start_vertex_id, &mut edge.end_vertex_id);
            }
        }

        let (sewn, report) = sew_solid(&detached, 1e-6).unwrap();
        assert_eq!(report.edges_sewn, 4, "{report:?}");
        assert_eq!(report.open_edges_after, 0);
        assert!(report.issues.is_empty(), "{:?}", report.issues);
        let volume = solid_mass_properties(&sewn).unwrap().volume;
        assert!((volume - 24.0).abs() < 1e-9, "{volume}");
    }

    #[test]
    fn sews_a_holed_face_with_ring_aware_genus() {
        // Block with a vertical through-hole: the top face carries an INNER
        // loop (the hole rim). Detaching and re-sewing it exercises closed-rim
        // pairing on a ringed face and the ring-aware (V−E+F−H) genus
        // recompute — an outer-loops-only Euler count would corrupt genus
        // here and trip validation.
        let block = make_box_brep(Vec3::default(), 4.0, 4.0, 2.0).unwrap();
        let drill = make_cylinder_brep(
            Vec3::new(2.0, 2.0, -1.0),
            Vec3::new(0.0, 0.0, 1.0),
            1.0,
            4.0,
        )
        .unwrap();
        let holed = crate::boolean_operation(
            &block,
            &drill,
            crate::BooleanOperation::Subtract,
            &crate::BooleanOptions::default(),
        )
        .unwrap();
        let top = holed
            .shells
            .iter()
            .flat_map(|shell| &shell.faces)
            .find(|face| {
                face.loops.len() > 1
                    && face
                        .surface
                        .evaluate(0.05, 0.05)
                        .map(|point| point.z > 1.9)
                        .unwrap_or(false)
            })
            .expect("holed top face")
            .id;
        let detached = detach_faces(&holed, &[top]);

        let (sewn, report) = sew_solid(&detached, 1e-6).unwrap();
        assert_eq!(report.open_edges_after, 0, "{report:?}");
        assert!(report.issues.is_empty(), "{:?}", report.issues);
        assert_eq!(sewn.shells.len(), 1);
        let expected = 4.0 * 4.0 * 2.0 - std::f64::consts::PI * 2.0;
        let volume = solid_mass_properties(&sewn).unwrap().volume;
        assert!((volume - expected).abs() < 1e-6, "{volume} vs {expected}");
    }

    #[test]
    fn watertight_input_is_a_no_op() {
        let solid = make_box_brep(Vec3::default(), 4.0, 3.0, 2.0).unwrap();
        let (sewn, report) = sew_solid(&solid, 1e-6).unwrap();
        assert_eq!(report.edges_sewn, 0);
        assert_eq!(report.shells_merged, 0);
        assert_eq!(report.open_edges_before, 0);
        assert_eq!(report.open_edges_after, 0);
        assert!(report.issues.is_empty(), "{:?}", report.issues);
        let volume = solid_mass_properties(&sewn).unwrap().volume;
        assert!((volume - 24.0).abs() < 1e-9, "{volume}");
    }

    #[test]
    fn reports_an_unsewable_gap_honestly() {
        let solid = make_box_brep(Vec3::default(), 4.0, 3.0, 2.0).unwrap();
        let moved = faces_beyond_x(&solid, 3.9);
        let mut detached = detach_faces(&solid, &moved);
        translate_shell(&mut detached, 1, Vec3::new(0.5, 0.0, 0.0));

        let open_before = one_use_count(&detached);
        let (sewn, report) = sew_solid(&detached, 1e-4).unwrap();
        assert_eq!(report.edges_sewn, 0, "{report:?}");
        assert_eq!(report.shells_merged, 0);
        assert_eq!(report.open_edges_after, open_before);
        assert_eq!(sewn.shells.len(), 2, "nothing force-welded");
    }

    /// The cap face of a z-up cylinder: the one whose surface centroid sits
    /// highest along the axis.
    fn top_cap_face_id(solid: &BrepSolid) -> u64 {
        solid.shells[0]
            .faces
            .iter()
            .max_by(|first, second| {
                let height = |face: &FaceRecord| {
                    face.surface
                        .evaluate(0.5, 0.5)
                        .map(|point| point.z)
                        .unwrap_or(f64::NEG_INFINITY)
                };
                height(first).total_cmp(&height(second))
            })
            .unwrap()
            .id
    }

    /// A cylinder whose detached cap carries a SEAM-ROTATED copy of the shared
    /// rim: the same circle, rebuilt as a full turn starting 0.7 rad round, so
    /// the two rims share a locus and share NO point.
    ///
    /// This is what a periodic wall meeting an independently-authored cap
    /// actually looks like, and it is the state `apply_rebind`'s fence was
    /// written for — "Only coincident endpoints weld — a rotated closed rim
    /// keeps its own seam vertex" — and which nothing drove, because every
    /// other fixture duplicates a rim TOGETHER WITH its seam, so the seams
    /// coincide and the weld always fires.
    struct SeamRotatedCap {
        solid: BrepSolid,
        /// The cap's copy of the rim — the edge whose seam was moved.
        rim_id: u64,
        /// Where the wall's rim still starts.
        original_seam: Vec3,
        /// Where the cap's rim now starts.
        rotated_seam: Vec3,
    }

    fn seam_rotated_cap() -> SeamRotatedCap {
        const SEAM_ROTATION: f64 = 0.7;
        let axis = Vec3::new(0.0, 0.0, 1.0);
        let solid = make_cylinder_brep(Vec3::default(), axis, 1.5, 3.0).unwrap();
        let top = top_cap_face_id(&solid);
        let mut detached = detach_faces(&solid, &[top]);
        assert_eq!(one_use_count(&detached), 2, "one closed rim duplicated");

        // The cap's own copy of the rim: the closed one-use edge inside the
        // detached shell.
        let cap_edges = detached.shells[1]
            .faces
            .iter()
            .flat_map(|face| &face.loops)
            .flat_map(|loop_record| &loop_record.coedges)
            .map(|coedge| coedge.edge_id)
            .collect::<HashSet<_>>();
        let rim_id = *cap_edges
            .iter()
            .find(|id| {
                detached
                    .edges
                    .iter()
                    .find(|edge| edge.id == **id)
                    .is_some_and(|edge| {
                        edge.start_vertex_id == edge.end_vertex_id && !edge.degenerate
                    })
            })
            .expect("the cap carries one closed rim");

        // Recover the circle's frame from the rim itself rather than assuming
        // where `make_cylinder_brep` put it: four evenly-spaced samples average
        // to the centre, and the quarter-turn sample gives an orthogonal axis.
        let (centre, x_axis, y_axis, radius) = {
            let rim = detached
                .edges
                .iter()
                .find(|edge| edge.id == rim_id)
                .unwrap();
            let sample = |fraction: f64| curve_point(rim, fraction).unwrap();
            let centre = sample(0.0)
                .add(sample(0.25))
                .add(sample(0.5))
                .add(sample(0.75))
                .scale(0.25);
            let radial = sample(0.0).sub(centre);
            let quarter = sample(0.25).sub(centre);
            (
                centre,
                radial.normalized().unwrap(),
                quarter.normalized().unwrap(),
                radial.length(),
            )
        };
        let rotated = make_arc(
            centre,
            x_axis,
            y_axis,
            radius,
            SEAM_ROTATION,
            SEAM_ROTATION + std::f64::consts::TAU,
        )
        .unwrap();
        let rotated_domain = rotated.domain().unwrap();
        let rotated_seam = rotated.evaluate(rotated_domain[0]).unwrap();
        let original_seam = {
            let rim = detached
                .edges
                .iter()
                .find(|edge| edge.id == rim_id)
                .unwrap();
            curve_point(rim, 0.0).unwrap()
        };
        assert!(
            rotated_seam.sub(original_seam).length() > 0.5,
            "the fixture must actually move the seam or it proves nothing"
        );

        // Re-seat the cap's rim: new curve, new domain, and its seam vertex
        // moved to the new start point (that vertex is the cap's own duplicate,
        // referenced by nothing else).
        let seam_vertex_id;
        {
            let rim = detached
                .edges
                .iter_mut()
                .find(|edge| edge.id == rim_id)
                .unwrap();
            seam_vertex_id = rim.start_vertex_id;
            rim.curve = rotated.clone();
            [rim.t0, rim.t1] = rotated_domain;
        }
        for vertex in &mut detached.vertices {
            if vertex.id == seam_vertex_id {
                vertex.point = rotated_seam;
            }
        }
        // The cap's pcurve must trace the NEW parametrization, exactly the way
        // `plan_rebind` derives one.
        for face in &mut detached.shells[1].faces {
            for loop_record in &mut face.loops {
                for coedge in &mut loop_record.coedges {
                    if coedge.edge_id != rim_id {
                        continue;
                    }
                    let traversal = if coedge.forward {
                        rotated.clone()
                    } else {
                        rotated.reversed().unwrap()
                    };
                    coedge.pcurve = build_pcurve_on_surface(&face.surface, &traversal).unwrap();
                }
            }
        }

        SeamRotatedCap {
            solid: detached,
            rim_id,
            original_seam,
            rotated_seam,
        }
    }

    #[test]
    fn seam_rotated_closed_rim_sews_without_dragging_a_seam_vertex() {
        // Three things must hold at once: closed-rim pairing (which reads the
        // locus alone, never a vertex) must still find the pair; `plan_rebind`
        // must refit the pcurve because `same_parametrization` is false; and
        // `apply_rebind` must leave the seam vertices unwelded rather than
        // dragging the wall's seam 0.7 rad round the circle onto the cap's.
        let SeamRotatedCap {
            solid: detached,
            rim_id,
            original_seam,
            rotated_seam,
        } = seam_rotated_cap();

        let (sewn, report) = sew_solid(&detached, 1e-6).unwrap();
        assert_eq!(report.edges_sewn, 1, "{report:?}");
        assert_eq!(report.open_edges_after, 0, "{report:?}");
        assert!(report.issues.is_empty(), "{:?}", report.issues);
        assert_eq!(sewn.shells.len(), 1);
        let volume = solid_mass_properties(&sewn).unwrap().volume;
        let expected = std::f64::consts::PI * 1.5 * 1.5 * 3.0;
        assert!((volume - expected).abs() < 1e-6, "{volume} vs {expected}");

        // The point of the fence: the surviving seam vertex is the WALL's,
        // still at angle 0, and nothing was dragged 0.7 rad round to meet the
        // cap's. The cap's rotated seam vertex is simply unreferenced now and
        // dropped with the removed duplicate edge.
        assert!(
            sewn.vertices
                .iter()
                .any(|vertex| vertex.point.sub(original_seam).length() < 1e-12),
            "the wall's seam vertex must be untouched"
        );
        assert!(
            !sewn.vertices
                .iter()
                .any(|vertex| vertex.point.sub(rotated_seam).length() < 1e-9),
            "the rotated seam vertex must be dropped, not welded onto"
        );

        // And the SAME result with the cap's rim first in edge order, which is
        // what makes `orient_pair` load-bearing rather than decorative. The
        // scan keeps `first` and removes `second`, so this ordering asks for
        // the wall's rim to be the one deleted — and the wall's seam vertex is
        // shared with the cylinder's vertical seam edge, so deleting it tears
        // the wall's loop open. Before `orient_pair` this returned
        // `edges_sewn: 1, open_edges_after: 0` over a solid carrying four
        // validation issues ("loop 100 of face 105 is open between coedges 102
        // and 103", a 1.028693-vs-0.079287 pcurve deviation, and an Euler
        // characteristic of 3 where 2 was expected). Now the pair is oriented
        // the only way that closes both loops, and the answer is
        // order-independent.
        let mut cap_rim_first = detached.clone();
        let at = cap_rim_first
            .edges
            .iter()
            .position(|edge| edge.id == rim_id)
            .unwrap();
        let rim = cap_rim_first.edges.remove(at);
        cap_rim_first.edges.insert(0, rim);
        let (flipped, flipped_report) = sew_solid(&cap_rim_first, 1e-6).unwrap();
        assert_eq!(flipped_report.edges_sewn, 1, "{flipped_report:?}");
        assert_eq!(flipped_report.open_edges_after, 0, "{flipped_report:?}");
        assert!(
            flipped_report.issues.is_empty(),
            "edge order must not decide whether the sew is valid: {:?}",
            flipped_report.issues
        );
        let flipped_volume = solid_mass_properties(&flipped).unwrap().volume;
        assert!(
            (flipped_volume - expected).abs() < 1e-6,
            "{flipped_volume} vs {expected}"
        );
    }

    #[test]
    fn a_seam_rotated_pair_whose_rebind_cannot_be_planned_is_blocked_once() {
        // `blocked` is keyed UNORDERED because `orient_pair` may swap which of
        // the two rims is kept. Keyed by (keep, remove) and looked up as
        // (first, second), a swapped pair whose rebind plan fails is never
        // recognised as blocked: nothing was mutated, so the next outer
        // iteration re-offers the same pair, re-swaps it, fails to plan it
        // again, and re-inserts the key it already holds — forever. That is
        // the one promise the `blocked` set exists to keep ("left open rather
        // than retried forever") and the swap quietly broke it.
        //
        // The fixture is the seam-rotated cap with the cap's rim FIRST in edge
        // order, so `orient_pair` swaps (the cap's seam vertex bounds nothing
        // else; the wall's anchors the vertical seam edge), plus an
        // unevaluable pcurve on the coedge `plan_rebind` must walk. An extra
        // knot breaks the knots-imply-control-points invariant; the solid is
        // then round-tripped through JSON because `NurbsCurve::validated` is
        // `#[serde(skip)]`, so a deserialized curve re-runs the full checks on
        // first use — which is also the realistic route in: a saved part file
        // carrying a malformed curve-on-surface.
        //
        // Run on a worker thread with a clock, because the regression is a
        // HANG, and a hang is not a test failure unless someone is timing it.
        let SeamRotatedCap {
            solid: mut fixture,
            rim_id,
            ..
        } = seam_rotated_cap();
        let at = fixture
            .edges
            .iter()
            .position(|edge| edge.id == rim_id)
            .unwrap();
        let rim = fixture.edges.remove(at);
        fixture.edges.insert(0, rim);
        // `orient_pair` keeps the wall's rim and removes the cap's, so
        // `plan_rebind` walks the CAP's coedge. Break that one.
        let mut broken = 0usize;
        for face in &mut fixture.shells[1].faces {
            for loop_record in &mut face.loops {
                for coedge in &mut loop_record.coedges {
                    if coedge.edge_id == rim_id {
                        coedge.pcurve.knots.push(2.0);
                        broken += 1;
                    }
                }
            }
        }
        assert_eq!(broken, 1, "exactly one cap coedge rides the rotated rim");
        let fixture: BrepSolid =
            serde_json::from_str(&serde_json::to_string(&fixture).unwrap()).unwrap();
        assert!(
            fixture.shells[1].faces[0].loops[0].coedges[0]
                .pcurve
                .domain()
                .is_err(),
            "the fixture must actually make plan_rebind fail"
        );

        let (sender, receiver) = std::sync::mpsc::channel();
        std::thread::spawn(move || {
            let _ = sender.send(sew_solid(&fixture, 1e-6).map(|(_solid, report)| report));
        });
        let report = receiver
            .recv_timeout(std::time::Duration::from_secs(30))
            .expect("sew_solid must terminate on a pair it cannot plan, not spin on it")
            .expect("sew_solid itself does not error here");
        assert_eq!(report.edges_sewn, 0, "{report:?}");
        assert_eq!(
            report.open_edges_after, report.open_edges_before,
            "an unplannable pair stays open and is counted: {report:?}"
        );
    }

    // -----------------------------------------------------------------
    // Spatial-hash prefilter: parity with the brute-force scan, and scale
    // -----------------------------------------------------------------

    /// Every id in `solid` shifted by `offset`, so two copies can live in one
    /// `BrepSolid` without colliding. Geometry is untouched.
    fn shift_ids(solid: &BrepSolid, offset: u64) -> BrepSolid {
        let mut copy = solid.clone();
        copy.id += offset;
        for vertex in &mut copy.vertices {
            vertex.id += offset;
        }
        for edge in &mut copy.edges {
            edge.id += offset;
            edge.start_vertex_id += offset;
            edge.end_vertex_id += offset;
        }
        for shell in &mut copy.shells {
            shell.id += offset;
            for face in &mut shell.faces {
                face.id += offset;
                for loop_record in &mut face.loops {
                    loop_record.id += offset;
                    for coedge in &mut loop_record.coedges {
                        coedge.id += offset;
                        coedge.edge_id += offset;
                    }
                }
            }
        }
        copy
    }

    /// A second, geometrically identical copy of `shell_index`'s faces with
    /// fresh ids — THREE coincident rims where there were two. With three
    /// mutually-sewable copies the pair search must pick a winner, so
    /// first-match-wins stops being an implementation detail and becomes an
    /// observable the parity comparison can catch.
    fn duplicate_shell(solid: &BrepSolid, shell_index: usize, offset: u64) -> BrepSolid {
        let mut result = solid.clone();
        let mut faces = result.shells[shell_index].faces.clone();
        let edge_ids = faces
            .iter()
            .flat_map(|face| &face.loops)
            .flat_map(|loop_record| &loop_record.coedges)
            .map(|coedge| coedge.edge_id)
            .collect::<HashSet<_>>();
        let mut edges = result
            .edges
            .iter()
            .filter(|edge| edge_ids.contains(&edge.id))
            .cloned()
            .collect::<Vec<_>>();
        let vertex_ids = edges
            .iter()
            .flat_map(|edge| [edge.start_vertex_id, edge.end_vertex_id])
            .collect::<HashSet<_>>();
        let mut vertices = result
            .vertices
            .iter()
            .filter(|vertex| vertex_ids.contains(&vertex.id))
            .cloned()
            .collect::<Vec<_>>();
        for vertex in &mut vertices {
            vertex.id += offset;
        }
        for edge in &mut edges {
            edge.id += offset;
            edge.start_vertex_id += offset;
            edge.end_vertex_id += offset;
        }
        for face in &mut faces {
            face.id += offset;
            for loop_record in &mut face.loops {
                loop_record.id += offset;
                for coedge in &mut loop_record.coedges {
                    coedge.id += offset;
                    coedge.edge_id += offset;
                }
            }
        }
        let shell_id = result.shells[shell_index].id + offset;
        result.vertices.extend(vertices);
        result.edges.extend(edges);
        result.shells.push(ShellRecord {
            id: shell_id,
            faces,
        });
        result
    }

    /// Deterministic Fisher-Yates over `solid.edges`. The candidate list is
    /// `edges` order filtered, and first-match-wins reads that order, so a
    /// permutation is the sharpest adversary a prefilter has: it must offer
    /// the same candidates in the same ASCENDING positions no matter how the
    /// edge vector is laid out.
    fn permute_edges(solid: &BrepSolid, seed: u64) -> BrepSolid {
        let mut result = solid.clone();
        let mut state = seed | 1;
        for index in (1..result.edges.len()).rev() {
            state = state
                .wrapping_mul(6364136223846793005)
                .wrapping_add(1442695040888963407);
            let pick = ((state >> 33) as usize) % (index + 1);
            result.edges.swap(index, pick);
        }
        result
    }

    /// The complete observable outcome of a sew: the serialized solid plus the
    /// report, or the error text. Bit-equality of this string between the two
    /// strategies is the parity claim.
    fn sew_outcome(solid: &BrepSolid, tolerance: f64, search: PairSearch) -> String {
        match sew_solid_with_search(solid, tolerance, search) {
            Ok((sewn, report, _examined)) => format!(
                "ok\n{}\n{report:?}",
                serde_json::to_string(&sewn).expect("a sewn solid serializes")
            ),
            Err(message) => format!("err {message}"),
        }
    }

    /// Every box of `count` fully detached — each face in its own shell with
    /// its own duplicate edges and vertices — the boxes 100 units apart so no
    /// pair can sew across boxes. Six faces per box, 24 one-use edges, 12
    /// sewable pairs.
    fn detached_box_soup(count: usize) -> BrepSolid {
        let unit = {
            let mut soup = make_box_brep(Vec3::default(), 4.0, 3.0, 2.0).unwrap();
            let face_ids = soup.shells[0]
                .faces
                .iter()
                .map(|face| face.id)
                .collect::<Vec<_>>();
            // All but one: an emptied shell is not a state any real importer
            // produces, and validate() would rightly object to it.
            for face_id in face_ids.iter().skip(1) {
                soup = detach_faces(&soup, &[*face_id]);
            }
            soup
        };
        let mut soup = BrepSolid {
            id: 1,
            vertices: Vec::new(),
            edges: Vec::new(),
            shells: Vec::new(),
            genus: 0,
        };
        for index in 0..count {
            let mut copy = shift_ids(&unit, 1_000 * (index as u64 + 1));
            for shell_index in 0..copy.shells.len() {
                translate_shell(&mut copy, shell_index, Vec3::new(index as f64 * 100.0, 0.0, 0.0));
            }
            soup.vertices.extend(copy.vertices);
            soup.edges.extend(copy.edges);
            soup.shells.extend(copy.shells);
        }
        soup
    }

    #[test]
    fn endpoint_grid_offers_a_partner_exactly_at_tolerance_across_a_cell_boundary() {
        // Two endpoints exactly `tolerance` apart, placed so that they fall in
        // DIFFERENT cells: `a` sits a hair below a cell boundary, `b` a full
        // tolerance beyond it. This is the case a cell-per-tolerance grid with
        // a naive same-cell lookup would drop.
        let tolerance = 1e-6;
        let cell = 2.0 * tolerance;
        let boundary = cell * 7.0;
        let a = Vec3::new(boundary - 1e-13, 0.0, 0.0);
        let b = Vec3::new(a.x + tolerance, 0.0, 0.0);
        assert_ne!(
            grid_cell(a, cell).unwrap(),
            grid_cell(b, cell).unwrap(),
            "the fixture must actually straddle a boundary or it proves nothing"
        );
        let grid = EndpointGrid::build(
            &[Some([a, Vec3::new(1.0, 0.0, 0.0)]), Some([b, b])],
            &[false, false],
            tolerance,
        )
        .expect("bucketable");
        let mut found = Vec::new();
        assert!(grid.near(a, &mut found));
        assert!(found.contains(&1), "partner across the boundary: {found:?}");
        assert!(grid.near(b, &mut found));
        assert!(found.contains(&0), "and symmetrically: {found:?}");
    }

    #[test]
    fn endpoint_grid_refuses_to_bucket_unusable_coordinates() {
        // A non-finite or astronomically distant coordinate cannot be keyed;
        // `build` must return None so the caller falls back to the exhaustive
        // scan instead of silently losing that candidate.
        let tolerance = 1e-6;
        assert!(grid_cell(Vec3::new(f64::NAN, 0.0, 0.0), 2e-6).is_none());
        assert!(grid_cell(Vec3::new(1e30, 0.0, 0.0), 2e-6).is_none());
        let far = Vec3::new(1e30, 0.0, 0.0);
        assert!(
            EndpointGrid::build(&[Some([far, far])], &[false], tolerance).is_none(),
            "one unbucketable candidate disables the whole prefilter"
        );
    }

    /// The eight shapes of dirt the parity comparison runs over: a single
    /// detached face, a three-face corner group whose duplicated rims share
    /// corners, a closed cylinder rim (the lane the grid deliberately does NOT
    /// key on a point), reversed duplicates that force a pcurve refit, three
    /// coincident copies where first-match-wins decides the winner, an
    /// unsewable gap, and a dangling vertex reference the exhaustive scan
    /// skips outright — plus the multi-box soup the prefilter exists for.
    fn parity_battery() -> [(&'static str, BrepSolid, f64); 8] {
        let box_solid = make_box_brep(Vec3::default(), 4.0, 3.0, 2.0).unwrap();
        let plus_x = faces_beyond_x(&box_solid, 3.9);
        let one_face = detach_faces(&box_solid, &plus_x);

        let corner_faces = box_solid.shells[0]
            .faces
            .iter()
            .filter(|face| {
                face.surface
                    .evaluate(0.5, 0.5)
                    .map(|point| point.x > 3.9 || point.y > 2.9 || point.z > 1.9)
                    .unwrap_or(false)
            })
            .map(|face| face.id)
            .collect::<Vec<_>>();
        let corner_group = detach_faces(&box_solid, &corner_faces);

        let cylinder = make_cylinder_brep(Vec3::default(), Vec3::new(0.0, 0.0, 1.0), 1.5, 3.0)
            .unwrap();
        let top_cap = cylinder.shells[0]
            .faces
            .iter()
            .max_by(|first, second| {
                let height = |face: &FaceRecord| {
                    face.surface
                        .evaluate(0.5, 0.5)
                        .map(|point| point.z)
                        .unwrap_or(f64::NEG_INFINITY)
                };
                height(first).total_cmp(&height(second))
            })
            .unwrap()
            .id;
        let closed_rim = detach_faces(&cylinder, &[top_cap]);

        let reversed = {
            let mut detached = one_face.clone();
            let duplicated = detached.shells[1]
                .faces
                .iter()
                .flat_map(|face| &face.loops)
                .flat_map(|loop_record| &loop_record.coedges)
                .map(|coedge| coedge.edge_id)
                .collect::<HashSet<_>>();
            for edge in &mut detached.edges {
                if duplicated.contains(&edge.id) {
                    edge.curve = edge.curve.reversed().unwrap();
                    [edge.t0, edge.t1] = edge.curve.domain().unwrap();
                    std::mem::swap(&mut edge.start_vertex_id, &mut edge.end_vertex_id);
                }
            }
            detached
        };

        let three_copies = duplicate_shell(&one_face, 1, 500_000);

        let gapped = {
            let mut detached = one_face.clone();
            translate_shell(&mut detached, 1, Vec3::new(0.5, 0.0, 0.0));
            detached
        };

        let dangling = {
            let mut detached = one_face.clone();
            let victim = detached.edges[0].id;
            for edge in &mut detached.edges {
                if edge.id == victim {
                    // A vertex id no record answers to: the exhaustive scan
                    // destructures `points.get(..)` and `continue`s, so the
                    // prefilter is free to omit it — but only if it does so on
                    // exactly the same pairs.
                    edge.start_vertex_id = 9_999_999;
                }
            }
            detached
        };

        // The soup is where the prefilter earns its keep, so it belongs in the
        // parity battery too, not only in the scaling smoke.
        let soup = detached_box_soup(3);

        [
            ("one detached face", one_face, 1e-6),
            ("three-face corner group", corner_group, 1e-6),
            ("closed cylinder rim", closed_rim, 1e-6),
            ("reversed duplicates", reversed, 1e-6),
            ("three coincident rims", three_copies, 1e-6),
            ("unsewable gap", gapped, 1e-4),
            ("dangling vertex reference", dangling, 1e-6),
            ("three-box soup", soup, 1e-6),
        ]
    }

    /// The hashed prefilter must reproduce the exhaustive scan EXACTLY, over
    /// every shape of dirt and every edge ordering.
    ///
    /// `Exhaustive` is the original pair loop, refactored only to walk
    /// candidate POSITIONS instead of cloned records. That refactor was itself
    /// checked against the pre-prefilter module verbatim from git (a throwaway
    /// `#[path]` copy of it compiled alongside): all 8 cases x 4 permutations
    /// produced byte-identical serialized solids and reports, so AT THAT
    /// COMMIT "hashed == exhaustive" was also "hashed == what shipped". The
    /// shared loop body has since changed on purpose (`orient_pair`), so what
    /// this test guards from here on is the prefilter alone — which is all it
    /// ever needed to guard.
    #[test]
    fn sew_prefilter_parity() {
        for (name, solid, tolerance) in parity_battery() {
            for seed in [0u64, 1, 12345, 0xdead_beef] {
                let permuted = permute_edges(&solid, seed);
                let hashed = sew_outcome(&permuted, tolerance, PairSearch::Hashed);
                let exhaustive = sew_outcome(&permuted, tolerance, PairSearch::Exhaustive);
                assert_eq!(
                    hashed, exhaustive,
                    "{name} (edge order seed {seed}): the prefilter changed the result"
                );
            }
        }
    }

    /// The soup a dirty import actually hands `sew_solid`: half the boxes sew
    /// cleanly, half have every face nudged off its neighbours by 1e-3 — gaps
    /// three orders of magnitude above the sew tolerance, so those rims can
    /// NEVER pair. They sit at the front of the edge order and are therefore
    /// rescanned from position 0 on every accepted pair. That prefix of
    /// permanently-open candidates is where the O(k·n²) scan actually burns,
    /// and it is the shape of a real torn shell or sheet body — not a
    /// synthetic worst case. The fully-sewable soup is the EASY case for the
    /// old scan (candidate 0 finds its partner on the first inner pass), which
    /// is exactly why it must not be the only thing measured.
    fn half_unsewable_box_soup(count: usize) -> BrepSolid {
        let mut soup = detached_box_soup(count);
        let nudged = soup.shells.len() / 2;
        for shell_index in 0..nudged {
            let step = (shell_index + 1) as f64 * 1e-3;
            translate_shell(&mut soup, shell_index, Vec3::new(step, -step, step));
        }
        soup
    }

    #[test]
    fn sew_prefilter_scales_sub_quadratically() {
        // `pairs_examined` is a deterministic function of the input, so a band
        // set from it is a real bound rather than a machine-speed lottery.
        //
        // Measured on this fixture (2026-09-03, cargo test --lib):
        //   fully sewable, 21 boxes / 126 faces /  504 one-use edges:
        //       hashed 546      exhaustive 1995     (3.7x)
        //   fully sewable, 84 boxes / 504 faces / 2016 one-use edges:
        //       hashed 2184     exhaustive 7980     (3.7x)
        //   half unsewable, 12 boxes /  72 faces /  288 one-use edges:
        //       hashed 5412     exhaustive 1509042  (278.8x)
        // The plan's 500-face soup is the 84-box row; the half-unsewable row
        // is where the prefilter earns its keep, because a candidate with no
        // partner anywhere is what makes the inner scan run to the end.
        //
        // Wall clock on the half-unsewable soup (same run, debug build) shows
        // the quadratic term leaving:
        //   12 boxes /  288 edges:  326 ms hashed vs   500 ms exhaustive (1.5x)
        //   24 boxes /  576 edges: 1.09 s  hashed vs  3.20 s  exhaustive (2.9x)
        //   36 boxes /  864 edges: 2.30 s  hashed vs 10.21 s  exhaustive (4.4x)
        // The ratio GROWS with n, which is the point; at the smallest size the
        // per-iteration candidate/vertex-map rebuild still dominates, so the
        // 278x cut in pair examinations only buys 1.5x there. Making that
        // rebuild incremental is a separate follow-up, deliberately not taken
        // here: it cannot be done without touching what the loop observes.
        let small = detached_box_soup(21);
        let large = detached_box_soup(84);
        let (_, small_report, small_examined) =
            sew_solid_with_search(&small, 1e-6, PairSearch::Hashed).unwrap();
        let (_, large_report, large_examined) =
            sew_solid_with_search(&large, 1e-6, PairSearch::Hashed).unwrap();
        assert_eq!(small_report.open_edges_after, 0, "{small_report:?}");
        assert_eq!(large_report.open_edges_after, 0, "{large_report:?}");
        // 4x the boxes is 4x the one-use edges AND 4x the accepted pairs, so a
        // per-iteration O(n) scan would grow 16x and the old O(n²) one 64x.
        // Measured growth is exactly 4.00x — the neighbourhood query is O(1)
        // per candidate, so the total tracks the number of accepted pairs
        // alone. Band 4.5x: the measured value plus room for a hash-order
        // reshuffle, well under the 16x that would signal a lost prefilter.
        assert!(
            (large_examined as f64) <= 4.5 * (small_examined as f64),
            "hashed pair examinations grew super-linearly with the soup: \
             {small_examined} -> {large_examined}"
        );

        // The case the prefilter exists for. Both strategies must still agree
        // (parity is the contract), and the hashed one must examine orders of
        // magnitude fewer pairs.
        let torn = half_unsewable_box_soup(12);
        let (torn_hashed, torn_report, torn_examined) =
            sew_solid_with_search(&torn, 1e-6, PairSearch::Hashed).unwrap();
        let (torn_reference, reference_report, reference_examined) =
            sew_solid_with_search(&torn, 1e-6, PairSearch::Exhaustive).unwrap();
        assert_eq!(
            serde_json::to_string(&torn_hashed).unwrap(),
            serde_json::to_string(&torn_reference).unwrap(),
            "the torn soup must sew identically under both strategies"
        );
        assert_eq!(format!("{torn_report:?}"), format!("{reference_report:?}"));
        assert!(
            torn_report.open_edges_after > 0,
            "the nudged half must stay honestly open: {torn_report:?}"
        );
        // Measured 1_509_042 -> 5_412, a 278.8x reduction. Band 50x: far
        // below the measured win, so it fails loudly if the prefilter degrades
        // towards a full scan, and does not re-break on a tie-order change
        // that shifts either count by a few percent.
        assert!(
            (reference_examined as f64) >= 50.0 * (torn_examined as f64),
            "prefilter did not cut the torn-soup scan: \
             exhaustive {reference_examined} vs hashed {torn_examined}"
        );
    }
}

/// Split PINCHED vertices — points where two (or more) umbrella fans of
/// faces meet at a single vertex record. Local manifold checks (edge use
/// counts, loop closure, orientation) cannot see a pinch; it surfaces only
/// as an odd Euler characteristic. The link of a manifold boundary vertex is
/// a single edge-connected fan: union incident edges through every loop
/// CORNER at the vertex (consecutive coedges meeting there inside one face);
/// more than one component means distinct fans sharing the record — give
/// each extra fan its own vertex at the same point and reassign that fan's
/// edge endpoints. Geometry is untouched; only identity is repaired.
pub fn split_pinched_vertices(solid: &mut BrepSolid) -> Result<usize, String> {
    let mut split_count = 0usize;
    let vertex_ids: Vec<u64> = solid.vertices.iter().map(|vertex| vertex.id).collect();
    let mut next_id = solid
        .vertices
        .iter()
        .map(|vertex| vertex.id)
        .chain(solid.edges.iter().map(|edge| edge.id))
        .max()
        .unwrap_or(0)
        + 1;
    // id -> index maps built once. This routine only mutates edge endpoints
    // in place and appends vertices, so indices never shift: every lookup
    // returns the identical live record the previous `.iter().find(id==)`
    // scans returned. Kills the O(V*coedges*E) nested edge find below.
    let edge_of_id: HashMap<u64, usize> = solid
        .edges
        .iter()
        .enumerate()
        .map(|(index, edge)| (edge.id, index))
        .collect();
    let vertex_of_id: HashMap<u64, usize> = solid
        .vertices
        .iter()
        .enumerate()
        .map(|(index, vertex)| (vertex.id, index))
        .collect();
    for vertex_id in vertex_ids {
        // Incident edges (either endpoint; closed and degenerate included so
        // pole/seam structures stay connected through their corners).
        let incident: Vec<u64> = solid
            .edges
            .iter()
            .filter(|edge| edge.start_vertex_id == vertex_id || edge.end_vertex_id == vertex_id)
            .map(|edge| edge.id)
            .collect();
        if incident.len() < 4 {
            // A pinch needs at least two fans of >= 2 edges each.
            continue;
        }
        let index_of: HashMap<u64, usize> = incident
            .iter()
            .enumerate()
            .map(|(index, id)| (*id, index))
            .collect();
        let mut parent: Vec<usize> = (0..incident.len()).collect();
        fn root(parent: &mut [usize], index: usize) -> usize {
            if parent[index] != index {
                parent[index] = root(parent, parent[index]);
            }
            parent[index]
        }
        let edge_end = |edge_id: u64, forward: bool| -> Option<u64> {
            edge_of_id.get(&edge_id).map(|&index| {
                let edge = &solid.edges[index];
                if forward {
                    edge.end_vertex_id
                } else {
                    edge.start_vertex_id
                }
            })
        };
        for shell in &solid.shells {
            for face in &shell.faces {
                for loop_record in &face.loops {
                    let count = loop_record.coedges.len();
                    for index in 0..count {
                        let current = &loop_record.coedges[index];
                        let next = &loop_record.coedges[(index + 1) % count];
                        // The corner between current and next sits at
                        // current's traversal END vertex.
                        let Some(junction) = edge_end(current.edge_id, current.forward) else {
                            continue;
                        };
                        if junction != vertex_id {
                            continue;
                        }
                        let (Some(&a), Some(&b)) =
                            (index_of.get(&current.edge_id), index_of.get(&next.edge_id))
                        else {
                            continue;
                        };
                        let ra = root(&mut parent, a);
                        let rb = root(&mut parent, b);
                        if ra != rb {
                            parent[rb] = ra;
                        }
                    }
                }
            }
        }
        let mut component_of: HashMap<usize, usize> = HashMap::default();
        let mut components = 0usize;
        let mut assignment: Vec<usize> = vec![0; incident.len()];
        for index in 0..incident.len() {
            let r = root(&mut parent, index);
            let component = *component_of.entry(r).or_insert_with(|| {
                components += 1;
                components - 1
            });
            assignment[index] = component;
        }
        if components < 2 {
            continue;
        }
        // Keep the original record for component 0; every further fan gets a
        // duplicate vertex at the same point.
        let point = vertex_of_id
            .get(&vertex_id)
            .map(|&index| solid.vertices[index].point)
            .ok_or("split_pinched_vertices: vertex vanished")?;
        let mut replacement_ids = vec![vertex_id];
        for _ in 1..components {
            let id = next_id;
            next_id += 1;
            solid.vertices.push(VertexRecord { id, point });
            replacement_ids.push(id);
        }
        for (offset, edge_id) in incident.iter().enumerate() {
            let replacement = replacement_ids[assignment[offset]];
            if replacement == vertex_id {
                continue;
            }
            if let Some(edge) = solid.edges.iter_mut().find(|edge| edge.id == *edge_id) {
                if edge.start_vertex_id == vertex_id {
                    edge.start_vertex_id = replacement;
                }
                if edge.end_vertex_id == vertex_id {
                    edge.end_vertex_id = replacement;
                }
            }
        }
        split_count += components - 1;
    }
    Ok(split_count)
}

#[cfg(test)]
mod pinch_tests {
    use super::*;
    use crate::make_box_brep;

    #[test]
    fn a_hand_welded_corner_pinch_splits_back_to_manifold() {
        // Two boxes meeting at exactly one corner point; welding that corner
        // into ONE vertex record creates the pinch (χ drops odd), and the
        // splitter must give each box its own vertex back.
        let first = make_box_brep(Vec3::default(), 2.0, 2.0, 2.0).unwrap();
        let second = make_box_brep(Vec3::new(2.0, 2.0, 2.0), 2.0, 2.0, 2.0).unwrap();
        // Merge into one two-shell solid with disjoint id spaces.
        let mut solid = first.clone();
        let offset = 1000u64;
        let mut moved = second.clone();
        for vertex in &mut moved.vertices {
            vertex.id += offset;
        }
        for edge in &mut moved.edges {
            edge.id += offset;
            edge.start_vertex_id += offset;
            edge.end_vertex_id += offset;
        }
        for shell in &mut moved.shells {
            shell.id += offset;
            for face in &mut shell.faces {
                face.id += offset;
                for loop_record in &mut face.loops {
                    loop_record.id += offset;
                    for coedge in &mut loop_record.coedges {
                        coedge.id += offset;
                        coedge.edge_id += offset;
                    }
                }
            }
        }
        solid.vertices.extend(moved.vertices);
        solid.edges.extend(moved.edges);
        solid.shells.extend(moved.shells);

        // Find the two coincident corner vertices at (2,2,2) and weld them.
        let corner = Vec3::new(2.0, 2.0, 2.0);
        let ids: Vec<u64> = solid
            .vertices
            .iter()
            .filter(|vertex| vertex.point.sub(corner).length() < 1e-9)
            .map(|vertex| vertex.id)
            .collect();
        assert_eq!(ids.len(), 2, "both corner vertices present");
        let (keep, remove) = (ids[0], ids[1]);
        for edge in &mut solid.edges {
            if edge.start_vertex_id == remove {
                edge.start_vertex_id = keep;
            }
            if edge.end_vertex_id == remove {
                edge.end_vertex_id = keep;
            }
        }
        solid.vertices.retain(|vertex| vertex.id != remove);

        // The welded solid is pinched: local checks pass but χ is odd.
        let split = split_pinched_vertices(&mut solid).unwrap();
        assert_eq!(split, 1, "exactly one fan duplicated");
        assert_eq!(
            solid
                .vertices
                .iter()
                .filter(|vertex| vertex.point.sub(corner).length() < 1e-9)
                .count(),
            2,
            "corner vertices restored"
        );
        // Each shell is a clean box again.
        assert!(solid.validate().is_empty(), "{:?}", solid.validate());
    }

    #[test]
    fn manifold_solids_are_untouched() {
        let mut solid = make_box_brep(Vec3::default(), 2.0, 2.0, 2.0).unwrap();
        let before = solid.vertices.len();
        assert_eq!(split_pinched_vertices(&mut solid).unwrap(), 0);
        assert_eq!(solid.vertices.len(), before);
    }
}