From 63176b49649dce2a50ce57ae096970d1cc4c2671 Mon Sep 17 00:00:00 2001 From: Bas Nijholt Date: Wed, 10 Jun 2026 11:57:04 -0700 Subject: [PATCH 1/3] refactor: let ? convert TriangulationError in the PyO3 layer Readability only; no behavior change (the full suite, stress sweeps, and benchmarks are identical). - impl From for PyErr replaces the inherent into_pyerr method, so the 31 call sites that wrapped every core call in .map_err(TriangulationError::into_pyerr) now just use ?; most of the awkward multi-line method chains in the bindings collapse back to single lines - given() collapses the repeated three-arm 'absent argument or explicit Python None' matches to one helper - the eight SimplicesProxy set operators delegate through a shared set_op/set_rop pair instead of repeating the snapshot dance --- src/py.rs | 204 ++++++++++++++++++++---------------------------------- 1 file changed, 75 insertions(+), 129 deletions(-) diff --git a/src/py.rs b/src/py.rs index 1f646f0..b2b743a 100644 --- a/src/py.rs +++ b/src/py.rs @@ -21,14 +21,14 @@ use crate::triangulation::{ index_out_of_range, reduced_simplex_positions, Simplex, Triangulation, TriangulationError, }; -impl TriangulationError { - pub(crate) fn into_pyerr(self) -> PyErr { - match self { - Self::Value(message) => PyValueError::new_err(message), - Self::Index(message) => PyIndexError::new_err(message), - Self::Runtime(message) => PyRuntimeError::new_err(message), - Self::Assertion(message) => PyAssertionError::new_err(message), - Self::Geometry(error) => PyValueError::new_err(error.to_string()), +impl From for PyErr { + fn from(error: TriangulationError) -> PyErr { + match error { + TriangulationError::Value(message) => PyValueError::new_err(message), + TriangulationError::Index(message) => PyIndexError::new_err(message), + TriangulationError::Runtime(message) => PyRuntimeError::new_err(message), + TriangulationError::Assertion(message) => PyAssertionError::new_err(message), + TriangulationError::Geometry(error) => PyValueError::new_err(error.to_string()), } } } @@ -137,12 +137,17 @@ pub(crate) fn parse_signed_index_set(obj: &Bound<'_, PyAny>) -> PyResult(obj: Option<&'a Bound<'py, PyAny>>) -> Option<&'a Bound<'py, PyAny>> { + obj.filter(|value| !value.is_none()) +} + pub(crate) fn parse_optional_transform( obj: Option<&Bound<'_, PyAny>>, ) -> PyResult>>> { - match obj { + match given(obj) { None => Ok(None), - Some(value) if value.is_none() => Ok(None), Some(value) => Ok(Some(parse_points_sized( value, "Expected an N x N transform matrix", @@ -191,25 +196,24 @@ pub(crate) fn normalize_index_set( } fn ordered_indices_from_py(obj: &Bound<'_, PyAny>, len: usize) -> PyResult> { - normalize_indices(&parse_signed_indices(obj)?, len).map_err(TriangulationError::into_pyerr) + normalize_indices(&parse_signed_indices(obj)?, len).map_err(PyErr::from) } fn canonical_simplex_from_py(obj: &Bound<'_, PyAny>, len: usize) -> PyResult { - canonicalize_simplex(&parse_signed_indices(obj)?, len).map_err(TriangulationError::into_pyerr) + canonicalize_simplex(&parse_signed_indices(obj)?, len).map_err(PyErr::from) } fn simplex_set_from_py(obj: &Bound<'_, PyAny>, len: usize) -> PyResult> { let simplices = parse_signed_simplex_set(obj)?; let mut normalized = FxHashSet::default(); for simplex in simplices { - normalized - .insert(normalize_indices(&simplex, len).map_err(TriangulationError::into_pyerr)?); + normalized.insert(normalize_indices(&simplex, len)?); } Ok(normalized) } fn vertex_index_set_from_py(obj: &Bound<'_, PyAny>, len: usize) -> PyResult> { - normalize_index_set(&parse_signed_index_set(obj)?, len).map_err(TriangulationError::into_pyerr) + normalize_index_set(&parse_signed_index_set(obj)?, len).map_err(PyErr::from) } pub(crate) fn point_tuple(py: Python<'_>, point: &[f64]) -> Py { PyTuple::new(py, point.iter().copied()).unwrap().into() @@ -330,6 +334,20 @@ impl PySimplicesProxy { None => simplex_set_py(py, triangulation.core.simplices()), } } + + /// `snapshot other`, delegated to Python's set type. + fn set_op(&self, py: Python<'_>, op: &str, other: &Bound<'_, PyAny>) -> PyResult> { + Ok(self + .snapshot_set(py)? + .bind(py) + .call_method1(op, (other,))? + .unbind()) + } + + /// `other snapshot`, delegated to Python's set type. + fn set_rop(&self, py: Python<'_>, op: &str, other: &Bound<'_, PyAny>) -> PyResult> { + Ok(other.call_method1(op, (self.snapshot_set(py)?,))?.unbind()) + } } /// Lazy, sequence-like view of the vertex coordinates (supports `__array__` @@ -427,59 +445,35 @@ impl PySimplicesProxy { // binary set operators to a snapshot set, in both operand orders. fn __sub__(&self, py: Python<'_>, other: &Bound<'_, PyAny>) -> PyResult> { - Ok(self - .snapshot_set(py)? - .bind(py) - .call_method1("__sub__", (other,))? - .unbind()) + self.set_op(py, "__sub__", other) } fn __rsub__(&self, py: Python<'_>, other: &Bound<'_, PyAny>) -> PyResult> { - Ok(other - .call_method1("__sub__", (self.snapshot_set(py)?,))? - .unbind()) + self.set_rop(py, "__sub__", other) } fn __and__(&self, py: Python<'_>, other: &Bound<'_, PyAny>) -> PyResult> { - Ok(self - .snapshot_set(py)? - .bind(py) - .call_method1("__and__", (other,))? - .unbind()) + self.set_op(py, "__and__", other) } fn __rand__(&self, py: Python<'_>, other: &Bound<'_, PyAny>) -> PyResult> { - Ok(other - .call_method1("__and__", (self.snapshot_set(py)?,))? - .unbind()) + self.set_rop(py, "__and__", other) } fn __or__(&self, py: Python<'_>, other: &Bound<'_, PyAny>) -> PyResult> { - Ok(self - .snapshot_set(py)? - .bind(py) - .call_method1("__or__", (other,))? - .unbind()) + self.set_op(py, "__or__", other) } fn __ror__(&self, py: Python<'_>, other: &Bound<'_, PyAny>) -> PyResult> { - Ok(other - .call_method1("__or__", (self.snapshot_set(py)?,))? - .unbind()) + self.set_rop(py, "__or__", other) } fn __xor__(&self, py: Python<'_>, other: &Bound<'_, PyAny>) -> PyResult> { - Ok(self - .snapshot_set(py)? - .bind(py) - .call_method1("__xor__", (other,))? - .unbind()) + self.set_op(py, "__xor__", other) } fn __rxor__(&self, py: Python<'_>, other: &Bound<'_, PyAny>) -> PyResult> { - Ok(other - .call_method1("__xor__", (self.snapshot_set(py)?,))? - .unbind()) + self.set_rop(py, "__xor__", other) } } @@ -518,8 +512,7 @@ impl PyVertexToSimplicesIter { impl PyVerticesProxy { fn __getitem__(&self, py: Python<'_>, index: isize) -> PyResult> { let triangulation = self.triangulation.bind(py).borrow(); - let index = normalize_index(index, triangulation.core.vertices.len()) - .map_err(TriangulationError::into_pyerr)?; + let index = normalize_index(index, triangulation.core.vertices.len())?; Ok(point_tuple(py, &triangulation.core.vertices[index])) } @@ -558,8 +551,7 @@ impl PyVerticesProxy { impl PyVertexToSimplicesProxy { fn __getitem__(&self, py: Python<'_>, index: isize) -> PyResult> { let triangulation = self.triangulation.bind(py).borrow(); - let index = normalize_index(index, triangulation.core.vertices.len()) - .map_err(TriangulationError::into_pyerr)?; + let index = normalize_index(index, triangulation.core.vertices.len())?; simplex_set_py(py, triangulation.core.simplices_of(index)) } @@ -581,14 +573,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")?; - let dim = Triangulation::validate_coords(&parsed_coords) - .map_err(TriangulationError::into_pyerr)?; + let dim = Triangulation::validate_coords(&parsed_coords)?; let core = match scipy_delaunay_simplices(py, &parsed_coords, dim)? { Some(initial) => Triangulation::from_simplices(parsed_coords, initial), None => Triangulation::new(parsed_coords), - } - .map_err(TriangulationError::into_pyerr)?; + }?; Ok(Self { core }) } @@ -598,14 +588,12 @@ impl PyTriangulation { simplex, self.core.vertices.len(), )?) - .map_err(TriangulationError::into_pyerr) + .map_err(PyErr::from) } fn delete_simplex(&mut self, simplex: &Bound<'_, PyAny>) -> PyResult<()> { let simplex = canonical_simplex_from_py(simplex, self.core.vertices.len())?; - self.core - .delete_simplex(&simplex) - .map_err(TriangulationError::into_pyerr) + self.core.delete_simplex(&simplex).map_err(PyErr::from) } #[pyo3(signature = (point, simplex=None, transform=None))] @@ -617,16 +605,12 @@ impl PyTriangulation { transform: Option<&Bound<'_, PyAny>>, ) -> PyResult<(Py, Py)> { let point = parse_point(point)?; - let simplex = match simplex { + let simplex = match given(simplex) { None => None, - Some(value) if value.is_none() => None, Some(value) => Some(canonical_simplex_from_py(value, self.core.vertices.len())?), }; let transform = parse_optional_transform(transform)?; - let (deleted, added) = self - .core - .add_point(point, simplex, transform) - .map_err(TriangulationError::into_pyerr)?; + let (deleted, added) = self.core.add_point(point, simplex, transform)?; Ok((simplex_set_py(py, &deleted)?, simplex_set_py(py, &added)?)) } @@ -661,14 +645,13 @@ impl PyTriangulation { py: Python<'_>, vertex: isize, ) -> PyResult> { - let vertex = normalize_index(vertex, slf.core.vertices.len()) - .map_err(TriangulationError::into_pyerr)?; + let vertex = normalize_index(vertex, slf.core.vertices.len())?; Py::new(py, PySimplicesProxy::for_vertex(slf.into(), vertex)) } #[getter(hull)] fn hull_property(&self, py: Python<'_>) -> PyResult> { - let hull = self.core.hull().map_err(TriangulationError::into_pyerr)?; + let hull = self.core.hull()?; let vertices: Vec = hull.into_iter().collect(); Ok(PySet::new(py, &vertices)?.into()) } @@ -699,8 +682,7 @@ impl PyTriangulation { items.push(py.None()); continue; } - let index = normalize_index(item.extract::()?, self.core.vertices.len()) - .map_err(TriangulationError::into_pyerr)?; + let index = normalize_index(item.extract::()?, self.core.vertices.len())?; items.push(point_tuple(py, &self.core.vertices[index]).into()); } Ok(PyList::new(py, items)?.into()) @@ -708,11 +690,7 @@ impl PyTriangulation { fn locate_point(&self, py: Python<'_>, point: &Bound<'_, PyAny>) -> PyResult> { let point = parse_point(point)?; - match self - .core - .locate_point(&point) - .map_err(TriangulationError::into_pyerr)? - { + match self.core.locate_point(&point)? { Some(simplex) => Ok(simplex_tuple(py, &simplex).into()), None => Ok(PyTuple::empty(py).into()), } @@ -731,17 +709,14 @@ impl PyTriangulation { // For a full simplex the reference echoes back the caller's indices // (including negative ones), so reduce positions rather than values. if simplex_signed.len() == self.core.dim + 1 { - self.core - .validate_point_dim(&point) - .map_err(TriangulationError::into_pyerr)?; - let simplex = normalize_indices(&simplex_signed, self.core.vertices.len()) - .map_err(TriangulationError::into_pyerr)?; + self.core.validate_point_dim(&point)?; + let simplex = normalize_indices(&simplex_signed, self.core.vertices.len())?; let alpha = match self.core.barycentric_alpha_for_simplex(&simplex, &point) { Ok(alpha) => alpha, Err(TriangulationError::Geometry(GeometryError::SingularMatrix)) => { return Err(numpy_linalg_error(py, "Singular matrix")); } - Err(other) => return Err(other.into_pyerr()), + Err(other) => return Err(other.into()), }; let reduced: Vec = reduced_simplex_positions(&alpha, eps) .into_iter() @@ -750,14 +725,13 @@ impl PyTriangulation { return signed_index_list_py(py, &reduced); } - let simplex = normalize_indices(&simplex_signed, self.core.vertices.len()) - .map_err(TriangulationError::into_pyerr)?; + let simplex = normalize_indices(&simplex_signed, self.core.vertices.len())?; let reduced = match self.core.get_reduced_simplex(&point, &simplex, eps) { Ok(reduced) => reduced, Err(TriangulationError::Geometry(GeometryError::SingularMatrix)) => { return Err(numpy_linalg_error(py, "Singular matrix")); } - Err(other) => return Err(other.into_pyerr()), + Err(other) => return Err(other.into()), }; index_list_py(py, &reduced) } @@ -771,10 +745,7 @@ impl PyTriangulation { ) -> PyResult<(Py, f64)> { let simplex = ordered_indices_from_py(simplex, self.core.vertices.len())?; let transform = parse_optional_transform(transform)?; - let (center, radius) = self - .core - .circumscribed_circle(&simplex, &transform) - .map_err(TriangulationError::into_pyerr)?; + let (center, radius) = self.core.circumscribed_circle(&simplex, &transform)?; Ok((point_tuple(py, ¢er), radius)) } @@ -785,13 +756,12 @@ impl PyTriangulation { simplex: &Bound<'_, PyAny>, transform: Option<&Bound<'_, PyAny>>, ) -> PyResult { - let pt_index = normalize_index(pt_index, self.core.vertices.len()) - .map_err(TriangulationError::into_pyerr)?; + let pt_index = normalize_index(pt_index, self.core.vertices.len())?; let simplex = ordered_indices_from_py(simplex, self.core.vertices.len())?; let transform = parse_optional_transform(transform)?; self.core .point_in_circumcircle(pt_index, &simplex, &transform) - .map_err(TriangulationError::into_pyerr) + .map_err(PyErr::from) } #[pyo3(name = "point_in_cicumcircle", signature = (pt_index, simplex, transform=None))] @@ -806,13 +776,11 @@ impl PyTriangulation { fn volume(&self, simplex: &Bound<'_, PyAny>) -> PyResult { let simplex = ordered_indices_from_py(simplex, self.core.vertices.len())?; - self.core - .volume(&simplex) - .map_err(TriangulationError::into_pyerr) + self.core.volume(&simplex).map_err(PyErr::from) } fn volumes(&self) -> PyResult> { - self.core.volumes().map_err(TriangulationError::into_pyerr) + self.core.volumes().map_err(PyErr::from) } #[pyo3(signature = (dim=None, simplices=None, vertices=None))] @@ -822,32 +790,23 @@ impl PyTriangulation { simplices: Option<&Bound<'_, PyAny>>, vertices: Option<&Bound<'_, PyAny>>, ) -> PyResult { - let simplices = match simplices { + let simplices = match given(simplices) { None => None, - Some(value) if value.is_none() => None, Some(value) => Some(simplex_set_from_py(value, self.core.vertices.len())?), }; - let vertices = match vertices { + let vertices = match given(vertices) { None => None, - Some(value) if value.is_none() => None, Some(value) => Some(vertex_index_set_from_py(value, self.core.vertices.len())?), }; let items = self .core - .faces(dim, simplices.as_ref(), vertices.as_ref()) - .map_err(TriangulationError::into_pyerr)?; + .faces(dim, simplices.as_ref(), vertices.as_ref())?; Ok(PyFacesIter { items, index: 0 }) } fn containing(&self, py: Python<'_>, face: &Bound<'_, PyAny>) -> PyResult> { let face = ordered_indices_from_py(face, self.core.vertices.len())?; - simplex_set_py( - py, - &self - .core - .containing(&face) - .map_err(TriangulationError::into_pyerr)?, - ) + simplex_set_py(py, &self.core.containing(&face)?) } fn reference_invariant(&self) -> bool { @@ -878,7 +837,7 @@ impl PyTriangulation { Err(TriangulationError::Geometry(GeometryError::SingularMatrix)) => { Err(numpy_linalg_error(py, "Singular matrix")) } - Err(other) => Err(other.into_pyerr()), + Err(other) => Err(other.into()), } } @@ -890,18 +849,15 @@ impl PyTriangulation { containing_simplex: Option<&Bound<'_, PyAny>>, transform: Option<&Bound<'_, PyAny>>, ) -> PyResult<(Py, Py)> { - let pt_index = normalize_index(pt_index, self.core.vertices.len()) - .map_err(TriangulationError::into_pyerr)?; - let containing_simplex = match containing_simplex { + let pt_index = normalize_index(pt_index, self.core.vertices.len())?; + let containing_simplex = match given(containing_simplex) { None => None, - Some(value) if value.is_none() => None, Some(value) => Some(canonical_simplex_from_py(value, self.core.vertices.len())?), }; let transform = parse_optional_transform(transform)?; let (deleted, added) = self .core - .bowyer_watson(pt_index, containing_simplex, &transform) - .map_err(TriangulationError::into_pyerr)?; + .bowyer_watson(pt_index, containing_simplex, &transform)?; Ok((simplex_set_py(py, &deleted)?, simplex_set_py(py, &added)?)) } @@ -910,8 +866,7 @@ impl PyTriangulation { match index { None => Ok(py.None()), Some(index) => { - let index = normalize_index(index, self.core.vertices.len()) - .map_err(TriangulationError::into_pyerr)?; + let index = normalize_index(index, self.core.vertices.len())?; Ok(point_tuple(py, &self.core.vertices[index]).into()) } } @@ -923,10 +878,7 @@ impl PyTriangulation { simplex: &Bound<'_, PyAny>, ) -> PyResult> { let indices = ordered_indices_from_py(simplex, self.core.vertices.len())?; - let neighbors = self - .core - .neighbors_from_vertices(&indices) - .map_err(TriangulationError::into_pyerr)?; + let neighbors = self.core.neighbors_from_vertices(&indices)?; simplex_set_py(py, &neighbors) } @@ -936,10 +888,7 @@ impl PyTriangulation { indices: &Bound<'_, PyAny>, ) -> PyResult> { let indices = ordered_indices_from_py(indices, self.core.vertices.len())?; - let attached = self - .core - .simplices_attached_to_points(&indices) - .map_err(TriangulationError::into_pyerr)?; + let attached = self.core.simplices_attached_to_points(&indices)?; simplex_set_py(py, &attached) } @@ -949,10 +898,7 @@ impl PyTriangulation { simplex: &Bound<'_, PyAny>, ) -> PyResult> { let simplex = ordered_indices_from_py(simplex, self.core.vertices.len())?; - let opposing = self - .core - .opposing_vertices(&simplex) - .map_err(TriangulationError::into_pyerr)?; + let opposing = self.core.opposing_vertices(&simplex)?; Ok(PyTuple::new(py, opposing)?.into()) } From f1ac88018fbda31d40e539c98c53bb01d4a395cf Mon Sep 17 00:00:00 2001 From: Bas Nijholt Date: Wed, 10 Jun 2026 11:57:04 -0700 Subject: [PATCH 2/3] refactor: shared helpers for the robust-predicate plumbing Readability only; no behavior change. - coord_2d/coord_3d replace the robust::Coord construction closures that were redefined inline in every predicate - sign() replaces the repeated three-way if/else blocks - validate_simplex() dedupes the 'dim + 1 points of equal dimension' check that circumsphere, volume, precise_volume, and robust_in_circumsphere each spelled out --- src/geometry.rs | 210 +++++++++++++++++++----------------------------- 1 file changed, 83 insertions(+), 127 deletions(-) diff --git a/src/geometry.rs b/src/geometry.rs index df71425..3dae551 100644 --- a/src/geometry.rs +++ b/src/geometry.rs @@ -303,14 +303,7 @@ pub fn fast_3d_circumsphere(points: &[[f64; 3]; 4]) -> ([f64; 3], f64) { /// `|x_i|^2 - |x_0|^2` form cancels catastrophically for small simplices far /// from the origin). Returns NaNs for degenerate input, matching numpy. pub fn circumsphere(pts: &[Vec]) -> Result<(Vec, f64), GeometryError> { - let dim = validate_points(pts)?; - if pts.len() != dim + 1 { - return Err(GeometryError::InvalidDimensions(format!( - "Expected {} points for a {}-dimensional simplex", - dim + 1, - dim - ))); - } + let dim = validate_simplex(pts)?; if dim == 2 { let points = [ @@ -441,6 +434,47 @@ pub fn orientation(face: &[Vec], origin: &[f64]) -> Result robust::Coord { + robust::Coord { + x: point[0], + y: point[1], + } +} + +/// `robust::Coord3D` from the first three coordinates of a point. +fn coord_3d(point: &[f64]) -> robust::Coord3D { + robust::Coord3D { + x: point[0], + y: point[1], + z: point[2], + } +} + +fn sign(value: f64) -> i32 { + if value > 0.0 { + 1 + } else if value < 0.0 { + -1 + } else { + 0 + } +} + +/// Validates that `points` are a full simplex (`dim + 1` points of equal +/// dimension) and returns `dim`. +fn validate_simplex(points: &[Vec]) -> Result { + let dim = validate_points(points)?; + if points.len() != dim + 1 { + return Err(GeometryError::InvalidDimensions(format!( + "Expected {} points for a {}-dimensional simplex", + dim + 1, + dim + ))); + } + Ok(dim) +} + /// Like [`orientation`], but exact in 2D and 3D: the sign is computed with /// Shewchuk's adaptive-precision predicates (via the `robust` crate), so it /// is correct even for slivers whose floating-point determinant rounds to @@ -460,54 +494,21 @@ pub fn robust_orientation(face: &[Vec], origin: &[f64]) -> Result { - 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 - }) - } + 2 => Ok(sign(robust::orient2d( + coord_2d(&face[0]), + coord_2d(&face[1]), + coord_2d(origin), + ))), + 3 => Ok(sign(robust::orient3d( + coord_3d(&face[0]), + coord_3d(&face[1]), + coord_3d(&face[2]), + coord_3d(origin), + ))), _ => orientation(face, origin), } } @@ -527,56 +528,36 @@ 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 - ))); + let dim = validate_simplex(simplex)?; + if point.len() != dim { + return Err(GeometryError::InvalidDimensions( + "Coordinates dimension mismatch".to_string(), + )); } + // incircle/insphere's sign follows the simplex's orientation; + // multiplying by it makes the test orientation-independent. 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])); + let [a, b, c] = [&simplex[0], &simplex[1], &simplex[2]].map(|p| coord_2d(p)); + let orient = robust::orient2d(a, b, c); 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)) + Ok(Some( + robust::incircle(a, b, c, coord_2d(point)) * 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]), - ); + let [a, b, c, d] = + [&simplex[0], &simplex[1], &simplex[2], &simplex[3]].map(|p| coord_3d(p)); + let orient = robust::orient3d(a, b, c, d); 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(Some( + robust::insphere(a, b, c, d, coord_3d(point)) * orient > 0.0, + )) } _ => Ok(None), } @@ -592,39 +573,21 @@ pub fn robust_in_circumsphere( /// 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 { + match validate_simplex(vertices)? { 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]), + coord_2d(&vertices[1]), + coord_2d(&vertices[2]), + coord_2d(&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]), + coord_3d(&vertices[1]), + coord_3d(&vertices[2]), + coord_3d(&vertices[3]), + coord_3d(&vertices[0]), ); Ok(det.abs() / 6.0) } @@ -635,14 +598,7 @@ pub fn precise_volume(vertices: &[Vec]) -> Result { /// Volume of a `dim`-simplex given its `dim + 1` vertices /// (`|det| / dim!` of the edge matrix). pub fn 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 - ))); - } + let dim = validate_simplex(vertices)?; let mut matrix = vec![vec![0.0; dim]; dim]; let x0 = &vertices[0]; From 6a3f6f6f7e6237f32fae193718ac00297625e700 Mon Sep 17 00:00:00 2001 From: Bas Nijholt Date: Wed, 10 Jun 2026 11:57:04 -0700 Subject: [PATCH 3/3] refactor: flatten the hard-to-read chains in the core Readability only; no behavior change. - TriangulationError::value() collapses the three-line Err(TriangulationError::Value("...".to_string())) returns - containing() uses plain loops instead of nested map/filter/collect chains wrapped across a dozen lines - facet_neighbour() names the 'other simplex across this facet' lookup shared by the walk and opposing_vertices - detach_vertex() dedupes the failed-insertion rollback that extend_hull and add_point both spelled out - summed_volume uses a for loop instead of a try_fold with a turbofish --- src/triangulation.rs | 196 ++++++++++++++++++++----------------------- 1 file changed, 93 insertions(+), 103 deletions(-) diff --git a/src/triangulation.rs b/src/triangulation.rs index c8e5c2b..7bdd40e 100644 --- a/src/triangulation.rs +++ b/src/triangulation.rs @@ -74,6 +74,12 @@ pub enum TriangulationError { Geometry(#[from] GeometryError), } +impl TriangulationError { + fn value(message: impl Into) -> Self { + Self::Value(message.into()) + } +} + pub(crate) fn index_out_of_range() -> TriangulationError { TriangulationError::Index("list index out of range".to_string()) } @@ -86,8 +92,8 @@ fn validate_transform( return Ok(()); }; if transform.len() != dim || transform.iter().any(|row| row.len() != dim) { - return Err(TriangulationError::Value( - "Transform must be an N x N matrix".to_string(), + return Err(TriangulationError::value( + "Transform must be an N x N matrix", )); } Ok(()) @@ -274,6 +280,13 @@ impl Triangulation { .is_some_and(|slot| slot.is_some()) } + /// The simplex on the other side of `facet` from `id`, or `None` when + /// the facet lies on the hull. + fn facet_neighbour(&self, facet: &[usize], id: SimplexId) -> Option { + let incident = self.facets.get(facet)?; + incident.iter().copied().find(|&other| other != id) + } + /// Intern `simplex` (which must already be sorted) and update every /// derived index. Returns `false` (and changes nothing) when it is /// already present. @@ -352,19 +365,17 @@ impl Triangulation { pub(crate) fn validate_coords(coords: &[Vec]) -> Result { if coords.is_empty() { - return Err(TriangulationError::Value( - "Please provide at least one simplex".to_string(), + return Err(TriangulationError::value( + "Please provide at least one simplex", )); } let dim = coords[0].len(); if coords.iter().any(|coord| coord.len() != dim) { - return Err(TriangulationError::Value( - "Coordinates dimension mismatch".to_string(), - )); + return Err(TriangulationError::value("Coordinates dimension mismatch")); } if coords.len() < dim + 1 { - return Err(TriangulationError::Value( - "Please provide at least one simplex".to_string(), + return Err(TriangulationError::value( + "Please provide at least one simplex", )); } @@ -373,8 +384,8 @@ impl Triangulation { .map(|coord| coord.iter().zip(&coords[0]).map(|(a, b)| a - b).collect()) .collect(); if geometry::numpy_matrix_rank(&vectors)? < dim { - return Err(TriangulationError::Value( - "Initial simplex has zero volumes (the points are linearly dependent)".to_string(), + return Err(TriangulationError::value( + "Initial simplex has zero volumes (the points are linearly dependent)", )); } @@ -383,9 +394,7 @@ impl Triangulation { pub(crate) fn validate_point_dim(&self, point: &[f64]) -> Result<(), TriangulationError> { if point.len() != self.dim { - return Err(TriangulationError::Value( - "Coordinates dimension mismatch".to_string(), - )); + return Err(TriangulationError::value("Coordinates dimension mismatch")); } Ok(()) } @@ -421,8 +430,8 @@ impl Triangulation { } } - Err(TriangulationError::Value( - "Initial simplex has zero volumes (the points are linearly dependent)".to_string(), + Err(TriangulationError::value( + "Initial simplex has zero volumes (the points are linearly dependent)", )) } @@ -500,8 +509,8 @@ impl Triangulation { let reduced_simplex = triangulation.get_reduced_simplex(&point, &actual_simplex, BARYCENTRIC_EPS)?; if reduced_simplex.is_empty() { - return Err(TriangulationError::Value( - "Point lies outside of the specified simplex.".to_string(), + return Err(TriangulationError::value( + "Point lies outside of the specified simplex.", )); } if reduced_simplex.len() == 1 { @@ -544,6 +553,16 @@ impl Triangulation { self.vertex_to_ids.pop(); } + /// Delete every simplex attached to `vertex` (rollback of a failed + /// insertion). + fn detach_vertex(&mut self, vertex: PointIndex) -> Result<(), TriangulationError> { + let attached: Vec = self.simplices_of(vertex).cloned().collect(); + for simplex in attached { + self.delete_simplex(&simplex)?; + } + Ok(()) + } + /// Number of simplices the given vertex belongs to. pub fn vertex_simplex_count(&self, vertex: PointIndex) -> usize { self.vertex_to_ids[vertex].len() @@ -584,7 +603,7 @@ impl Triangulation { let mut simplex = simplex.to_vec(); simplex.sort_unstable(); if !self.unlink_simplex(&simplex) { - return Err(TriangulationError::Value("Simplex not present".to_string())); + return Err(TriangulationError::value("Simplex not present")); } Ok(()) } @@ -649,11 +668,7 @@ impl Triangulation { } let face = facet_excluding(simplex, worst_idx); - let neighbour = self - .facets - .get(&face) - .and_then(|incident| incident.iter().copied().find(|&other| other != id)); - Ok(match neighbour { + Ok(match self.facet_neighbour(&face, id) { Some(next) => WalkStep::Neighbour(next), None => WalkStep::OutsideHull, }) @@ -741,8 +756,8 @@ impl Triangulation { vertices: Option<&FxHashSet>, ) -> Result, TriangulationError> { if simplices.is_some() && vertices.is_some() { - return Err(TriangulationError::Value( - "Only one of simplices and vertices is allowed.".to_string(), + return Err(TriangulationError::value( + "Only one of simplices and vertices is allowed.", )); } @@ -789,37 +804,36 @@ impl Triangulation { } self.validate_simplex_indices(face)?; + let mut result = FxHashSet::default(); + // A facet-sized face is a single lookup in the facet index. if face.len() == self.dim { let mut sorted = face.to_vec(); sorted.sort_unstable(); - return Ok(self - .facets - .get(&sorted) - .map(|incident| { - incident - .iter() - .map(|&id| self.simplex_by_id(id).clone()) - .collect() - }) - .unwrap_or_default()); + if let Some(incident) = self.facets.get(&sorted) { + for &id in incident { + result.insert(self.simplex_by_id(id).clone()); + } + } + return Ok(result); } + // Otherwise intersect the vertex incidence sets, scanning the vertex + // with the fewest simplices and testing each against the others. let first_vertex = face .iter() .copied() .min_by_key(|&vertex| self.vertex_to_ids[vertex].len()) .unwrap(); - Ok(self.vertex_to_ids[first_vertex] - .iter() - .copied() - .filter(|&id| { - face.iter().all(|&vertex| { - vertex == first_vertex || self.vertex_to_ids[vertex].contains(&id) - }) - }) - .map(|id| self.simplex_by_id(id).clone()) - .collect()) + for &id in &self.vertex_to_ids[first_vertex] { + let in_all = face + .iter() + .all(|&vertex| vertex == first_vertex || self.vertex_to_ids[vertex].contains(&id)); + if in_all { + result.insert(self.simplex_by_id(id).clone()); + } + } + Ok(result) } /// All simplices that share at least one vertex with the given vertex @@ -877,32 +891,27 @@ impl Triangulation { sorted.sort_unstable(); self.validate_simplex_indices(&sorted)?; let Some(&own_id) = self.ids.get(&sorted) else { - return Err(TriangulationError::Value( - "Provided simplex is not part of the triangulation".to_string(), + return Err(TriangulationError::value( + "Provided simplex is not part of the triangulation", )); }; - simplex - .iter() - .map(|&vertex| { - let position = sorted + let mut opposing = Vec::with_capacity(simplex.len()); + for &vertex in simplex { + let position = sorted + .iter() + .position(|&other| other == vertex) + .expect("vertex comes from the simplex itself"); + let facet = facet_excluding(&sorted, position); + opposing.push(self.facet_neighbour(&facet, own_id).map(|id| { + self.simplex_by_id(id) .iter() - .position(|&other| other == vertex) - .expect("vertex comes from the simplex itself"); - let facet = facet_excluding(&sorted, position); - let neighbour = self - .facets - .get(&facet) - .and_then(|incident| incident.iter().copied().find(|&id| id != own_id)); - Ok(neighbour.map(|id| { - self.simplex_by_id(id) - .iter() - .copied() - .find(|other| !sorted.contains(other)) - .expect("facet neighbour has exactly one vertex outside the simplex") - })) - }) - .collect() + .copied() + .find(|other| !sorted.contains(other)) + .expect("facet neighbour has exactly one vertex outside the simplex") + })); + } + Ok(opposing) } fn simplex_points( @@ -1112,9 +1121,7 @@ impl Triangulation { 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(), - )); + return Err(TriangulationError::value("Point already in triangulation.")); } } @@ -1132,9 +1139,7 @@ impl Triangulation { .iter() .all(|id| bad_ids.contains(id)) { - return Err(TriangulationError::Value( - "Point already in triangulation.".to_string(), - )); + return Err(TriangulationError::value("Point already in triangulation.")); } } @@ -1197,9 +1202,8 @@ impl Triangulation { 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)) + let mut new_triangles: FxHashSet = self + .simplices_of(pt_index) .filter(|simplex| !bad_set.contains(simplex)) .cloned() .collect(); @@ -1222,11 +1226,11 @@ impl Triangulation { /// (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)?)?, - ) - }) + let mut total = 0.0; + for simplex in simplices { + total += geometry::precise_volume(&self.get_vertices(simplex)?)?; + } + Ok(total) } /// Rebuild the cavity with exact in-circumsphere predicates (2D/3D @@ -1414,15 +1418,9 @@ impl Triangulation { } if new_simplices.is_empty() { - 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)?; - } - return Err(TriangulationError::Value( - "Candidate vertex is inside the hull.".to_string(), + self.detach_vertex(pt_index)?; + return Err(TriangulationError::value( + "Candidate vertex is inside the hull.", )); } @@ -1465,13 +1463,7 @@ impl Triangulation { // 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)?; - } + self.detach_vertex(pt_index)?; self.pop_vertex(); return Err(err); } @@ -1491,16 +1483,14 @@ impl Triangulation { let reduced_simplex = self.get_reduced_simplex(&point, &simplex, BARYCENTRIC_EPS)?; if reduced_simplex.is_empty() { - return Err(TriangulationError::Value( - "Point lies outside of the specified simplex.".to_string(), + return Err(TriangulationError::value( + "Point lies outside of the specified simplex.", )); } simplex = reduced_simplex; if simplex.len() == 1 { - return Err(TriangulationError::Value( - "Point already in triangulation.".to_string(), - )); + return Err(TriangulationError::value("Point already in triangulation.")); } let pt_index = self.add_vertex(point);