Constrained Delaunay tetrahedralization (boundary recovery) and mesh.tet.generate (M2) (#29)

* predicates: LNC points, generic insphere and its perturbed form (Diazzi et al. 2023 §4.3)

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

* mesh: constrained Delaunay tetrahedralization (segment recovery with LNC Steiner points, gift-wrapping facet recovery, labelling)

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

* cdt: rounding repair, mesh.tet.generate, corpus CDT column; review cleanups (shared perturbation rule, one facet scan, cached coordinates, worklist segment recovery)

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

* mesh.tet.generate is ready in the catalogue; corpus table with the CDT column

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

---------

Co-authored-by: jyviko <2501648+jyviko@users.noreply.github.com>
Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

crates/mesh/src/cdt/complex.rs

@@ -0,0 +1,590 @@
+//! The working tetrahedralization of the boundary recovery: tetrahedra over explicit and implicit
+//! ([`Lnc`](crate::predicates::Lnc)) vertices, with the operations the recovery needs: point
+//! location, Bowyer–Watson insertion of Steiner points, vertex stars, edge lookup, and the
+//! replacement of a set of tetrahedra by another filling the same region.
+//!
+//! Layout: per tetrahedron four vertex ids in tet4 order (`det[b − a, c − a, d − a] > 0`, decided
+//! exactly) and four neighbour ids, neighbour `i` across the face opposite vertex `i`; [`NONE`] on
+//! the hull. The hull is the bounding box the recovery adds around the surface, so every Steiner
+//! point is strictly inside it and no ghost tetrahedra are needed.
+
+use std::collections::HashMap;
+use std::hash::{BuildHasherDefault, Hasher};
+
+use crate::delaunay::Triangulation;
+use crate::predicates::{GenericPoint, Sign, generic_insphere_sos, generic_lex, generic_orient3d};
+use crate::{Error, Result};
+
+/// A multiplicative hasher for small integer keys (vertex triples). The maps it serves are only
+/// looked up, never iterated, so its order cannot leak into the output.
+#[derive(Default)]
+pub(super) struct IdHasher(u64);
+
+impl Hasher for IdHasher {
+    fn finish(&self) -> u64 {
+        self.0
+    }
+
+    fn write(&mut self, bytes: &[u8]) {
+        for &b in bytes {
+            self.write_u64(u64::from(b));
+        }
+    }
+
+    fn write_u32(&mut self, x: u32) {
+        self.write_u64(u64::from(x));
+    }
+
+    fn write_u64(&mut self, x: u64) {
+        self.0 = (self.0.rotate_left(5) ^ x).wrapping_mul(0x51_7c_c1_b7_27_22_0a_95);
+    }
+}
+
+/// A hash map keyed by vertex ids.
+pub(super) type IdMap<K, V> = HashMap<K, V, BuildHasherDefault<IdHasher>>;
+
+/// No tetrahedron / no vertex.
+pub(super) const NONE: u32 = u32::MAX;
+/// Marks a recycled tetrahedron slot.
+const DEAD: u32 = u32::MAX - 1;
+
+pub(super) fn corrupt(what: &str) -> Error {
+    Error::Failed(format!("constrained tetrahedralization is inconsistent: {what}"))
+}
+
+/// The face of `vs` opposite slot `k`, ordered so that vertex `k` is on its positive side
+/// ([`Complex::volume`] of the face plus that vertex is positive).
+pub(super) fn face_toward(vs: [u32; 4], k: usize) -> [u32; 3] {
+    let [a, b, c, d] = vs;
+    match k {
+        0 => [b, d, c],
+        1 => [a, c, d],
+        2 => [a, d, b],
+        _ => [a, b, c],
+    }
+}
+
+/// A triangle's vertices sorted: the key of an unoriented face.
+pub(super) fn sorted3([a, b, c]: [u32; 3]) -> [u32; 3] {
+    let mut f = [a, b, c];
+    f.sort_unstable();
+    f
+}
+
+/// The working tetrahedralization.
+#[derive(Debug, Clone)]
+pub(super) struct Complex {
+    points: Vec<GenericPoint>,
+    /// A box certainly containing each point (implicit points are not exactly representable).
+    lo: Vec<[f64; 3]>,
+    hi: Vec<[f64; 3]>,
+    /// Rounded coordinates (floating-point estimates only, never decisions).
+    approx: Vec<[f64; 3]>,
+    tets: Vec<[u32; 4]>,
+    adj: Vec<[u32; 4]>,
+    free: Vec<u32>,
+    vertex_tet: Vec<u32>,
+    /// Traversal marks: a tetrahedron is visited in the current traversal when its mark equals
+    /// `epoch`.
+    mark: Vec<u32>,
+    epoch: u32,
+    rng: u64,
+    /// Most tetrahedron slots allowed.
+    max_tets: usize,
+    live: usize,
+}
+
+/// Where a point lies ([`Complex::locate`]).
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub(super) enum Located {
+    /// In the closed tetrahedron and not one of its vertices.
+    Inside(u32),
+    /// Equal to this vertex.
+    Vertex(u32),
+}
+
+impl Complex {
+    /// The finite tetrahedra of a Delaunay triangulation of `points` (in the same vertex order),
+    /// whose hull faces get no neighbour.
+    pub(super) fn from_delaunay(tri: &Triangulation, points: Vec<GenericPoint>, max_tets: usize) -> Result<Self> {
+        let mut id = vec![NONE; tri.slots()];
+        let mut tets = Vec::new();
+        for (t, vs) in tri.finite_tets() {
+            if let Some(slot) = id.get_mut(t as usize) {
+                *slot = tets.len() as u32;
+            }
+            tets.push(vs);
+        }
+        if tets.len() > max_tets {
+            return Err(Error::Limit(format!("{} tetrahedra; the limit is {max_tets}", tets.len())));
+        }
+        let mut adj = Vec::with_capacity(tets.len());
+        for (t, _) in tri.finite_tets() {
+            let ns = tri.neighbors(t).ok_or_else(|| corrupt("Delaunay neighbours"))?;
+            adj.push(ns.map(|n| id.get(n as usize).copied().unwrap_or(NONE)));
+        }
+        let (lo, hi) = points.iter().map(GenericPoint::bounds).unzip();
+        let approx = points.iter().map(GenericPoint::approx).collect();
+        let mut cx = Complex {
+            vertex_tet: vec![NONE; points.len()],
+            points,
+            lo,
+            hi,
+            approx,
+            mark: vec![0; tets.len()],
+            live: tets.len(),
+            tets,
+            adj,
+            free: Vec::new(),
+            epoch: 0,
+            rng: 0x9E37_79B9_7F4A_7C15,
+            max_tets,
+        };
+        for (t, vs) in cx.tets.iter().enumerate() {
+            for &v in vs {
+                if let Some(s) = cx.vertex_tet.get_mut(v as usize) {
+                    *s = t as u32;
+                }
+            }
+        }
+        Ok(cx)
+    }
+
+    // ---- read access ----
+
+    pub(super) fn points(&self) -> &[GenericPoint] {
+        &self.points
+    }
+
+    pub(super) fn pt(&self, v: u32) -> Result<&GenericPoint> {
+        self.points.get(v as usize).ok_or_else(|| corrupt("vertex id"))
+    }
+
+    /// Vertex `v`'s coordinates rounded to doubles: for estimates and messages only.
+    pub(super) fn approx(&self, v: u32) -> Result<[f64; 3]> {
+        self.approx.get(v as usize).copied().ok_or_else(|| corrupt("vertex id"))
+    }
+
+    /// The box certainly containing vertex `v`.
+    pub(super) fn bounds(&self, v: u32) -> ([f64; 3], [f64; 3]) {
+        match (self.lo.get(v as usize), self.hi.get(v as usize)) {
+            (Some(&l), Some(&h)) => (l, h),
+            _ => ([f64::NEG_INFINITY; 3], [f64::INFINITY; 3]),
+        }
+    }
+
+    /// The box certainly containing all of `vs`.
+    pub(super) fn bounds_of(&self, vs: &[u32]) -> ([f64; 3], [f64; 3]) {
+        let mut lo = [f64::INFINITY; 3];
+        let mut hi = [f64::NEG_INFINITY; 3];
+        for &v in vs {
+            let (l, h) = self.bounds(v);
+            for k in 0..3 {
+                lo[k] = lo[k].min(l[k]);
+                hi[k] = hi[k].max(h[k]);
+            }
+        }
+        (lo, hi)
+    }
+
+    pub(super) fn slots(&self) -> usize {
+        self.tets.len()
+    }
+
+    pub(super) fn live_count(&self) -> usize {
+        self.live
+    }
+
+    /// The vertices of a live tetrahedron.
+    pub(super) fn tet(&self, t: u32) -> Result<[u32; 4]> {
+        self.tets.get(t as usize).copied().filter(|v| v[0] != DEAD).ok_or_else(|| corrupt("tetrahedron id"))
+    }
+
+    pub(super) fn is_live(&self, t: u32) -> bool {
+        self.tets.get(t as usize).is_some_and(|v| v[0] != DEAD)
+    }
+
+    pub(super) fn adj(&self, t: u32) -> Result<[u32; 4]> {
+        self.tet(t)?;
+        self.adj.get(t as usize).copied().ok_or_else(|| corrupt("tetrahedron id"))
+    }
+
+    /// The live tetrahedra in slot order.
+    pub(super) fn live_tets(&self) -> impl Iterator<Item = (u32, [u32; 4])> + '_ {
+        self.tets.iter().enumerate().filter(|(_, v)| v[0] != DEAD).map(|(t, v)| (t as u32, *v))
+    }
+
+    /// The sign of `det[b − a, c − a, d − a]` (tet4 orientation), exact.
+    pub(super) fn volume(&self, [a, b, c, d]: [u32; 4]) -> Result<Sign> {
+        Ok(-generic_orient3d(self.pt(a)?, self.pt(b)?, self.pt(c)?, self.pt(d)?))
+    }
+
+    /// Whether `x` is inside the perturbed circumsphere of the positive tetrahedron `vs`.
+    pub(super) fn in_sphere(&self, [a, b, c, d]: [u32; 4], x: u32) -> Result<bool> {
+        // `insphere` wants Shewchuk-positive orientation (tet4-negative): swap the first two.
+        Ok(generic_insphere_sos(self.pt(b)?, self.pt(a)?, self.pt(c)?, self.pt(d)?, self.pt(x)?) == Sign::Positive)
+    }
+
+    /// Whether two vertices are the same point (exact).
+    fn same_point(&self, u: u32, v: u32) -> Result<bool> {
+        Ok(generic_lex(self.pt(u)?, self.pt(v)?).is_eq())
+    }
+
+    // ---- traversal marks ----
+
+    /// A new traversal: marks `epoch` (visited) and `epoch + 1` (rejected) are fresh.
+    fn next_epoch(&mut self) {
+        if self.epoch >= u32::MAX - 4 {
+            self.mark.fill(0);
+            self.epoch = 0;
+        }
+        self.epoch += 2;
+    }
+
+    fn visited(&self, t: u32) -> bool {
+        self.mark.get(t as usize) == Some(&self.epoch)
+    }
+
+    fn visit(&mut self, t: u32) {
+        if let Some(m) = self.mark.get_mut(t as usize) {
+            *m = self.epoch;
+        }
+    }
+
+    fn rejected(&self, t: u32) -> bool {
+        self.mark.get(t as usize) == Some(&(self.epoch + 1))
+    }
+
+    fn reject(&mut self, t: u32) {
+        if let Some(m) = self.mark.get_mut(t as usize) {
+            *m = self.epoch + 1;
+        }
+    }
+
+    // ---- stars and edges ----
+
+    /// A live tetrahedron incident to `v`.
+    fn incident(&self, v: u32) -> Result<u32> {
+        let t = self.vertex_tet.get(v as usize).copied().unwrap_or(NONE);
+        if self.tets.get(t as usize).is_some_and(|vs| vs[0] != DEAD && vs.contains(&v)) {
+            return Ok(t);
+        }
+        // The hint is maintained on every allocation, so this scan is a safety net only.
+        self.live_tets().find(|(_, vs)| vs.contains(&v)).map(|(t, _)| t).ok_or_else(|| corrupt("a vertex is in no tetrahedron"))
+    }
+
+    /// The tetrahedra incident to `v` (its star), into `out`.
+    pub(super) fn star(&mut self, v: u32, out: &mut Vec<u32>) -> Result<()> {
+        out.clear();
+        let start = self.incident(v)?;
+        self.next_epoch();
+        self.visit(start);
+        out.push(start);
+        let mut i = 0;
+        while let Some(&t) = out.get(i) {
+            i += 1;
+            let vs = self.tet(t)?;
+            let ns = self.adj(t)?;
+            for (&u, &n) in vs.iter().zip(&ns) {
+                if u == v || n == NONE || self.visited(n) {
+                    continue;
+                }
+                self.visit(n);
+                out.push(n);
+            }
+            if out.len() > self.live {
+                return Err(corrupt("a vertex star larger than the mesh"));
+            }
+        }
+        Ok(())
+    }
+
+    /// Whether `u v` is an edge.
+    pub(super) fn has_edge(&mut self, u: u32, v: u32, scratch: &mut Vec<u32>) -> Result<bool> {
+        self.star(u, scratch)?;
+        for &t in scratch.iter() {
+            if self.tet(t)?.contains(&v) {
+                return Ok(true);
+            }
+        }
+        Ok(false)
+    }
+
+    /// The tetrahedra whose (conservative) box meets the box `lo..hi`, connected to the star of
+    /// `v` through such tetrahedra; at most `cap` of them.
+    pub(super) fn tets_near(&mut self, v: u32, lo: [f64; 3], hi: [f64; 3], cap: usize) -> Result<Vec<u32>> {
+        let mut out = Vec::new();
+        self.star(v, &mut out)?;
+        self.next_epoch();
+        for &t in &out {
+            self.visit(t);
+        }
+        let mut i = 0;
+        while let Some(&t) = out.get(i) {
+            i += 1;
+            if out.len() >= cap {
+                break;
+            }
+            for n in self.adj(t)? {
+                if n == NONE || self.visited(n) {
+                    continue;
+                }
+                self.visit(n);
+                let (l, h) = self.bounds_of(&self.tet(n)?);
+                if (0..3).all(|k| l[k] <= hi[k] && h[k] >= lo[k]) {
+                    out.push(n);
+                }
+            }
+        }
+        Ok(out)
+    }
+
+    // ---- point location ----
+
+    /// The visibility walk towards `p` from `start` (Devillers, Pion, Teillaud 2002), capped at
+    /// the number of tetrahedra, then a scan.
+    pub(super) fn locate(&mut self, p: u32, start: u32) -> Result<Located> {
+        let mut t = if self.is_live(start) { start } else { self.live_tets().next().map(|(t, _)| t).ok_or_else(|| corrupt("empty mesh"))? };
+        let mut prev = NONE;
+        for _ in 0..self.live.saturating_add(16) {
+            let vs = self.tet(t)?;
+            let ns = self.adj(t)?;
+            self.rng ^= self.rng << 13;
+            self.rng ^= self.rng >> 7;
+            self.rng ^= self.rng << 17;
+            let first = (self.rng >> 61) as usize;
+            let mut next = None;
+            for i in 0..4 {
+                let k = (first + i) % 4;
+                let n = ns.get(k).copied().unwrap_or(NONE);
+                if (n != prev || n == NONE) && self.beyond(vs, k, p)? {
+                    next = Some(n);
+                    break;
+                }
+            }
+            match next {
+                None => return self.inside(vs, t, p),
+                Some(NONE) => return Err(Error::Failed("a Steiner point lies outside the bounding box".into())),
+                Some(n) => {
+                    prev = t;
+                    t = n;
+                }
+            }
+        }
+        // The walk may cycle in a non-Delaunay mesh: scan.
+        let live: Vec<(u32, [u32; 4])> = self.live_tets().collect();
+        for (t, vs) in live {
+            let mut inside = true;
+            for k in 0..4 {
+                if self.beyond(vs, k, p)? {
+                    inside = false;
+                    break;
+                }
+            }
+            if inside {
+                return self.inside(vs, t, p);
+            }
+        }
+        Err(Error::Failed("a Steiner point is in no tetrahedron".into()))
+    }
+
+    /// Whether `p` is strictly beyond the face of `vs` opposite slot `k`.
+    fn beyond(&self, mut vs: [u32; 4], k: usize, p: u32) -> Result<bool> {
+        if let Some(s) = vs.get_mut(k) {
+            *s = p;
+        }
+        Ok(self.volume(vs)? == Sign::Negative)
+    }
+
+    fn inside(&self, vs: [u32; 4], t: u32, p: u32) -> Result<Located> {
+        for v in vs {
+            if self.same_point(v, p)? {
+                return Ok(Located::Vertex(v));
+            }
+        }
+        Ok(Located::Inside(t))
+    }
+
+    // ---- modification ----
+
+    /// Add a point (not yet in any tetrahedron); its vertex id.
+    pub(super) fn push_point(&mut self, p: GenericPoint) -> Result<u32> {
+        if self.points.len() >= (NONE - 2) as usize {
+            return Err(Error::Limit("too many vertices".into()));
+        }
+        let (l, h) = p.bounds();
+        self.approx.push(p.approx());
+        self.points.push(p);
+        self.lo.push(l);
+        self.hi.push(h);
+        self.vertex_tet.push(NONE);
+        Ok((self.points.len() - 1) as u32)
+    }
+
+    /// Insert vertex `v` (from [`Complex::push_point`]) by Bowyer–Watson with the perturbed
+    /// in-sphere test: the mesh must be Delaunay for that predicate, which it is until the facet
+    /// recovery starts. The walk starts at a tetrahedron incident to `near`. Returns the vertices
+    /// of the deleted tetrahedra (the only ones whose edges can have changed).
+    pub(super) fn insert(&mut self, v: u32, near: u32) -> Result<Vec<u32>> {
+        let start = self.incident(near)?;
+        let seed = match self.locate(v, start)? {
+            Located::Vertex(u) => {
+                let at = self.approx(u)?;
+                return Err(Error::Failed(format!("a Steiner point coincides with vertex {u} at {at:?}")));
+            }
+            Located::Inside(t) => t,
+        };
+        // The cavity: every tetrahedron whose perturbed circumsphere contains v, grown from the
+        // one containing it (which always conflicts: v is in the closed tetrahedron, so strictly
+        // inside its ball).
+        self.next_epoch();
+        let mut cavity = vec![seed];
+        self.visit(seed);
+        let mut boundary: Vec<(u32, usize)> = Vec::new();
+        let mut i = 0;
+        while let Some(&t) = cavity.get(i) {
+            i += 1;
+            for (j, n) in self.adj(t)?.into_iter().enumerate() {
+                if n == NONE || self.rejected(n) {
+                    boundary.push((t, j));
+                    continue;
+                }
+                if self.visited(n) {
+                    continue;
+                }
+                if self.in_sphere(self.tet(n)?, v)? {
+                    self.visit(n);
+                    cavity.push(n);
+                } else {
+                    self.reject(n);
+                    boundary.push((t, j));
+                }
+            }
+            if cavity.len() > self.live {
+                return Err(corrupt("a cavity larger than the mesh"));
+            }
+        }
+        let mut new = Vec::with_capacity(boundary.len());
+        for (t, j) in boundary {
+            let mut vs = self.tet(t)?;
+            if let Some(s) = vs.get_mut(j) {
+                *s = v;
+            }
+            if self.volume(vs)? != Sign::Positive {
+                return Err(corrupt("a Bowyer–Watson cavity is not star-shaped from the new point"));
+            }
+            new.push(vs);
+        }
+        let mut touched: Vec<u32> = Vec::with_capacity(4 * cavity.len());
+        for &t in &cavity {
+            touched.extend(self.tet(t)?);
+        }
+        self.replace(&cavity, &new)?;
+        touched.sort_unstable();
+        touched.dedup();
+        Ok(touched)
+    }
+
+    fn alloc(&mut self, vs: [u32; 4]) -> Result<u32> {
+        let t = match self.free.pop() {
+            Some(t) => {
+                *self.tets.get_mut(t as usize).ok_or_else(|| corrupt("free slot"))? = vs;
+                *self.adj.get_mut(t as usize).ok_or_else(|| corrupt("free slot"))? = [NONE; 4];
+                t
+            }
+            None => {
+                if self.tets.len() >= self.max_tets {
+                    return Err(Error::Limit(format!("more than {} tetrahedra", self.max_tets)));
+                }
+                self.tets.push(vs);
+                self.adj.push([NONE; 4]);
+                self.mark.push(0);
+                (self.tets.len() - 1) as u32
+            }
+        };
+        for v in vs {
+            if let Some(s) = self.vertex_tet.get_mut(v as usize) {
+                *s = t;
+            }
+        }
+        self.live += 1;
+        Ok(t)
+    }
+
+    fn kill(&mut self, t: u32) -> Result<()> {
+        *self.tets.get_mut(t as usize).ok_or_else(|| corrupt("tetrahedron id"))? = [DEAD; 4];
+        self.free.push(t);
+        self.live = self.live.saturating_sub(1);
+        Ok(())
+    }
+
+    fn set_adj(&mut self, t: u32, i: usize, n: u32) -> Result<()> {
+        *self.adj.get_mut(t as usize).and_then(|a| a.get_mut(i)).ok_or_else(|| corrupt("tetrahedron id"))? = n;
+        Ok(())
+    }
+
+    /// Replace the tetrahedra `old` by `new` (positive tet4 tuples) filling the same region: each
+    /// face of a new tetrahedron is either shared with exactly one other new tetrahedron or is a
+    /// face of the region's boundary, where it takes over the outside neighbour. Anything else is
+    /// reported, never patched. The new slot ids.
+    pub(super) fn replace(&mut self, old: &[u32], new: &[[u32; 4]]) -> Result<Vec<u32>> {
+        // Read the outside links before any slot is recycled.
+        // Boundary face → (outside neighbour, its slot back, covered yet).
+        let mut external: IdMap<[u32; 3], (u32, usize, bool)> = IdMap::with_capacity_and_hasher(old.len() * 2, Default::default());
+        self.next_epoch();
+        for &t in old {
+            self.visit(t);
+        }
+        for &t in old {
+            let vs = self.tet(t)?;
+            for (j, n) in self.adj(t)?.into_iter().enumerate() {
+                if n != NONE && self.visited(n) {
+                    continue;
+                }
+                let back = if n == NONE { 0 } else { self.adj(n)?.iter().position(|&m| m == t).ok_or_else(|| corrupt("one-sided adjacency"))? };
+                if external.insert(sorted3(face_toward(vs, j)), (n, back, false)).is_some() {
+                    return Err(corrupt("a region boundary face is listed twice"));
+                }
+            }
+        }
+        for &t in old {
+            self.kill(t)?;
+        }
+        let mut ids = Vec::with_capacity(new.len());
+        let mut open: IdMap<[u32; 3], (u32, usize)> = IdMap::with_capacity_and_hasher(new.len() * 2, Default::default());
+        let mut matched_external = 0usize;
+        for &vs in new {
+            let t = self.alloc(vs)?;
+            ids.push(t);
+            for k in 0..4 {
+                let key = sorted3(face_toward(vs, k));
+                if let Some(entry) = external.get_mut(&key) {
+                    let (n, back, covered) = *entry;
+                    if covered {
+                        return Err(corrupt("a region boundary face is covered twice"));
+                    }
+                    entry.2 = true;
+                    self.set_adj(t, k, n)?;
+                    if n != NONE {
+                        self.set_adj(n, back, t)?;
+                    }
+                    matched_external += 1;
+                } else if let Some((u, l)) = open.remove(&key) {
+                    self.set_adj(t, k, u)?;
+                    self.set_adj(u, l, t)?;
+                } else {
+                    open.insert(key, (t, k));
+                }
+            }
+        }
+        if !open.is_empty() || matched_external != external.len() {
+            return Err(corrupt(&format!(
+                "a retriangulated region does not close: {} unmatched inner faces, {} of {} boundary faces covered",
+                open.len(),
+                matched_external,
+                external.len()
+            )));
+        }
+        Ok(ids)
+    }
+}

crates/mesh/src/cdt/facets.rs

@@ -0,0 +1,414 @@
+//! Facet recovery (Diazzi et al. 2023, §3.4 and §4.4). Every input triangle is a facet; once all
+//! segments are edges, a facet is *missing* when an edge of the tetrahedralization pierces its
+//! interior. The tetrahedra with such edges form face-connected *cavities* (each meets the
+//! facet's plane inside the facet only, because the facet's edges are mesh edges). Each cavity is
+//! deleted and split by the facet's plane into two half-cavities, which are refilled one after
+//! the other by the modified gift-wrapping of §4.4 (after Shewchuk 2000, *Sweep algorithms for
+//! constructing higher-dimensional constrained Delaunay triangulations*): repeatedly take an
+//! open triangle `σ` of the half-cavity's front and join it to the vertex `w` of the half-cavity
+//! such that
+//! - (i) `w` is on the open side of `σ` (positive volume),
+//! - (ii) the tetrahedron meets every cavity, front and constraint triangle properly
+//!   ([`super::meet`]), and contains no other cavity vertex,
+//! - (iii) among those, the perturbed circumsphere of `σ w` contains no other candidate: the
+//!   candidate that is first in the pencil of spheres through `σ`, with the same symbolic
+//!   perturbation as the Delaunay tetrahedralization, so both half-cavities induce the same
+//!   triangulation on the facet.
+//!
+//! The first half-cavity leaves its triangles on the facet's plane open; the second closes them.
+//! Constraint triangles already recovered that lie inside a cavity stay constraints: both of
+//! their sides enter the front.
+//!
+//! Si's cavity expansion with local Delaunay tetrahedralizations (§3.4) is a speed-up the paper
+//! uses before falling back to gift-wrapping; gift-wrapping alone is the always-correct path
+//! (§4.4: it builds the constrained Delaunay tetrahedralization of any half-cavity that admits
+//! one), so it is the one implemented here.
+
+use std::collections::{HashMap, HashSet};
+
+use crate::predicates::{GenericPoint, Sign, generic_orient3d};
+use crate::{Error, Result};
+
+use super::Recovery;
+use super::complex::{NONE, corrupt, face_toward, sorted3};
+use super::meet::{in_closed_tet, meet_properly};
+
+/// Most rounds of finding and refilling the cavities of one facet (one suffices in exact
+/// arithmetic).
+const MAX_ROUNDS_PER_FACET: usize = 16;
+
+/// A box: `(lo, hi)`.
+type Bounds = ([f64; 3], [f64; 3]);
+
+/// The open front of a half-cavity: faces, their boxes, and whether each is still open.
+type Front<'a> = (&'a [[u32; 3]], &'a [Bounds], &'a [bool]);
+
+/// An oriented triangle rotated so its smallest id comes first: the key of an oriented face.
+fn canon([a, b, c]: [u32; 3]) -> [u32; 3] {
+    if a <= b && a <= c {
+        [a, b, c]
+    } else if b <= a && b <= c {
+        [b, c, a]
+    } else {
+        [c, a, b]
+    }
+}
+
+fn boxes_meet((l1, h1): Bounds, (l2, h2): Bounds) -> bool {
+    (0..3).all(|k| l1[k] <= h2[k] && l2[k] <= h1[k])
+}
+
+/// The side of the facet's plane each vertex is on, cached for one facet.
+struct Plane {
+    corners: [GenericPoint; 3],
+    side: HashMap<u32, Sign>,
+}
+
+impl Recovery {
+    fn plane(&self, f: usize) -> Result<Plane> {
+        let tri = *self.triangles.get(f).ok_or_else(|| corrupt("facet id"))?;
+        let corner = |v: u32| self.cx.pt(v).cloned();
+        Ok(Plane { corners: [corner(tri[0])?, corner(tri[1])?, corner(tri[2])?], side: HashMap::new() })
+    }
+
+    fn facet_side(&self, plane: &mut Plane, v: u32) -> Result<Sign> {
+        if let Some(&s) = plane.side.get(&v) {
+            return Ok(s);
+        }
+        let [a, b, c] = &plane.corners;
+        // tet4 volume sign of (a, b, c, v): positive on the side the facet's normal points to.
+        let s = -generic_orient3d(a, b, c, self.cx.pt(v)?);
+        plane.side.insert(v, s);
+        Ok(s)
+    }
+
+    /// Whether edge `p q` crosses the facet's open interior.
+    fn pierces(&self, plane: &mut Plane, p: u32, q: u32) -> Result<bool> {
+        if self.facet_side(plane, p)? * self.facet_side(plane, q)? != Sign::Negative {
+            return Ok(false);
+        }
+        let [a, b, c] = &plane.corners;
+        let (pp, pq) = (self.cx.pt(p)?, self.cx.pt(q)?);
+        let s0 = generic_orient3d(pp, pq, a, b);
+        if s0.is_zero() {
+            return Ok(false);
+        }
+        Ok(generic_orient3d(pp, pq, b, c) == s0 && generic_orient3d(pp, pq, c, a) == s0)
+    }
+
+    fn crosses(&self, plane: &mut Plane, t: u32) -> Result<bool> {
+        let vs = self.cx.tet(t)?;
+        for (i, j) in [(0, 1), (0, 2), (0, 3), (1, 2), (1, 3), (2, 3)] {
+            if let (Some(&p), Some(&q)) = (vs.get(i), vs.get(j))
+                && self.pierces(plane, p, q)?
+            {
+                return Ok(true);
+            }
+        }
+        Ok(false)
+    }
+
+    /// The vertices on facet `f`'s boundary: its corners and the Steiner points on its edges.
+    pub(super) fn facet_ring(&self, f: usize) -> Result<Vec<u32>> {
+        let tri = *self.triangles.get(f).ok_or_else(|| corrupt("facet id"))?;
+        let mut ring = Vec::new();
+        for k in 0..3 {
+            let (a, b) = (tri[k], tri[(k + 1) % 3]);
+            let s = *self.edge_segment.get(&(a.min(b), a.max(b))).ok_or_else(|| corrupt("facet edge without a segment"))?;
+            let chain = self.segment(s)?.chain_from(a);
+            ring.extend(chain.iter().take(chain.len().saturating_sub(1)));
+        }
+        Ok(ring)
+    }
+
+    /// The face-connected sets of tetrahedra crossing facet `f`'s interior. Each contains a
+    /// tetrahedron incident to a vertex of the facet's boundary `ring`.
+    fn crossing_cavities(&mut self, plane: &mut Plane, ring: &[u32]) -> Result<Vec<Vec<u32>>> {
+        let mut cavities = Vec::new();
+        let mut seen: HashSet<u32> = HashSet::new();
+        let mut star = Vec::new();
+        for &v in ring {
+            self.cx.star(v, &mut star)?;
+            for &t in &star {
+                if seen.contains(&t) || !self.crosses(plane, t)? {
+                    continue;
+                }
+                seen.insert(t);
+                let mut cavity = vec![t];
+                let mut i = 0;
+                while let Some(&c) = cavity.get(i) {
+                    i += 1;
+                    for n in self.cx.adj(c)? {
+                        if n != NONE && !seen.contains(&n) && self.crosses(plane, n)? {
+                            seen.insert(n);
+                            cavity.push(n);
+                        }
+                    }
+                }
+                cavities.push(cavity);
+            }
+        }
+        Ok(cavities)
+    }
+
+    /// Whether facet `f` is a union of mesh triangles: no tetrahedron crosses it.
+    pub(super) fn facet_present(&mut self, f: usize) -> Result<bool> {
+        let ring = self.facet_ring(f)?;
+        let mut plane = self.plane(f)?;
+        Ok(self.crossing_cavities(&mut plane, &ring)?.is_empty())
+    }
+
+    /// Recover facet `f`, then register its triangles as constraints. Filling one cavity leaves
+    /// the others untouched (they share no tetrahedron), so one scan finds them all; the scan is
+    /// repeated until it finds none.
+    pub(super) fn recover_facet(&mut self, f: usize) -> Result<()> {
+        let ring = self.facet_ring(f)?;
+        let mut plane = self.plane(f)?;
+        for _ in 0..MAX_ROUNDS_PER_FACET {
+            let cavities = self.crossing_cavities(&mut plane, &ring)?;
+            if cavities.is_empty() {
+                return self.register_facet(f, &ring);
+            }
+            for cavity in cavities {
+                self.stats.cavities += 1;
+                self.fill_cavity(f, &mut plane, &cavity)?;
+            }
+        }
+        Err(Error::Limit(format!("facet {f} was still crossed after {MAX_ROUNDS_PER_FACET} rounds of cavity filling")))
+    }
+
+    /// Record the mesh triangles that make up facet `f` (all three vertices on its boundary).
+    fn register_facet(&mut self, f: usize, ring: &[u32]) -> Result<()> {
+        let on: HashSet<u32> = ring.iter().copied().collect();
+        let mut star = Vec::new();
+        for &v in ring {
+            self.cx.star(v, &mut star)?;
+            for &t in &star {
+                let vs = self.cx.tet(t)?;
+                for k in 0..4 {
+                    let face = face_toward(vs, k);
+                    if face.contains(&v) && face.iter().all(|u| on.contains(u)) {
+                        let key = sorted3(face);
+                        match self.constraint.insert(key, f as u32) {
+                            Some(g) if g != f as u32 => {
+                                return Err(Error::Failed(format!("input triangles {g} and {f} overlap")));
+                            }
+                            _ => {}
+                        }
+                    }
+                }
+            }
+        }
+        Ok(())
+    }
+
+    /// Delete the cavity and refill its two halves by gift-wrapping.
+    fn fill_cavity(&mut self, f: usize, plane: &mut Plane, cavity: &[u32]) -> Result<()> {
+        let inside: HashSet<u32> = cavity.iter().copied().collect();
+        // Per half: the boundary faces (cavity on their positive side) and the constraint faces
+        // inside it, and the other sides of those constraints.
+        let (mut fixed_up, mut fixed_down) = (Vec::new(), Vec::new());
+        let (mut back_up, mut back_down) = (Vec::new(), Vec::new());
+        let mut vertices: Vec<u32> = Vec::new();
+        for &t in cavity {
+            let vs = self.cx.tet(t)?;
+            vertices.extend(vs);
+            for (j, n) in self.cx.adj(t)?.into_iter().enumerate() {
+                let face = face_toward(vs, j);
+                let interior = n != NONE && inside.contains(&n);
+                if interior && (n < t || !self.constraint.contains_key(&sorted3(face))) {
+                    continue;
+                }
+                let mut pos = false;
+                let mut neg = false;
+                for v in face {
+                    match self.facet_side(plane, v)? {
+                        Sign::Positive => pos = true,
+                        Sign::Negative => neg = true,
+                        Sign::Zero => {}
+                    }
+                }
+                let (fixed, back) = match (pos, neg) {
+                    (true, false) => (&mut fixed_up, &mut back_up),
+                    (false, true) => (&mut fixed_down, &mut back_down),
+                    _ => {
+                        let at = self.cx.approx(face[0])?;
+                        return Err(Error::Failed(format!("facet {f}: a cavity face near {at:?} is not on one side of the facet's plane")));
+                    }
+                };
+                fixed.push(face);
+                if interior {
+                    // A recovered constraint inside the cavity: both sides are open.
+                    let [a, b, c] = face;
+                    back.push([a, c, b]);
+                }
+            }
+        }
+        vertices.sort_unstable();
+        vertices.dedup();
+        let mut sides = Vec::with_capacity(vertices.len());
+        for &v in &vertices {
+            sides.push(self.facet_side(plane, v)?);
+        }
+        let above: Vec<u32> = vertices.iter().zip(&sides).filter(|(_, s)| **s != Sign::Negative).map(|(v, _)| *v).collect();
+        let below: Vec<u32> = vertices.iter().zip(&sides).filter(|(_, s)| **s != Sign::Positive).map(|(v, _)| *v).collect();
+        let on_plane = |v: u32| plane.side.get(&v) == Some(&Sign::Zero);
+
+        let front_up = fixed_up.iter().chain(&back_up).copied().collect();
+        let (mut new, open) = self.gift_wrap(f, front_up, &fixed_up, &above, Some(&on_plane))?;
+        if open.iter().any(|face| !face.iter().all(|&v| on_plane(v))) {
+            return Err(Error::Failed(format!("facet {f}: the upper half-cavity was left open off the facet")));
+        }
+        // The facet's triangles made by the upper half bound the lower one.
+        fixed_down.extend(open.iter().copied());
+        let front_down = fixed_down.iter().chain(&back_down).copied().collect();
+        let (lower, rest) = self.gift_wrap(f, front_down, &fixed_down, &below, None)?;
+        if !rest.is_empty() {
+            return Err(Error::Failed(format!("facet {f}: the lower half-cavity did not close")));
+        }
+        new.extend(lower);
+        let used: HashSet<u32> = new.iter().flatten().copied().collect();
+        if let Some(&lost) = vertices.iter().find(|v| !used.contains(v)) {
+            let at = self.cx.approx(lost)?;
+            return Err(Error::Failed(format!("facet {f}: refilling a cavity lost vertex {lost} at {at:?}")));
+        }
+        self.stats.gift_wrapped_tets += new.len();
+        self.cx.replace(cavity, &new)?;
+        Ok(())
+    }
+
+    /// Fill a half-cavity from its open `front` faces (each with the unfilled region on its
+    /// positive side), using only `vertices`. `fixed` are the triangles bounding it, which no new
+    /// tetrahedron may cross. Faces on the facet's plane are left open when `skip` says so. The new
+    /// tetrahedra and the faces left open.
+    fn gift_wrap(
+        &mut self,
+        f: usize,
+        front: Vec<[u32; 3]>,
+        fixed: &[[u32; 3]],
+        vertices: &[u32],
+        skip: Option<&dyn Fn(u32) -> bool>,
+    ) -> Result<(Vec<[u32; 4]>, Vec<[u32; 3]>)> {
+        let mut faces: Vec<[u32; 3]> = Vec::with_capacity(front.len() * 2);
+        let mut boxes: Vec<([f64; 3], [f64; 3])> = Vec::with_capacity(front.len() * 2);
+        let mut alive: Vec<bool> = Vec::new();
+        let mut index: HashMap<[u32; 3], usize> = HashMap::new();
+        let mut queue: Vec<usize> = Vec::new();
+        for face in front {
+            if index.insert(canon(face), faces.len()).is_some() {
+                return Err(corrupt("a cavity face is listed twice"));
+            }
+            queue.push(faces.len());
+            faces.push(face);
+            boxes.push(self.cx.bounds_of(&face));
+            alive.push(true);
+        }
+        queue.reverse();
+        let fixed_boxes: Vec<([u32; 3], Bounds)> = fixed.iter().map(|&t| (t, self.cx.bounds_of(&t))).collect();
+        let cap = 64 * (faces.len() + vertices.len()) + 64;
+        let mut tets: Vec<[u32; 4]> = Vec::new();
+        while let Some(i) = queue.pop() {
+            if !alive.get(i).copied().unwrap_or(false) {
+                continue;
+            }
+            let sigma = *faces.get(i).ok_or_else(|| corrupt("front index"))?;
+            if let Some(on) = skip
+                && sigma.iter().all(|&v| on(v))
+            {
+                continue;
+            }
+            let w = self.apex(f, sigma, (&faces, &boxes, &alive), &fixed_boxes, vertices)?;
+            let tet = [sigma[0], sigma[1], sigma[2], w];
+            tets.push(tet);
+            if tets.len() > cap {
+                return Err(Error::Limit(format!("facet {f}: gift-wrapping a cavity made more than {cap} tetrahedra")));
+            }
+            if let Some(a) = alive.get_mut(i) {
+                *a = false;
+            }
+            index.remove(&canon(sigma));
+            for k in 0..3 {
+                let g = face_toward(tet, k);
+                if let Some(j) = index.remove(&canon(g)) {
+                    if let Some(a) = alive.get_mut(j) {
+                        *a = false;
+                    }
+                    continue;
+                }
+                let r = [g[0], g[2], g[1]];
+                if index.insert(canon(r), faces.len()).is_some() {
+                    let at = self.cx.approx(r[0])?;
+                    return Err(Error::Failed(format!("facet {f}: gift-wrapping overlapped itself near {at:?}")));
+                }
+                queue.push(faces.len());
+                faces.push(r);
+                boxes.push(self.cx.bounds_of(&r));
+                alive.push(true);
+            }
+        }
+        let open: Vec<[u32; 3]> = faces.iter().zip(&alive).filter(|(_, a)| **a).map(|(f, _)| *f).collect();
+        Ok((tets, open))
+    }
+
+    /// The apex for front face `sigma`: conditions (i)–(iii) of the module docs.
+    fn apex(&mut self, f: usize, sigma: [u32; 3], front: Front<'_>, fixed: &[([u32; 3], Bounds)], vertices: &[u32]) -> Result<u32> {
+        let mut candidates = Vec::new();
+        for &w in vertices {
+            if !sigma.contains(&w) && self.cx.volume([sigma[0], sigma[1], sigma[2], w])? == Sign::Positive {
+                candidates.push(w);
+            }
+        }
+        let mut excluded = vec![false; candidates.len()];
+        for _ in 0..candidates.len() {
+            let mut best: Option<(usize, u32)> = None;
+            for (k, &w) in candidates.iter().enumerate() {
+                if excluded.get(k).copied().unwrap_or(true) {
+                    continue;
+                }
+                best = match best {
+                    None => Some((k, w)),
+                    Some((bk, b)) => {
+                        if self.cx.in_sphere([sigma[0], sigma[1], sigma[2], b], w)? {
+                            Some((k, w))
+                        } else {
+                            Some((bk, b))
+                        }
+                    }
+                };
+            }
+            let Some((k, w)) = best else { break };
+            if self.valid(sigma, w, front, fixed, vertices)? {
+                return Ok(w);
+            }
+            if let Some(e) = excluded.get_mut(k) {
+                *e = true;
+            }
+        }
+        let at = self.cx.approx(sigma[0])?;
+        Err(Error::Failed(format!("facet {f}: gift-wrapping found no valid apex for a cavity face near {at:?}")))
+    }
+
+    /// Condition (ii): the tetrahedron `sigma w` meets every front and fixed triangle properly
+    /// and contains no other cavity vertex.
+    fn valid(&self, sigma: [u32; 3], w: u32, (faces, boxes, alive): Front<'_>, fixed: &[([u32; 3], Bounds)], vertices: &[u32]) -> Result<bool> {
+        let tet = [sigma[0], sigma[1], sigma[2], w];
+        let bx = self.cx.bounds_of(&tet);
+        let points = self.cx.points();
+        for &u in vertices {
+            if !tet.contains(&u) && boxes_meet(bx, self.cx.bounds(u)) && in_closed_tet(points, tet, u)? {
+                return Ok(false);
+            }
+        }
+        for ((face, fb), _) in faces.iter().zip(boxes).zip(alive).filter(|(_, a)| **a) {
+            if boxes_meet(bx, *fb) && !meet_properly(points, &tet, face)? {
+                return Ok(false);
+            }
+        }
+        for (face, fb) in fixed {
+            if boxes_meet(bx, *fb) && !meet_properly(points, &tet, face)? {
+                return Ok(false);
+            }
+        }
+        Ok(true)
+    }
+}

crates/mesh/src/cdt/meet.rs

@@ -0,0 +1,159 @@
+//! Whether two closed simplices meet properly: their intersection is exactly the simplex spanned
+//! by their shared vertices (empty when they share none). The gift-wrapping step of the facet
+//! recovery accepts a new tetrahedron only if it meets every cavity and front triangle properly
+//! (Diazzi et al. 2023, §4.4 condition ii and App. B), and it must be exact on the coplanar
+//! configurations CAD boundaries are made of.
+//!
+//! # Method
+//! Two closed convex polytopes `A`, `B` whose relative interiors are disjoint are weakly separated
+//! by a hyperplane `H` (`A` on one closed side, `B` on the other), and then `A ∩ B = (A ∩ H) ∩
+//! (B ∩ H)`. The separating hyperplanes form a pointed polyhedral cone (in `(n, c)` space) whose
+//! extreme rays are hyperplanes through `d` affinely independent points of `A ∪ B` (`d` the
+//! dimension), so trying the hyperplanes through `d` of the vertices finds one whenever one
+//! exists. When none exists the relative interiors overlap and the intersection is not a common
+//! face. When one exists, every one gives the same `A ∩ B`, and the question recurses into `H`
+//! (one dimension lower) with the vertices on it; a point is the base case. Before that, if all
+//! of one simplex's vertices are shared, that simplex is a face of the other and they meet
+//! properly.
+//!
+//! Every test is an exact predicate: `orient3d` in space, `orient2d` projected along an axis the
+//! plane is not parallel to, and coordinate comparison along an axis the line is not
+//! perpendicular to.
+
+use crate::predicates::{Axis, GenericPoint, Sign, generic_compare, generic_orient2d, generic_orient3d};
+use crate::{Error, Result};
+
+use super::complex::face_toward;
+
+#[derive(Debug, Clone, Copy)]
+enum Level {
+    Space,
+    /// A plane, compared by orientation in its projection along this axis.
+    Plane(Axis),
+    /// A line, compared along this axis.
+    Line(Axis),
+    Point,
+}
+
+fn pt(p: &[GenericPoint], v: u32) -> Result<&GenericPoint> {
+    p.get(v as usize).ok_or_else(|| Error::Failed("intersection test: vertex id out of range".into()))
+}
+
+/// Whether the closed simplices with vertex ids `a` and `b` (each affinely independent, ids of
+/// distinct points) intersect exactly in the hull of their shared vertices.
+pub(super) fn meet_properly(p: &[GenericPoint], a: &[u32], b: &[u32]) -> Result<bool> {
+    proper(p, Level::Space, a, b)
+}
+
+/// Whether `x` is in the closed positive tetrahedron `t` (Diazzi et al. 2023, Alg. 3).
+pub(super) fn in_closed_tet(p: &[GenericPoint], t: [u32; 4], x: u32) -> Result<bool> {
+    for k in 0..4 {
+        let [u, v, w] = face_toward(t, k);
+        // The tetrahedron is on the face's positive side (tet4 volume, i.e. orient3d negative).
+        if generic_orient3d(pt(p, u)?, pt(p, v)?, pt(p, w)?, pt(p, x)?) == Sign::Positive {
+            return Ok(false);
+        }
+    }
+    Ok(true)
+}
+
+fn subset(x: &[u32], y: &[u32]) -> bool {
+    x.iter().all(|v| y.contains(v))
+}
+
+fn proper(p: &[GenericPoint], level: Level, a: &[u32], b: &[u32]) -> Result<bool> {
+    if subset(a, b) || subset(b, a) {
+        return Ok(true);
+    }
+    // Beyond the subset rule, a point holds one vertex, which both sides would share.
+    let k = match level {
+        Level::Space => 3,
+        Level::Plane(_) => 2,
+        Level::Line(_) => 1,
+        Level::Point => return Ok(true),
+    };
+    let mut all: Vec<u32> = a.to_vec();
+    all.extend(b.iter().filter(|v| !a.contains(v)));
+    // Candidate hyperplanes through k of the points, those spanned by one simplex first (they
+    // separate most pairs).
+    let mut candidates: Vec<[u32; 3]> = Vec::new();
+    combinations(&all, k, &mut [0; 3], 0, 0, &mut candidates);
+    candidates.sort_by_key(|h| !(subset(&h[..k], a) || subset(&h[..k], b)));
+    let mut spans = false;
+    for h in &candidates {
+        let h = &h[..k];
+        if let (Level::Space, &[u, v, w]) = (level, h)
+            && plane_axis(p, [u, v, w])?.is_none()
+        {
+            continue; // collinear: no plane
+        }
+        spans = true;
+        let mut signs: Vec<(u32, Sign)> = Vec::with_capacity(all.len());
+        for &x in &all {
+            signs.push((x, if h.contains(&x) { Sign::Zero } else { side(p, level, h, x)? }));
+        }
+        let sign_of = |v: &u32| signs.iter().find(|(x, _)| x == v).map_or(Sign::Zero, |(_, s)| *s);
+        let has = |set: &[u32], s: Sign| set.iter().any(|v| sign_of(v) == s);
+        let separated = (!has(a, Sign::Negative) && !has(b, Sign::Positive)) || (!has(a, Sign::Positive) && !has(b, Sign::Negative));
+        if !separated {
+            continue;
+        }
+        let on = |set: &[u32]| -> Vec<u32> { set.iter().copied().filter(|v| sign_of(v) == Sign::Zero).collect() };
+        let lower = match (level, h) {
+            (Level::Space, &[u, v, w]) => {
+                Level::Plane(plane_axis(p, [u, v, w])?.ok_or_else(|| Error::Failed("intersection test: degenerate plane".into()))?)
+            }
+            (Level::Plane(_), &[u, v]) => Level::Line(line_axis(p, u, v)?),
+            _ => Level::Point,
+        };
+        return proper(p, lower, &on(a), &on(b));
+    }
+    if !spans && let (Level::Space, &[u, v, ..]) = (level, all.as_slice()) {
+        // All points collinear: compare along their line.
+        return proper(p, Level::Line(line_axis(p, u, v)?), a, b);
+    }
+    Ok(false)
+}
+
+/// An axis the plane through `u v w` is not parallel to (`None`: the points are collinear).
+fn plane_axis(p: &[GenericPoint], [u, v, w]: [u32; 3]) -> Result<Option<Axis>> {
+    let (a, b, c) = (pt(p, u)?, pt(p, v)?, pt(p, w)?);
+    Ok(Axis::ALL.into_iter().find(|&ax| !generic_orient2d(a, b, c, ax).is_zero()))
+}
+
+/// An axis along which two distinct points differ.
+fn line_axis(p: &[GenericPoint], u: u32, v: u32) -> Result<Axis> {
+    let (a, b) = (pt(p, u)?, pt(p, v)?);
+    Axis::ALL
+        .into_iter()
+        .find(|&ax| generic_compare(a, b, ax).is_ne())
+        .ok_or_else(|| Error::Failed("intersection test: two vertices are the same point".into()))
+}
+
+fn side(p: &[GenericPoint], level: Level, h: &[u32], x: u32) -> Result<Sign> {
+    let px = pt(p, x)?;
+    Ok(match (level, h) {
+        (Level::Space, &[u, v, w]) => generic_orient3d(pt(p, u)?, pt(p, v)?, pt(p, w)?, px),
+        (Level::Plane(ax), &[u, v]) => generic_orient2d(pt(p, u)?, pt(p, v)?, px, ax),
+        (Level::Line(ax), &[u]) => match generic_compare(px, pt(p, u)?, ax) {
+            std::cmp::Ordering::Less => Sign::Negative,
+            std::cmp::Ordering::Equal => Sign::Zero,
+            std::cmp::Ordering::Greater => Sign::Positive,
+        },
+        _ => return Err(Error::Failed("intersection test: malformed hyperplane".into())),
+    })
+}
+
+/// The `k`-subsets of `all` (`k ≤ 3`), each in the first `k` slots.
+fn combinations(all: &[u32], k: usize, cur: &mut [u32; 3], len: usize, from: usize, out: &mut Vec<[u32; 3]>) {
+    if len == k {
+        out.push(*cur);
+        return;
+    }
+    for i in from..all.len() {
+        if let (Some(&v), Some(slot)) = (all.get(i), cur.get_mut(len)) {
+            *slot = v;
+            combinations(all, k, cur, len + 1, i + 1, out);
+        }
+    }
+}

crates/mesh/src/cdt/mod.rs

@@ -0,0 +1,361 @@
+//! Constrained Delaunay tetrahedralization of a closed boundary surface: every input triangle is
+//! a union of tetrahedron faces, with no tolerance anywhere in the decisions.
+//!
+//! # Method
+//! L. Diazzi, D. Panozzo, A. Vaxman, M. Attene, *Constrained Delaunay Tetrahedrization: A Robust
+//! and Practical Approach*, ACM TOG 42(6):181, 2023 (arXiv 2309.09805; paper only, the authors'
+//! code is GPL/LGPL and is not read), on the theory of H. Si, K. Gärtner, *Meshing Piecewise
+//! Linear Complexes by Constrained Delaunay Tetrahedralizations*, IMR 2005, and J. R. Shewchuk,
+//! *Constrained Delaunay Tetrahedralizations and Provably Good Boundary Recovery*, IMR 2002 (also
+//! papers only; TetGen is AGPL).
+//!
+//! 1. The Delaunay tetrahedralization ([`crate::delaunay`]) of the surface's vertices and the
+//!    eight corners of a box around them, so every Steiner point lies strictly inside the hull.
+//! 2. **Segment recovery** ([`segments`], §3.3): every edge of the surface is a segment; a missing
+//!    one is split at an implicit Steiner point (an [`Lnc`](crate::predicates::Lnc) point on the
+//!    original segment, §4.3) placed by Si's rule and inserted by Bowyer–Watson with the same
+//!    perturbed in-sphere test as the Delaunay step, until all segments are edges.
+//! 3. **Facet recovery** ([`facets`], §3.4, §4.4): every input triangle is a facet; the
+//!    tetrahedra crossing a missing facet are deleted and the two half-cavities refilled by the
+//!    modified gift-wrapping, which builds the constrained Delaunay tetrahedralization of each
+//!    half with exact predicates on implicit points.
+//! 4. **Inside/outside** (§4.5): a flood fill from the box corners that flips at every facet
+//!    triangle; the orientation of every facet triangle against its inside tetrahedron is then
+//!    checked exactly.
+//!
+//! Every decision is an exact predicate ([`crate::predicates`]): `orient3d`, `orient2d` and
+//! coordinate comparison on explicit and LNC points, and `insphere` with the Devillers–Teillaud
+//! symbolic perturbation. The result is checked before it is returned: every segment is an edge,
+//! no edge crosses any facet, every tetrahedron is positive (exactly), and the labelling is
+//! consistent with the surface's orientation.
+//!
+//! # Output
+//! [`TetMesh`]: the inside tetrahedra in tet4 order, the surface's vertices first (same ids) and
+//! then the Steiner points, each with its segment and parameter, rounded to doubles; and for each
+//! boundary tetrahedron face the input triangle it lies on (so the B-rep face tag carries over).
+//! Rounding a Steiner point to the nearest doubles moves it off its segment by under an ulp and
+//! can, on a nearly flat tetrahedron, flip it (§5.1). [`round`] then tries the doubles a few ulps
+//! around such a point for one that keeps every incident tetrahedron positive; [`TetMesh::gate`]
+//! runs the validity gate on the rounded mesh and any failure left is reported, never hidden.
+//!
+//! # Determinism and bounds
+//! Same surface → same mesh: the Delaunay step is deterministic, segments and facets are
+//! processed in a fixed order, and hash maps are only looked up. Steiner points, tetrahedra,
+//! passes and cavities are capped ([`CdtOptions`]); hitting a cap is an [`Error::Limit`].
+
+mod complex;
+mod facets;
+mod meet;
+mod round;
+mod segments;
+
+#[cfg(test)]
+mod tests;
+
+use std::collections::HashMap;
+
+use modelcraft_model::ElementKind;
+
+use crate::delaunay::Triangulation;
+use crate::predicates::{GenericPoint, Point3, Sign};
+use crate::quality::{GateOptions, GateReport, validity_gate_elements};
+use crate::surface::{BoundarySurface, Box3};
+use crate::{Error, Result};
+
+use complex::{Complex, NONE, corrupt, face_toward, sorted3};
+use segments::Segment;
+
+/// Caps on the work of [`tetrahedralize`].
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub struct CdtOptions {
+    /// Most Steiner points added on segments. A boundary whose triangles are long next to
+    /// features much closer than their length (needle triangles) needs Steiner points at that
+    /// spacing all along them; the default is well above what the corpus needs (30k at most) and
+    /// stops such input in seconds.
+    pub max_steiner: usize,
+    /// Most tetrahedra (inside and outside, while building).
+    pub max_tets: usize,
+}
+
+impl Default for CdtOptions {
+    fn default() -> Self {
+        CdtOptions { max_steiner: 50_000, max_tets: 20_000_000 }
+    }
+}
+
+/// A Steiner point: `points[a] + t (points[b] − points[a])` of the surface edge `[a, b]`.
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct SteinerPoint {
+    pub segment: [u32; 2],
+    pub t: f64,
+}
+
+/// A tetrahedron face on the boundary: tetrahedron `tet`'s face opposite its vertex
+/// `opposite` (0..4) lies on surface triangle `triangle`.
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub struct BoundaryFace {
+    pub tet: u32,
+    pub opposite: u8,
+    pub triangle: u32,
+}
+
+impl BoundaryFace {
+    /// The face's 1-based number in the model's tet4 face order ([`ElementKind::faces`]): the
+    /// face without the opposite vertex.
+    pub fn tet4_face(&self) -> Option<u8> {
+        let k = usize::from(self.opposite);
+        let n = ElementKind::Tet4.faces().iter().position(|f| !f.contains(&k))?;
+        u8::try_from(n + 1).ok()
+    }
+}
+
+/// What the recovery did.
+#[derive(Debug, Clone, Copy, Default, PartialEq, Eq)]
+pub struct CdtStats {
+    /// Tetrahedra of the initial Delaunay tetrahedralization (with the box).
+    pub delaunay_tets: usize,
+    /// Missing subsegments split.
+    pub missing_subsegments: usize,
+    /// Cavities refilled by gift-wrapping.
+    pub cavities: usize,
+    /// Tetrahedra made by gift-wrapping.
+    pub gift_wrapped_tets: usize,
+    /// Tetrahedra in the box at the end.
+    pub box_tets: usize,
+    /// Steiner points moved by a few ulps from their nearest doubles to keep every incident
+    /// tetrahedron positive ([`round`]).
+    pub rounding_adjusted: usize,
+}
+
+/// The constrained Delaunay tetrahedralization of a closed surface (see the module docs).
+#[derive(Debug, Clone, PartialEq)]
+pub struct TetMesh {
+    /// The surface's vertices (same ids as [`BoundarySurface::points`]), then the Steiner
+    /// points rounded to doubles.
+    pub points: Vec<[f64; 3]>,
+    /// The Steiner points, for `points[surface vertices..]`.
+    pub steiner: Vec<SteinerPoint>,
+    /// The inside tetrahedra, tet4 order (positive volume in exact arithmetic).
+    pub tets: Vec<[u32; 4]>,
+    /// Every tetrahedron face on the surface, with its input triangle.
+    pub boundary: Vec<BoundaryFace>,
+    pub stats: CdtStats,
+}
+
+impl TetMesh {
+    /// The validity gate on the rounded coordinates: every tetrahedron positive, and the volume
+    /// equal to `expected_volume` (the surface's, [`BoundarySurface::volume`]) within 1e-9.
+    pub fn gate(&self, expected_volume: f64) -> Result<GateReport> {
+        let coords = |t: &[u32; 4]| -> Result<[[f64; 3]; 4]> {
+            let p = |v: u32| self.points.get(v as usize).copied().ok_or_else(|| corrupt("output vertex id"));
+            Ok([p(t[0])?, p(t[1])?, p(t[2])?, p(t[3])?])
+        };
+        let elements = self.tets.iter().enumerate().map(|(i, t)| Ok((i as u32 + 1, ElementKind::Tet4, coords(t)?))).collect::<Result<Vec<_>>>()?;
+        validity_gate_elements(elements, &GateOptions { expected_volume: Some(expected_volume), volume_tolerance: 1e-9 })
+    }
+}
+
+/// The state of one recovery.
+pub(crate) struct Recovery {
+    cx: Complex,
+    /// The surface's triangles (the facets), vertex ids as in the complex.
+    triangles: Vec<[u32; 3]>,
+    segments: Vec<Segment>,
+    /// Sorted edge → segment.
+    edge_segment: HashMap<(u32, u32), usize>,
+    /// Recovered facet triangles (sorted vertex ids) → facet.
+    constraint: HashMap<[u32; 3], u32>,
+    /// The surface's vertices are the complex's first `surface` vertices, then 8 box corners,
+    /// then the Steiner points.
+    surface: usize,
+    /// Steiner points in creation order: (segment, t).
+    steiner: Vec<(usize, f64)>,
+    options: CdtOptions,
+    stats: CdtStats,
+}
+
+/// The constrained Delaunay tetrahedralization of `surface` (see the module docs).
+pub fn tetrahedralize(surface: &BoundarySurface, options: &CdtOptions) -> Result<TetMesh> {
+    let mut rec = Recovery::new(surface, options)?;
+    rec.classify_acute()?;
+    rec.recover_segments()?;
+    for f in 0..rec.triangles.len() {
+        rec.recover_facet(f)?;
+    }
+    rec.verify()?;
+    let inside = rec.label()?;
+    rec.output(&inside)
+}
+
+impl Recovery {
+    fn new(surface: &BoundarySurface, options: &CdtOptions) -> Result<Self> {
+        let n = surface.points().len();
+        let triangles = surface.triangles().to_vec();
+        if triangles.is_empty() {
+            return Err(Error::Invalid("the surface has no triangles".into()));
+        }
+        let mut points: Vec<[f64; 3]> = surface.points().iter().map(|p| p.map(|x| x + 0.0)).collect();
+        let mut bbox = Box3::EMPTY;
+        for p in &points {
+            bbox.add_point(*p);
+        }
+        let size = (0..3).map(|k| bbox.hi[k] - bbox.lo[k]).fold(0.0, f64::max).max(1.0);
+        for corner in 0..8 {
+            points.push(std::array::from_fn(|k| {
+                let x = if corner >> k & 1 == 0 { bbox.lo[k] - size } else { bbox.hi[k] + size };
+                crate::predicates::snap_to_domain(x) + 0.0
+            }));
+        }
+        let tri = Triangulation::new(&points)?;
+        if tri.dimension() != 3 || (0..points.len() as u32).any(|v| tri.duplicate_of(v).is_some()) {
+            return Err(Error::Invalid("the surface's vertices must be distinct".into()));
+        }
+        let generic = points.iter().map(|&p| Point3::try_from(p).map(GenericPoint::Explicit)).collect::<Result<Vec<_>>>()?;
+        let cx = Complex::from_delaunay(&tri, generic, options.max_tets)?;
+        let mut edges: Vec<(u32, u32)> =
+            triangles.iter().flat_map(|t| (0..3).map(move |k| (t[k].min(t[(k + 1) % 3]), t[k].max(t[(k + 1) % 3])))).collect();
+        edges.sort_unstable();
+        edges.dedup();
+        if let Some(&(a, b)) = edges.iter().find(|(a, b)| a == b || *b as usize >= n) {
+            return Err(Error::Invalid(format!("surface edge {a}–{b} is degenerate or has no vertex")));
+        }
+        let segments: Vec<Segment> = edges.iter().map(|&(a, b)| Segment::new(a, b)).collect();
+        let edge_segment = edges.iter().enumerate().map(|(i, &e)| (e, i)).collect();
+        let stats = CdtStats { delaunay_tets: cx.live_count(), ..CdtStats::default() };
+        Ok(Recovery { cx, triangles, segments, edge_segment, constraint: HashMap::new(), surface: n, steiner: Vec::new(), options: *options, stats })
+    }
+
+    /// Every segment is an edge, no edge crosses a facet, every facet has triangles, every
+    /// tetrahedron is positive.
+    fn verify(&mut self) -> Result<()> {
+        let mut scratch = Vec::new();
+        for s in 0..self.segments.len() {
+            let chain: Vec<u32> = self.segment(s)?.splits.iter().map(|&(_, v)| v).collect();
+            for w in chain.windows(2) {
+                if let &[u, v] = w
+                    && !self.cx.has_edge(u, v, &mut scratch)?
+                {
+                    let at = self.cx.approx(u)?;
+                    return Err(Error::Failed(format!("a boundary edge near {at:?} was lost during facet recovery")));
+                }
+            }
+        }
+        let mut has_faces = vec![false; self.triangles.len()];
+        for &f in self.constraint.values() {
+            if let Some(h) = has_faces.get_mut(f as usize) {
+                *h = true;
+            }
+        }
+        for f in 0..self.triangles.len() {
+            if !has_faces.get(f).copied().unwrap_or(false) || !self.facet_present(f)? {
+                return Err(Error::Failed(format!("input triangle {f} is not a union of tetrahedron faces")));
+            }
+        }
+        let live: Vec<[u32; 4]> = self.cx.live_tets().map(|(_, vs)| vs).collect();
+        for vs in live {
+            if self.cx.volume(vs)? != Sign::Positive {
+                let at = self.cx.approx(vs[0])?;
+                return Err(corrupt(&format!("a tetrahedron near {at:?} is not positive")));
+            }
+        }
+        Ok(())
+    }
+
+    /// Inside (`true`) or outside, per tetrahedron slot: a flood fill from the box corners that
+    /// flips across every facet triangle (Diazzi et al. §4.5), then the exact check that each
+    /// facet triangle has its inside tetrahedron on the side its normal points away from.
+    fn label(&mut self) -> Result<Vec<bool>> {
+        let mut label: Vec<Option<bool>> = vec![None; self.cx.slots()];
+        let corners = self.surface..self.surface + 8;
+        let corner = self.cx.live_tets().find(|(_, vs)| vs.iter().any(|&v| corners.contains(&(v as usize))));
+        let (start, _) = corner.ok_or_else(|| corrupt("no tetrahedron touches the bounding box"))?;
+        if let Some(l) = label.get_mut(start as usize) {
+            *l = Some(false);
+        }
+        let mut queue = vec![start];
+        while let Some(t) = queue.pop() {
+            let inside = label.get(t as usize).copied().flatten().unwrap_or(false);
+            let vs = self.cx.tet(t)?;
+            for (k, nb) in self.cx.adj(t)?.into_iter().enumerate() {
+                let want = inside != self.constraint.contains_key(&sorted3(face_toward(vs, k)));
+                match label.get_mut(nb as usize) {
+                    None => {} // the hull
+                    Some(l @ None) => {
+                        *l = Some(want);
+                        queue.push(nb);
+                    }
+                    Some(Some(got)) if *got != want => {
+                        let at = self.cx.approx(vs[k])?;
+                        return Err(Error::Failed(format!(
+                            "inside/outside labelling is inconsistent near {at:?}: the surface does not bound a volume"
+                        )));
+                    }
+                    Some(_) => {}
+                }
+            }
+        }
+        let live: Vec<(u32, [u32; 4])> = self.cx.live_tets().collect();
+        for (t, vs) in live {
+            match label.get(t as usize).copied().flatten() {
+                None => return Err(corrupt("a tetrahedron the flood fill did not reach")),
+                Some(false) => {}
+                Some(true) => {
+                    for k in 0..4 {
+                        let Some(&f) = self.constraint.get(&sorted3(face_toward(vs, k))) else { continue };
+                        let [a, b, c] = *self.triangles.get(f as usize).ok_or_else(|| corrupt("facet id"))?;
+                        let x = vs.get(k).copied().unwrap_or(NONE);
+                        if self.cx.volume([a, b, c, x])? != Sign::Negative {
+                            let at = self.cx.approx(a)?;
+                            return Err(Error::Failed(format!("input triangle {f} near {at:?} faces into the solid: the surface is inverted")));
+                        }
+                    }
+                }
+            }
+        }
+        Ok(label.into_iter().map(|l| l == Some(true)).collect())
+    }
+
+    fn output(&self, inside: &[bool]) -> Result<TetMesh> {
+        let n = self.surface;
+        // Surface vertices keep their ids; box corners disappear; Steiner points follow.
+        let total = self.cx.points().len();
+        let mut id = vec![NONE; total];
+        let mut points = Vec::with_capacity(n + self.steiner.len());
+        for v in 0..n {
+            if let Some(s) = id.get_mut(v) {
+                *s = v as u32;
+            }
+            points.push(self.cx.approx(v as u32)?);
+        }
+        let mut steiner = Vec::with_capacity(self.steiner.len());
+        for (k, &(s, t)) in self.steiner.iter().enumerate() {
+            let v = n + 8 + k;
+            if let Some(slot) = id.get_mut(v) {
+                *slot = points.len() as u32;
+            }
+            points.push(self.cx.approx(v as u32)?);
+            steiner.push(SteinerPoint { segment: self.segment(s)?.ends, t });
+        }
+        let mut tets = Vec::new();
+        let mut boundary = Vec::new();
+        for (t, vs) in self.cx.live_tets() {
+            if !inside.get(t as usize).copied().unwrap_or(false) {
+                continue;
+            }
+            let mapped = vs.map(|v| id.get(v as usize).copied().unwrap_or(NONE));
+            if mapped.contains(&NONE) {
+                return Err(corrupt("an inside tetrahedron uses a bounding-box corner"));
+            }
+            for k in 0..4 {
+                if let Some(&f) = self.constraint.get(&sorted3(face_toward(vs, k))) {
+                    boundary.push(BoundaryFace { tet: tets.len() as u32, opposite: k as u8, triangle: f });
+                }
+            }
+            tets.push(mapped);
+        }
+        let mut stats = self.stats;
+        stats.box_tets = self.cx.live_count();
+        stats.rounding_adjusted = round::repair(&mut points, &tets, n)?;
+        Ok(TetMesh { points, steiner, tets, boundary, stats })
+    }
+}

crates/mesh/src/cdt/round.rs

@@ -0,0 +1,126 @@
+//! Rounding the Steiner points to doubles without inverting tetrahedra.
+//!
+//! The exact tetrahedralization is valid, but a Steiner point rounded to its nearest doubles is
+//! off its segment by up to half an ulp per coordinate, which flips a tetrahedron whose volume is
+//! below that scale (Diazzi et al. 2023, §5.1, §6.2: the nearly flat ones, typically four Steiner
+//! points on nearly parallel boundary edges). For each tetrahedron that is not positive after
+//! rounding, each of its Steiner points tries the doubles within [`REACH`] ulps of its rounded
+//! coordinates, nearest first, and takes the first with which every tetrahedron around it is
+//! positive (exact `orient3d` on the doubles). Surface vertices never move. The moves are a few
+//! ulps, so boundary faces stay on their triangles to rounding precision. The number of passes is
+//! capped; what is left is for the validity gate to report.
+//!
+//! This is the interim answer. The paper's own (§5.1) improves the exact mesh by flips that never
+//! cross a facet before rounding, which removes most nearly flat tetrahedra instead of nudging
+//! their vertices; it belongs with the flips of the quality stage.
+
+use crate::predicates::{Point3, Sign, tet_volume_sign};
+use crate::{Error, Result};
+
+/// How far, in ulps per coordinate, a rounded Steiner point may move.
+const REACH: i32 = 2;
+/// Passes over the tetrahedra that are not positive.
+const MAX_PASSES: usize = 4;
+
+fn positive(points: &[[f64; 3]], t: &[u32; 4]) -> Result<bool> {
+    let p = |v: u32| -> Result<Point3> {
+        let c = points.get(v as usize).copied().ok_or_else(|| Error::Failed("rounding: vertex id out of range".into()))?;
+        Point3::try_from(c)
+    };
+    Ok(tet_volume_sign(&p(t[0])?, &p(t[1])?, &p(t[2])?, &p(t[3])?) == Sign::Positive)
+}
+
+/// `x` moved by `k` ulps.
+fn step(mut x: f64, k: i32) -> f64 {
+    for _ in 0..k.unsigned_abs() {
+        x = if k > 0 { x.next_up() } else { x.next_down() };
+    }
+    x
+}
+
+/// Offsets within [`REACH`], nearest (smallest total) first, in a fixed order.
+fn offsets() -> Vec<[i32; 3]> {
+    let r = -REACH..=REACH;
+    let mut o: Vec<[i32; 3]> = r.clone().flat_map(|a| r.clone().flat_map(move |b| (-REACH..=REACH).map(move |c| [a, b, c]))).collect();
+    o.retain(|d| d != &[0, 0, 0]);
+    o.sort_by_key(|d| (d.iter().map(|k| k.abs()).sum::<i32>(), *d));
+    o
+}
+
+/// Move Steiner points (ids `>= surface`) so every tetrahedron is positive where a few-ulp move
+/// does it. The number of points moved.
+pub(super) fn repair(points: &mut [[f64; 3]], tets: &[[u32; 4]], surface: usize) -> Result<usize> {
+    let mut bad = Vec::new();
+    for (i, t) in tets.iter().enumerate() {
+        if !positive(points, t)? {
+            bad.push(i);
+        }
+    }
+    if bad.is_empty() {
+        return Ok(0);
+    }
+    let mut around: Vec<Vec<u32>> = vec![Vec::new(); points.len()];
+    for (i, t) in tets.iter().enumerate() {
+        for &v in t {
+            if (v as usize) >= surface
+                && let Some(a) = around.get_mut(v as usize)
+            {
+                a.push(i as u32);
+            }
+        }
+    }
+    let offsets = offsets();
+    let mut moved = vec![false; points.len()];
+    for _ in 0..MAX_PASSES {
+        let mut still = Vec::new();
+        for &i in &bad {
+            let t = *tets.get(i).ok_or_else(|| Error::Failed("rounding: tetrahedron id".into()))?;
+            if positive(points, &t)? {
+                continue;
+            }
+            let mut fixed = false;
+            for v in t.into_iter().filter(|&v| (v as usize) >= surface) {
+                let Some(origin) = points.get(v as usize).copied() else { continue };
+                let incident = around.get(v as usize).cloned().unwrap_or_default();
+                for d in &offsets {
+                    let candidate: [f64; 3] = std::array::from_fn(|k| step(origin[k], d[k]));
+                    if candidate.iter().any(|x| !crate::predicates::is_admissible(*x)) {
+                        continue;
+                    }
+                    if let Some(p) = points.get_mut(v as usize) {
+                        *p = candidate;
+                    }
+                    let mut ok = true;
+                    for &j in &incident {
+                        let u = tets.get(j as usize).ok_or_else(|| Error::Failed("rounding: tetrahedron id".into()))?;
+                        if !positive(points, u)? {
+                            ok = false;
+                            break;
+                        }
+                    }
+                    if ok {
+                        fixed = true;
+                        if let Some(m) = moved.get_mut(v as usize) {
+                            *m = true;
+                        }
+                        break;
+                    }
+                    if let Some(p) = points.get_mut(v as usize) {
+                        *p = origin;
+                    }
+                }
+                if fixed {
+                    break;
+                }
+            }
+            if !fixed {
+                still.push(i);
+            }
+        }
+        if still.is_empty() || still.len() == bad.len() {
+            break;
+        }
+        bad = still;
+    }
+    Ok(moved.iter().filter(|m| **m).count())
+}

crates/mesh/src/cdt/segments.rs

@@ -0,0 +1,279 @@
+//! Segment recovery (Si & Gärtner 2005, as restated by Diazzi et al. 2023, §3.3): while an edge
+//! of the input surface (a *segment*) is not an edge of the Delaunay tetrahedralization, split
+//! it at a Steiner point and insert the point (Bowyer–Watson). Each Steiner point is the LNC
+//! point `a + t (b − a)` of its original segment `a b` (Diazzi et al. §4.2–§4.3): it lies on the
+//! segment exactly, and only `t` is computed in floating point, so the placement rule may round
+//! but never moves a point off its segment.
+//!
+//! The placement follows Si's reference-point rule. The reference point `r` of a missing
+//! subsegment `v₁ v₂` is the vertex inside its diametral sphere for which the circle through
+//! `v₁, v₂, r` is largest. A vertex is *acute* when two segments meet there at less than 90°;
+//! the rule depends on how many of the original segment's endpoints are acute:
+//! - none: split where the smaller of the spheres centred at `v₁`, `v₂` through `r` cuts the
+//!   segment, or at the midpoint when both radii exceed half the length;
+//! - one (`w`): `p` is where the sphere centred at `w` through `r` cuts the segment; keep `p` if
+//!   it is closer to `r` than to the far endpoint `v_o`, else shift it towards the near endpoint
+//!   `v_w` by `|r p|` when that is under half of `|v_w p|`, else take the midpoint of `v_w p`;
+//! - two: split at the midpoint first, after which each half has one acute endpoint.
+//!
+//! Any value that falls outside the open subsegment becomes the midpoint. Termination is Si's
+//! (the spheres around acute vertices keep splits from cascading); the number of Steiner points
+//! and of passes is still capped.
+
+use crate::predicates::{Axis, GenericPoint, Lnc, Point3, generic_orient2d};
+use crate::vec3::{cross, dot, norm, sub};
+use crate::{Error, Result};
+
+use super::Recovery;
+
+/// Most tetrahedra searched for a reference point.
+const NEAR_CAP: usize = 20_000;
+/// Most passes over the segments.
+const MAX_PASSES: usize = 10_000;
+
+/// An input segment and its split points.
+#[derive(Debug, Clone)]
+pub(super) struct Segment {
+    /// The original endpoints (explicit vertices), `ends[0] < ends[1]`.
+    pub ends: [u32; 2],
+    /// `(t, vertex)` along the segment from `ends[0]` (`t = 0`) to `ends[1]` (`t = 1`).
+    pub splits: Vec<(f64, u32)>,
+    /// Whether each original endpoint is acute.
+    pub acute: [bool; 2],
+}
+
+impl Segment {
+    pub(super) fn new(a: u32, b: u32) -> Self {
+        let ends = [a.min(b), a.max(b)];
+        Segment { ends, splits: vec![(0.0, ends[0]), (1.0, ends[1])], acute: [false; 2] }
+    }
+
+    /// The vertices from `from` to the other end, in order.
+    pub(super) fn chain_from(&self, from: u32) -> Vec<u32> {
+        let mut v: Vec<u32> = self.splits.iter().map(|&(_, v)| v).collect();
+        if from != self.ends[0] {
+            v.reverse();
+        }
+        v
+    }
+}
+
+fn dist(a: [f64; 3], b: [f64; 3]) -> f64 {
+    norm(sub(a, b))
+}
+
+impl Recovery {
+    /// Mark acute vertices: two segments at a vertex make an angle under 90° (floating point:
+    /// this only steers where Steiner points go).
+    pub(super) fn classify_acute(&mut self) -> Result<()> {
+        let n = self.cx.points().len();
+        let mut dirs: Vec<Vec<[f64; 3]>> = vec![Vec::new(); n];
+        for s in &self.segments {
+            let [a, b] = s.ends;
+            let (pa, pb) = (self.cx.approx(a)?, self.cx.approx(b)?);
+            if let Some(d) = dirs.get_mut(a as usize) {
+                d.push(sub(pb, pa));
+            }
+            if let Some(d) = dirs.get_mut(b as usize) {
+                d.push(sub(pa, pb));
+            }
+        }
+        let acute: Vec<bool> = dirs.iter().map(|d| d.iter().enumerate().any(|(i, x)| d.iter().skip(i + 1).any(|y| dot(*x, *y) > 0.0))).collect();
+        for s in &mut self.segments {
+            s.acute = s.ends.map(|v| acute.get(v as usize).copied().unwrap_or(false));
+        }
+        Ok(())
+    }
+
+    /// Split missing segments until every one is an edge of the (Delaunay) tetrahedralization.
+    /// Alg. 2 sweeps all segments until none is missing; an insertion only changes edges between
+    /// vertices of the tetrahedra it deletes, so after the first sweep only the segments through
+    /// those vertices are checked again.
+    pub(super) fn recover_segments(&mut self) -> Result<()> {
+        let mut scratch = Vec::new();
+        let mut through: Vec<Vec<u32>> = vec![Vec::new(); self.cx.points().len()];
+        for (s, seg) in self.segments.iter().enumerate() {
+            for v in seg.ends {
+                if let Some(l) = through.get_mut(v as usize) {
+                    l.push(s as u32);
+                }
+            }
+        }
+        let mut queued = vec![false; self.segments.len()];
+        let mut work: Vec<usize> = (0..self.segments.len()).collect();
+        for _ in 0..MAX_PASSES {
+            if work.is_empty() {
+                return Ok(());
+            }
+            queued.fill(false);
+            let mut next: Vec<usize> = Vec::new();
+            let mut requeue = |s: usize, next: &mut Vec<usize>| {
+                if let Some(q) = queued.get_mut(s)
+                    && !*q
+                {
+                    *q = true;
+                    next.push(s);
+                }
+            };
+            for &s in &work {
+                let mut i = 0;
+                while let (Some(&(tu, u)), Some(&(tv, v))) = (self.segment(s)?.splits.get(i), self.segment(s)?.splits.get(i + 1)) {
+                    if self.cx.has_edge(u, v, &mut scratch)? {
+                        i += 1;
+                        continue;
+                    }
+                    self.stats.missing_subsegments += 1;
+                    let t = self.split_parameter(s, (tu, u), (tv, v))?;
+                    let [a, b] = self.segment(s)?.ends;
+                    let (pa, pb) = (explicit(self.cx.pt(a)?)?, explicit(self.cx.pt(b)?)?);
+                    let w = self.cx.push_point(GenericPoint::Lnc(Lnc::new(&pa, &pb, t)?))?;
+                    through.push(vec![s as u32]);
+                    self.steiner.push((s, t));
+                    if self.steiner.len() > self.options.max_steiner {
+                        let at = self.cx.approx(u)?;
+                        return Err(Error::Limit(format!(
+                            "more than {} Steiner points on boundary edges (the last near {at:?}): the boundary has features much thinner than its triangles are long; raise the cap or tessellate coarser",
+                            self.options.max_steiner
+                        )));
+                    }
+                    let touched = self.cx.insert(w, u)?;
+                    if let Some(seg) = self.segments.get_mut(s) {
+                        seg.splits.insert(i + 1, (t, w));
+                    }
+                    requeue(s, &mut next);
+                    for x in touched {
+                        for &r in through.get(x as usize).map(Vec::as_slice).unwrap_or_default() {
+                            requeue(r as usize, &mut next);
+                        }
+                    }
+                    // Both halves wait for the next sweep (Alg. 2 splits each missing segment
+                    // once per sweep).
+                    i += 2;
+                }
+            }
+            work = next;
+        }
+        Err(Error::Limit(format!("segment recovery did not finish in {MAX_PASSES} passes")))
+    }
+
+    pub(super) fn segment(&self, s: usize) -> Result<&Segment> {
+        self.segments.get(s).ok_or_else(|| super::complex::corrupt("segment id"))
+    }
+
+    /// Si's split parameter for the missing subsegment `(tu, u)`–`(tv, v)` of segment `s`.
+    fn split_parameter(&mut self, s: usize, (tu, u): (f64, u32), (tv, v): (f64, u32)) -> Result<f64> {
+        let (ends, acute) = {
+            let seg = self.segment(s)?;
+            (seg.ends, seg.acute)
+        };
+        let [a, b] = ends;
+        let (pa, pb) = (self.cx.approx(a)?, self.cx.approx(b)?);
+        let length = dist(pa, pb);
+        let mid = 0.5 * (tu + tv);
+        let r = self.reference_point(u, v)?;
+        let chosen = match (r, acute) {
+            (None, _) => mid,
+            // Two acute ends: the first split is the midpoint.
+            (_, [true, true]) if tu < 0.5 && tv > 0.5 => 0.5,
+            (Some(r), acute) => {
+                let pr = self.cx.approx(r)?;
+                let (pu, pv) = (self.cx.approx(u)?, self.cx.approx(v)?);
+                match acute {
+                    [false, false] => {
+                        let (ru, rv) = (dist(pu, pr), dist(pv, pr));
+                        let half = 0.5 * dist(pu, pv);
+                        if ru > half && rv > half {
+                            mid
+                        } else if ru <= rv {
+                            tu + ru / length
+                        } else {
+                            tv - rv / length
+                        }
+                    }
+                    _ => {
+                        // The acute end w: the only one, or the end of the half this
+                        // subsegment lies in.
+                        let w_is_a = match acute {
+                            [true, false] => true,
+                            [false, true] => false,
+                            _ => tv <= 0.5,
+                        };
+                        let pw = if w_is_a { pa } else { pb };
+                        let radius = dist(pw, pr);
+                        let tp = if w_is_a { radius / length } else { 1.0 - radius / length };
+                        let (t_near, t_far) = if w_is_a { (tu, tv) } else { (tv, tu) };
+                        let at = |t: f64| -> [f64; 3] { std::array::from_fn(|k| pa[k] + t * (pb[k] - pa[k])) };
+                        let p = at(tp);
+                        if dist(p, pr) < dist(p, at(t_far)) {
+                            tp
+                        } else {
+                            let d = dist(pr, p);
+                            let near = dist(at(t_near), p);
+                            if d < 0.5 * near { if t_near < tp { tp - d / length } else { tp + d / length } } else { 0.5 * (t_near + tp) }
+                        }
+                    }
+                }
+            }
+        };
+        let t = if chosen.is_finite() && chosen > tu && chosen < tv { chosen } else { mid };
+        if !(t > tu && t < tv) {
+            let at = self.cx.approx(u)?;
+            return Err(Error::Failed(format!("a missing segment near {at:?} is too short to split further in double precision")));
+        }
+        Ok(t)
+    }
+
+    /// The reference point of the missing subsegment `u v`: the vertex strictly inside its
+    /// diametral sphere with the largest circumcircle through `u, v` (floating point; `None` if
+    /// rounding hides every encroaching vertex). Candidates are the vertices of the tetrahedra
+    /// whose boxes meet the subsegment's box, connected to `u`'s star: they include every
+    /// tetrahedron the subsegment crosses, whose vertices are the ones keeping it from being an
+    /// edge. (On the corpus this picks the same point as searching the whole diametral ball, at a
+    /// fraction of the cost.) A vertex exactly on the open subsegment makes the
+    /// input invalid (a segment through a vertex), which is reported.
+    fn reference_point(&mut self, u: u32, v: u32) -> Result<Option<u32>> {
+        let (pu, pv) = (self.cx.approx(u)?, self.cx.approx(v)?);
+        let c: [f64; 3] = std::array::from_fn(|k| 0.5 * (pu[k] + pv[k]));
+        let r2 = 0.25 * dot(sub(pu, pv), sub(pu, pv));
+        let lo: [f64; 3] = std::array::from_fn(|k| pu[k].min(pv[k]));
+        let hi: [f64; 3] = std::array::from_fn(|k| pu[k].max(pv[k]));
+        let near = self.cx.tets_near(u, lo, hi, NEAR_CAP)?;
+        let mut seen: std::collections::HashSet<u32> = std::collections::HashSet::new();
+        let mut best: Option<(f64, u32)> = None;
+        for t in near {
+            for x in self.cx.tet(t)? {
+                if x == u || x == v || !seen.insert(x) {
+                    continue;
+                }
+                let px = self.cx.approx(x)?;
+                let d = sub(px, c);
+                if dot(d, d) >= r2 {
+                    continue;
+                }
+                let (ea, eb) = (sub(pu, px), sub(pv, px));
+                let area2 = norm(cross(ea, eb));
+                let radius = if area2 > 0.0 { norm(ea) * norm(eb) * norm(sub(ea, eb)) / (2.0 * area2) } else { f64::INFINITY };
+                if !radius.is_finite() && self.on_segment(u, v, x)? {
+                    let at = self.cx.approx(x)?;
+                    return Err(Error::Invalid(format!("vertex {x} at {at:?} lies on a boundary edge: the surface is not a valid complex")));
+                }
+                if best.is_none_or(|(b, _)| radius > b) {
+                    best = Some((radius, x));
+                }
+            }
+        }
+        Ok(best.map(|(_, x)| x))
+    }
+
+    /// Whether `x` is exactly collinear with `u v` (it is inside their diametral ball, so then on
+    /// the open segment).
+    fn on_segment(&self, u: u32, v: u32, x: u32) -> Result<bool> {
+        let (a, b, c) = (self.cx.pt(u)?, self.cx.pt(v)?, self.cx.pt(x)?);
+        Ok(Axis::ALL.into_iter().all(|ax| generic_orient2d(a, b, c, ax).is_zero()))
+    }
+}
+
+/// The explicit point of an original segment endpoint.
+fn explicit(p: &GenericPoint) -> Result<Point3> {
+    p.explicit().copied().ok_or_else(|| super::complex::corrupt("a segment endpoint is not an input vertex"))
+}

crates/mesh/src/cdt/tests.rs

@@ -0,0 +1,320 @@
+//! Tests of the constrained Delaunay tetrahedralization: exact conformity to every input
+//! triangle, the enclosed volume, the validity gate, and the classic polyhedra that have no
+//! tetrahedralization without Steiner points.
+
+use super::meet::{in_closed_tet, meet_properly};
+use super::*;
+use crate::predicates::{GenericPoint, Lnc, Point3};
+use crate::surface::{BoundarySurface, SurfaceOptions};
+use crate::vec3::{cross, dot, norm, sub};
+
+/// A closed surface from points and triangles, each triangle oriented outward from `centre`
+/// (the shapes here are star-shaped about it).
+fn surface(points: &[[f64; 3]], triangles: &[[u32; 3]], centre: [f64; 3]) -> BoundarySurface {
+    let oriented: Vec<[u32; 3]> = triangles
+        .iter()
+        .map(|&[a, b, c]| {
+            let (pa, pb, pc) = (points[a as usize], points[b as usize], points[c as usize]);
+            let n = cross(sub(pb, pa), sub(pc, pa));
+            if dot(n, sub(pa, centre)) > 0.0 { [a, b, c] } else { [a, c, b] }
+        })
+        .collect();
+    let faces: Vec<u32> = (0..oriented.len() as u32).collect();
+    BoundarySurface::new(points, &oriented, &faces, &SurfaceOptions::default()).unwrap()
+}
+
+/// The checks every result must pass: the gate (validity, volume), every input triangle covered
+/// by boundary faces of exactly its area, each boundary face on its triangle's plane.
+fn check(s: &BoundarySurface, m: &TetMesh) {
+    let report = m.gate(s.volume()).unwrap();
+    assert!(report.passed(), "{report:?}");
+    let mut area = vec![0.0; s.triangles().len()];
+    for b in &m.boundary {
+        let t = m.tets[b.tet as usize];
+        let face = face_toward(t, b.opposite as usize).map(|v| m.points[v as usize]);
+        area[b.triangle as usize] += 0.5 * norm(cross(sub(face[1], face[0]), sub(face[2], face[0])));
+        let tri = s.triangles()[b.triangle as usize].map(|v| s.points()[v as usize]);
+        let tn = cross(sub(tri[1], tri[0]), sub(tri[2], tri[0]));
+        for p in face {
+            let d = dot(tn, sub(p, tri[0])) / norm(tn);
+            assert!(d.abs() < 1e-9, "{d}");
+        }
+    }
+    for (i, t) in s.triangles().iter().enumerate() {
+        let p = t.map(|v| s.points()[v as usize]);
+        let want = 0.5 * norm(cross(sub(p[1], p[0]), sub(p[2], p[0])));
+        assert!((area[i] - want).abs() <= 1e-9 * want.max(1e-300), "triangle {i}: {} vs {want}", area[i]);
+    }
+}
+
+fn cube() -> (Vec<[f64; 3]>, Vec<[u32; 3]>) {
+    let p: Vec<[f64; 3]> = (0..8).map(|i| [(i & 1) as f64, (i >> 1 & 1) as f64, (i >> 2 & 1) as f64]).collect();
+    let quads = [[0, 1, 3, 2], [4, 5, 7, 6], [0, 1, 5, 4], [2, 3, 7, 6], [0, 2, 6, 4], [1, 3, 7, 5]];
+    let t = quads.iter().flat_map(|q| [[q[0], q[1], q[2]], [q[0], q[2], q[3]]]).collect();
+    (p, t)
+}
+
+#[test]
+fn a_cube_conforms() {
+    let (p, t) = cube();
+    let s = surface(&p, &t, [0.5; 3]);
+    let m = tetrahedralize(&s, &CdtOptions::default()).unwrap();
+    check(&s, &m);
+    // A square's diagonal is cocircular with the other two corners, so it is not strongly
+    // Delaunay; at most one split per face (its midpoint) recovers it.
+    assert!(m.steiner.len() <= 6, "{:?}", m.steiner);
+    assert_eq!(m.boundary.len(), 12 + 2 * m.steiner.len());
+}
+
+/// Schönhardt's twisted prism (Schönhardt 1928): no tetrahedralization without Steiner points.
+fn schonhardt(twist_deg: f64) -> (Vec<[f64; 3]>, Vec<[u32; 3]>) {
+    let mut p = Vec::new();
+    for z in 0..2 {
+        for i in 0..3 {
+            let a = (i as f64) * 2.0 * std::f64::consts::PI / 3.0 + if z == 1 { twist_deg.to_radians() } else { 0.0 };
+            p.push([a.cos(), a.sin(), z as f64 * 1.5]);
+        }
+    }
+    let mut t = vec![[0, 1, 2], [3, 4, 5]];
+    for i in 0..3u32 {
+        let j = (i + 1) % 3;
+        // The reflex diagonal: from bottom i to top j (the top is turned towards j).
+        t.push([i, j, 3 + j]);
+        t.push([i, 3 + j, 3 + i]);
+    }
+    (p, t)
+}
+
+#[test]
+fn schonhardt_prism_gets_steiner_points_and_conforms() {
+    let (p, t) = schonhardt(30.0);
+    let s = surface(&p, &t, [0.0, 0.0, 0.75]);
+    let m = tetrahedralize(&s, &CdtOptions::default()).unwrap();
+    check(&s, &m);
+    assert!(!m.steiner.is_empty(), "Schönhardt's prism has no tetrahedralization without Steiner points");
+    for sp in &m.steiner {
+        assert!(sp.t > 0.0 && sp.t < 1.0);
+    }
+}
+
+/// A UV sphere: many cospherical vertices.
+fn uv_sphere(rings: usize, sectors: usize) -> (Vec<[f64; 3]>, Vec<[u32; 3]>) {
+    let mut p = vec![[0.0, 0.0, 1.0], [0.0, 0.0, -1.0]];
+    for r in 1..rings {
+        let th = std::f64::consts::PI * r as f64 / rings as f64;
+        for s in 0..sectors {
+            let ph = 2.0 * std::f64::consts::PI * s as f64 / sectors as f64;
+            p.push([th.sin() * ph.cos(), th.sin() * ph.sin(), th.cos()]);
+        }
+    }
+    let id = |r: usize, s: usize| (2 + (r - 1) * sectors + s % sectors) as u32;
+    let mut t = Vec::new();
+    for s in 0..sectors {
+        t.push([0, id(1, s), id(1, s + 1)]);
+        t.push([1, id(rings - 1, s + 1), id(rings - 1, s)]);
+        for r in 1..rings - 1 {
+            t.push([id(r, s), id(r + 1, s), id(r + 1, s + 1)]);
+            t.push([id(r, s), id(r + 1, s + 1), id(r, s + 1)]);
+        }
+    }
+    (p, t)
+}
+
+#[test]
+fn a_sphere_conforms_and_is_deterministic() {
+    let (p, t) = uv_sphere(12, 24);
+    let s = surface(&p, &t, [0.0; 3]);
+    let a = tetrahedralize(&s, &CdtOptions::default()).unwrap();
+    check(&s, &a);
+    let b = tetrahedralize(&s, &CdtOptions::default()).unwrap();
+    assert_eq!(a, b);
+}
+
+#[test]
+fn caps_are_reported() {
+    let (p, t) = schonhardt(30.0);
+    let s = surface(&p, &t, [0.0, 0.0, 0.75]);
+    let none = CdtOptions { max_steiner: 0, max_tets: 1_000_000 };
+    assert!(matches!(tetrahedralize(&s, &none), Err(Error::Limit(_))));
+    let few = CdtOptions { max_steiner: 1000, max_tets: 10 };
+    assert!(matches!(tetrahedralize(&s, &few), Err(Error::Limit(_))));
+}
+
+fn e(x: f64, y: f64, z: f64) -> GenericPoint {
+    GenericPoint::Explicit(Point3::new(x, y, z).unwrap())
+}
+
+#[test]
+fn proper_meeting_handles_coplanar_and_shared_cases() {
+    let p = vec![
+        e(0.0, 0.0, 0.0),  // 0
+        e(1.0, 0.0, 0.0),  // 1
+        e(0.0, 1.0, 0.0),  // 2
+        e(0.0, 0.0, 1.0),  // 3: tet 0 1 2 3
+        e(1.0, 1.0, 0.0),  // 4: beside face 0 1 2, in its plane
+        e(0.2, 0.2, 0.0),  // 5: inside face 0 1 2
+        e(0.0, 0.0, -1.0), // 6: below
+        e(0.5, 0.5, 2.0),  // 7: above, off the tet
+        e(2.0, 2.0, 2.0),  // 8: far
+        e(0.1, 0.1, 0.1),  // 9: inside the tet
+    ];
+    let tet = [0, 1, 2, 3];
+    let yes = |b: &[u32]| assert!(meet_properly(&p, &tet, b).unwrap(), "{b:?}");
+    let no = |b: &[u32]| assert!(!meet_properly(&p, &tet, b).unwrap(), "{b:?}");
+    yes(&[0, 1, 2]); // a face
+    yes(&[1, 2, 4]); // coplanar, sharing edge 1 2, on the other side
+    yes(&[0, 1, 6]); // shares an edge, below
+    yes(&[3, 7, 8]); // shares a vertex, away
+    yes(&[6, 8, 4]); // disjoint
+    no(&[0, 1, 5]); // coplanar, overlapping the face 0 1 2
+    no(&[1, 2, 5]); // coplanar, overlapping
+    no(&[6, 7, 4]); // passes through the tet
+    no(&[3, 9, 7]); // shares a vertex, enters the tet
+    no(&[5, 6, 8]); // touches the face 0 1 2 at an interior point
+    assert!(in_closed_tet(&p, tet, 9).unwrap());
+    assert!(in_closed_tet(&p, tet, 5).unwrap());
+    assert!(!in_closed_tet(&p, tet, 8).unwrap());
+    // An LNC point on edge 0 3 is on the tet's boundary.
+    let mut q = p.clone();
+    q.push(GenericPoint::Lnc(Lnc::new(&Point3::new(0.0, 0.0, 0.0).unwrap(), &Point3::new(0.0, 0.0, 1.0).unwrap(), 0.3).unwrap()));
+    assert!(in_closed_tet(&q, tet, 10).unwrap());
+    assert!(!meet_properly(&q, &tet, &[10, 6, 8]).unwrap(), "touches edge 0 3 at a non-shared point");
+    // Tetrahedra sharing a face meet properly; overlapping ones do not.
+    assert!(meet_properly(&p, &tet, &[0, 1, 2, 6]).unwrap());
+    assert!(!meet_properly(&p, &tet, &[0, 1, 2, 9]).unwrap());
+}
+
+/// Two solids: a wide tetrahedron standing on the big triangle z = 0, and a small one whose top
+/// vertex is just under that triangle's centre. The Delaunay edge from the big apex to the small
+/// top pierces the big triangle while all its edges are present: only facet recovery (the
+/// gift-wrapping of a cavity) restores it.
+fn pierced() -> (Vec<[f64; 3]>, Vec<[u32; 3]>, Vec<[f64; 3]>) {
+    let p = vec![
+        [-10.0, -10.0, 0.0],
+        [10.0, -10.0, 0.0],
+        [0.0, 10.0, 0.0],
+        [0.0, 0.0, 5.0],
+        [-1.0, -1.0, -1.0],
+        [1.0, -1.0, -1.0],
+        [0.0, 1.0, -1.0],
+        [0.0, 0.0, -0.01],
+    ];
+    let t = vec![[0, 1, 2], [0, 1, 3], [1, 2, 3], [2, 0, 3], [4, 5, 6], [4, 5, 7], [5, 6, 7], [6, 4, 7]];
+    (p, t, vec![[0.0, -3.0, 1.0], [0.0, -0.3, -0.7]])
+}
+
+/// Orients each triangle outward from the centre of the shell it belongs to.
+fn shells(points: &[[f64; 3]], triangles: &[[u32; 3]], centres: &[[f64; 3]], per_shell: usize) -> BoundarySurface {
+    let oriented: Vec<[u32; 3]> = triangles
+        .iter()
+        .enumerate()
+        .map(|(i, &[a, b, c])| {
+            let centre = centres[i / per_shell];
+            let (pa, pb, pc) = (points[a as usize], points[b as usize], points[c as usize]);
+            let n = cross(sub(pb, pa), sub(pc, pa));
+            if dot(n, sub(pa, centre)) > 0.0 { [a, b, c] } else { [a, c, b] }
+        })
+        .collect();
+    let faces: Vec<u32> = (0..oriented.len() as u32).collect();
+    BoundarySurface::new(points, &oriented, &faces, &SurfaceOptions::default()).unwrap()
+}
+
+#[test]
+fn a_pierced_facet_is_recovered_by_gift_wrapping() {
+    let (p, t, centres) = pierced();
+    let s = shells(&p, &t, &centres, 4);
+    assert_eq!(s.shells().len(), 2);
+    let m = tetrahedralize(&s, &CdtOptions::default()).unwrap();
+    check(&s, &m);
+    assert!(m.stats.cavities > 0, "{:?}", m.stats);
+    assert!(m.stats.gift_wrapped_tets > 0);
+}
+
+/// A Chazelle-style polyhedron (Chazelle 1984): a block whose top is cut by V-notches along `x`
+/// and whose bottom is cut by V-notches along `y`, the notch edges almost touching. Its
+/// tetrahedralizations need many Steiner points.
+fn chazelle(notches: usize, gap: f64) -> (Vec<[f64; 3]>, Vec<[u32; 3]>) {
+    let n = 2 * notches;
+    let at = |i: usize| i as f64 / n as f64;
+    // Top height by y index: ridges at even, valleys at odd. Bottom height by x index.
+    let top = |j: usize| if j.is_multiple_of(2) { 1.0 } else { 0.5 + gap / 2.0 };
+    let bottom = |i: usize| if i.is_multiple_of(2) { 0.0 } else { 0.5 - gap / 2.0 };
+    let mut p = Vec::new();
+    let mut id = std::collections::HashMap::new();
+    let mut vertex = |x: usize, y: usize, up: bool, p: &mut Vec<[f64; 3]>| -> u32 {
+        *id.entry((x, y, up)).or_insert_with(|| {
+            p.push([at(x), at(y), if up { top(y) } else { bottom(x) }]);
+            (p.len() - 1) as u32
+        })
+    };
+    let mut t = Vec::new();
+    let quad = |q: [u32; 4], out: [f64; 3], p: &Vec<[f64; 3]>, t: &mut Vec<[u32; 3]>| {
+        for tri in [[q[0], q[1], q[2]], [q[0], q[2], q[3]]] {
+            let (a, b, c) = (p[tri[0] as usize], p[tri[1] as usize], p[tri[2] as usize]);
+            let nrm = cross(sub(b, a), sub(c, a));
+            t.push(if dot(nrm, out) > 0.0 { tri } else { [tri[0], tri[2], tri[1]] });
+        }
+    };
+    for i in 0..n {
+        for j in 0..n {
+            let q = [vertex(i, j, true, &mut p), vertex(i + 1, j, true, &mut p), vertex(i + 1, j + 1, true, &mut p), vertex(i, j + 1, true, &mut p)];
+            quad(q, [0.0, 0.0, 1.0], &p, &mut t);
+            let q =
+                [vertex(i, j, false, &mut p), vertex(i + 1, j, false, &mut p), vertex(i + 1, j + 1, false, &mut p), vertex(i, j + 1, false, &mut p)];
+            quad(q, [0.0, 0.0, -1.0], &p, &mut t);
+        }
+    }
+    for k in 0..n {
+        for (fixed, out) in [(0, -1.0), (n, 1.0)] {
+            // Walls x = fixed (along y) and y = fixed (along x).
+            let q = [
+                vertex(fixed, k, false, &mut p),
+                vertex(fixed, k + 1, false, &mut p),
+                vertex(fixed, k + 1, true, &mut p),
+                vertex(fixed, k, true, &mut p),
+            ];
+            quad(q, [out, 0.0, 0.0], &p, &mut t);
+            let q = [
+                vertex(k, fixed, false, &mut p),
+                vertex(k + 1, fixed, false, &mut p),
+                vertex(k + 1, fixed, true, &mut p),
+                vertex(k, fixed, true, &mut p),
+            ];
+            quad(q, [0.0, out, 0.0], &p, &mut t);
+        }
+    }
+    (p, t)
+}
+
+#[test]
+fn chazelle_polyhedron_conforms() {
+    let (p, t) = chazelle(3, 0.02);
+    let faces: Vec<u32> = (0..t.len() as u32).collect();
+    let s = BoundarySurface::new(&p, &t, &faces, &SurfaceOptions::default()).unwrap();
+    assert_eq!(s.shells().len(), 1);
+    let m = tetrahedralize(&s, &CdtOptions::default()).unwrap();
+    check(&s, &m);
+    assert!(!m.steiner.is_empty());
+}
+
+#[test]
+fn hostile_and_degenerate_surfaces_are_errors() {
+    // A flat "solid" (all points coplanar) bounds no volume.
+    let p = vec![[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [1.0, 1.0, 0.0]];
+    let t = vec![[0, 1, 2], [1, 3, 2], [0, 2, 1], [1, 2, 3]];
+    let faces = vec![0; 4];
+    if let Ok(s) = BoundarySurface::new(&p, &t, &faces, &SurfaceOptions::default()) {
+        assert!(tetrahedralize(&s, &CdtOptions::default()).is_err());
+    }
+}
+
+#[test]
+fn the_gate_reports_an_inverted_tetrahedron() {
+    let (p, t) = cube();
+    let s = surface(&p, &t, [0.5; 3]);
+    let mut m = tetrahedralize(&s, &CdtOptions::default()).unwrap();
+    m.tets[0].swap(0, 1);
+    let report = m.gate(s.volume()).unwrap();
+    assert_eq!(report.invalid, vec![1]);
+    assert!(report.into_result().is_err());
+}

crates/mesh/src/lib.rs

@@ -7,6 +7,7 @@
 #![deny(clippy::unwrap_used, clippy::expect_used, clippy::panic, clippy::unimplemented, clippy::todo, clippy::unreachable)]
 #![forbid(unsafe_code)]
 
+pub mod cdt;
 pub mod delaunay;
 pub mod quad;
 pub mod quality;

crates/mesh/src/predicates/indirect.rs

@@ -1,7 +1,10 @@
 //! Indirect predicates on implicit points (Attene 2020, §4–§5), from the paper only.
 //!
 //! An implicit point is kept as the explicit points that define it and is never rounded: a
-//! line–plane intersection ([`Lpi`], §4.2) or a three-plane intersection ([`Tpi`]). Each has
+//! line–plane intersection ([`Lpi`], §4.2), a three-plane intersection ([`Tpi`]), or a linear
+//! combination of two points ([`Lnc`], the Steiner points of Diazzi, Panozzo, Vaxman, Attene,
+//! *Constrained Delaunay Tetrahedrization: A Robust and Practical Approach*, ACM TOG 42(6), 2023,
+//! §4.3). Each has
 //! homogeneous coordinates `(X, W)` that are polynomials in the defining coordinates, so a
 //! predicate on implicit points is the sign of a polynomial times the signs of the `W`s
 //! (Attene §4.3–§4.6, §5.1: "count how many of them have a negative sign").
@@ -19,9 +22,10 @@
 use std::cmp::Ordering;
 
 use super::dyadic::Dyadic;
+use super::eft::{two_diff, two_product, two_sum};
 use super::interval::Interval;
-use super::ring::{Ring, Signed, det2, det3};
-use super::{Point2, Point3, Sign, orient2d, orient3d};
+use super::ring::{Ring, Signed, det2, det3, sq3};
+use super::{Point2, Point3, Sign, insphere, orient2d, orient3d};
 use crate::{Error, Result};
 
 /// A coordinate axis: the axis compared by [`generic_compare`], or dropped by
@@ -34,6 +38,9 @@ pub enum Axis {
 }
 
 impl Axis {
+    /// The three axes, in order.
+    pub const ALL: [Axis; 3] = [Axis::X, Axis::Y, Axis::Z];
+
     pub(crate) fn index(self) -> usize {
         match self {
             Axis::X => 0,
@@ -127,6 +134,15 @@ fn tpi_hom<T: Ring>(planes: &[[[f64; 3]; 3]; 3]) -> Hom<T> {
     Hom { x, w: Some(det3(&[n1, n2, n3])) }
 }
 
+/// `a + t (b - a)` (Diazzi et al. 2023, §4.3 write `t v₁ + (1 - t) v₂`; this form keeps every
+/// operand a double, so it is exact in the ring). `W = 1`.
+fn lnc_hom<T: Ring>(a: &[f64; 3], b: &[f64; 3], t: f64) -> Hom<T> {
+    let (a, b, t) = (lift::<T>(a), lift::<T>(b), T::from_f64(t));
+    let d = sub3(&b, &a);
+    let x = std::array::from_fn(|k| coord(&a, k) + t.clone() * coord(&d, k));
+    Hom { x, w: None }
+}
+
 /// Decides the sign of `W` exactly: filter first, integers if needed.
 fn exact_w_sign(filter: &Hom<Interval>, exact: impl FnOnce() -> Hom<Dyadic>) -> Sign {
     w_sign(&filter.w).or_else(|| w_sign(&exact().w)).unwrap_or(Sign::Zero)
@@ -172,12 +188,61 @@ impl Tpi {
     }
 }
 
+/// The point `a + t (b - a)` of segment `a b`, `0 < t < 1`: an LNC point (LiNear Combination,
+/// Diazzi et al. 2023, §4.3), the Steiner points of constrained Delaunay boundary recovery. It
+/// lies exactly on the segment, and points on one segment are ordered by `t`.
+#[derive(Debug, Clone)]
+pub struct Lnc {
+    a: [f64; 3],
+    b: [f64; 3],
+    t: f64,
+    filter: [Interval; 3],
+}
+
+impl Lnc {
+    /// Fails unless `a ≠ b` and `t` is an admissible double strictly between 0 and 1.
+    pub fn new(a: &Point3, b: &Point3, t: f64) -> Result<Self> {
+        if !(t > 0.0 && t < 1.0 && super::is_admissible(t)) {
+            return Err(Error::Invalid(format!("segment parameter {t:e} is not an admissible number strictly between 0 and 1")));
+        }
+        if a == b {
+            return Err(Error::Invalid("a point on a segment needs two distinct endpoints".into()));
+        }
+        let filter = lnc_hom::<Interval>(&a.0, &b.0, t).x;
+        Ok(Self { a: a.0, b: b.0, t, filter })
+    }
+
+    /// The segment parameter `t`.
+    pub fn t(&self) -> f64 {
+        self.t
+    }
+
+    /// The segment's endpoints `a` (at `t = 0`) and `b` (at `t = 1`).
+    pub fn endpoints(&self) -> [[f64; 3]; 2] {
+        [self.a, self.b]
+    }
+
+    /// The coordinates rounded to doubles, within about one unit in the last place of each
+    /// (error-free transforms: `b - a` and `t (b - a)` are kept as two-term sums).
+    pub fn rounded(&self) -> [f64; 3] {
+        std::array::from_fn(|k| {
+            let (a, b) = (self.a[k], self.b[k]);
+            let (d, d_lo) = two_diff(b, a);
+            let (p, p_lo) = two_product(self.t, d);
+            let tail = p_lo + self.t * d_lo;
+            let (s, s_lo) = two_sum(a, p);
+            s + (s_lo + tail)
+        })
+    }
+}
+
 /// An explicit point or an implicit one, for the `generic_*` predicates.
 #[derive(Debug, Clone)]
 pub enum GenericPoint {
     Explicit(Point3),
     Lpi(Lpi),
     Tpi(Tpi),
+    Lnc(Lnc),
 }
 
 impl From<Point3> for GenericPoint {
@@ -187,7 +252,8 @@ impl From<Point3> for GenericPoint {
 }
 
 impl GenericPoint {
-    fn explicit(&self) -> Option<&Point3> {
+    /// The point, when it is explicit.
+    pub fn explicit(&self) -> Option<&Point3> {
         match self {
             GenericPoint::Explicit(p) => Some(p),
             _ => None,
@@ -199,6 +265,7 @@ impl GenericPoint {
             GenericPoint::Explicit(p) => Hom { x: lift(&p.0), w: None },
             GenericPoint::Lpi(l) => l.filter.clone(),
             GenericPoint::Tpi(t) => t.filter.clone(),
+            GenericPoint::Lnc(l) => Hom { x: l.filter, w: None },
         }
     }
 
@@ -207,6 +274,7 @@ impl GenericPoint {
             GenericPoint::Explicit(p) => Hom { x: lift(&p.0), w: None },
             GenericPoint::Lpi(l) => lpi_hom(&l.points),
             GenericPoint::Tpi(t) => tpi_hom(&t.planes),
+            GenericPoint::Lnc(l) => lnc_hom(&l.a, &l.b, l.t),
         }
     }
 
@@ -215,6 +283,7 @@ impl GenericPoint {
     pub fn approx(&self) -> [f64; 3] {
         match self {
             GenericPoint::Explicit(p) => p.0,
+            GenericPoint::Lnc(l) => l.rounded(),
             _ => {
                 let h = self.filter_hom();
                 let w = h.w.map_or(1.0, |w| w.mid());
@@ -222,6 +291,30 @@ impl GenericPoint {
             }
         }
     }
+
+    /// A box that certainly contains the point: `(lo, hi)` per axis. Unbounded (`±∞`) when the
+    /// filter cannot bound the division of an LPI or TPI.
+    pub fn bounds(&self) -> ([f64; 3], [f64; 3]) {
+        match self {
+            GenericPoint::Explicit(p) => (p.0, p.0),
+            GenericPoint::Lnc(l) => (l.filter.map(|x| x.lo()), l.filter.map(|x| x.hi())),
+            _ => {
+                let h = self.filter_hom();
+                let mut lo = [f64::NEG_INFINITY; 3];
+                let mut hi = [f64::INFINITY; 3];
+                if let Some(w) = h.w.filter(|w| w.lo() > 0.0 || w.hi() < 0.0) {
+                    for (k, x) in h.x.iter().enumerate() {
+                        let q = [x.lo() / w.lo(), x.lo() / w.hi(), x.hi() / w.lo(), x.hi() / w.hi()];
+                        if q.iter().all(|v| v.is_finite()) {
+                            lo[k] = q.iter().copied().fold(f64::INFINITY, f64::min).next_down();
+                            hi[k] = q.iter().copied().fold(f64::NEG_INFINITY, f64::max).next_up();
+                        }
+                    }
+                }
+                (lo, hi)
+            }
+        }
+    }
 }
 
 /// Runs a formula through the interval filter, then exactly.
@@ -264,6 +357,41 @@ pub fn generic_orient3d(a: &GenericPoint, b: &GenericPoint, c: &GenericPoint, d:
     if flipped { -sign } else { sign }
 }
 
+/// `det[p_i - e, |p_i - e|²]` (rows `a, b, c, d`), Shewchuk's in-sphere determinant, in
+/// homogeneous form: `R_i = W_e X_i - W_i X_e = W_i W_e (p_i - e)`, and row `i` scaled by the
+/// positive `(W_i W_e)²` becomes `[W_i W_e R_i, |R_i|²]`, so the sign needs no division.
+fn insphere_formula<T: Ring + Signed>([a, b, c, d, e]: [Hom<T>; 5]) -> Option<Sign> {
+    let row = |p: &Hom<T>| -> ([T; 3], T) {
+        let r: [T; 3] = std::array::from_fn(|k| mulw(coord(&p.x, k), &e.w) - mulw(coord(&e.x, k), &p.w));
+        let [x, y, z] = &r;
+        let lift = sq3(x, y, z);
+        let scale = match (&p.w, &e.w) {
+            (None, None) => None,
+            (Some(w), None) | (None, Some(w)) => Some(w.clone()),
+            (Some(w), Some(v)) => Some(w.clone() * v.clone()),
+        };
+        (r.map(|v| mulw(v, &scale)), lift)
+    };
+    let [(ra, la), (rb, lb), (rc, lc), (rd, ld)] = [row(&a), row(&b), row(&c), row(&d)];
+    // Expansion along the lift column: -la·M_a + lb·M_b - lc·M_c + ld·M_d.
+    let m = |p: &[T; 3], q: &[T; 3], r: &[T; 3]| det3(&[p.clone(), q.clone(), r.clone()]);
+    let det = ld * m(&ra, &rb, &rc) - lc * m(&ra, &rb, &rd) + lb * m(&ra, &rc, &rd) - la * m(&rb, &rc, &rd);
+    det.sign()
+}
+
+/// [`insphere`] on explicit or implicit points, same sign convention (positive when `e` is
+/// inside the sphere through `a, b, c, d`, for `orient3d(a, b, c, d) > 0`). Exact.
+pub fn generic_insphere(a: &GenericPoint, b: &GenericPoint, c: &GenericPoint, d: &GenericPoint, e: &GenericPoint) -> Sign {
+    if let (Some(a), Some(b), Some(c), Some(d), Some(e)) = (a.explicit(), b.explicit(), c.explicit(), d.explicit(), e.explicit()) {
+        return insphere(a, b, c, d, e);
+    }
+    // The 5×5 form with rows `[p, |p|², 1]` is alternating in its arguments, so moving an
+    // explicit point to the pivot only flips the sign.
+    let (points, flipped) = pivot_last([a, b, c, d, e]);
+    let sign = decide(points, insphere_formula, insphere_formula);
+    if flipped { -sign } else { sign }
+}
+
 /// `det[a' - c'; b' - c']` of the projections: `det[R_a; R_b] / (W_a W_b W_c²)`.
 fn orient2d_formula<T: Ring + Signed>([a, b, c]: [Hom<T>; 3], (u, v): (usize, usize)) -> Option<Sign> {
     let row = |p: &Hom<T>, k: usize| mulw(coord(&p.x, k), &c.w) - mulw(coord(&c.x, k), &p.w);

crates/mesh/src/predicates/interval.rs

@@ -27,6 +27,16 @@ impl Interval {
         if lo.is_nan() || hi.is_nan() { Self::NAN } else { Self { lo: lo.next_down(), hi: hi.next_up() } }
     }
 
+    /// The lower bound.
+    pub(crate) fn lo(&self) -> f64 {
+        self.lo
+    }
+
+    /// The upper bound.
+    pub(crate) fn hi(&self) -> f64 {
+        self.hi
+    }
+
     /// The midpoint, an approximation of the value.
     pub(crate) fn mid(&self) -> f64 {
         0.5 * self.lo + 0.5 * self.hi

crates/mesh/src/predicates/mod.rs

@@ -13,6 +13,11 @@
 //! - M. Attene, *Indirect Predicates for Geometric Constructions*, CAD 126:102856, 2020 (paper
 //!   only; his code is LGPL and is not read): [`Lpi`], [`Tpi`], [`GenericPoint`] and the
 //!   `generic_*` predicates.
+//! - L. Diazzi, D. Panozzo, A. Vaxman, M. Attene, *Constrained Delaunay Tetrahedrization: A
+//!   Robust and Practical Approach*, ACM TOG 42(6):181, 2023, §4.3 (paper only; the authors'
+//!   code is GPL/LGPL and is not read): the [`Lnc`] point on a segment, [`generic_insphere`],
+//!   and its perturbed form [`generic_insphere_sos`] with the same lexicographic lift as
+//!   [`insphere_sos`].
 //!
 //! # Arithmetic assumptions
 //! IEEE 754 binary64, round-to-nearest-even, no extended intermediate precision, no FMA
@@ -40,8 +45,8 @@ mod ring;
 mod sos;
 
 pub use adaptive::{incircle, incircle_filter, insphere, insphere_filter, orient2d, orient2d_filter, orient3d, orient3d_filter};
-pub use indirect::{Axis, GenericPoint, Lpi, Tpi, generic_compare, generic_orient2d, generic_orient3d, orient2d_projected};
-pub use sos::{incircle_sos, insphere_sos};
+pub use indirect::{Axis, GenericPoint, Lnc, Lpi, Tpi, generic_compare, generic_insphere, generic_orient2d, generic_orient3d, orient2d_projected};
+pub use sos::{generic_insphere_sos, generic_lex, incircle_sos, insphere_sos};
 
 use crate::{Error, Result};
 

crates/mesh/src/predicates/sos.rs

@@ -20,6 +20,7 @@
 
 use std::cmp::Ordering;
 
+use super::indirect::{Axis, GenericPoint, generic_compare, generic_insphere, generic_orient3d};
 use super::{Point2, Point3, Sign, incircle, insphere, orient2d, orient3d};
 
 /// Lexicographic order on coordinates (total: `f64::total_cmp`).
@@ -50,7 +51,14 @@ pub fn insphere_sos(a: &Point3, b: &Point3, c: &Point3, d: &Point3, e: &Point3)
 
 #[cold]
 fn insphere_perturbed(a: &Point3, b: &Point3, c: &Point3, d: &Point3, e: &Point3) -> Sign {
-    for i in ranked(&[a.0, b.0, c.0, d.0, e.0]) {
+    perturbed_cofactors([a, b, c, d, e], ranked(&[a.0, b.0, c.0, d.0, e.0]), orient3d)
+}
+
+/// The sign of the first nonzero lift cofactor, in rank order (`order`: argument positions by
+/// decreasing rank). The single statement of the perturbation rule for every point type.
+fn perturbed_cofactors<P>(points: [&P; 5], order: [usize; 5], orient: impl Fn(&P, &P, &P, &P) -> Sign) -> Sign {
+    let [a, b, c, d, e] = points;
+    for i in order {
         let [p, q, r, s] = match i {
             0 => [b, c, d, e],
             1 => [a, c, d, e],
@@ -58,7 +66,7 @@ fn insphere_perturbed(a: &Point3, b: &Point3, c: &Point3, d: &Point3, e: &Point3
             3 => [a, b, c, e],
             _ => [a, b, c, d],
         };
-        let minor = orient3d(p, q, r, s);
+        let minor = orient(p, q, r, s);
         if !minor.is_zero() {
             // Row i + 1 (1-based), lift column 4: cofactor sign (-1)^(i + 5).
             return if i % 2 == 1 { minor } else { -minor };
@@ -67,6 +75,46 @@ fn insphere_perturbed(a: &Point3, b: &Point3, c: &Point3, d: &Point3, e: &Point3
     Sign::Zero
 }
 
+/// The lexicographic order of two points' exact coordinates (`x`, then `y`, then `z`): the rank
+/// order of the perturbation, for explicit and implicit points alike. On explicit points it is
+/// the order [`insphere_sos`] uses (callers admit `+0.0`, never `-0.0`).
+pub fn generic_lex(p: &GenericPoint, q: &GenericPoint) -> Ordering {
+    if let (Some(p), Some(q)) = (p.explicit(), q.explicit()) {
+        return lex(&p.coords(), &q.coords());
+    }
+    Axis::ALL.into_iter().map(|a| generic_compare(p, q, a)).find(|o| o.is_ne()).unwrap_or(Ordering::Equal)
+}
+
+/// [`insphere_sos`] on explicit or implicit points: the same lift by lexicographic rank, so a
+/// triangulation built with [`insphere_sos`] stays consistent when implicit points are added.
+/// Never [`Sign::Zero`] when `a, b, c, d` is not coplanar; points must be pairwise distinct.
+pub fn generic_insphere_sos(a: &GenericPoint, b: &GenericPoint, c: &GenericPoint, d: &GenericPoint, e: &GenericPoint) -> Sign {
+    if let (Some(a), Some(b), Some(c), Some(d), Some(e)) = (a.explicit(), b.explicit(), c.explicit(), d.explicit(), e.explicit()) {
+        return insphere_sos(a, b, c, d, e);
+    }
+    let exact = generic_insphere(a, b, c, d, e);
+    if !exact.is_zero() {
+        return exact;
+    }
+    let points = [a, b, c, d, e];
+    // Positions by decreasing rank, ties by position: an insertion sort, which needs no total
+    // order guarantee from the comparator to stay in bounds.
+    let mut order = [0usize, 1, 2, 3, 4];
+    for i in 1..order.len() {
+        let mut j = i;
+        while j > 0 {
+            let (Some(&x), Some(&y)) = (order.get(j - 1), order.get(j)) else { break };
+            let (Some(px), Some(py)) = (points.get(x), points.get(y)) else { break };
+            if generic_lex(py, px) != Ordering::Greater {
+                break;
+            }
+            order.swap(j - 1, j);
+            j -= 1;
+        }
+    }
+    perturbed_cofactors(points, order, generic_orient3d)
+}
+
 /// [`incircle`] with symbolic perturbation: never [`Sign::Zero`] when `a, b, c` is not
 /// collinear. Same conventions as [`insphere_sos`].
 #[inline]

crates/mesh/src/predicates/tests.rs

@@ -605,3 +605,138 @@ fn undefined_constructions_are_rejected() {
         assert_eq!(generic_compare(&t, &GenericPoint::Explicit(o), axis), Ordering::Equal);
     }
 }
+
+// ---------------------------------------------------------------------------------------------
+// LNC points and the generic in-sphere (Diazzi et al. 2023, §4.3)
+// ---------------------------------------------------------------------------------------------
+
+fn r_lnc(a: &Point3, b: &Point3, t: f64) -> RPoint {
+    let (a, b) = (rp(a), rp(b));
+    std::array::from_fn(|k| a[k].clone() + q(t) * (b[k].clone() - a[k].clone()))
+}
+
+fn r_insphere(p: &[RPoint; 5]) -> Sign {
+    let rows: Vec<Vec<BigRational>> = (0..4).map(|i| lifted(rsub(&p[i], &p[4]).to_vec())).collect();
+    rsign(&rdet(&rows))
+}
+
+#[test]
+fn lnc_points_lie_on_their_segment_in_parameter_order() {
+    let mut r = Rng(11);
+    for _ in 0..200 {
+        let (a, b) = (rand_pt(&mut r, 10.0), rand_pt(&mut r, 10.0));
+        let (t0, t1) = (r.unit().clamp(1e-9, 0.5), 0.5 + r.unit().clamp(1e-9, 0.49));
+        let (p0, p1) = (GenericPoint::Lnc(Lnc::new(&a, &b, t0).unwrap()), GenericPoint::Lnc(Lnc::new(&a, &b, t1).unwrap()));
+        let (ga, gb) = (GenericPoint::Explicit(a), GenericPoint::Explicit(b));
+        for axis in [Axis::X, Axis::Y, Axis::Z] {
+            assert_eq!(generic_orient2d(&ga, &gb, &p0, axis), Sign::Zero);
+            assert_eq!(generic_orient2d(&p0, &p1, &gb, axis), Sign::Zero);
+            let k = axis.index();
+            let want = r_lnc(&a, &b, t0)[k].cmp(&r_lnc(&a, &b, t1)[k]);
+            assert_eq!(generic_compare(&p0, &p1, axis), want);
+        }
+        let exact = r_lnc(&a, &b, t0);
+        let rounded = p0.approx();
+        let (lo, hi) = p0.bounds();
+        for k in 0..3 {
+            let x = exact[k].to_f64().unwrap();
+            assert!((rounded[k] - x).abs() <= 2.0 * f64::EPSILON * x.abs().max(1e-300), "{} vs {x}", rounded[k]);
+            assert!(q(lo[k]) <= exact[k] && exact[k] <= q(hi[k]));
+        }
+    }
+    let (a, b) = (p3(0.0, 0.0, 0.0), p3(1.0, 1.0, 1.0));
+    for bad in [0.0, 1.0, -0.5, 1.5, f64::NAN, f64::INFINITY, 1e-300] {
+        assert!(Lnc::new(&a, &b, bad).is_err(), "{bad}");
+    }
+    assert!(Lnc::new(&a, &a, 0.5).is_err());
+}
+
+#[test]
+fn generic_orient3d_and_insphere_match_rationals_on_lnc_points() {
+    let mut r = Rng(12);
+    let mut checked = 0;
+    for i in 0..120 {
+        let pts: Vec<Point3> = (0..9).map(|_| rand_pt(&mut r, 1.0)).collect();
+        let t = r.unit().clamp(1e-6, 1.0 - 1e-6);
+        let lnc = Lnc::new(&pts[0], &pts[1], t).unwrap();
+        let Ok(lpi) = Lpi::new(&pts[2], &pts[3], &pts[4], &pts[5], &pts[6]) else { continue };
+        let g = [GenericPoint::Lnc(lnc), GenericPoint::Lpi(lpi)];
+        let rg = [r_lnc(&pts[0], &pts[1], t), r_lpi(&pts[2], &pts[3], &pts[4], &pts[5], &pts[6])];
+        let ex: Vec<Point3> = (0..5).map(|_| rand_pt(&mut r, 1.0)).collect();
+        for mask in 0u32..243 {
+            if !(i + mask as usize).is_multiple_of(5) {
+                continue;
+            }
+            let (mut gp, mut rpv) = (Vec::new(), Vec::new());
+            let mut m = mask;
+            for e in &ex {
+                match m % 3 {
+                    0 => {
+                        gp.push(GenericPoint::Explicit(*e));
+                        rpv.push(rp(e));
+                    }
+                    k => {
+                        gp.push(g[k as usize - 1].clone());
+                        rpv.push(rg[k as usize - 1].clone());
+                    }
+                }
+                m /= 3;
+            }
+            let want = r_orient3d(&[rpv[0].clone(), rpv[1].clone(), rpv[2].clone(), rpv[3].clone()]);
+            assert_eq!(generic_orient3d(&gp[0], &gp[1], &gp[2], &gp[3]), want, "orient mask {mask}");
+            let want = r_insphere(&[rpv[0].clone(), rpv[1].clone(), rpv[2].clone(), rpv[3].clone(), rpv[4].clone()]);
+            assert_eq!(generic_insphere(&gp[0], &gp[1], &gp[2], &gp[3], &gp[4]), want, "insphere mask {mask}");
+            checked += 1;
+        }
+    }
+    assert!(checked > 3000, "{checked}");
+    // The explicit fast path and the formula agree on the sign convention.
+    let (a, b, c, d) = (p3(0.0, 0.0, 0.0), p3(1.0, 0.0, 0.0), p3(0.0, 1.0, 0.0), p3(0.0, 0.0, 1.0));
+    let inside = GenericPoint::Lnc(Lnc::new(&a, &p3(1.0, 1.0, 1.0), 0.25).unwrap());
+    let e = |p: &Point3| GenericPoint::Explicit(*p);
+    let shewchuk_positive = orient3d(&b, &a, &c, &d) == Sign::Positive;
+    let (x, y) = if shewchuk_positive { (b, a) } else { (a, b) };
+    assert_eq!(insphere(&x, &y, &c, &d, &p3(0.25, 0.25, 0.25)), Sign::Positive);
+    assert_eq!(generic_insphere(&e(&x), &e(&y), &e(&c), &e(&d), &inside), Sign::Positive);
+}
+
+/// On a grid, a segment's midpoint (an LNC point) is a grid point: the perturbed in-sphere of
+/// the implicit form must give the explicit form's answer, including every cospherical tie.
+#[test]
+fn perturbed_generic_insphere_agrees_with_explicit_points() {
+    let mut r = Rng(13);
+    let mut ties = 0;
+    for _ in 0..800 {
+        let v: Vec<[f64; 3]> = (0..5).map(|_| [r.int(1) as f64, r.int(1) as f64, r.int(1) as f64]).collect();
+        if (0..5).any(|i| (0..i).any(|j| v[i] == v[j])) {
+            continue;
+        }
+        let p: Vec<Point3> = v.iter().map(|c| p3(c[0], c[1], c[2])).collect();
+        if orient3d(&p[0], &p[1], &p[2], &p[3]).is_zero() {
+            continue;
+        }
+        // Each point as the midpoint of a segment through it (direction from the generator).
+        let g: Vec<GenericPoint> = p
+            .iter()
+            .map(|c| {
+                if r.int(1) == 0 {
+                    return GenericPoint::Explicit(*c);
+                }
+                let d = [r.int(2) as f64, r.int(2) as f64, (r.int(2) as f64) + 7.0];
+                let x = c.coords();
+                let a = p3(x[0] - d[0], x[1] - d[1], x[2] - d[2]);
+                let b = p3(x[0] + d[0], x[1] + d[1], x[2] + d[2]);
+                GenericPoint::Lnc(Lnc::new(&a, &b, 0.5).unwrap())
+            })
+            .collect();
+        let want = insphere_sos(&p[0], &p[1], &p[2], &p[3], &p[4]);
+        assert_eq!(generic_insphere_sos(&g[0], &g[1], &g[2], &g[3], &g[4]), want, "{v:?}");
+        ties += insphere(&p[0], &p[1], &p[2], &p[3], &p[4]).is_zero() as usize;
+        for i in 0..5 {
+            for j in 0..5 {
+                assert_eq!(generic_lex(&g[i], &g[j]), v[i].partial_cmp(&v[j]).unwrap());
+            }
+        }
+    }
+    assert!(ties > 20, "{ties}");
+}