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