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;
            }
        }
    }
}

/// 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()
}

/// 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> {
    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;
    // Pairs whose rebind plan failed (unprojectable pcurve) — left open
    // rather than retried forever.
    let mut blocked = HashSet::<(u64, u64)>::default();
    loop {
        let one_use = one_use_edge_ids(&result);
        let candidates = result
            .edges
            .iter()
            .filter(|edge| !edge.degenerate && one_use.contains(&edge.id))
            .cloned()
            .collect::<Vec<_>>();
        let points = result
            .vertices
            .iter()
            .map(|vertex| (vertex.id, vertex.point))
            .collect::<HashMap<_, _>>();
        let mut chosen = None;
        'pairs: for (index, first) in candidates.iter().enumerate() {
            for second in candidates.iter().skip(index + 1) {
                if blocked.contains(&(first.id, second.id)) {
                    continue;
                }
                let first_closed = first.start_vertex_id == first.end_vertex_id;
                let second_closed = second.start_vertex_id == second.end_vertex_id;
                if first_closed != second_closed {
                    continue;
                }
                if !first_closed {
                    let (Some(&fs), Some(&fe), Some(&ss), Some(&se)) = (
                        points.get(&first.start_vertex_id),
                        points.get(&first.end_vertex_id),
                        points.get(&second.start_vertex_id),
                        points.get(&second.end_vertex_id),
                    ) 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;
                }
                chosen = Some((first.clone(), second.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((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,
        },
    ))
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::topology::VertexRecord;
    use crate::{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");
    }
}

/// 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);
    }
}