kernel: stitch splits an edge at another sheet's vertex lying on it

Before sewing, each edge of one input is split exactly (Brep::split_edges) at every vertex
of another input that lies on it within the stitch tolerance away from its own ends, so a
patch whose corner meets a body's edge part way (the five-plus-sided fill's midpoints)
sews with no T-junction. The pieces carry their edge's history. The search is the nearest
of 64 samples refined by golden section, over the vertices within a chord's reach of the
edge's samples. Nothing else in stitch changes.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_017itGh5fpMEousnZ91cYcVt

crates/kernel/src/surfaces.rs

@@ -658,6 +658,74 @@ fn carried(h: &History, c: &Canonical, offsets: &Offsets, operands: usize, face:
     Some(History { operands, faces, edges })
 }
 
+/// The joined record `c` with each edge of one input split exactly at every vertex of another
+/// input that lies on it within `tol` (away from the edge's own ends), so that a sheet's corner
+/// on another's edge is a vertex of both and the sheets sew there (`Brep::split_edges`); `c`'s
+/// maps carried through.
+fn without_t_junctions(c: Canonical, offsets: &Offsets, tol: f64) -> Result<Canonical> {
+    let r = &c.brep;
+    let input_of_edge: Vec<Option<usize>> = c.edges.iter().map(|p| owner(offsets, *p, 1).map(|(k, _)| k)).collect();
+    let mut input_of_vertex: Vec<Option<usize>> = vec![None; r.vertices.len()];
+    for (e, k) in r.edges.iter().zip(&input_of_edge) {
+        for v in e.vertices {
+            if let Some(slot) = input_of_vertex.get_mut(v as usize) {
+                *slot = slot.or(*k);
+            }
+        }
+    }
+    let point = |v: u32| r.vertices.get(v as usize).map(|x| Vec3::new(x.point[0], x.point[1], x.point[2]));
+    let mut splits: Vec<(u32, f64)> = Vec::new();
+    for (e, (edge, ke)) in r.edges.iter().zip(&input_of_edge).enumerate() {
+        let Ok(curve) = CurveEval::new(&edge.curve) else { continue };
+        let ends: Vec<Vec3> = edge.vertices.iter().filter_map(|v| point(*v)).collect();
+        // Only vertices within a chord's reach of the edge's samples can lie on it.
+        let [a, z] = edge.interval;
+        let samples: Vec<Vec3> = (0..=SAMPLES).filter_map(|k| curve.point(a + (z - a) * k as f64 / SAMPLES as f64)).collect();
+        let chord = samples.windows(2).map(|w| w[0].dist(w[1])).fold(0.0, f64::max);
+        for (v, kv) in input_of_vertex.iter().enumerate() {
+            let Some(p) = point(index(v)?) else { continue };
+            if kv == ke || ends.iter().any(|q| q.dist(p) <= tol) || samples.iter().all(|q| q.dist(p) > chord + tol) {
+                continue;
+            }
+            if let Some(t) = nearest_parameter(&curve, edge.interval, p).filter(|(_, d)| *d <= tol).map(|(t, _)| t) {
+                splits.push((index(e)?, t));
+            }
+        }
+    }
+    if splits.is_empty() {
+        return Ok(c);
+    }
+    let split = r.split_edges(&splits).map_err(|e| failed(format!("stitch: splitting edges at the other sheets' vertices: {e}")))?;
+    let faces =
+        split.faces.iter().map(|f| c.faces.get(*f as usize).copied().ok_or_else(|| failed("stitch: a face reference"))).collect::<Result<_>>()?;
+    let edges =
+        split.edges.iter().map(|e| c.edges.get(*e as usize).copied().ok_or_else(|| failed("stitch: an edge reference"))).collect::<Result<_>>()?;
+    Ok(Canonical { brep: split.brep, faces, edges })
+}
+
+/// Samples along an edge for finding the points of it near a vertex.
+const SAMPLES: usize = 64;
+
+/// The parameter of the curve's point nearest `p` inside `interval`, and the distance: the
+/// nearest of evenly spaced samples, then a golden-section search between its neighbours.
+fn nearest_parameter(c: &CurveEval, [a, b]: [f64; 2], p: Vec3) -> Option<(f64, f64)> {
+    let dist = |t: f64| c.point(t).map_or(f64::INFINITY, |q| q.dist(p));
+    let at = |k: usize| a + (b - a) * k as f64 / SAMPLES as f64;
+    let best = (0..=SAMPLES).min_by(|&i, &j| dist(at(i)).total_cmp(&dist(at(j))))?;
+    let (mut lo, mut hi) = (at(best.saturating_sub(1)), at((best + 1).min(SAMPLES)));
+    let ratio = 0.5 * (5f64.sqrt() - 1.0);
+    for _ in 0..80 {
+        let (x, y) = (hi - ratio * (hi - lo), lo + ratio * (hi - lo));
+        if dist(x) < dist(y) {
+            hi = y;
+        } else {
+            lo = x;
+        }
+    }
+    let t = 0.5 * (lo + hi);
+    (t > a && t < b).then(|| (t, dist(t)))
+}
+
 /// Sew surface bodies along edges that coincide (within `tol`; the bodies' working tolerance when
 /// not positive): one body, a solid when the faces close up. SMLib turns the faces so that
 /// neighbours agree (outward on a solid). The result's history is over `bodies`: each face the one
@@ -673,7 +741,7 @@ pub fn stitch(bodies: &[&Body], tol: f64) -> Result<Body> {
     let tol = if tol.is_finite() && tol > 0.0 { tol } else { size * 1e-6 + 1e-6 };
     let records: Vec<&Brep> = bodies.iter().map(|b| b.record()).collect();
     let (parts, offsets) = Parts::joined(&records)?;
-    let c = parts.canonical()?;
+    let c = without_t_junctions(parts.canonical()?, &offsets, tol)?;
     let (body, h) = imp::stitched(&c.brep, tol)?;
     let h = h.and_then(|h| carried(&h, &c, &offsets, bodies.len(), |operand, face| FaceOrigin::Face { operand, face }));
     Ok(body.with_checked_history(h, &record::counts(bodies)))