tess: a face's sampled boundary never crosses itself

Each edge is sampled to the chord tolerance on its own, so where two curves of one face run
closer together than their chords sag (the two arcs of a thin lune, the sliver caps of a
joined ring of cylinder segments) the chords crossed although the curves do not, and no
triangulation follows such a boundary ("boundary crosses itself").

After the edges are sampled, every face's boundary is checked for chords that cross another
of the face, in its loop or another, and on a closed surface across whole periods (chords
ending on a pole are left out). The edge segments that cross are halved at their curve
parameter's middle, on the edge, so every face of the edge gets the same new nodes, and the
faces of the edges split are checked again until no chords cross. A boundary still crossing
after 16 rounds, or an edge past its node cap, is an error naming the face and the edges.

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

(cherry picked from commit d2b2bdb8e2720690d79123a16bc18c768b1598e8)

crates/kernel/src/tess/edges.rs

@@ -118,6 +118,50 @@ pub(crate) fn sample_edge(b: &brep::Brep, ei: usize, rules: &Rules, h: f64) -> R
     Ok(Nodes { points: out.iter().map(|p| p.1).collect(), params: out.iter().map(|p| p.0).collect() })
 }
 
+/// Edge `ei`'s nodes with each of its segments `segs` (counted from its start) halved at the
+/// middle of its curve parameters, the new node on the curve (on a line, halfway between the two,
+/// as [`sample_edge`] spaces a line's nodes). A segment too short to halve in floating point keeps
+/// its two nodes. More than [`MAX_EDGE_NODES`] nodes, or a curve that does not evaluate there, is
+/// an error naming the edge.
+pub(crate) fn split_segments(b: &brep::Brep, ei: usize, nodes: &mut Nodes, segs: &[u32]) -> Result<()> {
+    let fail = |why: String| KernelError::Failed(format!("tessellation: edge {ei}: {why}"));
+    let edge = b.edges.get(ei).ok_or_else(|| fail("no such edge".into()))?;
+    let mut segs = segs.to_vec();
+    segs.sort_unstable();
+    segs.dedup();
+    let n = nodes.points.len().min(nodes.params.len());
+    if n + segs.len() > MAX_EDGE_NODES {
+        return Err(fail(format!("more than {MAX_EDGE_NODES} nodes")));
+    }
+    let curve = match edge.curve {
+        Curve::Line { .. } => None,
+        _ => Some(CurveEval::new(edge.parameters()).map_err(|e| fail(e.to_string()))?),
+    };
+    let (mut points, mut params) = (Vec::with_capacity(n + segs.len()), Vec::with_capacity(n + segs.len()));
+    let mut next = segs.iter().copied().peekable();
+    for k in 0..n {
+        let (Some(&p), Some(&t)) = (nodes.points.get(k), nodes.params.get(k)) else { continue };
+        points.push(p);
+        params.push(t);
+        if next.next_if_eq(&(k as u32)).is_none() {
+            continue;
+        }
+        let (Some(&q), Some(&u)) = (nodes.points.get(k + 1), nodes.params.get(k + 1)) else { continue };
+        let m = 0.5 * (t + u);
+        if !(m > t.min(u) && m < t.max(u)) {
+            continue;
+        }
+        let pm = match &curve {
+            None => (p + q) * 0.5,
+            Some(c) => c.point(m).filter(|p| p.is_finite()).ok_or_else(|| fail(format!("the curve does not evaluate at t = {m}")))?,
+        };
+        points.push(pm);
+        params.push(m);
+    }
+    *nodes = Nodes { points, params };
+    Ok(())
+}
+
 /// Angle (radians) between two directions; zero when either is empty.
 fn turn(a: Vec3, b: Vec3) -> f64 {
     let (la, lb) = (a.len(), b.len());

crates/kernel/src/tess/mod.rs

@@ -30,6 +30,7 @@ pub(crate) mod sag;
 pub(crate) mod seam;
 pub(crate) mod size;
 pub(crate) mod surface;
+pub(crate) mod tangle;
 pub(crate) mod uv;
 
 use std::collections::HashMap;
@@ -64,6 +65,9 @@ const PARALLEL_EDGES: usize = 256;
 const MAX_LINES: usize = 8192;
 /// Most points Ruppert refinement may add to one face, beyond its seeds.
 const MAX_REFINE: usize = 20_000;
+/// Most rounds of splitting the edge segments where a face's sampled boundary crosses itself:
+/// each round halves them, so their chords sag a quarter as far.
+const MAX_UNTANGLE: u32 = 16;
 
 fn unsupported(what: impl std::fmt::Display) -> KernelError {
     KernelError::Failed(format!("not supported yet: tessellation: {what}"))
@@ -154,13 +158,21 @@ pub(crate) fn plan(b: &brep::Brep, tol: f64, turn: Option<size::Turn>) -> Result
     .collect::<Result<_>>()?;
     relax_crowded_faces(b, &jobs, &rules, &mut nodes)?;
     let indexed: Vec<usize> = (0..faces.len()).collect();
-    let planned = crate::par::map("face tessellation", &indexed, curved >= PARALLEL_FACES, |&fi| {
-        let (Some(face), Some(eval), Some(map), Some(size)) = (faces.get(fi), evals.get(fi), params.get(fi), sizes.get(fi)) else {
+    let face_boundary = |fi: usize, nodes: &[edges::Nodes]| {
+        let (Some(face), Some(eval), Some(map)) = (faces.get(fi), evals.get(fi), params.get(fi)) else {
+            return Err(unsupported(format!("face {fi}: missing")));
+        };
+        boundary(b, fi, face, eval, map, nodes, &rules)
+    };
+    let mut boundaries: Vec<Boundary> =
+        crate::par::map("face boundaries", &indexed, curved >= PARALLEL_FACES, |&fi| face_boundary(fi, &nodes)).into_iter().collect::<Result<_>>()?;
+    untangle(b, &evals, &edge_faces, &mut nodes, &mut boundaries, &face_boundary)?;
+    let planned = crate::par::map("face tessellation", &indexed, curved >= PARALLEL_FACES, |&fi| {
+        let (Some(face), Some(eval), Some(size), Some(fb)) = (faces.get(fi), evals.get(fi), sizes.get(fi), boundaries.get(fi)) else {
             return Err(unsupported(format!("face {fi}: missing")));
         };
-        let loops = boundary(b, fi, face, eval, map, &nodes, &rules)?;
         let mut stats = PlanStats::default();
-        let fp = plan_face(fi, kind(face.parameters()), eval, &loops, &rules, *size, &mut stats)?;
+        let fp = plan_face(fi, kind(face.parameters()), eval, &fb.loops, &rules, *size, &mut stats)?;
         let fwd = forward.get(fi).copied().unwrap_or(true);
         let fp =
             FacePlan { flip: !fwd, triangles: if fwd { fp.triangles } else { fp.triangles.into_iter().map(|[a, b, c]| [a, c, b]).collect() }, ..fp };
@@ -204,13 +216,14 @@ fn boundary(
     map: &surface::Params,
     nodes: &[edges::Nodes],
     rules: &size::Rules,
-) -> Result<Vec<seam::Loop>> {
+) -> Result<Boundary> {
     let dom = uv::domain(s);
     let fail = |e: u32, why: String| unsupported(format!("face {fi} ({}): edge {e}: {why}", kind(&face.surface)));
-    let mut loops = Vec::with_capacity(face.loops.len());
+    let mut out = Boundary { loops: Vec::with_capacity(face.loops.len()), segments: Vec::with_capacity(face.loops.len()) };
     for l in &face.loops {
         let Loop::Edges { edges } = l else { continue };
         let (mut pts, mut uvs): (Vec<Vec3>, Vec<[f64; 2]>) = (Vec::new(), Vec::new());
+        let mut segs: Vec<(u32, u32)> = Vec::new();
         // Where the previous edge ended, and the periods its pcurve was moved by.
         let (mut end, mut shift): (Option<[f64; 2]>, [f64; 2]) = (None, [0.0; 2]);
         for oe in edges {
@@ -255,6 +268,9 @@ fn boundary(
             end = eq.last().copied();
             ep.pop();
             eq.pop();
+            // The segment from each node on: the edge's own, counted from its start.
+            let m = ep.len() as u32;
+            segs.extend((0..m).map(|j| (oe.edge, if oe.forward { j } else { m - 1 - j })));
             pts.extend(ep);
             uvs.extend(eq);
         }
@@ -263,12 +279,88 @@ fn boundary(
         if let Some(k) = (0..n).find(|&k| uvs.get(k).is_some_and(|q| !pole(q)) && uvs.get((k + n - 1) % n).is_some_and(|q| !pole(q))) {
             uvs.rotate_left(k);
             pts.rotate_left(k);
+            segs.rotate_left(k);
         }
         if pts.len() >= 2 {
-            loops.push(seam::Loop { uv: uvs, pts });
+            out.loops.push(seam::Loop { uv: uvs, pts });
+            out.segments.push(segs);
         }
     }
-    Ok(loops)
+    Ok(out)
+}
+
+/// A face's boundary loops, and for each loop node the edge segment that starts there: the
+/// edge, and the segment's place along it from its start vertex.
+struct Boundary {
+    loops: Vec<seam::Loop>,
+    segments: Vec<Vec<(u32, u32)>>,
+}
+
+/// Edges are sampled to the chord tolerance each on its own, so where two curves of one face
+/// run closer than their chords sag (the two arcs of a thin lune, a sliver cap between two
+/// cylinders), the chords may cross although the curves do not, and no triangulation follows
+/// such a boundary. The edge segments whose chords cross another of the face (in its loop or
+/// another, see [`tangle::crossing`]) are halved at their curve parameter's middle, on the edge,
+/// so every face of the edge gets the new nodes; the faces of every edge split are checked again,
+/// until no chords cross. Boundaries that still cross after [`MAX_UNTANGLE`] rounds, or would
+/// put more than [`edges::MAX_EDGE_NODES`] on an edge, are an error naming the face and edges.
+fn untangle(
+    b: &brep::Brep,
+    evals: &[EvalSurface],
+    edge_faces: &[Vec<usize>],
+    nodes: &mut [edges::Nodes],
+    boundaries: &mut [Boundary],
+    face_boundary: &(dyn Fn(usize, &[edges::Nodes]) -> Result<Boundary> + Sync),
+) -> Result<()> {
+    let crossing = |fi: usize, fb: &Boundary| -> Vec<(u32, u32)> {
+        let Some(s) = evals.get(fi) else { return Vec::new() };
+        tangle::crossing(s, &fb.loops).into_iter().filter_map(|(l, k)| fb.segments.get(l)?.get(k).copied()).collect()
+    };
+    let indexed: Vec<usize> = (0..boundaries.len()).collect();
+    let curved = evals.iter().filter(|e| !matches!(e, EvalSurface::Plane { .. })).count();
+    let found: Vec<Vec<(u32, u32)>> = crate::par::map("boundary crossings", &indexed, curved >= PARALLEL_FACES, |&fi| {
+        Ok(boundaries.get(fi).map_or_else(Vec::new, |fb| crossing(fi, fb)))
+    })
+    .into_iter()
+    .collect::<Result<_>>()?;
+    let mut tangled: Vec<(usize, Vec<(u32, u32)>)> = found.into_iter().enumerate().filter(|(_, c)| !c.is_empty()).collect();
+    for round in 0..=MAX_UNTANGLE {
+        let Some((fi, first)) = tangled.first() else { return Ok(()) };
+        if round == MAX_UNTANGLE {
+            let mut es: Vec<u32> = first.iter().map(|(e, _)| *e).collect();
+            es.sort_unstable();
+            es.dedup();
+            let kind = b.faces.get(*fi).map_or("?", |f| kind(&f.surface));
+            return Err(unsupported(format!(
+                "face {fi} ({kind}): boundary crosses itself: edges {es:?} run closer together than their chords can follow after {MAX_UNTANGLE} halvings"
+            )));
+        }
+        let mut split: Vec<(u32, u32)> = tangled.iter().flat_map(|(_, c)| c.iter().copied()).collect();
+        split.sort_unstable();
+        split.dedup();
+        let mut redo: Vec<usize> = Vec::new();
+        for run in split.chunk_by(|a, c| a.0 == c.0) {
+            let Some(&(e, _)) = run.first() else { continue };
+            let segs: Vec<u32> = run.iter().map(|(_, k)| *k).collect();
+            let Some(slot) = nodes.get_mut(e as usize) else { continue };
+            edges::split_segments(b, e as usize, slot, &segs)?;
+            redo.extend(edge_faces.get(e as usize).into_iter().flatten().copied());
+        }
+        redo.sort_unstable();
+        redo.dedup();
+        tangled.clear();
+        for fi in redo {
+            let fb = face_boundary(fi, nodes)?;
+            let c = crossing(fi, &fb);
+            if let Some(slot) = boundaries.get_mut(fi) {
+                *slot = fb;
+            }
+            if !c.is_empty() {
+                tangled.push((fi, c));
+            }
+        }
+    }
+    Ok(())
 }
 
 /// Bounded work for the facet angle: a face whose boundary would get more than

crates/kernel/src/tess/tangle.rs

@@ -0,0 +1,135 @@
+//! Where a face's sampled boundary crosses itself: chords of its loops, in the face's
+//! parameters, that cross another chord of the face although the curves they stand for do not
+//! (two curves running closer together than their chords sag).
+
+use solidcraft_geom::tess::EvalSurface;
+
+use super::seam;
+use super::uv;
+
+/// Two chords cross only where each one's ends lie strictly either side of the other, by more
+/// than this fraction of the product of their lengths (a chord ending on another, as at a
+/// shared vertex or along a seam's two uses a period apart, is not a crossing).
+const STRADDLE: f64 = 1e-10;
+
+/// One chord of a loop, or a copy of it whole periods along.
+struct Chord {
+    /// The chord's place in the face: its loop and its first node.
+    at: (usize, usize),
+    a: [f64; 2],
+    b: [f64; 2],
+}
+
+impl Chord {
+    fn lo(&self, d: usize) -> f64 {
+        self.a[d].min(self.b[d])
+    }
+    fn hi(&self, d: usize) -> f64 {
+        self.a[d].max(self.b[d])
+    }
+}
+
+/// The chords (loop, first node) of `loops` that cross another chord of the face, in its loop or
+/// another, each once and in order. On a closed surface the boundary crosses itself where any
+/// copies of two chords, whole periods apart, cross; chords ending on a pole, where the
+/// parameter along the pole is free, are left out. Adjacent chords of a loop share a node and
+/// never cross.
+pub(crate) fn crossing(s: &EvalSurface, loops: &[seam::Loop]) -> Vec<(usize, usize)> {
+    let dom = uv::domain(s);
+    let pole = |q: [f64; 2]| (0..2).any(|d| uv::degenerate(s, q, d));
+    let mut chords: Vec<Chord> = Vec::new();
+    for (li, l) in loops.iter().enumerate() {
+        let n = l.uv.len();
+        if n < 3 {
+            continue;
+        }
+        for k in 0..n {
+            let (Some(&a), Some(&b)) = (l.uv.get(k), l.uv.get((k + 1) % n)) else { continue };
+            if !(a.iter().chain(&b).all(|x| x.is_finite())) || pole(a) || pole(b) {
+                continue;
+            }
+            // The loop's closing chord of a loop that winds round the surface ends at the copy of
+            // its first node nearest its last; every chord then starts in one period.
+            let (mut a, mut b) = (a, b);
+            for d in 0..2 {
+                let Some(per) = dom.period[d] else { continue };
+                b[d] = uv::unwrap(b[d], a[d], Some(per));
+                let back = (a[d].min(b[d]) / per).floor() * per;
+                a[d] -= back;
+                b[d] -= back;
+            }
+            chords.push(Chord { at: (li, k), a, b });
+        }
+    }
+    // With every chord starting in one period, two chords' copies can only cross where one of
+    // them is moved a period along at most: each chord also as its copies a period along.
+    for d in 0..2 {
+        let Some(per) = dom.period[d] else { continue };
+        let along: Vec<Chord> = chords
+            .iter()
+            .map(|c| {
+                let (mut a, mut b) = (c.a, c.b);
+                a[d] += per;
+                b[d] += per;
+                Chord { at: c.at, a, b }
+            })
+            .collect();
+        chords.extend(along);
+    }
+    let adjacent = |p: (usize, usize), q: (usize, usize)| {
+        let n = loops.get(p.0).map_or(0, |l| l.uv.len());
+        p == q || (p.0 == q.0 && ((p.1 + 1) % n == q.1 || (q.1 + 1) % n == p.1))
+    };
+    // Sweep along u: only chords whose u ranges overlap can cross.
+    chords.sort_by(|x, y| x.lo(0).total_cmp(&y.lo(0)));
+    let mut out: Vec<(usize, usize)> = Vec::new();
+    for (i, x) in chords.iter().enumerate() {
+        for y in chords.iter().skip(i + 1) {
+            if y.lo(0) > x.hi(0) {
+                break;
+            }
+            if y.lo(1) > x.hi(1) || x.lo(1) > y.hi(1) || adjacent(x.at, y.at) {
+                continue;
+            }
+            if cross(x.a, x.b, y.a, y.b) {
+                out.push(x.at);
+                out.push(y.at);
+            }
+        }
+    }
+    out.sort_unstable();
+    out.dedup();
+    out
+}
+
+/// Whether the segments `ab` and `cd` cross: each one's ends strictly either side of the other.
+fn cross(a: [f64; 2], b: [f64; 2], c: [f64; 2], d: [f64; 2]) -> bool {
+    let orient = |p: [f64; 2], q: [f64; 2], r: [f64; 2]| (q[0] - p[0]) * (r[1] - p[1]) - (q[1] - p[1]) * (r[0] - p[0]);
+    let len = |p: [f64; 2], q: [f64; 2]| (q[0] - p[0]).hypot(q[1] - p[1]);
+    let eps = STRADDLE * len(a, b) * len(c, d);
+    let side = |o: f64| {
+        if o > eps {
+            1
+        } else if o < -eps {
+            -1
+        } else {
+            0
+        }
+    };
+    let (s1, s2, s3, s4) = (side(orient(a, b, c)), side(orient(a, b, d)), side(orient(c, d, a)), side(orient(c, d, b)));
+    s1 * s2 < 0 && s3 * s4 < 0
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn crossing_chords_are_found_and_touching_ones_are_not() {
+        assert!(cross([0.0, 0.0], [2.0, 2.0], [0.0, 2.0], [2.0, 0.0]));
+        // Ending on the other, sharing an end, collinear: no crossing.
+        assert!(!cross([0.0, 0.0], [2.0, 0.0], [1.0, 0.0], [1.0, 1.0]));
+        assert!(!cross([0.0, 0.0], [2.0, 0.0], [2.0, 0.0], [3.0, 1.0]));
+        assert!(!cross([0.0, 0.0], [2.0, 0.0], [1.0, 0.0], [3.0, 0.0]));
+    }
+}

crates/kernel/src/tess/tests.rs

@@ -636,3 +636,37 @@ fn loops_touching_at_a_shared_vertex_bound_the_face() {
     let Some((_, triangles)) = cdt::triangulate(&uvs, &ranges, fill) else { panic!("triangulates") };
     assert!(closed_boundary(&ranges, &fixed, &uvs, &triangles), "the triangles follow both holes through their shared vertex");
 }
+
+/// A face whose loop is a thin lune meshes: a half disc of radius 2.5 less a circle through the
+/// same ends 0.03 mm inside it at the middle, its inner side split at 11.25°. Each edge is
+/// sampled to the 0.1 mm chord on its own: the outer arc's chords (every 22.5°) sag 0.048 mm,
+/// more than the curves are apart, so the one from 0° to 22.5° passes inside the inner side's
+/// node at 11.25° and crosses its chords; no triangulation follows a boundary that crosses
+/// itself (the sliver caps of a joined ring of cylinder segments). The crossing chords are
+/// halved on their edges until none cross.
+#[test]
+fn a_thin_lune_whose_chords_cross_meshes() {
+    use solidcraft_geom::{Loop2, Plane, Region2, Seg2, Vec2};
+    let (r, gap, h, tol) = (2.5, 0.03, 5.0, 0.1);
+    // The inner circle: centre (−d, 0) and radius r + d·…, through (0, ±r), its middle at r − gap.
+    let d = (r * r - (r - gap) * (r - gap)) / (2.0 * (r - gap));
+    let (centre, r2) = (Vec2::new(-d, 0.0), r - gap + d);
+    let top = r.atan2(d);
+    let split = 11.25f64.to_radians();
+    let segs = vec![
+        Seg2::Arc { center: Vec2::ZERO, radius: r, start: -0.5 * PI, sweep: PI },
+        Seg2::Arc { center: centre, radius: r2, start: top, sweep: split - top },
+        Seg2::Arc { center: centre, radius: r2, start: split, sweep: -top - split },
+    ];
+    let lune = Region2 { outer: Loop2 { segs }, holes: vec![] };
+    let Some(body) = ok(crate::extrude(&Plane::XY, &[lune], 0.0, h)).pop() else { panic!("one body") };
+    let rec = record(&body);
+    let (p, stats) = match plan(&rec, tol, None) {
+        Ok(r) => r,
+        Err(e) => panic!("lune: {e}"),
+    };
+    assert_eq!(stats.planar + stats.structured + stats.delaunay, body.face_count(), "every face planned");
+    let m = p.evaluate();
+    assert_watertight(&m, "lune");
+    assert!(m.measure().volume > 0.0, "inside out");
+}