diff --git a/src/geometry.rs b/src/geometry.rs index fca6cba..9cc5d98 100644 --- a/src/geometry.rs +++ b/src/geometry.rs @@ -424,10 +424,12 @@ pub fn simplex_volume_in_embedding(vertices: &[Vec]) -> Result]) -> Result, + coords: &[Vec], + dim: usize, +) -> PyResult>> { + if dim == 1 { + return Ok(None); + } + + let Ok(spatial) = PyModule::import(py, "scipy.spatial") else { + return Ok(None); + }; + let coords_array = PyArray2::from_vec2(py, coords)?; + let Ok(delaunay) = spatial.getattr("Delaunay")?.call1((coords_array,)) else { + return Ok(None); + }; + + let simplices = delaunay.getattr("simplices")?; + let mut initial = Vec::new(); + for simplex in simplices.try_iter()? { + let simplex = simplex?; + let mut indices = Vec::new(); + for item in simplex.try_iter()? { + indices.push(item?.extract::()?); + } + indices.sort_unstable(); + initial.push(indices); + } + Ok(Some(initial)) +} + fn is_close(a: f64, b: f64) -> bool { - (a - b).abs() <= 1e-8 + 1e-5 * b.abs() + let scale = a.abs().max(b.abs()); + (a - b).abs() <= DEFAULT_EPS + 1e-5 * scale } #[derive(Debug)] @@ -383,11 +415,6 @@ impl Triangulation { "Coordinates dimension mismatch".to_string(), )); } - if dim == 1 { - return Err(TriangulationError::Value( - "Triangulation class only supports dim >= 2".to_string(), - )); - } if coords.len() < dim + 1 { return Err(TriangulationError::Value( "Please provide at least one simplex".to_string(), @@ -908,7 +935,7 @@ impl Triangulation { simplex.push(pt_index); simplex.sort_unstable(); - if self.volume(&simplex)? < 1e-8 { + if self.simplex_is_numerically_degenerate(&simplex)? { continue; } self.add_simplex(simplex)?; @@ -1062,6 +1089,40 @@ impl Triangulation { Ok(geometry::volume(&self.get_vertices(simplex)?)?) } + fn normalized_volume(&self, simplex: &[usize]) -> Result { + let vertices = self.get_vertices(simplex)?; + let base = &vertices[0]; + let mut total_abs_coordinate_delta = 0.0; + let mut delta_count = 0usize; + for vertex in vertices.iter().skip(1) { + for (coord, origin) in vertex.iter().zip(base) { + total_abs_coordinate_delta += (coord - origin).abs(); + delta_count += 1; + } + } + + if delta_count == 0 { + return Ok(0.0); + } + let characteristic_length = total_abs_coordinate_delta / delta_count as f64; + if characteristic_length == 0.0 { + return Ok(0.0); + } + + Ok(self.volume(simplex)? / characteristic_length.powi(self.dim as i32)) + } + + fn simplex_is_numerically_degenerate( + &self, + simplex: &[usize], + ) -> Result { + if self.dim == 1 { + // In 1D we only want to reject coincident endpoints, not tiny but valid intervals. + return Ok(self.normalized_volume(simplex)? < DEFAULT_EPS); + } + Ok(self.volume(simplex)? < DEFAULT_EPS) + } + pub fn has_simplex(&self, simplex: &[usize]) -> Result { let mut simplex = simplex.to_vec(); simplex.sort_unstable(); @@ -1350,30 +1411,12 @@ impl PyTriangulation { fn new(py: Python<'_>, coords: &Bound<'_, PyAny>) -> PyResult { let parsed_coords = parse_points_sized(coords, "Please provide a 2-dimensional list of points")?; - Triangulation::validate_coords(&parsed_coords).map_err(TriangulationError::into_pyerr)?; - - let core = match PyModule::import(py, "scipy.spatial") { - Ok(spatial) => { - let coords_array = PyArray2::from_vec2(py, &parsed_coords)?; - match spatial.getattr("Delaunay")?.call1((coords_array,)) { - Ok(delaunay) => { - let simplices = delaunay.getattr("simplices")?; - let mut initial = Vec::new(); - for simplex in simplices.try_iter()? { - let simplex = simplex?; - let mut indices = Vec::new(); - for item in simplex.try_iter()? { - indices.push(item?.extract::()?); - } - indices.sort_unstable(); - initial.push(indices); - } - Triangulation::from_simplices(parsed_coords.clone(), initial) - } - Err(_) => Triangulation::new(parsed_coords.clone()), - } - } - Err(_) => Triangulation::new(parsed_coords.clone()), + let dim = Triangulation::validate_coords(&parsed_coords) + .map_err(TriangulationError::into_pyerr)?; + + let core = match scipy_delaunay_simplices(py, &parsed_coords, dim)? { + Some(initial) => Triangulation::from_simplices(parsed_coords.clone(), initial), + None => Triangulation::new(parsed_coords.clone()), } .map_err(TriangulationError::into_pyerr)?; Ok(Self { core }) @@ -1736,4 +1779,30 @@ mod tests { assert!(simplices.contains(&vec![0, 1, 3])); assert!(simplices.contains(&vec![0, 2, 3])); } + + #[test] + fn one_dimensional_triangulation_connects_adjacent_points() { + let tri = Triangulation::new(vec![vec![2.0], vec![0.0], vec![1.0], vec![3.0]]).unwrap(); + + assert_eq!(tri.dim, 1); + assert_eq!( + tri.simplices, + FxHashSet::from_iter([vec![0, 2], vec![1, 2], vec![0, 3]]) + ); + assert_eq!(tri.hull().unwrap(), FxHashSet::from_iter([1, 3])); + assert!(tri.reference_invariant()); + } + + #[test] + fn one_dimensional_tiny_interval_survives_bowyer_watson() { + let mut tri = Triangulation::new(vec![vec![0.0], vec![1e-12]]).unwrap(); + let (deleted, added) = tri.add_point(vec![5e-13], None, None).unwrap(); + + assert_eq!(deleted, FxHashSet::from_iter([vec![0, 1]])); + assert_eq!(added, FxHashSet::from_iter([vec![0, 2], vec![1, 2]])); + assert_eq!( + tri.simplices, + FxHashSet::from_iter([vec![0, 2], vec![1, 2]]) + ); + } } diff --git a/tests/test_triangulation.py b/tests/test_triangulation.py index 3cea919..aca0b15 100644 --- a/tests/test_triangulation.py +++ b/tests/test_triangulation.py @@ -1,5 +1,6 @@ from __future__ import annotations +import itertools import math from collections import Counter @@ -27,6 +28,21 @@ def face_counter(iterator) -> Counter[tuple[int, ...]]: return Counter(tuple(face) for face in iterator) +def expected_1d_order(vertices) -> np.ndarray: + coords = np.asarray(vertices, dtype=float).reshape(-1) + return np.argsort(coords, kind="mergesort") + + +def expected_1d_simplices(vertices) -> set[tuple[int, int]]: + order = expected_1d_order(vertices) + return {tuple(sorted((int(left), int(right)))) for left, right in itertools.pairwise(order)} + + +def expected_1d_hull(vertices) -> set[int]: + order = expected_1d_order(vertices) + return {int(order[0]), int(order[-1])} + + def assert_points_close(lhs, rhs, atol: float = 1e-8) -> None: lhs = np.asarray(lhs, dtype=float) rhs = np.asarray(rhs, dtype=float) @@ -55,6 +71,20 @@ def assert_same_exception_type_name(rust_callable, reference_callable) -> None: assert type(rust_exc.value).__name__ == type(ref_exc.value).__name__ +def assert_1d_triangulation_state(tri) -> None: + assert as_simplex_set(tri.simplices) == expected_1d_simplices(tri.vertices) + assert set(tri.hull) == expected_1d_hull(tri.vertices) + assert tri.reference_invariant() is True + + +def assert_1d_midpoints_locate_adjacent_segments(tri) -> None: + vertices = np.asarray(tri.vertices, dtype=float).reshape(-1) + order = expected_1d_order(vertices) + for left, right in itertools.pairwise(order): + midpoint = 0.5 * (vertices[left] + vertices[right]) + assert locate_result(tri.locate_point([midpoint])) == tuple(sorted((int(left), int(right)))) + + class Seq: def __init__(self, items): self._items = items @@ -171,6 +201,45 @@ def test_geometry_functions_match_reference(): assert fast3_alias_radius == pytest.approx(ref_fast3_radius) +def test_geometry_functions_support_1d_segments(): + segment = np.array([[0.0], [2.0]]) + embedded_segment = np.array([[0.0, 0.0, 0.0], [3.0, 4.0, 0.0]]) + + center, radius = rust_tri.circumsphere(segment) + assert_points_close(center, [1.0]) + assert radius == pytest.approx(1.0) + + assert rust_tri.point_in_simplex([0.5], segment) is True + assert rust_tri.point_in_simplex([2.5], segment) is False + assert rust_tri.volume(segment) == pytest.approx(2.0) + assert rust_tri.simplex_volume_in_embedding(embedded_segment) == pytest.approx(5.0) + + with pytest.raises(ValueError, match="Provided vertices do not form a simplex"): + rust_tri.simplex_volume_in_embedding([[1.0, 2.0], [1.0, 2.0]]) + + +def test_circumscribed_circle_supports_1d(): + tri = rust_tri.Triangulation([[0.0], [2.0], [1.0]]) + + center, radius = tri.circumscribed_circle((0, 1)) + assert_points_close(center, [1.0]) + assert radius == pytest.approx(1.0) + assert tri.point_in_circumcircle(2, (0, 1)) is True + + scaled_center, scaled_radius = tri.circumscribed_circle((0, 1), transform=np.array([[2.0]])) + assert_points_close(scaled_center, [2.0]) + assert scaled_radius == pytest.approx(2.0) + + +def test_construction_1d_creates_adjacent_segments(): + coords = np.array([[2.0], [0.0], [1.0], [3.0]]) + tri = rust_tri.Triangulation(coords) + + assert tri.dim == 1 + assert_points_close(tri.vertices, coords) + assert_1d_triangulation_state(tri) + + def test_faces_containing_and_hull_match_reference(): coords = [[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0], [0.4, 0.2]] rust = rust_tri.Triangulation(coords) @@ -188,6 +257,15 @@ def test_faces_containing_and_hull_match_reference(): assert set(rust.hull) == set(reference.hull) +def test_faces_containing_and_hull_work_in_1d(): + tri = rust_tri.Triangulation([[0.0], [1.0], [2.0]]) + + assert face_counter(tri.faces()) == Counter({(0,): 1, (1,): 2, (2,): 1}) + assert face_counter(tri.faces(dim=1)) == Counter({(0,): 1, (1,): 2, (2,): 1}) + assert as_simplex_set(tri.containing((1,))) == {(0, 1), (1, 2)} + assert_1d_triangulation_state(tri) + + def test_add_point_inside_hull_matches_reference(): coords = [[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]] point = [0.3, 0.4] @@ -216,6 +294,51 @@ def test_add_point_outside_hull_matches_reference(): assert_triangulation_equal(rust, reference) +@pytest.mark.parametrize( + ("coords", "point", "expected_deleted", "expected_added", "expected_volumes"), + [ + ( + [[0.0], [2.0]], + [1.0], + {(0, 1)}, + {(0, 2), (1, 2)}, + {}, + ), + ( + [[0.0], [1.0]], + [2.0], + set(), + {(1, 2)}, + {}, + ), + ( + [[0.0], [1e-12]], + [5e-13], + {(0, 1)}, + {(0, 2), (1, 2)}, + {(0, 2): 5e-13, (1, 2): 5e-13}, + ), + ], + ids=["inside-hull", "outside-hull", "tiny-interval"], +) +def test_add_point_works_in_1d( + coords, + point, + expected_deleted, + expected_added, + expected_volumes, +): + tri = rust_tri.Triangulation(coords) + + deleted, added = tri.add_point(point) + + assert as_simplex_set(deleted) == expected_deleted + assert as_simplex_set(added) == expected_added + assert_1d_triangulation_state(tri) + for simplex, expected_volume in expected_volumes.items(): + assert tri.volume(simplex) == pytest.approx(expected_volume) + + def test_add_point_with_transform_matches_reference(): coords = [[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]] point = [0.4, 0.25] @@ -274,6 +397,13 @@ def test_duplicate_point_is_rejected(): tri.add_point([0.0, 0.0], simplex=simplex) +def test_duplicate_point_is_rejected_in_1d(): + tri = rust_tri.Triangulation([[0.0], [1.0]]) + + with pytest.raises(ValueError, match="Point already in triangulation"): + tri.add_point([1.0]) + + def test_random_cross_validation_2d(): rng = np.random.default_rng(1234) coords = rng.random((6, 2)) @@ -526,15 +656,10 @@ def test_tiny_triangle_point_in_simplex_matches_reference(): def test_simplex_volume_in_embedding_matches_reference_edge_case(): - assert_same_exception_type_name( - lambda: rust_tri.simplex_volume_in_embedding([[0.0, 0.0], [1.0, 0.0]]), - lambda: reference_module.simplex_volume_in_embedding([[0.0, 0.0], [1.0, 0.0]]), - ) + assert rust_tri.simplex_volume_in_embedding([[0.0, 0.0], [1.0, 0.0]]) == pytest.approx(1.0) assert rust_tri.simplex_volume_in_embedding( [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]] - ) == pytest.approx( - reference_module.simplex_volume_in_embedding([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]]) - ) + ) == pytest.approx(1.0) def test_orientation_matches_reference_at_large_scale(): @@ -626,6 +751,18 @@ def test_random_cross_validation_4d(): assert_triangulation_equal(rust, reference) +def test_random_cross_validation_1d(): + rng = np.random.default_rng(8642) + values = rng.random(10) + np.arange(10) * 1e-6 + coords = values.reshape(-1, 1)[rng.permutation(10)] + tri = rust_tri.Triangulation(coords[:2]) + + for point in coords[2:]: + tri.add_point(point) + assert_1d_triangulation_state(tri) + assert_1d_midpoints_locate_adjacent_segments(tri) + + def test_public_bowyer_watson_is_exposed(): tri = rust_tri.Triangulation([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]]) assert hasattr(tri, "bowyer_watson")