From aaa721132e951fc2ac4836839add63a120e74082 Mon Sep 17 00:00:00 2001 From: Bas Nijholt Date: Wed, 10 Jun 2026 11:30:10 -0700 Subject: [PATCH 1/3] feat: exact geometric predicates (robust crate) and precise volumes Add three primitives to geometry, all backed by Shewchuk's adaptive-precision predicates via the robust crate in 2D/3D: - robust_orientation: like orientation() but exact, with no log-det cutoff; resolves the true sign of determinants that slogdet rounds to zero on slivers (falls back to orientation() in other dimensions) - robust_in_circumsphere: strict exact in-circumsphere test (incircle/insphere normalized by the simplex orientation so the answer is orientation-independent); None outside 2D/3D - precise_volume: simplex volume from the adaptive-precision determinant value, immune to the cancellation that corrupts plain f64 volumes of high-aspect slivers These deliberately do NOT replace the reference-compatible tolerance based tests; they exist for the cavity validation/repair path added in the next commit, where geometric truth matters more than bug compatibility. Sign conventions are pinned by unit tests against orientation() on well-conditioned inputs plus a below-cutoff sliver case. --- Cargo.lock | 7 ++ Cargo.toml | 1 + src/geometry.rs | 241 ++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 249 insertions(+) diff --git a/Cargo.lock b/Cargo.lock index 19c3830..2c6dcea 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -10,6 +10,7 @@ dependencies = [ "nalgebra", "numpy", "pyo3", + "robust", "rustc-hash", "smallvec", "thiserror", @@ -601,6 +602,12 @@ version = "0.8.10" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "dc897dd8d9e8bd1ed8cdad82b5966c3e0ecae09fb1907d58efaa013543185d0a" +[[package]] +name = "robust" +version = "1.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4e27ee8bb91ca0adcf0ecb116293afa12d393f9c2b9b9cd54d33e8078fe19839" + [[package]] name = "rustc-hash" version = "2.1.2" diff --git a/Cargo.toml b/Cargo.toml index 90003d6..5cbfa7c 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -23,6 +23,7 @@ nalgebra = "0.33" thiserror = "2" rustc-hash = "2" smallvec = "1" +robust = "1.2.0" [profile.release] lto = true diff --git a/src/geometry.rs b/src/geometry.rs index ae8f7c1..df71425 100644 --- a/src/geometry.rs +++ b/src/geometry.rs @@ -441,6 +441,197 @@ pub fn orientation(face: &[Vec], origin: &[f64]) -> Result], origin: &[f64]) -> Result { + let dim = validate_points(face)?; + if face.len() != dim || origin.len() != dim { + return Err(GeometryError::InvalidDimensions( + "Face and origin dimensions do not match".to_string(), + )); + } + + match dim { + // orient2d(a, b, c) = det [[a - c], [b - c]], exactly the + // determinant `orientation` takes the sign of. + 2 => { + let det = robust::orient2d( + robust::Coord { + x: face[0][0], + y: face[0][1], + }, + robust::Coord { + x: face[1][0], + y: face[1][1], + }, + robust::Coord { + x: origin[0], + y: origin[1], + }, + ); + Ok(if det > 0.0 { + 1 + } else if det < 0.0 { + -1 + } else { + 0 + }) + } + // orient3d(a, b, c, d) = det [[a - d], [b - d], [c - d]], exactly the + // determinant `orientation` takes the sign of. + 3 => { + let coord = |point: &[f64]| robust::Coord3D { + x: point[0], + y: point[1], + z: point[2], + }; + let det = robust::orient3d( + coord(&face[0]), + coord(&face[1]), + coord(&face[2]), + coord(origin), + ); + Ok(if det > 0.0 { + 1 + } else if det < 0.0 { + -1 + } else { + 0 + }) + } + _ => orientation(face, origin), + } +} + +/// Whether `point` lies strictly inside the circumsphere of `simplex`, +/// decided exactly with Shewchuk's adaptive-precision incircle/insphere +/// predicates. Returns `None` for dimensions other than 2 and 3 (no exact +/// predicate available). A degenerate simplex (zero orientation) has no +/// finite circumsphere and reports `false`, matching the NaN-radius +/// convention of [`circumsphere`]. +/// +/// Note the deliberate difference from the reference's tolerance test +/// (`distance < radius * (1 + CIRCUMCIRCLE_RTOL)`): this is strict +/// containment with no slack. Use it where geometric truth matters more +/// than bug compatibility (e.g. cavity repair). +pub fn robust_in_circumsphere( + simplex: &[Vec], + point: &[f64], +) -> Result, GeometryError> { + let dim = validate_points(simplex)?; + if simplex.len() != dim + 1 || point.len() != dim { + return Err(GeometryError::InvalidDimensions(format!( + "Expected {} points for a {}-dimensional simplex", + dim + 1, + dim + ))); + } + + match dim { + // incircle's sign follows the triangle's orientation; multiplying + // by orient2d makes the test orientation-independent. + 2 => { + let coord = |p: &[f64]| robust::Coord { x: p[0], y: p[1] }; + let orient = + robust::orient2d(coord(&simplex[0]), coord(&simplex[1]), coord(&simplex[2])); + if orient == 0.0 { + return Ok(Some(false)); + } + let incircle = robust::incircle( + coord(&simplex[0]), + coord(&simplex[1]), + coord(&simplex[2]), + coord(point), + ); + Ok(Some(incircle * orient > 0.0)) + } + 3 => { + let coord = |p: &[f64]| robust::Coord3D { + x: p[0], + y: p[1], + z: p[2], + }; + let orient = robust::orient3d( + coord(&simplex[0]), + coord(&simplex[1]), + coord(&simplex[2]), + coord(&simplex[3]), + ); + if orient == 0.0 { + return Ok(Some(false)); + } + let insphere = robust::insphere( + coord(&simplex[0]), + coord(&simplex[1]), + coord(&simplex[2]), + coord(&simplex[3]), + coord(point), + ); + Ok(Some(insphere * orient > 0.0)) + } + _ => Ok(None), + } +} + +/// Like [`volume`], but computed with adaptive-precision determinants in 2D +/// and 3D (the value of Shewchuk's orientation predicate is the exact edge +/// determinant, evaluated to machine relative accuracy), so cancellation +/// cannot corrupt the result for sliver simplices. Higher dimensions fall +/// back to [`volume`]. +/// +/// Not used for reference-compatible answers (the reference computes plain +/// floating-point determinants); use it where an accurate value matters, +/// e.g. the volume-conservation check in Bowyer-Watson. +pub fn precise_volume(vertices: &[Vec]) -> Result { + let dim = validate_points(vertices)?; + if vertices.len() != dim + 1 { + return Err(GeometryError::InvalidDimensions(format!( + "Expected {} points for a {}-dimensional simplex", + dim + 1, + dim + ))); + } + + match dim { + 2 => { + let coord = |point: &[f64]| robust::Coord { + x: point[0], + y: point[1], + }; + let det = robust::orient2d( + coord(&vertices[1]), + coord(&vertices[2]), + coord(&vertices[0]), + ); + Ok(det.abs() / 2.0) + } + 3 => { + let coord = |point: &[f64]| robust::Coord3D { + x: point[0], + y: point[1], + z: point[2], + }; + let det = robust::orient3d( + coord(&vertices[1]), + coord(&vertices[2]), + coord(&vertices[3]), + coord(&vertices[0]), + ); + Ok(det.abs() / 6.0) + } + _ => volume(vertices), + } +} + /// Volume of a `dim`-simplex given its `dim + 1` vertices /// (`|det| / dim!` of the edge matrix). pub fn volume(vertices: &[Vec]) -> Result { @@ -543,6 +734,56 @@ mod tests { assert!((radius - 1.0).abs() < 1e-12); } + #[test] + fn robust_orientation_matches_orientation_on_well_conditioned_input() { + let faces_2d = [ + (vec![vec![1.0, 0.0], vec![0.0, 1.0]], vec![0.1, 0.1]), + (vec![vec![0.0, 1.0], vec![1.0, 0.0]], vec![0.1, 0.1]), + (vec![vec![3.0, -2.0], vec![-1.5, 4.0]], vec![10.0, -7.0]), + ]; + for (face, origin) in &faces_2d { + assert_eq!( + robust_orientation(face, origin).unwrap(), + orientation(face, origin).unwrap() + ); + } + + let faces_3d = [ + ( + vec![ + vec![1.0, 0.0, 0.0], + vec![0.0, 1.0, 0.0], + vec![0.0, 0.0, 1.0], + ], + vec![0.0, 0.0, 0.0], + ), + ( + vec![ + vec![2.0, 1.0, -1.0], + vec![-1.0, 3.0, 0.5], + vec![0.5, -2.0, 4.0], + ], + vec![5.0, 5.0, 5.0], + ), + ]; + for (face, origin) in &faces_3d { + assert_eq!( + robust_orientation(face, origin).unwrap(), + orientation(face, origin).unwrap() + ); + } + } + + #[test] + fn robust_orientation_resolves_signs_below_the_log_det_cutoff() { + // A sliver whose determinant is far below exp(-50): the slogdet-based + // test gives up and reports 0, the exact predicate knows the sign. + let face = vec![vec![1e-12, 0.0], vec![0.0, 1e-12]]; + let origin = vec![1e-13, 1e-13]; + assert_eq!(orientation(&face, &origin).unwrap(), 0); + assert_eq!(robust_orientation(&face, &origin).unwrap(), 1); + } + #[test] fn circumsphere_is_accurate_for_small_simplices_far_from_origin() { let (center, radius) = circumsphere(&[vec![0.5], vec![0.5 + 1e-6]]).unwrap(); From 59246dad93b8d24107dcdee4f2a586cf3f4f7f2e Mon Sep 17 00:00:00 2001 From: Bas Nijholt Date: Wed, 10 Jun 2026 11:30:10 -0700 Subject: [PATCH 2/3] fix: validate Bowyer-Watson cavities before mutating; repair or reject Sliver simplices (mixed coordinate scales force aspect ratios ~1e7) feed the floating-point circumsphere test cancellation noise, so the cascade can assemble a cavity that is not re-triangulable: voids or protrusions (volume not conserved), a falsely-empty cavity that leaves the inserted point unconnected, or a cavity that swallows every simplex of the point while its total volume sits below the conservation check's absolute tolerance (making that check vacuous). The Python reference asserts volume conservation only AFTER mutating - a failing insertion corrupts its state - and silently orphans vertices in the cases its check cannot see. bowyer_watson now computes the would-be (deleted, created) outcome and validates it read-only, BEFORE any mutation: - volume conservation is checked with adaptive-precision determinant volumes (geometry::precise_volume), so the check compares true volumes instead of cancellation noise - the point must end up connected (catches both the empty-cavity and the swallowed-hull-simplices failure modes) When validation fails, the cavity is repaired and re-validated: in 2D/3D it is rebuilt from the seeds with exact insphere predicates (geometry::robust_in_circumsphere) - the exact Delaunay cavity is connected, void-free, and star-shaped around the point - and in higher dimensions it is pruned to the seed component and shrunk until every hole face is visible from the point. If even the repaired cavity fails validation, the insertion is rejected with the triangulation untouched (same AssertionError as before, or ValueError when the point cannot be connected), and add_point rolls back its vertex on these rejections, so a failed insertion is fully atomic and callers can skip the point and continue. Well-conditioned insertions never enter the repair path and behave bit-identically to the reference (the 64 cross-validation tests pass unchanged). Stress sweeps over 50 seeds x {mixed-scale, off-origin, near-duplicate, random} x {2D,3D}: corrupted end states drop from 99/400 runs on main to 0/400, with every former corruption now either repaired (e.g. all 22 previously-failing insertions in the hull-heavy 3D benchmark now succeed) or atomically rejected. A/B benchmark shows no performance change on well-conditioned workloads. Policy documented in src/tolerances.rs. --- src/tolerances.rs | 33 ++-- src/triangulation.rs | 348 +++++++++++++++++++++++++++++++++++-------- 2 files changed, 313 insertions(+), 68 deletions(-) diff --git a/src/tolerances.rs b/src/tolerances.rs index 50b9e7b..961be75 100644 --- a/src/tolerances.rs +++ b/src/tolerances.rs @@ -12,19 +12,34 @@ //! simplices far from the origin get misclassified — see the regression tests //! added for exactly that failure in PR #3. //! -//! # Known limitation +//! # Degenerate-input policy (deliberate divergence from the reference) //! //! Point sets that mix widely separated coordinate scales in one //! triangulation (e.g. a unit-sized cloud plus a 1e-5-sized cluster 100 //! away) force sliver simplices with aspect ratios around 1e7 into the mesh. -//! Floating-point circumsphere tests on such slivers are unreliable, and the -//! Bowyer-Watson cascade can then assemble an inconsistent cavity; this is -//! caught loudly by the volume-conservation assertion rather than corrupting -//! state silently. The Python reference fails the same way (more often, in -//! fact). Fixing it would require exact geometric predicates, i.e. -//! adaptive-precision orientation/insphere tests, which is out of scope for -//! a tolerance constant. Homogeneous-scale inputs — at any magnitude and any -//! offset from the origin — are unaffected; see +//! Floating-point circumsphere tests on such slivers are unreliable, so the +//! Bowyer-Watson cascade can assemble a cavity that is not re-triangulable +//! (it has voids or protrusions, misses the simplex containing the point, +//! or would disconnect the point). The Python reference detects part of +//! this with a volume-conservation `assert` only *after* mutating, so a +//! failing insertion corrupts its state; failures it cannot see (cavities +//! whose total volume sits below the check's absolute tolerance) silently +//! orphan vertices. +//! +//! This implementation instead validates the cavity *before* mutating +//! (volume conservation, computed with adaptive-precision determinants, and +//! point connectivity). When validation fails it repairs the cavity — in 2D +//! and 3D by rebuilding it with Shewchuk's exact insphere predicates (the +//! `robust` crate; see `geometry::robust_in_circumsphere`), in higher +//! dimensions by shrinking it until star-shaped — and re-validates. If even +//! the repaired cavity fails, the insertion is rejected (the same +//! `AssertionError`, or `ValueError` for a point that cannot be connected) +//! with the triangulation untouched, so callers can skip the point and +//! continue. Well-conditioned insertions never enter the repair path and +//! behave bit-identically to the reference. +//! +//! Homogeneous-scale inputs — at any magnitude and any offset from the +//! origin — never trigger any of this; see //! `test_random_insertions_small_scale_far_from_origin`. /// Relative slack on barycentric coordinates when deciding whether a point diff --git a/src/triangulation.rs b/src/triangulation.rs index 68cecc8..9f8e06e 100644 --- a/src/triangulation.rs +++ b/src/triangulation.rs @@ -953,6 +953,7 @@ impl Triangulation { } } } + let seed_ids: Vec = queue.iter().copied().collect(); while let Some(id) = queue.pop_front() { let simplex = self.simplex_by_id(id); @@ -980,29 +981,58 @@ impl Triangulation { } } - let bad_simplices: Vec = bad_ids - .iter() - .map(|&id| self.simplex_by_id(id).clone()) - .collect(); - - let hole_faces: Vec = Self::facet_multiplicities(bad_simplices.iter()) - .into_iter() - .filter_map(|(face, count)| (count == 1).then_some(face)) - .collect(); - - let mut candidates: Vec = Vec::new(); - for face in hole_faces { - if face.contains(&pt_index) { - continue; + let (mut bad_simplices, mut candidates) = self.cavity_candidates(pt_index, &bad_ids)?; + let (mut deleted_simplices, mut new_simplices, mut point_connected) = + self.cavity_outcome(pt_index, &bad_simplices, &candidates); + + let mut old_vol = self.summed_volume(&deleted_simplices)?; + let mut new_vol = self.summed_volume(&new_simplices)?; + // The cavity is unusable when re-triangulating it would not cover + // the volume it frees, or when it would leave the point connected + // to nothing: either an empty cavity (the circumsphere test falsely + // excluded even the simplex the point lies in) or a cavity that + // swallows every simplex of the point while producing no candidates + // (its total volume can sit below the conservation check's absolute + // tolerance, making that check vacuous). Both happen when sliver + // simplices feed the floating-point circumsphere test cancellation + // noise; see src/tolerances.rs. + if !is_close(old_vol, new_vol) || !point_connected { + // Repair: in 2D/3D rebuild the cavity with exact insphere + // predicates — the exact Delaunay cavity is connected, + // void-free, and star-shaped around the point, so filling it + // conserves volume by construction (provided the surrounding + // mesh itself is still consistent). Where no exact predicate + // exists (4D+), fall back to shrinking the cavity until every + // hole face is visible from the point. Nothing has been mutated + // yet, so when even the repaired cavity fails validation the + // insertion is rejected with the triangulation untouched; the + // Python reference instead mutates first and corrupts its + // state. + if let Some(exact_cavity) = + self.rebuild_cavity_exact(pt_index, &seed_ids, transform.as_deref())? + { + bad_ids = exact_cavity; + } else { + self.shrink_cavity_to_star(pt_index, &mut bad_ids, &seed_ids)?; } - let mut simplex = face; - simplex.push(pt_index); - simplex.sort_unstable(); - - if self.simplex_is_numerically_degenerate(&simplex)? { - continue; + (bad_simplices, candidates) = self.cavity_candidates(pt_index, &bad_ids)?; + (deleted_simplices, new_simplices, point_connected) = + self.cavity_outcome(pt_index, &bad_simplices, &candidates); + old_vol = self.summed_volume(&deleted_simplices)?; + new_vol = self.summed_volume(&new_simplices)?; + + if !is_close(old_vol, new_vol) { + return Err(TriangulationError::Assertion(format!( + "{old_vol} !== {new_vol}" + ))); + } + if !point_connected { + // Even the repaired cavity cannot connect the point: it is + // numerically indistinguishable from existing structure. + return Err(TriangulationError::Value( + "Point already in triangulation.".to_string(), + )); } - candidates.push(simplex); } // A cavity vertex that would lose all of its simplices without @@ -1032,29 +1062,227 @@ impl Triangulation { self.add_simplex(simplex)?; } - let new_triangles: FxHashSet = self.vertex_to_ids[pt_index] + Ok((deleted_simplices, new_simplices)) + } + + /// The cavity simplices (cloned out of the slab) and the candidate + /// simplices that would fill the cavity's hole: one per hole face (a + /// facet belonging to exactly one cavity simplex) that does not contain + /// the point, joined with the point, minus numerically degenerate ones. + fn cavity_candidates( + &self, + pt_index: usize, + bad_ids: &FxHashSet, + ) -> Result<(Vec, Vec), TriangulationError> { + let bad_simplices: Vec = bad_ids .iter() .map(|&id| self.simplex_by_id(id).clone()) .collect(); - let bad_triangles: FxHashSet = bad_simplices.into_iter().collect(); - let deleted_simplices: FxHashSet = - bad_triangles.difference(&new_triangles).cloned().collect(); - let new_simplices: FxHashSet = - new_triangles.difference(&bad_triangles).cloned().collect(); - - let old_vol = deleted_simplices.iter().try_fold(0.0, |acc, simplex| { - Ok::(acc + self.volume(simplex)?) - })?; - let new_vol = new_simplices.iter().try_fold(0.0, |acc, simplex| { - Ok::(acc + self.volume(simplex)?) - })?; - if !is_close(old_vol, new_vol) { - return Err(TriangulationError::Assertion(format!( - "{old_vol} !== {new_vol}" - ))); + + let hole_faces: Vec = Self::facet_multiplicities(bad_simplices.iter()) + .into_iter() + .filter_map(|(face, count)| (count == 1).then_some(face)) + .collect(); + + let mut candidates: Vec = Vec::new(); + for face in hole_faces { + if face.contains(&pt_index) { + continue; + } + let mut simplex = face; + simplex.push(pt_index); + simplex.sort_unstable(); + + if self.simplex_is_numerically_degenerate(&simplex)? { + continue; + } + candidates.push(simplex); } + Ok((bad_simplices, candidates)) + } - Ok((deleted_simplices, new_simplices)) + /// The `(deleted, created, point_connected)` outcome that deleting the + /// cavity and adding the candidates *would* produce, computed without + /// mutating anything: the point's simplices afterwards are its current + /// ones minus the cavity plus the candidates (`point_connected` reports + /// whether that set is non-empty), and a simplex both deleted and + /// recreated cancels out of the result. + fn cavity_outcome( + &self, + pt_index: usize, + bad_simplices: &[Simplex], + candidates: &[Simplex], + ) -> (FxHashSet, FxHashSet, bool) { + let bad_set: FxHashSet<&Simplex> = bad_simplices.iter().collect(); + let mut new_triangles: FxHashSet = self.vertex_to_ids[pt_index] + .iter() + .map(|&id| self.simplex_by_id(id)) + .filter(|simplex| !bad_set.contains(simplex)) + .cloned() + .collect(); + new_triangles.extend(candidates.iter().cloned()); + let point_connected = !new_triangles.is_empty(); + + let deleted: FxHashSet = bad_simplices + .iter() + .filter(|simplex| !new_triangles.contains(*simplex)) + .cloned() + .collect(); + let created: FxHashSet = new_triangles + .into_iter() + .filter(|simplex| !bad_set.contains(simplex)) + .collect(); + (deleted, created, point_connected) + } + + /// Total volume of `simplices`, using adaptive-precision determinants + /// (in 2D/3D) so the conservation check compares true volumes instead of + /// cancellation noise when the cavity contains slivers. + fn summed_volume(&self, simplices: &FxHashSet) -> Result { + simplices.iter().try_fold(0.0, |acc, simplex| { + Ok::( + acc + geometry::precise_volume(&self.get_vertices(simplex)?)?, + ) + }) + } + + /// Rebuild the cavity with exact in-circumsphere predicates (2D/3D + /// only; `None` elsewhere): cascade through facet neighbours from the + /// seeds, keeping the seeds themselves (the point lies in or on them) + /// plus every simplex whose circumsphere strictly contains the point. + /// The exact Delaunay cavity is connected, void-free, and star-shaped + /// around the point, so filling it from the point conserves volume even + /// when slivers fooled the floating-point tests that assembled the + /// original cavity. + fn rebuild_cavity_exact( + &self, + pt_index: usize, + seed_ids: &[SimplexId], + transform: Option<&[Vec]>, + ) -> Result>, TriangulationError> { + if self.dim != 2 && self.dim != 3 { + return Ok(None); + } + let point = match transform { + Some(matrix) => apply_transform(&self.vertices[pt_index], matrix), + None => self.vertices[pt_index].clone(), + }; + + let mut cavity: FxHashSet = FxHashSet::default(); + let mut examined: FxHashSet = FxHashSet::default(); + let mut queue: VecDeque = VecDeque::new(); + for &seed in seed_ids { + if self.id_is_live(seed) && examined.insert(seed) { + cavity.insert(seed); + queue.push_back(seed); + } + } + + while let Some(id) = queue.pop_front() { + let simplex = self.simplex_by_id(id); + for skip in 0..simplex.len() { + let Some(incident) = self.facets.get(&facet_excluding(simplex, skip)) else { + continue; + }; + for &neighbour in incident { + if neighbour == id || !examined.insert(neighbour) { + continue; + } + let points = self.simplex_points(self.simplex_by_id(neighbour), transform)?; + if geometry::robust_in_circumsphere(&points, &point)?.unwrap_or(false) { + cavity.insert(neighbour); + queue.push_back(neighbour); + } + } + } + } + Ok(Some(cavity)) + } + + /// Shrink the cavity until it is star-shaped as seen from `pt_index`: + /// every hole face not containing the point must be strictly visible + /// from it, i.e. the point lies on the same side of the face as the + /// cavity simplex owning it. Owners of invisible faces are evicted and + /// the hole boundary recomputed until no violation remains. Visibility + /// uses exact orientation predicates in 2D/3D + /// ([`geometry::robust_orientation`]), so a sliver cannot fool it. + fn shrink_cavity_to_star( + &self, + pt_index: usize, + bad_ids: &mut FxHashSet, + seed_ids: &[SimplexId], + ) -> Result<(), TriangulationError> { + let point = &self.vertices[pt_index]; + self.prune_to_seed_component(bad_ids, seed_ids); + loop { + // Hole face -> (multiplicity within the cavity, owning cavity + // simplex, its vertex opposite the face). + let mut faces: FxHashMap = FxHashMap::default(); + for &id in bad_ids.iter() { + let simplex = self.simplex_by_id(id); + for skip in 0..simplex.len() { + faces + .entry(facet_excluding(simplex, skip)) + .and_modify(|entry| entry.0 += 1) + .or_insert((1, id, simplex[skip])); + } + } + + let mut evicted = false; + for (face, &(count, owner, opposite)) in &faces { + if count != 1 || face.contains(&pt_index) || !bad_ids.contains(&owner) { + continue; + } + let face_points = self.get_vertices(face)?; + let side_of_point = geometry::robust_orientation(&face_points, point)?; + let side_of_cavity = + geometry::robust_orientation(&face_points, &self.vertices[opposite])?; + if side_of_point == 0 || side_of_point != side_of_cavity { + bad_ids.remove(&owner); + evicted = true; + } + } + if !evicted || bad_ids.is_empty() { + return Ok(()); + } + // Evictions can split the cavity; fragments no longer reachable + // from the seeds cannot be filled from the point. + self.prune_to_seed_component(bad_ids, seed_ids); + if bad_ids.is_empty() { + return Ok(()); + } + } + } + + /// Restrict the cavity to the component reachable from the seed + /// simplices through shared facets. Disconnected fragments are + /// circumsphere false positives on slivers; keeping them makes the hole + /// unfillable from the point. + fn prune_to_seed_component(&self, bad_ids: &mut FxHashSet, seed_ids: &[SimplexId]) { + let mut reachable: FxHashSet = FxHashSet::default(); + let mut queue: VecDeque = VecDeque::new(); + for &seed in seed_ids { + if bad_ids.contains(&seed) && reachable.insert(seed) { + queue.push_back(seed); + } + } + while let Some(id) = queue.pop_front() { + let simplex = self.simplex_by_id(id); + for skip in 0..simplex.len() { + let Some(incident) = self.facets.get(&facet_excluding(simplex, skip)) else { + continue; + }; + for &neighbour in incident { + if neighbour != id + && bad_ids.contains(&neighbour) + && reachable.insert(neighbour) + { + queue.push_back(neighbour); + } + } + } + } + *bad_ids = reachable; } /// Connect vertex `pt_index`, which must lie outside the convex hull, to @@ -1146,24 +1374,26 @@ impl Triangulation { return Err(err); } }; - let (deleted_simplices, added_simplices) = - match self.bowyer_watson(pt_index, None, &transform) { - Ok(result) => result, - // Value errors are raised before bowyer_watson mutates, - // so only the hull extension needs rolling back. - Err(err @ TriangulationError::Value(_)) => { - let attached: Vec = self.vertex_to_ids[pt_index] - .iter() - .map(|&id| self.simplex_by_id(id).clone()) - .collect(); - for simplex in attached { - self.delete_simplex(&simplex)?; - } - self.pop_vertex(); - return Err(err); + let (deleted_simplices, added_simplices) = match self + .bowyer_watson(pt_index, None, &transform) + { + Ok(result) => result, + // Value and assertion errors are raised before + // bowyer_watson mutates, so only the hull extension + // needs rolling back. + Err(err @ (TriangulationError::Value(_) | TriangulationError::Assertion(_))) => { + let attached: Vec = self.vertex_to_ids[pt_index] + .iter() + .map(|&id| self.simplex_by_id(id).clone()) + .collect(); + for simplex in attached { + self.delete_simplex(&simplex)?; } - Err(err) => return Err(err), - }; + self.pop_vertex(); + return Err(err); + } + Err(err) => return Err(err), + }; let deleted: FxHashSet = deleted_simplices .difference(&temporary_simplices) @@ -1192,9 +1422,9 @@ impl Triangulation { let pt_index = self.add_vertex(point); match self.bowyer_watson(pt_index, Some(actual_simplex), &transform) { - // Value errors are raised before bowyer_watson mutates, so only - // the freshly pushed vertex needs rolling back. - Err(err @ TriangulationError::Value(_)) => { + // Value and assertion errors are raised before bowyer_watson + // mutates, so only the freshly pushed vertex needs rolling back. + Err(err @ (TriangulationError::Value(_) | TriangulationError::Assertion(_))) => { self.pop_vertex(); Err(err) } From 2adb5c187f98c2e4495fb7412fccb5f0749e6e36 Mon Sep 17 00:00:00 2001 From: Bas Nijholt Date: Wed, 10 Jun 2026 11:30:10 -0700 Subject: [PATCH 3/3] test: pin the degenerate-input robustness contract Seeded mixed-scale and near-duplicate sweeps (2D/3D, 10 seeds each) assert the contract from src/tolerances.rs: insertions either succeed and connect the new vertex, or raise leaving vertices and simplices exactly as they were; no vertex ever loses its last simplex; the internal indexes stay consistent. Plus a regression test for the empty-cavity orphan (mixed-scale seed 2) and a guard that random unit-cube sweeps never reject (the repair path must not fire on well-conditioned input). --- tests/test_robustness.py | 109 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 109 insertions(+) create mode 100644 tests/test_robustness.py diff --git a/tests/test_robustness.py b/tests/test_robustness.py new file mode 100644 index 0000000..e86a2eb --- /dev/null +++ b/tests/test_robustness.py @@ -0,0 +1,109 @@ +"""Degenerate-input robustness: failed insertions must never corrupt state. + +The Python reference corrupts itself on these inputs (it mutates before its +volume-conservation assert, and orphans vertices when the assert cannot see +the problem), so unlike test_triangulation.py nothing here cross-validates +against it. The contract pinned here is documented in src/tolerances.rs: + +- every insertion either succeeds, or raises (ValueError/AssertionError) + leaving the triangulation exactly as it was; +- a successful insertion always connects the new vertex; +- vertices never lose their last simplex; +- the internal indexes stay consistent (reference_invariant). +""" + +from __future__ import annotations + +import adaptive_triangulation as at +import numpy as np +import pytest + + +def mixed_scale_points(seed: int, dim: int) -> np.ndarray: + """A unit-sized cloud plus a tiny far-away cluster: forces ~1e7 aspect + ratio slivers, the worst known case for floating-point predicates.""" + rng = np.random.default_rng(seed) + big = rng.random((20, dim)) + tiny = 100.0 + 1e-5 * rng.random((20, dim)) + pts = np.vstack([big, tiny]) + rng.shuffle(pts) + return pts + + +def near_duplicate_points(seed: int, dim: int) -> np.ndarray: + rng = np.random.default_rng(seed) + base = rng.random((15, dim)) + dup = base[rng.integers(0, 15, 15)] + 1e-12 * rng.random((15, dim)) + pts = np.vstack([base, dup]) + rng.shuffle(pts) + return pts + + +def simplex_count_per_vertex(tri) -> list[int]: + return [len(tri.vertex_to_simplices[i]) for i in range(len(tri.vertices))] + + +@pytest.mark.parametrize("dim", [2, 3]) +@pytest.mark.parametrize("seed", range(10)) +@pytest.mark.parametrize("generate", [mixed_scale_points, near_duplicate_points]) +def test_degenerate_insertions_never_corrupt_state(generate, seed, dim): + pts = generate(seed, dim) + tri = at.Triangulation(pts[: dim + 2].tolist()) + # scipy/QHull may merge near-degenerate construction points (the + # reference does the same); those stay unconnected by design. + construction_orphans = { + i for i in range(len(tri.vertices)) if len(tri.vertex_to_simplices[i]) == 0 + } + + for p in pts[dim + 2 :]: + n_vertices = len(tri.vertices) + simplices_before = {tuple(s) for s in tri.simplices} + try: + tri.add_point(p.tolist()) + except (ValueError, AssertionError): + # Rejections must be atomic: no vertex leaked, no simplex changed. + assert len(tri.vertices) == n_vertices + assert {tuple(s) for s in tri.simplices} == simplices_before + else: + # Successful insertions must connect the new vertex. + assert len(tri.vertices) == n_vertices + 1 + assert len(tri.vertex_to_simplices[n_vertices]) > 0 + + assert tri.reference_invariant() + counts = simplex_count_per_vertex(tri) + orphans = {i for i, count in enumerate(counts) if count == 0} + assert orphans == construction_orphans + + +def test_empty_cavity_insertion_is_repaired_not_orphaned(): + # Regression (mixed-scale seed 2): the point lies inside a huge bridge + # sliver whose floating-point circumsphere test reports it outside, so + # the cavity came back empty and add_point "succeeded" with (set(), set()), + # leaving the new vertex unconnected. The exact-predicate rebuild must + # either connect it or reject the insertion outright. + pts = mixed_scale_points(2, 2) + dim = 2 + tri = at.Triangulation(pts[: dim + 2].tolist()) + for p in pts[dim + 2 :]: + n_vertices = len(tri.vertices) + try: + deleted, added = tri.add_point(p.tolist()) + except (ValueError, AssertionError): + assert len(tri.vertices) == n_vertices + continue + assert added, "successful insertion must create at least one simplex" + assert len(tri.vertex_to_simplices[len(tri.vertices) - 1]) > 0 + assert tri.reference_invariant() + + +def test_well_conditioned_insertions_never_reject(): + # The repair path must not change behavior for ordinary inputs: random + # unit-cube sweeps go through the reference-compatible fast path only. + for dim in (2, 3): + rng = np.random.default_rng(99) + pts = rng.random((60, dim)) + tri = at.Triangulation(pts[: dim + 2].tolist()) + for p in pts[dim + 2 :]: + tri.add_point(p.tolist()) # must not raise + assert tri.reference_invariant() + assert all(len(tri.vertex_to_simplices[i]) > 0 for i in range(len(tri.vertices)))