diff --git a/CHANGELOG.md b/CHANGELOG.md index 46ed617..9fad682 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,10 +5,28 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](http://keepachangelog.com/) and this project adheres to [Semantic Versioning](http://semver.org/). -## [2.13.1] - 2025-04-03 +## [2.15.0] - 2025-08-16 + +### Added + - Implements `ConstrainedDelaunayTriangulation::try_bulk_load_cdt` for loading CDTs without panicking + if any constraints intersect (#137) + +### Changed + - `ConstrainedDelaunayTriangulation::bulk_load_cdt` does load all vertices in order by default ("stable"). + - `ConstrainedDelaunayTriangulation::bulk_load_cdt_stable` is deprecated. Use `...::bulk_load_cdt` instead which + behaves exactly the same. + - Optimized implementation of now stable `ConstrainedDelaunayTriangulation::bulk_load_cdt`. It is nearly as fast + as its unstable alternative, `Triangulation::bulk_load`. + - Bulk loading a triangulation with a `HierarchyHintGenerator` does not load the hierarchy sequentially but uses + bulk loading. This should speed up creating such triangulations. + +## Fixed + - Fixes various panics related to bulk loading very exotic vertex sets. + +## [2.14.0] - 2025-04-03 ### Fix - - Fixes #116 and #119 + - Fixes #116 and #130 ## [2.13.0] - 2025-02-20 @@ -486,7 +504,9 @@ A lot has changed for the 1.0. release, only larger changes are shown. Initial commit -[2.13.1]: https://github.com/Stoeoef/spade/compare/v2.13.0...v2.13.1 +[2.15.0]: https://github.com/Stoeoef/spade/compare/v2.14.0...v2.15.0 + +[2.14.0]: https://github.com/Stoeoef/spade/compare/v2.13.0...v2.14.0 [2.13.0]: https://github.com/Stoeoef/spade/compare/v2.12.1...v2.13.0 diff --git a/Cargo.toml b/Cargo.toml index 918f613..5092afb 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -8,13 +8,7 @@ description = "Delaunay triangulations for the rust ecosystem" repository = "https://github.com/Stoeoef/spade" license = "MIT OR Apache-2.0" autobenches = false -categories = [ - "algorithms", - "data-structures", - "graphics", - "mathematics", - "science::geo", -] +categories = ["algorithms", "data-structures", "graphics", "mathematics"] keywords = ["Delaunay", "CDT", "geometry", "triangulation", "voronoi"] [lib] @@ -25,15 +19,15 @@ default = ["std"] std = [] [dependencies] -smallvec = "1.13" -robust = "1.1.0" -num-traits = "0.2" -hashbrown = "0.15.2" +smallvec = "1.15.1" +robust = "1.2.0" +num-traits = "0.2.19" +hashbrown = "0.15.5" [dependencies.serde] package = "serde" optional = true -version = "1.0.218" +version = "1.0.219" default-features = false features = ["derive", "alloc"] @@ -47,18 +41,18 @@ version = "0.5.9" members = ["delaunay_compare"] [dev-dependencies] -approx = "0.5" -rand = "0.9.1" +approx = "0.5.1" +rand = "0.9.2" cgmath = "0.18.0" svg = "0.18.0" -float_next_after = "1" -image = "0.25.1" -tiny-skia = "0.11.3" -criterion = { version = "0.5.1", features = ["html_reports"] } +float_next_after = "1.0.0" +image = "0.25.6" +tiny-skia = "0.11.4" +criterion = { version = "0.7.0", features = ["html_reports"] } base64 = "0.22.1" -anyhow = "1.0.97" -shapefile = "0.6.0" -proptest = "1.5.0" +anyhow = "1.0.99" +shapefile = "0.7.0" +proptest = "1.7.0" [[bench]] name = "benchmarks" diff --git a/delaunay_compare/examples/europe.png b/delaunay_compare/examples/europe.png index c86cbb2..e98fa8c 100644 Binary files a/delaunay_compare/examples/europe.png and b/delaunay_compare/examples/europe.png differ diff --git a/delaunay_compare/examples/real_data_benchmark.rs b/delaunay_compare/examples/real_data_benchmark.rs index 7bcc289..f40d851 100644 --- a/delaunay_compare/examples/real_data_benchmark.rs +++ b/delaunay_compare/examples/real_data_benchmark.rs @@ -67,7 +67,7 @@ fn load_with_spade(vertices: Vec>, edges: Vec<[usize; 2]>) -> anyhow let edges_clone = edges.clone(); let now = Instant::now(); - ConstrainedDelaunayTriangulation::<_>::bulk_load_cdt_stable(vertices_clone, edges_clone)?; + ConstrainedDelaunayTriangulation::<_>::bulk_load_cdt(vertices_clone, edges_clone)?; println!( "loading time (spade cdt bulk load stable): {}ms", diff --git a/examples/svg_renderer/scenario_list.rs b/examples/svg_renderer/scenario_list.rs index 26781ef..e8d2cd6 100644 --- a/examples/svg_renderer/scenario_list.rs +++ b/examples/svg_renderer/scenario_list.rs @@ -1576,6 +1576,6 @@ fn get_cdt_for_try_add_constraint() -> Result { VertexType::new(50.0, -34.0), ]; - let result = Cdt::bulk_load_cdt_stable(vertices, vec![[3, 2], [5, 4], [7, 6]])?; + let result = Cdt::bulk_load_cdt(vertices, vec![[3, 2], [5, 4], [7, 6]])?; Ok(result) } diff --git a/fuzz/fuzz_targets/bulk_load_cdt_fuzz.rs b/fuzz/fuzz_targets/bulk_load_cdt_fuzz.rs index aacb7b1..d8c3405 100644 --- a/fuzz/fuzz_targets/bulk_load_cdt_fuzz.rs +++ b/fuzz/fuzz_targets/bulk_load_cdt_fuzz.rs @@ -47,7 +47,7 @@ fuzz_target!(|input: (Vec, Vec<[usize; 2]>)| { edges.truncate(last_index); let bulk_loaded = - ConstrainedDelaunayTriangulation::::bulk_load_cdt(data, edges).unwrap(); + ConstrainedDelaunayTriangulation::::bulk_load_cdt(data.clone(), edges).unwrap(); bulk_loaded.cdt_sanity_check(); }); diff --git a/src/cdt.rs b/src/cdt.rs index 05b2655..8c897db 100644 --- a/src/cdt.rs +++ b/src/cdt.rs @@ -6,8 +6,10 @@ use num_traits::{Float, NumCast}; use serde::{Deserialize, Serialize}; use crate::cdt::ConflictRegionEnd::{EdgeOverlap, Existing}; -use crate::delaunay_core::dcel_operations::flip_cw; -use crate::delaunay_core::{bulk_load_cdt, bulk_load_stable, try_bulk_load_cdt}; +use crate::delaunay_core::dcel_operations::{ + append_unconnected_vertex, flip_cw, new_with_fixed_vertices, +}; +use crate::delaunay_core::{bulk_load_cdt, try_bulk_load_cdt}; use crate::{ delaunay_core::Dcel, intersection_iterator::LineIntersectionIterator, PositionInTriangulation, SpadeNum, @@ -61,6 +63,10 @@ impl CdtEdge { pub fn data_mut(&mut self) -> &mut UE { &mut self.1 } + + pub(crate) fn deconstruct(self) -> (bool, UE) { + (self.0, self.1) + } } impl Default for CdtEdge { @@ -311,13 +317,12 @@ where /// /// The edges are given as pairs of vertex indices. /// - /// Note that the vertex order is not preserved by this function - iterating through all vertices will not result in - /// the same sequence as the input vertices. Use [ConstrainedDelaunayTriangulation::bulk_load_cdt_stable] for a - /// slower but order preserving variant. + /// The vertex order is preserved by this function. /// /// Input vertices may have the same position. However, only one vertex for each position will be kept. Edges /// that go to a discarded vertex are rerouted and still inserted. - /// It is arbitrary which duplicated vertex remains. + /// For each set of duplicate vertices, the vertex with the smallest index is kept. Any resulting gap is filled + /// by shifting the remaining vertices down. /// /// # Example /// ``` @@ -328,15 +333,18 @@ where /// Point2::new(1.0, 2.0), /// Point2::new(3.0, -3.0), /// Point2::new(-1.0, -2.0), + /// Point2::new(3.0, -3.0), // Duplicate /// Point2::new(-4.0, -5.0), /// ]; - /// let mut edges = vec![[0, 1], [1, 2], [2, 3], [3, 4]]; + /// let mut edges = vec![[0, 1], [1, 2], [2, 3], [4, 5]]; /// let cdt = ConstrainedDelaunayTriangulation::<_>::bulk_load_cdt(vertices.clone(), edges)?; /// - /// assert_eq!(cdt.num_vertices(), 5); + /// assert_ne!(cdt.num_vertices(), vertices.len()); // Duplicates are removed /// assert_eq!(cdt.num_constraints(), 4); - /// // The order will usually change - /// assert_ne!(cdt.vertices().map(|v| v.position()).collect::>(), vertices); + /// + /// // The order is preserved (excluding duplicates): + /// vertices.remove(4); + /// assert_eq!(cdt.vertices().map(|v| v.position()).collect::>(), vertices); /// # Ok(()) /// # } /// ``` @@ -345,13 +353,17 @@ where /// /// Panics if any constraint edges overlap. Panics if the edges contain an invalid index (out of range). pub fn bulk_load_cdt(vertices: Vec, edges: Vec<[usize; 2]>) -> Result { - let mut result = bulk_load_cdt(vertices, edges).unwrap(); - *result.hint_generator_mut() = L::initialize_from_triangulation(&result); - Ok(result) + bulk_load_cdt(vertices, edges) } - /// Same behaviour as [bulk_load_cdt], but rather than panicking, - /// ignores and calls the parameter function `on_conflict_found` for each conflicting edges encountered. + /// Same behavior as [bulk_load_cdt], but rather than panicking, + /// skips and calls the parameter function `on_conflict_found` whenever a conflict occurs. + /// + /// For any conflicting pair of constraint edges it is unspecified which constraint is reported. + /// + /// `on_conflict_found` is called with the edge indices _after_ potential duplicates were removed. + /// Consider checking for removed duplicates by comparing the size of the triangulation with the + /// input vertices. /// /// # Example /// ``` @@ -361,111 +373,45 @@ where /// Point2::new(0.0, 1.0), /// Point2::new(1.0, 2.0), /// Point2::new(3.0, -3.0), - /// Point2::new(-1.0, -2.0), + /// Point2::new(-4.0, -2.0), + /// Point2::new(3.0, -3.0), // Duplicate /// Point2::new(-4.0, -5.0), /// ]; /// let mut conflicting_edges = Vec::new(); - /// let mut edges = vec![[0, 1], [1, 2], [2, 3], [3, 4]]; + /// let mut edges = vec![[0, 1], [1, 2], [2, 3], [4, 5], [5, 0]]; /// let cdt = ConstrainedDelaunayTriangulation::<_>::try_bulk_load_cdt(vertices.clone(), edges, |e| conflicting_edges.push(e))?; /// /// assert_eq!(cdt.num_vertices(), 5); - /// assert_eq!(cdt.num_constraints(), 4); - /// // The order will usually change - /// assert_ne!(cdt.vertices().map(|v| v.position()).collect::>(), vertices); - /// # Ok(()) - /// # } - /// ``` - pub fn try_bulk_load_cdt( - vertices: Vec, - edges: Vec<[usize; 2]>, - on_conflict_found: impl FnMut([usize; 2]), - ) -> Result { - let mut result = try_bulk_load_cdt(vertices, edges, on_conflict_found)?; - *result.hint_generator_mut() = L::initialize_from_triangulation(&result); - Ok(result) - } - - /// Stable bulk load variant that preserves the input vertex order + /// assert_eq!(cdt.num_constraints(), 4); // One constraint was not generated /// - /// The resulting vertex set will be equal to the input vertex set if their positions are all distinct. - /// See [ConstrainedDelaunayTriangulation::bulk_load_cdt] for additional details like panic behavior and duplicate - /// handling. - /// - /// # Example - /// ``` - /// # fn main() -> Result<(), spade::InsertionError> { - /// use spade::{ConstrainedDelaunayTriangulation, Point2, Triangulation}; - /// let mut vertices = vec![ - /// Point2::new(0.0, 1.0), - /// Point2::new(1.0, 2.0), - /// Point2::new(3.0, -3.0), - /// Point2::new(-1.0, -2.0), - /// Point2::new(-4.0, -5.0), - /// ]; - /// let mut edges = vec![[0, 1], [1, 2], [2, 3], [3, 4]]; - /// let cdt = ConstrainedDelaunayTriangulation::<_>::bulk_load_cdt_stable(vertices.clone(), edges)?; + /// // The third and last edge generate a conflict. One of them will be reported via `on_conflict_found`. + /// // In this case, the last edge is returned. Note the index change from v5 to v4 due to the duplicate removal. + /// assert_eq!(conflicting_edges, [[4, 0]]); /// - /// // The ordered will be preserved: + /// // The order is preserved (excluding duplicates): + /// vertices.remove(4); /// assert_eq!(cdt.vertices().map(|v| v.position()).collect::>(), vertices); /// # Ok(()) /// # } /// ``` - /// - /// It is fine to include vertex positions multiple times. The resulting order will be the same as if - /// the duplicates were removed prior to insertion. However, it is unclear *which* duplicates are - /// removed - e.g. do not assume that always the first duplicated vertex remains. - /// - /// ``` - /// # fn main() -> Result<(), spade::InsertionError> { - /// use spade::{ConstrainedDelaunayTriangulation, Point2, Triangulation}; - /// let mut vertices = vec![ - /// Point2::new(0.0, 1.0), - /// Point2::new(1.0, 2.0), // Duplicate - /// Point2::new(1.0, 2.0), - /// Point2::new(3.0, -3.0), - /// Point2::new(3.0, -3.0), // Duplicate - /// Point2::new(-4.0, -5.0), - /// ]; - /// let mut edges = vec![[0, 1], [2, 3], [4, 5]]; - /// let cdt = ConstrainedDelaunayTriangulation::<_>::bulk_load_cdt_stable(vertices.clone(), edges)?; - /// - /// // The choice of deduplicated vertices is arbitrary. In this example, dedup[1] and dedup[2] could - /// // have been swapped - /// let dedup = [ - /// Point2::new(0.0, 1.0), - /// Point2::new(1.0, 2.0), - /// Point2::new(3.0, -3.0), - /// Point2::new(-4.0, -5.0), - /// ]; - /// assert_eq!(cdt.vertices().map(|v| v.position()).collect::>(), dedup); - /// # Ok(()) - /// # } - /// ``` - /// - /// # See also - /// - /// See also [Self::try_bulk_load_cdt_stable] - pub fn bulk_load_cdt_stable( + pub fn try_bulk_load_cdt( vertices: Vec, edges: Vec<[usize; 2]>, + on_conflict_found: impl FnMut([usize; 2]), ) -> Result { - Self::try_bulk_load_cdt_stable(vertices, edges, |e| { - panic!("Conflicting edge encountered: {};{}", e[0], e[1]) - }) + try_bulk_load_cdt(vertices, edges, on_conflict_found) } - /// See [Self::bulk_load_cdt_stable] - pub fn try_bulk_load_cdt_stable( + /// Deprecated. Use [Self::bulk_load_cdt] instead which is now stable by default. + #[deprecated( + since = "1.15.0", + note = "Use bulk_load_cdt instead. It is now stable by default." + )] + pub fn bulk_load_cdt_stable( vertices: Vec, edges: Vec<[usize; 2]>, - on_conflict_found: impl FnMut([usize; 2]), ) -> Result { - let mut result: Self = bulk_load_stable( - |vertices| try_bulk_load_cdt(vertices, edges, on_conflict_found), - vertices, - )?; - *result.hint_generator_mut() = L::initialize_from_triangulation(&result); - Ok(result) + Self::bulk_load_cdt(vertices, edges) } /// # Handle invalidation @@ -746,7 +692,7 @@ where &self, from: Point2<::Scalar>, to: Point2<::Scalar>, - ) -> impl Iterator, F>> { + ) -> impl Iterator, F>> { LineIntersectionIterator::new(self, from, to) .flat_map(|intersection| intersection.as_edge_intersection()) .filter(|e| e.is_constraint_edge()) @@ -762,7 +708,7 @@ where &self, from: FixedVertexHandle, to: FixedVertexHandle, - ) -> impl Iterator, F>> { + ) -> impl Iterator, F>> { LineIntersectionIterator::new_from_handles(self, from, to) .flat_map(|intersection| intersection.as_edge_intersection()) .filter(|e| e.is_constraint_edge()) @@ -1015,7 +961,8 @@ where alternative_vertex } else { let edge = edge.fix(); - let (new_vertex, [e0, e1]) = self.insert_on_edge(edge, new_vertex); + let new_vertex = append_unconnected_vertex(self.s_mut(), new_vertex); + let [e0, e1] = self.insert_on_edge(edge, new_vertex); self.handle_legal_edge_split([e0, e1]); self.legalize_vertex(new_vertex); new_vertex @@ -1167,6 +1114,32 @@ where constraint_edges } + + pub(crate) fn create_with_fixed_vertices( + elements: Vec, + first_vertex: FixedVertexHandle, + second_vertex: FixedVertexHandle, + ) -> Self { + Self { + dcel: new_with_fixed_vertices(elements, first_vertex, second_vertex), + num_constraints: 0, + // Doesn't make sense to initialize with fixed vertices. Bulk loading will initialize + // and replace this afterwards anyway. + hint_generator: Default::default(), + } + } + + /// Swaps out and re-initializes the hint generator. + pub(crate) fn adjust_hint_generator( + self, + ) -> ConstrainedDelaunayTriangulation + where + L2: HintGenerator<::Scalar>, + { + let hint_generator = L2::initialize_from_triangulation(&self); + let (dcel, _, num_constraints) = self.into_parts(); + ConstrainedDelaunayTriangulation::from_parts(dcel, hint_generator, num_constraints) + } } impl ConstrainedDelaunayTriangulation @@ -1349,6 +1322,8 @@ pub fn get_edge_intersections( #[cfg(test)] mod test { + use std::collections::HashSet; + use approx::assert_abs_diff_eq; use proptest::prelude::*; @@ -1546,7 +1521,7 @@ mod test { for p in delaunay_points { d.insert(p)?; } - let mut used_vertices = hashbrown::HashSet::new(); + let mut used_vertices = HashSet::new(); let mut inserted_constraints = Vec::new(); for v in d.vertices() { @@ -2044,7 +2019,7 @@ mod test { Point2::new(50.0, -34.0), ]; - Cdt::bulk_load_cdt_stable(vertices, vec![[3, 2], [5, 4], [7, 6]]) + Cdt::bulk_load_cdt(vertices, vec![[3, 2], [5, 4], [7, 6]]) } #[test] @@ -2061,11 +2036,11 @@ mod test { Point2::new(50.0, -34.0), ]; let mut conflicting_edges = Vec::new(); - let cdt = Cdt::try_bulk_load_cdt_stable(vertices, vec![[3, 2], [5, 4], [7, 6]], |e| { + let cdt = Cdt::try_bulk_load_cdt(vertices, vec![[3, 2], [5, 4], [7, 6]], |e| { conflicting_edges.push(e) })?; // Hardcoded values, may change if CDT algorithm change - assert_eq!(&conflicting_edges, &[[6, 7,], [4, 5,]]); + assert_eq!(&conflicting_edges, &[[5, 6,], [3, 4,]]); assert_eq!(cdt.num_constraints, 1); Ok(()) } @@ -2079,7 +2054,7 @@ mod test { Point2::new(0.0, 1.0), ]; - let mut cdt = Cdt::bulk_load_cdt_stable(vertices, vec![[2, 3]])?; + let mut cdt = Cdt::bulk_load_cdt(vertices, vec![[2, 3]])?; let initial_num_vertices = cdt.num_vertices(); let from = FixedVertexHandle::from_index(0); @@ -2139,16 +2114,15 @@ mod test { #[test] fn test_remove_vertex_respects_constraints() -> Result<(), InsertionError> { - let mut triangulation = - ConstrainedDelaunayTriangulation::>::bulk_load_cdt_stable( - vec![ - Point2::new(-1.0, 0.0), - Point2::new(0.0, 0.0), - Point2::new(0.0, 1.0), - Point2::new(1.0, 1.0), - ], - vec![[0, 1], [1, 2], [2, 3]], - )?; + let mut triangulation = ConstrainedDelaunayTriangulation::>::bulk_load_cdt( + vec![ + Point2::new(-1.0, 0.0), + Point2::new(0.0, 0.0), + Point2::new(0.0, 1.0), + Point2::new(1.0, 1.0), + ], + vec![[0, 1], [1, 2], [2, 3]], + )?; let mut copy = triangulation.clone(); @@ -2287,7 +2261,7 @@ mod test { (68.78571, 71.541046), ]; - run_proptest(&points) + run_bulk_load_proptest(&points) } #[test] @@ -2299,7 +2273,7 @@ mod test { (6.3033395, 32.111538), (89.09958, 1.0720834), ]; - run_proptest(&points) + run_bulk_load_proptest(&points) } #[test] @@ -2313,7 +2287,7 @@ mod test { (6.3033395, 32.111538), (89.09958, 1.0720834), ]; - run_proptest(&points) + run_bulk_load_proptest(&points) } #[test] @@ -2371,7 +2345,7 @@ mod test { } #[test] - fn foo() -> Result<(), InsertionError> { + fn cdt_test_inserting_few_points() -> Result<(), InsertionError> { let vertices = [ Point2::new(43.44877, 78.04408), Point2::new(43.498154, 81.00147), @@ -2418,10 +2392,10 @@ mod test { (75.59119, 42.91933), (25.808325, 97.32154), ]; - run_proptest(&points) + run_bulk_load_proptest(&points) } - fn run_proptest(points: &[(f32, f32)]) -> Result<(), InsertionError> { + fn run_bulk_load_proptest(points: &[(f32, f32)]) -> Result<(), InsertionError> { let mut cdt = ConstrainedDelaunayTriangulation::>::new(); let mut last_handle = None; @@ -2451,7 +2425,7 @@ mod test { (38.64656, 24.732662), ]; - run_proptest(&points) + run_bulk_load_proptest(&points) } #[test] @@ -2473,7 +2447,167 @@ mod test { (27.416088, 88.37326), ]; - run_proptest(&points) + run_bulk_load_proptest(&points) + } + + proptest! { + #![proptest_config(ProptestConfig::with_cases(200))] + #[test] + #[ignore] + fn + prop_test_bulk_load_cdt_stable( + points in proptest::collection::vec((1.0f32..100.0f32, 1.0f32..100.0f32), 1..100), + edges in proptest::collection::vec((0usize..100, 0usize..100), 1..100)) { + prop_test_bulk_load_cdt_stable_impl(points, edges); + } + } + + #[test] + fn fuzz_test_bulk_load_cdt() { + let points = vec![ + Point2::new(-2.7049442424493675e-11, -2.704772938494617e-11), + Point2::new(-2.7049442424493675e-11, -2.2374774353910703e-35), + Point2::new(-2.704944241771741e-11, -2.704922037988875e-11), + Point2::new(-2.7049442424493675e-11, -2.7049442424493572e-11), + Point2::new(-2.7049442424493675e-11, 0.0), + ]; + + let cdt = + ConstrainedDelaunayTriangulation::>::bulk_load_cdt(points, Vec::new()) + .unwrap(); + + cdt.cdt_sanity_check(); + } + + #[test] + fn fuzz_test_bulk_load_cdt_2() { + { + let points = vec![ + Point2::new(0.11617647291654172, -2.7049442424493675e-11), + Point2::new(-3.684373062283352e-14, -2.7049442424493675e-11), + Point2::new(-2.7049442424493268e-11, -2.7049442424493675e-11), + Point2::new(-2.7049442424493675e-11, -2.7049442424493675e-11), + Point2::new(-2.704943949713427e-11, -2.7049442424493856e-11), + ]; + let edges = vec![[3, 0]]; + + let cdt = ConstrainedDelaunayTriangulation::>::bulk_load_cdt(points, edges) + .unwrap(); + + cdt.cdt_sanity_check(); + } + } + + #[test] + fn fuzz_test_bulk_load_cdt_3() { + { + let points = vec![ + Point2::new(-2.7049442424378697e-11, -2.7049442424493675e-11), + Point2::new(-2.7049434889184558e-11, -2.7049442424493675e-11), + Point2::new(-2.7049442395058873e-11, -2.7049442424493675e-11), + Point2::new(-2.7049442424546208e-11, -2.7049442424493675e-11), + Point2::new(0.0004538143382352941, -2.705036194996328e-11), + Point2::new(-2.704944239484691e-11, -2.7049442424493675e-11), + Point2::new(-2.7049442424546615e-11, -2.7049442424493675e-11), + Point2::new(-2.7049442424493882e-11, -2.70494406897702e-11), + Point2::new(0.11617551691391888, -2.7049442424493675e-11), + Point2::new(-2.7049442424493572e-11, -2.7049442424493675e-11), + Point2::new(-2.704944241771741e-11, -2.704944242438945e-11), + Point2::new(-2.7049442424493678e-11, -2.7049442424493662e-11), + Point2::new(-2.704944241771741e-11, -2.7049442424493675e-11), + Point2::new(-2.7049442424467205e-11, -2.7049442424493675e-11), + Point2::new(-2.7049442424493675e-11, -2.7049442424493675e-11), + Point2::new(-2.7049442424493675e-11, -2.7049441611342046e-11), + Point2::new(-2.7049442424493675e-11, -2.7049442424493675e-11), + Point2::new(-2.7049442424493675e-11, -2.7049442424493675e-11), + Point2::new(-2.704944241771741e-11, -2.7049442424493675e-11), + Point2::new(-2.7049442424493675e-11, -2.7049442424493675e-11), + Point2::new(-2.7049442424493675e-11, -2.7049442424493675e-11), + Point2::new(-2.704944242448623e-11, -2.7049442424493675e-11), + Point2::new(-2.7049442424493675e-11, -2.7049442424493675e-11), + Point2::new(-2.7049442424493675e-11, -2.7049442424493675e-11), + Point2::new(-2.7049442424493675e-11, -2.6502324517958398e-11), + Point2::new(-2.7049442424493675e-11, -2.7049442424493675e-11), + Point2::new(-2.7049442424493675e-11, -2.7049442423646642e-11), + Point2::new(-2.704944242437787e-11, -2.7049442424493675e-11), + Point2::new(-2.7049442424493675e-11, -2.7049442424493675e-11), + Point2::new(-2.7049442424493675e-11, -2.7049442424493675e-11), + Point2::new(-2.7049442424493646e-11, -2.6375537048546205e-11), + ]; + let edges = vec![]; + + let cdt = ConstrainedDelaunayTriangulation::>::bulk_load_cdt(points, edges) + .unwrap(); + + cdt.cdt_sanity_check(); + } + } + + #[test] + fn fuzz_test_bulk_load_cdt_4() { + { + let points = vec![ + Point2::new(5.532608426405843e-5, -2.705000112790462e-11), + Point2::new(-2.7049442384471368e-11, -2.6555608406219246e-11), + Point2::new(-4.6189490530823523e-10, -2.7049442424337338e-11), + Point2::new(-2.688601759526906e-11, -2.7049442424517663e-11), + Point2::new(-2.704944242448623e-11, -2.7049442424493675e-11), + Point2::new(-2.704944240840005e-11, -2.7049442424493675e-11), + Point2::new(-1.1276731182827942e-13, -2.7049442424493675e-11), + ]; + + let edges = vec![]; + + ConstrainedDelaunayTriangulation::>::bulk_load(points.clone()).unwrap(); + + let cdt = ConstrainedDelaunayTriangulation::>::bulk_load_cdt(points, edges) + .unwrap(); + + cdt.cdt_sanity_check(); + } + } + + fn prop_test_bulk_load_cdt_stable_impl(points: Vec<(f32, f32)>, edges: Vec<(usize, usize)>) { + let max_index = points.len() - 1; + let edges = edges + .iter() + .copied() + .map(|(from, to)| [from.min(max_index), to.min(max_index)]) + .collect::>(); + + let vertices = points + .iter() + .copied() + .map(|(x, y)| Point2::new(x, y)) + .collect::>(); + let result = ConstrainedDelaunayTriangulation::>::try_bulk_load_cdt( + vertices, + edges, + |_| {}, + ) + .unwrap(); + + result.cdt_sanity_check(); + } + + #[test] + fn bulk_load_cdt_stable_proptest_1() { + let points = vec![(1.0, 1.0), (61.541332, 1.0), (1.0, 69.54), (1.0, 19.13391)]; + let edges = vec![(3, 3)]; + prop_test_bulk_load_cdt_stable_impl(points, edges) + } + + #[test] + fn bulk_load_cdt_stable_proptest_2() { + let points = vec![ + (1.0, 98.26258), + (61.541332, 59.135162), + (1.0, 1.0), + (57.31029, 1.0), + ]; + let edges = vec![(0, 3)]; + + prop_test_bulk_load_cdt_stable_impl(points, edges) } fn check_returned_edges( diff --git a/src/delaunay_core/bulk_load.rs b/src/delaunay_core/bulk_load.rs index 1de2d41..e4dded4 100644 --- a/src/delaunay_core/bulk_load.rs +++ b/src/delaunay_core/bulk_load.rs @@ -1,16 +1,26 @@ +use core::f64; +use core::panic; + +use crate::delaunay_core::math; +use crate::delaunay_core::triangulation_ext::VertexToInsert; use crate::{ - ConstrainedDelaunayTriangulation, HasPosition, HintGenerator, InsertionError, Point2, - Triangulation, TriangulationExt, + ConstrainedDelaunayTriangulation, HasPosition, HintGenerator, InsertionError, + LastUsedVertexHintGenerator, Point2, Triangulation, TriangulationExt, }; use core::cmp::{Ordering, Reverse}; + +#[cfg(not(feature = "std"))] +use hashbrown::HashSet; +#[cfg(feature = "std")] +use std::collections::HashSet; + use num_traits::Zero; use super::{ dcel_operations, FixedDirectedEdgeHandle, FixedUndirectedEdgeHandle, FixedVertexHandle, }; -use alloc::vec; -use alloc::vec::Vec; +use alloc::{vec, vec::Vec}; /// An `f64` wrapper implementing `Ord` and `Eq`. /// @@ -20,6 +30,14 @@ use alloc::vec::Vec; #[derive(Debug, PartialEq, PartialOrd, Clone, Copy)] struct FloatOrd(S); +fn bulk_load_distance_fn(center: Point2, position: Point2) -> impl Ord { + ( + Reverse(FloatOrd(center.distance_2(position))), + FloatOrd(position.x), + FloatOrd(position.y), + ) +} + #[allow(clippy::derive_ord_xor_partial_ord)] impl Ord for FloatOrd where @@ -104,9 +122,7 @@ where let mut result = T::with_capacity(elements.len(), elements.len() * 3, elements.len() * 2); // Sort by distance, smallest values last. This allows to pop values depending on their distance. - elements.sort_unstable_by_key(|e| { - Reverse(FloatOrd(initial_center.distance_2(e.position().to_f64()))) - }); + elements.sort_unstable_by_key(|e| bulk_load_distance_fn(initial_center, e.position().to_f64())); let mut hull = loop { let Some(next) = elements.pop() else { @@ -115,11 +131,20 @@ where result.insert(next)?; - if let Some(hull) = try_get_hull_center(&result) + if result.num_vertices() < 3 { + continue; + } + + let thee_vertices = [ + FixedVertexHandle::new(result.num_vertices() - 3), + FixedVertexHandle::new(result.num_vertices() - 2), + FixedVertexHandle::new(result.num_vertices() - 1), + ]; + + if let Some(hull) = try_get_hull_center(&result, thee_vertices) .and_then(|center| Hull::from_triangulation(&result, center)) { hull_sanity_check(&result, &hull); - break hull; } }; @@ -129,12 +154,16 @@ where } let mut buffer = Vec::new(); - let mut skipped_elements = Vec::::new(); + let mut skipped_elements = Vec::new(); while let Some(next) = elements.pop() { - skipped_elements.extend( - single_bulk_insertion_step(&mut result, false, &mut hull, next, &mut buffer).err(), - ); + skipped_elements.extend(single_bulk_insertion_step( + &mut result, + false, + &mut hull, + VertexToInsert::NewVertex(next), + &mut buffer, + )) } if cfg!(any(fuzzing, test)) { @@ -143,8 +172,8 @@ where fix_convexity(&mut result); - for element in skipped_elements { - result.insert(element)?; + for skipped in skipped_elements { + result.insert_with_hint_option_impl(skipped, None); } Ok(result) @@ -169,7 +198,7 @@ where L: HintGenerator<::Scalar>, { try_bulk_load_cdt(elements, edges, |e| { - panic!("Conflicting edge encountered: {};{}", e[0], e[1]) + panic!("Conflicting edge encountered: [{}, {}]", e[0], e[1]) }) } @@ -179,7 +208,7 @@ where /// This function does not panic if any two constraints intersect. /// It will instead call `on_conflict_found` on all edges that could not be added as they intersect a previously added constraint. pub fn try_bulk_load_cdt( - elements: Vec, + mut elements: Vec, mut edges: Vec<[usize; 2]>, mut on_conflict_found: impl FnMut([usize; 2]), ) -> Result, InsertionError> @@ -194,327 +223,341 @@ where return Ok(ConstrainedDelaunayTriangulation::new()); } - if edges.is_empty() { - return bulk_load(elements); - } - - let mut point_sum = Point2::::new(0.0, 0.0); - - for element in &elements { - crate::validate_vertex(element)?; - let position = element.position(); - - point_sum = point_sum.add(position.to_f64()); - } + let element_indices = sort_and_dedup(&mut elements, &mut edges)?; - // Set the initial center to the average of all positions. This should be a good choice for most triangulations. - // - // The research paper uses a different approach by taking the center of the points' bounding box. - // However, this position might be far off the center off mass if the triangulation has just a few outliers. - // This could lead to a very uneven angle distribution as nearly all points are might be in a very small angle segment - // around the center. This degrades the hull-structure's lookup and insertion performance. - // For this reason, taking the average appears to be a safer option as most vertices should be distributed around the - // initial center. - let initial_center = point_sum.mul(1.0 / (elements.len() as f64)); - - let mut result = ConstrainedDelaunayTriangulation::with_capacity( - elements.len(), - elements.len() * 3, - elements.len() * 2, - ); - - let distance_fn = |position: Point2<::Scalar>| { - ( - Reverse(FloatOrd(initial_center.distance_2(position.to_f64()))), - FloatOrd(position.x), - FloatOrd(position.y), - ) - }; - - for edge in &mut edges { - let [d1, d2] = edge.map(|vertex| distance_fn(elements[vertex].position())); - if d1 > d2 { - edge.reverse(); + if elements.len() < 3 { + // Special case for fewer than 3 vertices. Starting the hull calculation requires at least 3 vertices. + let mut result = ConstrainedDelaunayTriangulation::new(); + for elem in elements { + result.insert(elem)?; } - } - - edges.sort_by_cached_key(|[from, _]| distance_fn(elements[*from].position())); - - let mut elements = elements.into_iter().enumerate().collect::>(); - - // Sort by distance, smallest values last. This allows to pop values depending on their distance. - elements.sort_unstable_by_key(|(_, e)| distance_fn(e.position())); - let mut old_to_new = vec![usize::MAX; elements.len()]; - let mut last_position = None; - let mut last_index = 0; - for (old_index, e) in elements.iter().rev() { - let position = e.position(); - if last_position.is_some() && Some(position) != last_position { - last_index += 1; + for &[from, to] in &edges { + result.add_constraint(FixedVertexHandle::new(from), FixedVertexHandle::new(to)); } - old_to_new[*old_index] = last_index; - last_position = Some(position); + return Ok(result); } - let mut next_constraint = edges.pop(); - let mut buffer = Vec::new(); + + let mut skipped_vertices = HashSet::new(); + let mut skipped_edges = HashSet::new(); + + let fix = FixedVertexHandle::new; + let mut add_constraints_for_new_vertex = - |result: &mut ConstrainedDelaunayTriangulation, index| { - while let Some([from, to]) = next_constraint { - // Check if next creates any constraint edge - if old_to_new[from] == old_to_new[index] { - let [new_from, new_to] = - [from, to].map(|v| FixedVertexHandle::new(old_to_new[v])); + |result: &mut ConstrainedDelaunayTriangulation, + edges: &mut Vec<[usize; 2]>, + skipped_vertices: &HashSet, + index| { + while let Some([from, to]) = edges.last().copied() { + if skipped_vertices.contains(&from) || skipped_vertices.contains(&to) { + skipped_edges.insert([from, to]); + edges.pop(); + continue; + } + + // Check if the next edge creates any constraint edge + if from == index { + let [from, to] = [from, to].map(fix); // Insert constraint edge - if result.try_add_constraint(new_from, new_to).is_empty() { - on_conflict_found([from, to]); + if result.try_add_constraint(from, to).is_empty() { + on_conflict_found([from.index(), to.index()]); } - next_constraint = edges.pop(); + edges.pop(); } else { break; } } }; - let mut hull = loop { - let Some((old_index, next)) = elements.pop() else { - return Ok(result); - }; - result.insert(next)?; - add_constraints_for_new_vertex(&mut result, old_index); + // Goes through all element indices, closest to farthest away + let mut elements_iter = element_indices.iter().rev().copied(); + + let first = fix(elements_iter.next().unwrap()); + let second = fix(elements_iter.next().unwrap()); - if let Some(hull) = try_get_hull_center(&result) + let mut result = ConstrainedDelaunayTriangulation::::create_with_fixed_vertices( + elements, + first, + second, + ); + + let third = result.insert_with_hint_option_impl( + VertexToInsert::ExistingVertex(fix(elements_iter.next().unwrap())), + Some(first), + ); + + add_constraints_for_new_vertex(&mut result, &mut edges, &skipped_vertices, first.index()); + add_constraints_for_new_vertex(&mut result, &mut edges, &skipped_vertices, second.index()); + add_constraints_for_new_vertex(&mut result, &mut edges, &skipped_vertices, third.index()); + + let mut last_3_vertices = [first, second, third]; + let mut next_handle_to_drop = 0; + + let mut hull = loop { + if let Some(hull) = try_get_hull_center(&result, last_3_vertices) .and_then(|center| Hull::from_triangulation(&result, center)) { break hull; } + + let Some(index) = elements_iter.next() else { + return Ok(result.adjust_hint_generator()); + }; + + let next_handle = result.insert_with_hint_option_impl( + VertexToInsert::ExistingVertex(FixedVertexHandle::new(index)), + Some(first), + ); + + add_constraints_for_new_vertex(&mut result, &mut edges, &skipped_vertices, index); + + last_3_vertices[next_handle_to_drop] = next_handle; + next_handle_to_drop = (next_handle_to_drop + 1) % 3; }; - while let Some((old_index, next)) = elements.pop() { - if let Err(skipped) = - single_bulk_insertion_step(&mut result, true, &mut hull, next, &mut buffer) + let mut require_convexity = true; + + for next_index in elements_iter { + // Optimization: Once all edges are inserted, convexity isn't required anymore. + require_convexity &= !edges.is_empty(); + + if single_bulk_insertion_step( + &mut result, + require_convexity, + &mut hull, + VertexToInsert::ExistingVertex(FixedVertexHandle::new(next_index)), + &mut buffer, + ) + .is_some() { - // Sometimes the bulk insertion step fails due to floating point inaccuracies. - // The easiest way to handle these rare occurrences is by skipping them. However, this doesn't - // work as CDT vertices **must** be inserted in their predefined order (after sorting for distance) - // to keep `old_to_new` lookup accurate. - // Instead, this code leverages that the triangulation for CDTs is always convex: This - // means that `result.insert` should work. Unfortunately, using `insert` will invalidate - // the hull structure. We'll recreate it with a loop similar to the initial hull creation. - // - // This process is certainly confusing and inefficient but, luckily, rarely required for real inputs. - - // Push the element again, it will be popped directly. This seems to be somewhat simpler than - // the alternatives. - elements.push((old_index, skipped)); - hull = loop { - let Some((old_index, next)) = elements.pop() else { - return Ok(result); - }; - result.insert(next)?; - add_constraints_for_new_vertex(&mut result, old_index); - - if let Some(hull) = Hull::from_triangulation(&result, hull.center) { - break hull; - }; - }; - } else { - add_constraints_for_new_vertex(&mut result, old_index); + skipped_vertices.insert(next_index); } - } - assert_eq!(edges.len(), 0); + add_constraints_for_new_vertex(&mut result, &mut edges, &skipped_vertices, next_index); + } if cfg!(any(fuzzing, test)) { hull_sanity_check(&result, &hull); } - Ok(result) -} + if !require_convexity { + fix_convexity(&mut result); + } -fn try_get_hull_center(result: &T) -> Option> -where - V: HasPosition, - T: Triangulation, -{ - let zero = ::Scalar::zero(); - if !result.all_vertices_on_line() && result.num_vertices() >= 4 { - // We'll need 4 vertices to calculate a center position with good precision. - // Otherwise, dividing by 3.0 can introduce precision loss and errors. - - // Get new center that is usually within the convex hull - let center_positions = || result.vertices().rev().take(4).map(|v| v.position()); - - let sum_x = center_positions() - .map(|p| p.x) - .fold(zero, |num, acc| num + acc); - let sum_y = center_positions() - .map(|p| p.y) - .fold(zero, |num, acc| num + acc); - - // Note that we don't re-sort the elements according to their distance to the newest center. This doesn't seem to - // be required for the algorithms performance, probably due to the `center` being close to `initial_center`. - // As of now, it is unclear how to construct point sets that result in a `center` being farther off - // `initial center` and what the impact of this would be. - let center = Point2::new(sum_x, sum_y).mul(::Scalar::from(0.25f32)); - - if let crate::PositionInTriangulation::OnFace(_) = result.locate(center) { - return Some(center.to_f64()); + for element in skipped_vertices { + result.insert_with_hint_option_impl( + VertexToInsert::ExistingVertex(FixedVertexHandle::new(element)), + Some(first), + ); + } + + for [from, to] in skipped_edges { + let [from, to] = [from, to].map(FixedVertexHandle::new); + if result.try_add_constraint(from, to).is_empty() { + on_conflict_found([from.index(), to.index()]); } } - None -} + assert_eq!(edges, Vec::<[usize; 2]>::new()); -pub(crate) struct PointWithIndex { - data: V, - index: usize, + Ok(result.adjust_hint_generator()) } -impl HasPosition for PointWithIndex +fn sort_and_dedup( + elements: &mut Vec, + edges: &mut [[usize; 2]], +) -> Result, InsertionError> where V: HasPosition, { - type Scalar = V::Scalar; - fn position(&self) -> Point2 { - self.data.position() + let mut point_sum = Point2::::new(0.0, 0.0); + for element in &*elements { + crate::validate_vertex(element)?; + let position = element.position(); + + point_sum = point_sum.add(position.to_f64()); } -} + let initial_center = point_sum.mul(1.0 / (elements.len() as f64)); -pub fn bulk_load_stable( - constructor: Constructor, - elements: Vec, -) -> Result -where - V: HasPosition, - T: Triangulation, - T2: Triangulation< - Vertex = PointWithIndex, - DirectedEdge = T::DirectedEdge, - UndirectedEdge = T::UndirectedEdge, - Face = T::Face, - HintGenerator = T::HintGenerator, - >, - Constructor: FnOnce(Vec>) -> Result, -{ - let elements = elements - .into_iter() - .enumerate() - .map(|(index, data)| PointWithIndex { index, data }) - .collect::>(); - - let num_original_elements = elements.len(); - - let mut with_indices = constructor(elements)?; - - if with_indices.num_vertices() != num_original_elements { - // Handling duplicates is more complicated - we cannot simply swap the elements into - // their target position indices as these indices may contain gaps. The following code - // fills those gaps. - // - // Running example: The original indices in with_indices could look like - // - // [3, 0, 1, 4, 6] - // - // which indicates that the original elements at indices 2 and 5 were duplicates. - let mut no_gap = (0usize..with_indices.num_vertices()).collect::>(); - - // This will be sorted by their original index: - // no_gap (before sorting): [0, 1, 2, 3, 4] - // keys for sorting : [3, 0, 1, 4, 6] - // no_gap (after sorting) : [1, 2, 0, 3, 4] - // sorted keys : [0, 1, 3, 4, 6] - no_gap.sort_unstable_by_key(|elem| { - with_indices - .vertex(FixedVertexHandle::new(*elem)) - .data() - .index + for edge in edges.iter_mut() { + let [d1, d2] = edge.map(|vertex| { + bulk_load_distance_fn(initial_center, elements[vertex].position().to_f64()) }); - - // Now, the sequential target index for FixedVertexHandle(no_gap[i]) is i - // - // Example: - // Vertex index in with_indices: [0, 1, 2, 3, 4] - // Original target indices : [3, 0, 1, 4, 6] - // Sequential target index : [2, 0, 1, 3, 4] - for (sequential_index, vertex) in no_gap.into_iter().enumerate() { - with_indices - .vertex_data_mut(FixedVertexHandle::new(vertex)) - .index = sequential_index; + if d1 > d2 { + edge.reverse(); } } - // Swap elements until the target order is restored. - // The attached indices for each vertex are guaranteed to form a permutation over all index - // since gaps are eliminated in the step above. - let mut current_index = 0; - loop { - if current_index >= with_indices.num_vertices() { - break; - } + edges.sort_unstable_by_key(|[from, _]| { + bulk_load_distance_fn(initial_center, elements[*from].position().to_f64()) + }); + let mut element_indices = (0..elements.len()).collect::>(); + element_indices.sort_unstable_by_key(|index| { + bulk_load_distance_fn(initial_center, elements[*index].position().to_f64()) + }); + let mut removed_indices = HashSet::new(); + let mut last_position = None; + let mut duplicate_sets = Vec::new(); + let mut current_duplicates = Vec::new(); + let mut current_min_duplicate = element_indices[0]; + element_indices.retain(|&index| { + let position = elements[index].position(); + + let mut retain = true; + if Some(position) == last_position { + retain = false; // Duplicate found - only the first index is kept + + // Check if the index is the smallest index of the current duplicate set + let new_duplicate = if index < current_min_duplicate { + let previous_min = current_min_duplicate; + current_min_duplicate = index; + previous_min + } else { + index + }; - // Example: The permutation [0 -> 2, 1 -> 0, 2 -> 1, 3 -> 3, 4 -> 4] - // (written as [2, 0, 1, 3, 4]) will lead to the following swaps: - // Swap 2, 0 (leading to [1, 0, 2, 3, 4]) - // Swap 1, 0 (leading to [0, 1, 2, 3, 4]) - // Done - let new_index = FixedVertexHandle::new(current_index); - let old_index = with_indices.vertex(new_index).data().index; - if current_index == old_index { - current_index += 1; + removed_indices.insert(new_duplicate); + current_duplicates.push(new_duplicate); } else { - with_indices - .s_mut() - .swap_vertices(FixedVertexHandle::new(old_index), new_index); + if !current_duplicates.is_empty() { + if duplicate_sets.is_empty() { + // Lazy initialize to avoid allocation. Each element is initially mapped to itself. + duplicate_sets = (0..elements.len()).collect::>(); + } + + for dupe in current_duplicates.iter() { + duplicate_sets[*dupe] = current_min_duplicate; + } + } + + // Reset duplicate information + current_duplicates.clear(); + last_position = Some(position); + current_min_duplicate = index; + }; + + if !current_duplicates.is_empty() { + if duplicate_sets.is_empty() { + duplicate_sets = (0..elements.len()).collect::>(); + } + + for dupe in current_duplicates.iter() { + duplicate_sets[*dupe] = current_min_duplicate; + } } + retain + }); + + if !removed_indices.is_empty() { + let mut old_to_new = vec![usize::MAX; elements.len()]; + let mut index_counter = 0; + let mut i = 0; + elements.retain(|_| { + old_to_new[i] = index_counter; + if removed_indices.contains(&i) { + old_to_new[i] = old_to_new[duplicate_sets[i]]; + i += 1; + return false; + } + old_to_new[i] = index_counter; + + i += 1; + index_counter += 1; + true + }); + + for elem in &mut element_indices { + *elem = old_to_new[*elem]; + } + + for &mut [ref mut from, ref mut to] in edges { + *from = old_to_new[*from]; + *to = old_to_new[*to]; + } + } + Ok(element_indices) +} + +fn try_get_hull_center( + result: &T, + three_vertices: [FixedVertexHandle; 3], +) -> Option> +where + V: HasPosition, + T: Triangulation, +{ + if result.all_vertices_on_line() { + return None; + } + + let zero = ::Scalar::zero(); + // Get new center that is usually within the convex hull + let mut sum_x = zero; + let mut sum_y = zero; + + for handle in three_vertices.iter() { + let position = result.vertex(*handle).position(); + sum_x = sum_x + position.x; + sum_y = sum_y + position.y; } - let (dcel, hint_generator, num_constraints) = with_indices.into_parts(); - let dcel = dcel.map_vertices(|point_with_index| point_with_index.data); + // Note that we don't re-sort the elements according to their distance to the newest center. This doesn't seem to + // be required for the algorithms performance, probably due to the `center` being close to `initial_center`. + // As of now, it is unclear how to construct point sets that result in a `center` being farther off + // `initial center` and what the impact of this would be. + let center = Point2::new(sum_x, sum_y).mul(::Scalar::from(1f32 / 3f32)); - Ok(T::from_parts(dcel, hint_generator, num_constraints)) + if let crate::PositionInTriangulation::OnFace(_) = + result.locate_with_hint(center, three_vertices[0]) + { + return Some(center.to_f64()); + } + None } +/// Returns `handle` if it couldn't be inserted due to precision issues. #[inline(never)] // Prevent inlining for better profiling data fn single_bulk_insertion_step( result: &mut TR, require_convexity: bool, hull: &mut Hull, - element: T, + handle: VertexToInsert, buffer_for_edge_legalization: &mut Vec, -) -> Result<(), T> +) -> Option> where T: HasPosition, TR: Triangulation, { - let next_position = element.position(); - let current_angle = pseudo_angle(next_position.to_f64(), hull.center); + let next_position = handle.position(result.s()); - let edge_hint = hull.get(current_angle); + let current_angle = pseudo_angle(next_position.to_f64(), hull.center); - let edge = result.directed_edge(edge_hint); + let edge = result.directed_edge(hull.get(current_angle)); let [from, to] = edge.positions(); if next_position == from || next_position == to { - return Ok(()); + return None; } if edge.side_query(next_position).is_on_right_side_or_on_line() { // The position is, for some reason, not on the left side of the edge. This can e.g. happen // if the vertices have the same angle. The safest way to include these elements appears to // skip them and insert them individually at the end (albeit that's very slow) - return Err(element); + return Some(handle); } assert!(edge.is_outer_edge()); let edge = edge.fix(); + let handle = handle.resolve(result.s_mut()); + let new_vertex = - dcel_operations::create_new_face_adjacent_to_edge(result.s_mut(), edge, element); + dcel_operations::create_new_face_adjacent_to_edge(result.s_mut(), edge, handle); + let ccw_walk_start = result.directed_edge(edge).prev().rev().fix(); let cw_walk_start = result.directed_edge(edge).next().rev().fix(); @@ -552,22 +595,17 @@ where let prev = handle.prev(); let handle = handle.fix(); - let [prev_from, prev_to] = prev.positions(); + let [prev_from, prev_to] = prev.positions().map(|p| p.to_f64()); // `!point_projection.is_behind_edge` is used to identify if the new face's angle will be less // than 90° let angle_condition = require_convexity - || !super::math::project_point(next_position, prev_to, prev_from).is_behind_edge(); + || !math::is_behind_edge(next_position.to_f64(), prev_to, prev_from) + // Prevents inserting the same angle twice into the hull which can cause issues + || current_angle == pseudo_angle(prev_to, hull.center); current_edge = prev.fix(); if angle_condition && prev.side_query(next_position).is_on_left_side() { - let prev_prev = prev.prev(); - if prev - .side_query(prev_prev.from().position()) - .is_on_left_side_or_on_line() - { - assert!(prev_prev.side_query(next_position).is_on_left_side()); - } let new_edge = dcel_operations::create_single_face_between_edge_and_next( result.s_mut(), current_edge, @@ -592,17 +630,18 @@ where let next = handle.next(); let handle = handle.fix(); + let middle_position = next.from().position().to_f64(); + let angle_condition = require_convexity - || !super::math::project_point( - next.from().position(), - next_position, - next.to().position(), + || !super::math::is_behind_edge( + middle_position, + next_position.to_f64(), + next.to().position().to_f64(), ) - .is_behind_edge(); + || current_angle == pseudo_angle(middle_position, hull.center); let next_fix = next.fix(); let is_on_left_side = next.side_query(next_position).is_on_left_side(); - if angle_condition && is_on_left_side { let new_edge = dcel_operations::create_single_face_between_edge_and_next( result.s_mut(), @@ -637,8 +676,13 @@ where first_edge.fix(), second_edge.fix(), ); + + if cfg!(any(fuzzing, test)) { + hull_sanity_check(result, hull); + } } - Ok(()) + + None } /// Makes the outer hull convex. Similar to a graham scan. @@ -884,27 +928,22 @@ impl Hull { left_index = current_node.right; } - let mut right_index; - if left_angle == right_angle { - right_index = left_index; - } else { - right_index = self.data[left_index].right; - loop { - let current_node = self.data[right_index]; - if current_node.angle == right_angle { - break; - } - - if cfg!(any(fuzzing, test)) { - assert!(!self.empty.contains(&right_index)); - } + let mut right_index = self.data[left_index].right; + loop { + let current_node = self.data[right_index]; + if current_node.angle == right_angle { + break; + } - // Remove current_node - it is completely overlapped by the new segment - self.empty.push(right_index); - self.data[current_node.left].right = current_node.right; - self.data[current_node.right].left = current_node.left; - right_index = current_node.right; + if cfg!(any(fuzzing, test)) { + assert!(!self.empty.contains(&right_index)); } + + // Remove current_node - it is completely overlapped by the new segment + self.empty.push(right_index); + self.data[current_node.left].right = current_node.right; + self.data[current_node.right].left = current_node.left; + right_index = current_node.right; } let new_index = self.get_next_index(); @@ -1048,9 +1087,7 @@ impl Hull { #[inline] fn pseudo_angle(a: Point2, center: Point2) -> FloatOrd { let diff = a.sub(center); - let p = diff.x / (diff.x.abs() + diff.y.abs()); - FloatOrd(1.0 - (if diff.y > 0.0 { 3.0 - p } else { 1.0 + p }) * 0.25) } @@ -1096,6 +1133,7 @@ mod test { use float_next_after::NextAfter; use rand::{seq::SliceRandom, SeedableRng}; + use crate::delaunay_core::triangulation_ext::VertexToInsert; use crate::handles::FixedVertexHandle; use crate::test_utilities::{random_points_with_seed, SEED2}; @@ -1170,8 +1208,7 @@ mod test { #[test] fn test_bulk_load_on_epsilon_grid() -> Result<(), InsertionError> { - // TODO: Setting this to 20 currently generates an inexplicably failing test case. Investigate! - const GRID_SIZE: usize = 18; + const GRID_SIZE: usize = 20; let mut rng = rand::rngs::StdRng::from_seed(*SEED2); const TEST_REPETITIONS: usize = 30; @@ -1197,8 +1234,7 @@ mod test { // Creates a zig zag pattern let edges = (0..vertices.len() - 1).map(|i| [i, i + 1]).collect(); - let triangulation = - ConstrainedDelaunayTriangulation::<_>::bulk_load_cdt_stable(vertices, edges)?; + let triangulation = ConstrainedDelaunayTriangulation::<_>::bulk_load_cdt(vertices, edges)?; triangulation.cdt_sanity_check(); assert_eq!(triangulation.num_vertices(), GRID_SIZE * GRID_SIZE); @@ -1259,8 +1295,7 @@ mod test { let num_constraints = edges.len(); let num_vertices = vertices.len(); - let cdt = - ConstrainedDelaunayTriangulation::<_>::bulk_load_cdt_stable(vertices, edges.clone())?; + let cdt = ConstrainedDelaunayTriangulation::<_>::bulk_load_cdt(vertices, edges.clone())?; cdt.cdt_sanity_check(); assert_eq!(cdt.num_vertices(), num_vertices); @@ -1299,7 +1334,7 @@ mod test { let edges = vec![[0, 1], [9, 3], [11, 2], [10, 0]]; let num_constraints = edges.len(); - let cdt = ConstrainedDelaunayTriangulation::<_>::bulk_load_cdt_stable(vertices, edges)?; + let cdt = ConstrainedDelaunayTriangulation::<_>::bulk_load_cdt(vertices, edges)?; assert_eq!(cdt.num_constraints(), num_constraints); Ok(()) } @@ -1355,10 +1390,8 @@ mod test { #[test] fn test_empty() -> Result<(), InsertionError> { - let cdt = ConstrainedDelaunayTriangulation::>::bulk_load_cdt_stable( - Vec::new(), - Vec::new(), - )?; + let cdt = + ConstrainedDelaunayTriangulation::>::bulk_load_cdt(Vec::new(), Vec::new())?; assert_eq!(cdt.num_vertices(), 0); assert_eq!(cdt.num_constraints(), 0); @@ -1429,14 +1462,15 @@ mod test { ]; for (index, element) in additional_elements.iter().enumerate() { - super::single_bulk_insertion_step( + assert!(super::single_bulk_insertion_step( &mut triangulation, false, &mut hull, - *element, + VertexToInsert::NewVertex(*element), &mut Vec::new(), ) - .unwrap(); + .is_none()); + if index != 0 { super::hull_sanity_check(&triangulation, &hull) } diff --git a/src/delaunay_core/bulk_load_fuzz_tests.rs b/src/delaunay_core/bulk_load_fuzz_tests.rs index 40e363e..436102a 100644 --- a/src/delaunay_core/bulk_load_fuzz_tests.rs +++ b/src/delaunay_core/bulk_load_fuzz_tests.rs @@ -8,6 +8,10 @@ fn fuzz_test(vertices: Vec>) { clone.dedup(); let expected_size = clone.len(); + let triangulation = DelaunayTriangulation::<_>::bulk_load_stable(vertices.clone()).unwrap(); + triangulation.sanity_check(); + assert_eq!(triangulation.num_vertices(), expected_size); + let triangulation = DelaunayTriangulation::<_>::bulk_load(vertices).unwrap(); triangulation.sanity_check(); assert_eq!(triangulation.num_vertices(), expected_size); @@ -354,12 +358,7 @@ fn test_bulk_load_fuzz_21() { Point2::new(-2.7049442424493675e-11, -2.7049442384314453e-11), Point2::new(-2.7049442424493675e-11, -2.7049442424493675e-11), Point2::new(-2.7049442424493675e-11, -2.6992599005632867e-11), - Point2::new(-2.7049442424493675e-11, -2.7049442424493675e-11), - Point2::new(-2.7049442424493675e-11, -2.7049442424493675e-11), - Point2::new(-2.7049442424493675e-11, -2.7049442424493675e-11), Point2::new(-2.693575558677206e-11, -2.7049442424493675e-11), - Point2::new(-2.7049442424493675e-11, -2.7049442424493675e-11), - Point2::new(-2.7049442424493675e-11, -2.7049442424493675e-11), Point2::new(-2.7049442424493675e-11, -2.704943548559977e-11), Point2::new(-2.7049442424493662e-11, -2.7049442424493675e-11), Point2::new(-2.6811408873290566e-11, -2.7049442424493675e-11), @@ -381,11 +380,35 @@ fn test_bulk_load_fuzz_21() { Point2::new(-2.7049442424493675e-11, -2.7049442424493672e-11), Point2::new(-2.7049442424467205e-11, -2.7049441936602697e-11), Point2::new(-2.7049442424493675e-11, -2.704944239865917e-11), - Point2::new(-2.7049442424493675e-11, -2.7049442424493675e-11), Point2::new(-2.704944242258785e-11, -2.7049442424493675e-11), Point2::new(-2.7049442398980546e-11, -2.7049442395376918e-11), - Point2::new(-2.7049442424493675e-11, -2.7049442424493675e-11), Point2::new(-1.600710678326356e-10, -2.7049442424493675e-11), Point2::new(-2.7049442424909747e-11, -2.7049442424493675e-11), ]); } + +#[test] +fn test_bulk_load_fuzz_22() { + fuzz_test(vec![ + Point2::new(-0.0003849186150368648, 4.8578744028535106e-33), + Point2::new(4.860884325263086e-33, 4.8578750597249806e-33), + Point2::new(4.8578750597249806e-33, 4.8578750597249806e-33), + Point2::new(4.8578750597249806e-33, 4.8579279569707704e-33), + Point2::new(4.857875059725681e-33, 4.8578750597249806e-33), + Point2::new(8.961197796880322e-14, 4.8578750597249806e-33), + Point2::new(8.955648037009912e-14, 4.884958449567926e-33), + ]); +} + +#[test] +fn test_bulk_load_fuzz_23() { + fuzz_test(vec![ + Point2::new(-2.1571074338520272e-11, -2.1132376894135366e-13), + Point2::new(-2.7226523566846033e-40, 0.11608455882352732), + Point2::new(-6.806630891709812e-41, 0.1238970588235294), + Point2::new(-2.70478881122592e-11, -2.722652356683964e-40), + Point2::new(-2.722652356831814e-40, 0.11608453241362017), + Point2::new(-2.6601800500862655e-11, -1.127672851467415e-13), + Point2::new(0.0, 0.0), + ]); +} diff --git a/src/delaunay_core/dcel.rs b/src/delaunay_core/dcel.rs index 1a58647..c57db8e 100644 --- a/src/delaunay_core/dcel.rs +++ b/src/delaunay_core/dcel.rs @@ -7,14 +7,6 @@ use serde::{Deserialize, Serialize}; use alloc::vec::Vec; -#[derive(Default, PartialEq, Eq, PartialOrd, Ord, Debug, Clone, Copy, Hash)] -pub struct EdgeData { - directed_data: [DE; 2], - pub undirected_data: UE, -} - -impl EdgeData {} - #[derive(Clone, Copy, Debug, PartialEq, Eq, PartialOrd, Ord, Hash)] #[cfg_attr( feature = "serde", @@ -133,7 +125,7 @@ impl Dcel { self.faces.truncate(1); // Keep outer face } - pub fn get_vertex(&self, handle: FixedVertexHandle) -> Option> { + pub fn get_vertex(&self, handle: FixedVertexHandle) -> Option> { (handle.index() < self.vertices.len()).then(|| self.vertex(handle)) } @@ -190,7 +182,7 @@ impl Dcel { self.faces.len() } - pub fn vertex(&self, handle: FixedVertexHandle) -> VertexHandle { + pub fn vertex(&self, handle: FixedVertexHandle) -> VertexHandle<'_, V, DE, UE, F> { DynamicHandleImpl::new(self, handle) } @@ -203,18 +195,18 @@ impl Dcel { pub fn directed_edge( &self, handle: FixedDirectedEdgeHandle, - ) -> DirectedEdgeHandle { + ) -> DirectedEdgeHandle<'_, V, DE, UE, F> { DirectedEdgeHandle::new(self, handle) } pub fn undirected_edge( &self, handle: FixedUndirectedEdgeHandle, - ) -> UndirectedEdgeHandle { + ) -> UndirectedEdgeHandle<'_, V, DE, UE, F> { UndirectedEdgeHandle::new(self, handle) } - pub fn outer_face(&self) -> FaceHandle { + pub fn outer_face(&self) -> FaceHandle<'_, PossiblyOuterTag, V, DE, UE, F> { let outer_face = super::dcel_operations::OUTER_FACE_HANDLE; self.face(outer_face) } @@ -267,7 +259,7 @@ impl Dcel { pub fn face( &self, handle: FixedHandleImpl, - ) -> DynamicHandleImpl { + ) -> DynamicHandleImpl<'_, V, DE, UE, F, FaceTag, InnerOuter> { DynamicHandleImpl::new(self, handle) } @@ -310,7 +302,7 @@ impl Dcel { &self, from: FixedVertexHandle, to: FixedVertexHandle, - ) -> Option> { + ) -> Option> { let vertex = self.vertex(from); for edge in vertex.out_edges() { if edge.to().fix() == to.adjust_inner_outer() { @@ -328,15 +320,15 @@ impl Dcel { self.vertices[handle.index()].data = data; } - pub fn directed_edges(&self) -> DirectedEdgeIterator { + pub fn directed_edges(&self) -> DirectedEdgeIterator<'_, V, DE, UE, F> { DirectedEdgeIterator::new(self) } - pub fn undirected_edges(&self) -> UndirectedEdgeIterator { + pub fn undirected_edges(&self) -> UndirectedEdgeIterator<'_, V, DE, UE, F> { UndirectedEdgeIterator::new(self) } - pub fn vertices(&self) -> VertexIterator { + pub fn vertices(&self) -> VertexIterator<'_, V, DE, UE, F> { VertexIterator::new(self) } @@ -344,11 +336,11 @@ impl Dcel { FixedVertexIterator::new(self.num_vertices()) } - pub fn faces(&self) -> FaceIterator { + pub fn faces(&self) -> FaceIterator<'_, V, DE, UE, F> { FaceIterator::new(self) } - pub fn inner_faces(&self) -> InnerFaceIterator { + pub fn inner_faces(&self) -> InnerFaceIterator<'_, V, DE, UE, F> { let mut iterator = InnerFaceIterator::new(self); iterator.next(); // Skip the outer face iterator @@ -390,7 +382,6 @@ impl Dcel { } } - #[cfg(any(test, fuzzing))] pub fn sanity_check(&self) { if self.num_vertices() <= 1 { assert_eq!(self.num_faces(), 1); diff --git a/src/delaunay_core/dcel_operations.rs b/src/delaunay_core/dcel_operations.rs index b671509..a535f86 100644 --- a/src/delaunay_core/dcel_operations.rs +++ b/src/delaunay_core/dcel_operations.rs @@ -297,7 +297,7 @@ where pub fn create_new_face_adjacent_to_edge( dcel: &mut Dcel, edge: FixedDirectedEdgeHandle, - new_vertex: V, + new_vertex_handle: FixedVertexHandle, ) -> FixedVertexHandle where DE: Default, @@ -327,7 +327,6 @@ where let new_prev_handle = FixedUndirectedEdgeHandle::new(dcel.num_undirected_edges() + 1); let new_face_handle = FixedFaceHandle::::new(dcel.num_faces()); - let new_vertex_handle = FixedVertexHandle::new(dcel.num_vertices()); let new_next = HalfEdgeEntry { next: new_prev_handle.normalized(), @@ -374,10 +373,7 @@ where data: F::default(), }); - dcel.vertices.push(VertexEntry { - data: new_vertex, - out_edge: Some(new_prev_handle.normalized()), - }); + dcel.vertices[new_vertex_handle.index()].out_edge = Some(new_prev_handle.normalized()); *dcel.half_edge_mut(edge) = HalfEdgeEntry { prev: new_prev_handle.normalized(), @@ -397,7 +393,7 @@ where pub fn extend_line( dcel: &mut Dcel, end_vertex: FixedVertexHandle, - new_vertex: V, + new_vertex_handle: FixedVertexHandle, ) -> FixedVertexHandle where DE: Default, @@ -423,7 +419,6 @@ where let new_edge = FixedDirectedEdgeHandle::new_normalized(dcel.num_undirected_edges()); let new_edge_rev = new_edge.rev(); - let new_vertex_handle = FixedVertexHandle::new(dcel.num_vertices()); let face = out_edge.face().fix(); let out_edge = out_edge.fix(); @@ -448,10 +443,7 @@ where }, )); - dcel.vertices.push(VertexEntry { - data: new_vertex, - out_edge: Some(new_edge), - }); + dcel.vertices[new_vertex_handle.index()].out_edge = Some(new_edge); new_vertex_handle } @@ -459,8 +451,8 @@ where pub fn split_edge_when_all_vertices_on_line( dcel: &mut Dcel, edge: FixedDirectedEdgeHandle, - new_vertex: V, -) -> ([FixedDirectedEdgeHandle; 2], FixedVertexHandle) + new_vertex_handle: FixedVertexHandle, +) -> [FixedDirectedEdgeHandle; 2] where DE: Default, UE: Default, @@ -476,7 +468,6 @@ where let rev_prev = rev.prev().fix(); let to = edge.to().fix(); - let new_vertex_handle = FixedVertexHandle::new(dcel.vertices.len()); let face = edge.face().fix(); let edge = edge.fix(); @@ -522,12 +513,9 @@ where }, )); - dcel.vertices.push(VertexEntry { - data: new_vertex, - out_edge: Some(new_edge), - }); + dcel.vertices[new_vertex_handle.index()].out_edge = Some(new_edge); - ([edge, new_edge], new_vertex_handle) + [edge, new_edge] } /// Splits `edge_handle` only one side. Used to split edges on the convex hull. @@ -540,8 +528,8 @@ where pub fn split_half_edge( dcel: &mut Dcel, edge_handle: FixedDirectedEdgeHandle, - new_vertex_data: V, -) -> (FixedVertexHandle, [FixedDirectedEdgeHandle; 2]) + new_vertex_handle: FixedVertexHandle, +) -> [FixedDirectedEdgeHandle; 2] where DE: Default, UE: Default, @@ -607,8 +595,6 @@ where let e2 = FixedDirectedEdgeHandle::new_normalized(dcel.edges.len() + 1); let t2 = e2.rev(); - let nv = FixedVertexHandle::new(dcel.vertices.len()); - let edge1 = HalfEdgeEntry { next: e2, prev: edge_next, @@ -620,13 +606,13 @@ where next: edge_prev, prev: edge_handle, face: f1, - origin: nv, + origin: new_vertex_handle, }; let edge2 = HalfEdgeEntry { next: edge_next, prev: e1, - origin: nv, + origin: new_vertex_handle, face: nf, }; @@ -642,15 +628,10 @@ where data: Default::default(), }; - let new_vertex_entry = VertexEntry { - data: new_vertex_data, - out_edge: Some(e2), - }; - dcel.edges.push(EdgeEntry::new(edge1, twin1)); dcel.edges.push(EdgeEntry::new(edge2, twin2)); dcel.faces.push(new_face); - dcel.vertices.push(new_vertex_entry); + dcel.vertices[new_vertex_handle.index()].out_edge = Some(e2); dcel.half_edge_mut(edge_twin_prev).next = t2; @@ -661,12 +642,12 @@ where dcel.half_edge_mut(edge_handle).next = t1; dcel.half_edge_mut(edge_next).face = nf; - dcel.half_edge_mut(edge_twin).origin = nv; + dcel.half_edge_mut(edge_twin).origin = new_vertex_handle; dcel.vertices[to.index()].out_edge = Some(e2.rev()); dcel.faces[f1.index()].adjacent_edge = Some(edge_handle); - (nv, [edge_handle, e2]) + [edge_handle, e2] } /// Splits `edge_handle`, introducing 6 new half edges, two new faces and one @@ -677,8 +658,8 @@ where pub fn split_edge( dcel: &mut Dcel, edge_handle: FixedDirectedEdgeHandle, - new_vertex: V, -) -> (FixedVertexHandle, [FixedDirectedEdgeHandle; 2]) + new_vertex_handle: FixedVertexHandle, +) -> [FixedDirectedEdgeHandle; 2] where DE: Default, UE: Default, @@ -742,7 +723,7 @@ where let tn = twin.next; let tp = twin.prev; - let v0 = FixedVertexHandle::new(dcel.vertices.len()); + let v0 = new_vertex_handle; let v1 = edge.origin; let v2 = dcel.half_edge(tp).origin; let v3 = twin.origin; @@ -804,11 +785,6 @@ where face: f0, }; - let new_vertex_entry = VertexEntry { - out_edge: Some(t0), - data: new_vertex, - }; - let face2 = FaceEntry { adjacent_edge: Some(e2), data: F::default(), @@ -836,7 +812,7 @@ where dcel.half_edge_mut(tn).next = e1; dcel.half_edge_mut(ep).prev = t3; - dcel.vertices.push(new_vertex_entry); + dcel.vertices[new_vertex_handle.index()].out_edge = Some(t0); dcel.vertices[v3.index()].out_edge = Some(e2); dcel.faces[f0.index()].adjacent_edge = Some(e0); @@ -844,31 +820,30 @@ where dcel.faces.push(face2); dcel.faces.push(face3); - (v0.adjust_inner_outer(), [e0, e2.rev()]) + [e0, e2.rev()] } -pub fn insert_first_vertex( +pub fn append_unconnected_vertex( dcel: &mut Dcel, vertex: V, ) -> FixedVertexHandle { - assert!(dcel.vertices.is_empty()); + let result = FixedVertexHandle::new(dcel.vertices.len()); dcel.vertices.push(VertexEntry { data: vertex, out_edge: None, }); - FixedVertexHandle::new(0) + + result } -pub fn insert_second_vertex( +pub fn setup_initial_two_vertices( dcel: &mut Dcel, - vertex: V, -) -> FixedVertexHandle -where + first_vertex: FixedVertexHandle, + second_vertex: FixedVertexHandle, +) where DE: Default, UE: Default, { - let first_vertex = FixedVertexHandle::new(0); - let second_vertex = FixedVertexHandle::new(1); let normalized = FixedDirectedEdgeHandle::new_normalized(0); let not_normalized = normalized.rev(); @@ -887,18 +862,14 @@ where }, )); - dcel.vertices.push(VertexEntry { - data: vertex, - out_edge: Some(not_normalized), - }); - dcel.vertices[0].out_edge = Some(normalized); + dcel.vertices[first_vertex.index()].out_edge = Some(normalized); + dcel.vertices[second_vertex.index()].out_edge = Some(not_normalized); dcel.faces[OUTER_FACE_HANDLE.index()].adjacent_edge = Some(normalized); - second_vertex } pub fn insert_into_triangle( dcel: &mut Dcel, - vertex: V, + v: FixedVertexHandle, f0: FixedFaceHandle, ) -> FixedVertexHandle where @@ -943,7 +914,6 @@ where let e7 = FixedDirectedEdgeHandle::new_normalized(dcel.edges.len() + 2); let e8 = e7.rev(); - let v = FixedVertexHandle::new(dcel.vertices.len()); let v0 = dcel.half_edge(e0).origin; let v1 = dcel.half_edge(e1).origin; let v2 = dcel.half_edge(e2).origin; @@ -964,11 +934,7 @@ where dcel.faces.push(face1); dcel.faces.push(face2); - let vertex = VertexEntry { - out_edge: Some(e4), - data: vertex, - }; - dcel.vertices.push(vertex); + dcel.vertices[v.index()].out_edge = Some(e4); dcel.half_edge_mut(e0).prev = e8; dcel.half_edge_mut(e0).next = e3; @@ -1045,6 +1011,41 @@ where } } +pub fn new_with_fixed_vertices( + vertices: Vec, + first_vertex: FixedVertexHandle, + second_vertex: FixedVertexHandle, +) -> Dcel +where + DE: Default, + UE: Default, + F: Default, +{ + let vertices_len = vertices.len(); + + let mut result = Dcel { + vertices: vertices + .into_iter() + .map(|v| VertexEntry { + data: v, + out_edge: None, + }) + .collect(), + faces: Vec::with_capacity(vertices_len * 2), + edges: Vec::with_capacity(vertices_len * 3), + }; + + let outer_face = FaceEntry { + adjacent_edge: None, + data: F::default(), + }; + result.faces.push(outer_face); + + setup_initial_two_vertices(&mut result, first_vertex, second_vertex); + + result +} + /// Flip an edge in cw direction pub fn flip_cw(dcel: &mut Dcel, e: FixedUndirectedEdgeHandle) { let e = e.as_directed(); @@ -1318,7 +1319,10 @@ fn remove_when_all_vertices_on_line( #[cfg(test)] mod test { - use crate::handles::{InnerTag, VertexHandle}; + use crate::{ + delaunay_core::dcel_operations::append_unconnected_vertex, + handles::{InnerTag, VertexHandle}, + }; use super::{Dcel, FixedDirectedEdgeHandle, FixedFaceHandle, FixedVertexHandle}; @@ -1430,7 +1434,8 @@ mod test { #[test] fn test_insert_into_triangle() { let mut dcel = default_triangle(); - super::insert_into_triangle(&mut dcel, 3, FixedFaceHandle::new(1)); + let new_vertex = super::append_unconnected_vertex(&mut dcel, 3); + super::insert_into_triangle(&mut dcel, new_vertex, FixedFaceHandle::new(1)); assert_eq!(dcel.faces.len(), 4); assert_eq!(dcel.num_directed_edges(), 12); assert_eq!(dcel.vertices.len(), 4); @@ -1448,8 +1453,9 @@ mod test { let mut dcel = default_triangle(); let face_to_insert = FixedFaceHandle::::new(1); + let v3 = super::append_unconnected_vertex(&mut dcel, 3); let new_vertex = - super::insert_into_triangle(&mut dcel, 3, face_to_insert.adjust_inner_outer()); + super::insert_into_triangle(&mut dcel, v3, face_to_insert.adjust_inner_outer()); let vertex_to_remove = FixedVertexHandle::new(3); let border_loop = get_border_loop(dcel.vertex(new_vertex)); @@ -1468,14 +1474,16 @@ mod test { fn test_insert_into_and_remove_from_quad() { let mut dcel = default_triangle(); - super::insert_into_triangle(&mut dcel, 3, FixedFaceHandle::new(1)); + let v3 = super::append_unconnected_vertex(&mut dcel, 3); + super::insert_into_triangle(&mut dcel, v3, FixedFaceHandle::new(1)); let e_split = dcel .get_edge_from_neighbors(FixedVertexHandle::new(0), FixedVertexHandle::new(3)) .unwrap() .fix(); - let (vertex_to_remove, _) = super::split_edge(&mut dcel, e_split, 3); + let vertex_to_remove = super::append_unconnected_vertex(&mut dcel, 4); + super::split_edge(&mut dcel, e_split, vertex_to_remove); let border_loop = get_border_loop(dcel.vertex(vertex_to_remove)); let mut result = @@ -1491,7 +1499,8 @@ mod test { #[test] fn test_cw_iterator() { let mut dcel = default_triangle(); - super::insert_into_triangle(&mut dcel, 3, FixedFaceHandle::new(1)); + let v3 = super::append_unconnected_vertex(&mut dcel, 3); + super::insert_into_triangle(&mut dcel, v3, FixedFaceHandle::new(1)); let vertex = dcel.vertices().next().unwrap(); assert_eq!(vertex.out_edges().count(), 3); assert_eq!(vertex.out_edges().rev().count(), 3); @@ -1507,7 +1516,8 @@ mod test { #[test] fn test_swap() { let mut dcel = default_triangle(); - super::insert_into_triangle(&mut dcel, 3, FixedFaceHandle::new(1)); + let v3 = super::append_unconnected_vertex(&mut dcel, 3); + super::insert_into_triangle(&mut dcel, v3, FixedFaceHandle::new(1)); let v = FixedVertexHandle::new; for (i0, i1) in [(0, 1), (0, 2), (1, 2), (3, 2), (3, 0), (0, 0), (1, 1)] { dcel.swap_vertices(v(i0), v(i1)); @@ -1518,7 +1528,8 @@ mod test { #[test] fn test_flip() { let mut dcel = default_triangle(); - super::insert_into_triangle(&mut dcel, 3, FixedFaceHandle::new(1)); + let v3 = super::append_unconnected_vertex(&mut dcel, 3); + super::insert_into_triangle(&mut dcel, v3, FixedFaceHandle::new(1)); let e_flip = dcel .get_edge_from_neighbors(FixedVertexHandle::new(0), FixedVertexHandle::new(3)) @@ -1536,7 +1547,8 @@ mod test { fn test_half_edge_split() { let mut dcel = default_triangle(); let edge_handle = FixedDirectedEdgeHandle::new(0); - super::split_half_edge(&mut dcel, edge_handle, 3); + let v3 = super::append_unconnected_vertex(&mut dcel, 3); + super::split_half_edge(&mut dcel, edge_handle, v3); assert_eq!(dcel.num_directed_edges(), 10); assert_eq!(dcel.faces.len(), 3); @@ -1547,14 +1559,17 @@ mod test { #[test] fn test_split() { let mut dcel = default_triangle(); - super::insert_into_triangle(&mut dcel, 3, FixedFaceHandle::new(1)); + let v3 = super::append_unconnected_vertex(&mut dcel, 3); + super::insert_into_triangle(&mut dcel, v3, FixedFaceHandle::new(1)); let e_split = dcel .get_edge_from_neighbors(FixedVertexHandle::new(0), FixedVertexHandle::new(3)) .unwrap() .fix(); dcel.sanity_check(); - super::split_edge(&mut dcel, e_split, 3); + + let v4 = super::append_unconnected_vertex(&mut dcel, 4); + super::split_edge(&mut dcel, e_split, v4); dcel.sanity_check(); assert!(dcel @@ -1568,8 +1583,14 @@ mod test { #[test] fn test_insert_first_and_second() { let mut dcel = Dcel::<_>::default(); - super::insert_first_vertex(&mut dcel, 0); - super::insert_second_vertex(&mut dcel, 1); + super::append_unconnected_vertex(&mut dcel, 0); + super::append_unconnected_vertex(&mut dcel, 1); + + super::setup_initial_two_vertices( + &mut dcel, + FixedVertexHandle::new(0), + FixedVertexHandle::new(1), + ); dcel.sanity_check(); } @@ -1577,14 +1598,16 @@ mod test { #[test] fn test_split_non_triangle_edge() { let mut dcel = Dcel::<_>::default(); - let v0 = super::insert_first_vertex(&mut dcel, 0); - let v1 = super::insert_second_vertex(&mut dcel, 1); + let v0 = super::append_unconnected_vertex(&mut dcel, 0); + let v1 = append_unconnected_vertex(&mut dcel, 1); + super::setup_initial_two_vertices(&mut dcel, v0, v1); dcel.sanity_check(); let edge = dcel.directed_edges().next().unwrap().fix(); // Test split isolated edge - let (_, v2) = super::split_edge_when_all_vertices_on_line(&mut dcel, edge, 2); + let v2 = super::append_unconnected_vertex(&mut dcel, 2); + super::split_edge_when_all_vertices_on_line(&mut dcel, edge, v2); dcel.sanity_check(); @@ -1594,7 +1617,8 @@ mod test { assert!(dcel.get_edge_from_neighbors(v1, v2).is_some()); let non_isolated_edge = dcel.get_edge_from_neighbors(v0, v2).unwrap().fix(); - super::split_edge_when_all_vertices_on_line(&mut dcel, non_isolated_edge, 3); + let v3 = super::append_unconnected_vertex(&mut dcel, 3); + super::split_edge_when_all_vertices_on_line(&mut dcel, non_isolated_edge, v3); dcel.sanity_check(); assert_eq!(dcel.num_directed_edges(), 6); @@ -1606,13 +1630,16 @@ mod test { fn test_extend_line() { let mut dcel = Dcel::<_>::default(); - let v0 = super::insert_first_vertex(&mut dcel, 0); - let v1 = super::insert_second_vertex(&mut dcel, 1); + let v0 = super::append_unconnected_vertex(&mut dcel, 0); + let v1 = append_unconnected_vertex(&mut dcel, 1); + super::setup_initial_two_vertices(&mut dcel, v0, v1); - super::extend_line(&mut dcel, v0, 2); + let v2 = super::append_unconnected_vertex(&mut dcel, 2); + super::extend_line(&mut dcel, v0, v2); dcel.sanity_check(); - super::extend_line(&mut dcel, v1, 3); + let v3 = super::append_unconnected_vertex(&mut dcel, 3); + super::extend_line(&mut dcel, v1, v3); dcel.sanity_check(); assert_eq!(dcel.num_vertices(), 4); @@ -1623,9 +1650,11 @@ mod test { #[test] fn test_create_single_face_between_edge_and_next() { let mut dcel = Dcel::<_>::default(); - super::insert_first_vertex(&mut dcel, 0); - super::insert_second_vertex(&mut dcel, 1); - super::split_edge_when_all_vertices_on_line(&mut dcel, FixedDirectedEdgeHandle::new(0), 2); + let v0 = super::append_unconnected_vertex(&mut dcel, 0); + let v1 = super::append_unconnected_vertex(&mut dcel, 1); + super::setup_initial_two_vertices(&mut dcel, v0, v1); + let v2 = super::append_unconnected_vertex(&mut dcel, 2); + super::split_edge_when_all_vertices_on_line(&mut dcel, FixedDirectedEdgeHandle::new(0), v2); super::create_single_face_between_edge_and_next(&mut dcel, FixedDirectedEdgeHandle::new(0)); @@ -1636,7 +1665,8 @@ mod test { fn test_create_new_face_adjacent_to_edge() { let mut dcel = default_triangle(); - super::create_new_face_adjacent_to_edge(&mut dcel, FixedDirectedEdgeHandle::new(0), 3); + let v3 = super::append_unconnected_vertex(&mut dcel, 3); + super::create_new_face_adjacent_to_edge(&mut dcel, FixedDirectedEdgeHandle::new(0), v3); dcel.sanity_check(); } } diff --git a/src/delaunay_core/handles/handle_defs.rs b/src/delaunay_core/handles/handle_defs.rs index 4fb6df0..346a7b6 100644 --- a/src/delaunay_core/handles/handle_defs.rs +++ b/src/delaunay_core/handles/handle_defs.rs @@ -173,8 +173,7 @@ pub struct DirectedVoronoiEdgeTag; )] #[derive(Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Debug, Default, Hash)] pub struct UndirectedVoronoiEdgeTag; -#[derive(Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Debug, Default, Hash)] -pub struct VoronoiVertexTag; + #[derive(Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Debug, Default, Hash)] pub struct VoronoiFaceTag; diff --git a/src/delaunay_core/hint_generator.rs b/src/delaunay_core/hint_generator.rs index e710b02..57f826d 100644 --- a/src/delaunay_core/hint_generator.rs +++ b/src/delaunay_core/hint_generator.rs @@ -224,6 +224,12 @@ impl HintGenerator vertex: FixedVertexHandle, _vertex_position: Point2, ) { + if self.num_elements_of_base_triangulation == 1 { + // Special case: Triangulation has become empty + *self = Default::default(); + return; + } + let index = vertex.index() as u32; let mut current_divisor = BRANCH_FACTOR; @@ -275,9 +281,32 @@ impl HintGenerator V: HasPosition, { let mut result = Self::default(); - for vertex in triangulation.vertices() { - result.notify_vertex_inserted(vertex.fix(), vertex.position()); + + if triangulation.num_vertices() == 0 { + return result; } + + let mut current_step = BRANCH_FACTOR as usize; + loop { + let vertices = triangulation + .vertices() + .map(|v| v.position()) + .step_by(current_step) + .collect::>(); + + result.hierarchy.push( + DelaunayTriangulation::bulk_load_stable(vertices) + .expect("Failed to load hierarchy layer."), + ); + + if current_step >= triangulation.num_vertices() { + break; + } + + current_step *= BRANCH_FACTOR as usize; + } + + result.num_elements_of_base_triangulation = triangulation.num_vertices(); result } } @@ -306,13 +335,30 @@ mod test { #[test] fn hierarchy_hint_generator_test() -> Result<(), InsertionError> { let vertices = test_utilities::random_points_with_seed(1025, test_utilities::SEED); - let triangulation = HierarchyTriangulation::bulk_load(vertices)?; + let triangulation = HierarchyTriangulation::bulk_load(vertices.clone())?; + hierarchy_sanity_check(&triangulation); + + // Test sequential load + let mut triangulation = HierarchyTriangulation::new(); + for vertex in vertices { + triangulation.insert(vertex)?; + } hierarchy_sanity_check(&triangulation); Ok(()) } fn hierarchy_sanity_check(triangulation: &HierarchyTriangulation) { + let expected_len = if triangulation.num_vertices() <= 1 { + triangulation.num_vertices() + } else { + (triangulation.num_vertices() as f64) + .log(BRANCH_FACTOR as f64) + .ceil() as usize + }; + + assert_eq!(triangulation.hint_generator.hierarchy.len(), expected_len); + for vertex in triangulation.vertices() { let position = vertex.position(); let base_index = vertex.fix().index() as u32; @@ -349,7 +395,7 @@ mod test { for size in 1..20 { let vertices = test_utilities::random_points_with_seed(1 + size * 26, &seed_fn()); - let triangulation = HierarchyTriangulation::bulk_load(vertices)?; + let triangulation = HierarchyTriangulation::bulk_load_stable(vertices)?; hierarchy_sanity_check(&triangulation); } Ok(()) diff --git a/src/delaunay_core/math.rs b/src/delaunay_core/math.rs index 1d82ae3..4f3cab6 100644 --- a/src/delaunay_core/math.rs +++ b/src/delaunay_core/math.rs @@ -385,6 +385,18 @@ where let c = v2.sub(v0); (b.x * c.y - b.y * c.x).abs() * 0.5.into() } +/// Same as `project_point(e1, e2, query_point).is_behind_edge()` but with an added factor against +/// floating point inaccuracies. This doesn't use precise floating point arithmetics but rather errs +/// by incorrectly returning `false` if precision issues are detected. +pub(crate) fn is_behind_edge(e1: Point2, e2: Point2, query_point: Point2) -> bool { + let projection = project_point(e1, e2, query_point); + // There is probably an exact computational method for this. However, I'm not smart enough to + // figure that out. It shouldn't be an issue as this method is allowed to err on the safe side + // (returning false). + // But if you, dear reader, are into exact floating point computation and are looking for an + // exercise: This is your chance! + projection.factor > projection.length_2 + f64::EPSILON +} #[cfg(test)] mod test { diff --git a/src/delaunay_core/mod.rs b/src/delaunay_core/mod.rs index 016de02..56bfdfc 100644 --- a/src/delaunay_core/mod.rs +++ b/src/delaunay_core/mod.rs @@ -15,7 +15,7 @@ pub mod refinement; pub mod interpolation; pub mod math; -pub use bulk_load::{bulk_load, bulk_load_cdt, bulk_load_stable, try_bulk_load_cdt}; +pub use bulk_load::{bulk_load, bulk_load_cdt, try_bulk_load_cdt}; pub use triangulation_ext::{RemovalResult, TriangulationExt}; diff --git a/src/delaunay_core/refinement.rs b/src/delaunay_core/refinement.rs index abbb92e..ac8f67a 100644 --- a/src/delaunay_core/refinement.rs +++ b/src/delaunay_core/refinement.rs @@ -9,8 +9,9 @@ use alloc::vec::Vec; use num_traits::Float; use crate::{ - delaunay_core::math, CdtEdge, ConstrainedDelaunayTriangulation, HasPosition, HintGenerator, - Point2, PositionInTriangulation, SpadeNum, Triangulation, + delaunay_core::{dcel_operations::append_unconnected_vertex, math}, + CdtEdge, ConstrainedDelaunayTriangulation, HasPosition, HintGenerator, Point2, + PositionInTriangulation, SpadeNum, Triangulation, }; use super::{ @@ -829,7 +830,10 @@ where let segment = segment.fix(); let (v0, v1) = (v0.fix(), v1.fix()); - let (new_vertex, [e1, e2]) = self.insert_on_edge(segment, final_position.into()); + let new_vertex = append_unconnected_vertex(self.s_mut(), final_position.into()); + let [e1, e2] = self.insert_on_edge(segment, new_vertex); + self.hint_generator_mut() + .notify_vertex_inserted(new_vertex, final_position); let original_vertices = v0_constraint_vertex .or(v1_constraint_vertex) diff --git a/src/delaunay_core/triangulation_ext.rs b/src/delaunay_core/triangulation_ext.rs index ec82ada..a459398 100644 --- a/src/delaunay_core/triangulation_ext.rs +++ b/src/delaunay_core/triangulation_ext.rs @@ -5,6 +5,8 @@ use super::dcel_operations::IsolateVertexResult; use super::handles::*; use super::math; +use crate::delaunay_core::dcel_operations::append_unconnected_vertex; +use crate::delaunay_core::Dcel; use crate::HintGenerator; use crate::Point2; use crate::{HasPosition, InsertionError, PositionInTriangulation, Triangulation}; @@ -21,6 +23,41 @@ pub enum PositionWhenAllVerticesOnLine { ExtendingLine(FixedVertexHandle), } +#[derive(Copy, Clone, Debug, Ord, PartialOrd, Hash, Eq, PartialEq)] +pub enum VertexToInsert { + NewVertex(V), + ExistingVertex(FixedVertexHandle), +} + +impl VertexToInsert { + pub fn into_vertex(self) -> V { + match self { + VertexToInsert::NewVertex(v) => v, + VertexToInsert::ExistingVertex(_) => panic!("Cannot convert existing vertex to vertex"), + } + } + + pub fn position( + &self, + dcel: &Dcel, + ) -> Point2<::Scalar> + where + V: HasPosition, + { + match self { + VertexToInsert::NewVertex(v) => v.position(), + VertexToInsert::ExistingVertex(handle) => dcel.vertex(*handle).position(), + } + } + + pub fn resolve(self, dcel: &mut Dcel) -> FixedVertexHandle { + match self { + VertexToInsert::NewVertex(v) => append_unconnected_vertex(dcel, v), + VertexToInsert::ExistingVertex(handle) => handle, + } + } +} + impl PositionWhenAllVerticesOnLine { fn to_regular_position_in_triangulation( self, @@ -64,7 +101,7 @@ pub trait TriangulationExt: Triangulation { ) -> Result { math::validate_vertex(&t)?; let position = t.position(); - let result = self.insert_with_hint_option_impl(t, hint); + let result = self.insert_with_hint_option_impl(VertexToInsert::NewVertex(t), hint); self.hint_generator_mut() .notify_vertex_inserted(result, position); @@ -73,16 +110,16 @@ pub trait TriangulationExt: Triangulation { fn insert_with_hint_option_impl( &mut self, - t: Self::Vertex, + t: VertexToInsert, hint: Option, ) -> FixedVertexHandle { use PositionInTriangulation::*; let insertion_result = match self.num_vertices() { - 0 => return dcel_operations::insert_first_vertex(self.s_mut(), t), - 1 => self.insert_second_vertex(t), + 0 => return dcel_operations::append_unconnected_vertex(self.s_mut(), t.into_vertex()), + 1 => self.insert_second_vertex(t.into_vertex()), _ => { - let pos = t.position(); + let pos = t.position(self.s()); if self.all_vertices_on_line() { let location = self.locate_when_all_vertices_on_line(pos); @@ -90,14 +127,19 @@ pub trait TriangulationExt: Triangulation { } else { let position_in_triangulation = self.locate_with_hint_option_core(pos, hint); match position_in_triangulation { - OutsideOfConvexHull(edge) => InsertionResult::NewlyInserted( - self.insert_outside_of_convex_hull(edge, t), - ), + OutsideOfConvexHull(edge) => { + let resolved = t.resolve(self.s_mut()); + InsertionResult::NewlyInserted( + self.insert_outside_of_convex_hull(edge, resolved), + ) + } OnFace(face) => { - InsertionResult::NewlyInserted(self.insert_into_face(face, t)) + let resolved = t.resolve(self.s_mut()); + InsertionResult::NewlyInserted(self.insert_into_face(face, resolved)) } OnEdge(edge) => { - let (new_handle, split_parts) = self.insert_on_edge(edge, t); + let new_handle = t.resolve(self.s_mut()); + let split_parts = self.insert_on_edge(edge, new_handle); if self.is_defined_legal(edge.as_undirected()) { // If the edge is defined as legal the resulting edges must @@ -109,7 +151,7 @@ pub trait TriangulationExt: Triangulation { InsertionResult::NewlyInserted(new_handle) } OnVertex(vertex) => { - self.s_mut().update_vertex(vertex, t); + self.s_mut().update_vertex(vertex, t.into_vertex()); InsertionResult::Updated(vertex) } NoTriangulation => panic!("Error during vertex lookup. This is a bug."), @@ -142,7 +184,7 @@ pub trait TriangulationExt: Triangulation { return NotOnLine(edge.fix().rev()); } - let mut vertices: Vec<_> = self.vertices().collect(); + let mut vertices: Vec<_> = self.vertices().filter(|v| v.out_edge().is_some()).collect(); vertices.sort_by(|left, right| left.position().partial_cmp(&right.position()).unwrap()); let index_to_insert = @@ -176,7 +218,8 @@ pub trait TriangulationExt: Triangulation { return InsertionResult::Updated(first_vertex); } - let second_vertex = dcel_operations::insert_second_vertex(self.s_mut(), vertex); + let second_vertex = dcel_operations::append_unconnected_vertex(self.s_mut(), vertex); + dcel_operations::setup_initial_two_vertices(self.s_mut(), first_vertex, second_vertex); InsertionResult::NewlyInserted(second_vertex) } @@ -184,12 +227,13 @@ pub trait TriangulationExt: Triangulation { fn insert_when_all_vertices_on_line( &mut self, location: PositionWhenAllVerticesOnLine, - new_vertex: Self::Vertex, + new_vertex: VertexToInsert, ) -> InsertionResult { match location { PositionWhenAllVerticesOnLine::OnEdge(edge) => { + let new_vertex = new_vertex.resolve(self.s_mut()); let is_constraint_edge = self.is_defined_legal(edge.as_undirected()); - let (new_edges, new_vertex) = dcel_operations::split_edge_when_all_vertices_on_line( + let new_edges = dcel_operations::split_edge_when_all_vertices_on_line( self.s_mut(), edge, new_vertex, @@ -202,10 +246,12 @@ pub trait TriangulationExt: Triangulation { } PositionWhenAllVerticesOnLine::OnVertex(vertex) => InsertionResult::Updated(vertex), PositionWhenAllVerticesOnLine::NotOnLine(edge) => { + let new_vertex = new_vertex.resolve(self.s_mut()); let result = self.insert_outside_of_convex_hull(edge, new_vertex); InsertionResult::NewlyInserted(result) } PositionWhenAllVerticesOnLine::ExtendingLine(end_vertex) => { + let new_vertex = new_vertex.resolve(self.s_mut()); let result = dcel_operations::extend_line(self.s_mut(), end_vertex, new_vertex); InsertionResult::NewlyInserted(result) } @@ -224,9 +270,9 @@ pub trait TriangulationExt: Triangulation { fn insert_outside_of_convex_hull( &mut self, convex_hull_edge: FixedDirectedEdgeHandle, - new_vertex: Self::Vertex, + new_vertex: FixedVertexHandle, ) -> FixedVertexHandle { - let position = new_vertex.position(); + let position = self.vertex(new_vertex).position(); assert!(self .directed_edge(convex_hull_edge) @@ -284,7 +330,7 @@ pub trait TriangulationExt: Triangulation { fn insert_into_face( &mut self, face: FixedFaceHandle, - t: Self::Vertex, + t: FixedVertexHandle, ) -> FixedVertexHandle { let new_handle = dcel_operations::insert_into_triangle(self.s_mut(), t, face); self.legalize_vertex(new_handle); @@ -294,13 +340,12 @@ pub trait TriangulationExt: Triangulation { fn insert_on_edge( &mut self, edge: FixedDirectedEdgeHandle, - new_vertex: Self::Vertex, - ) -> (FixedVertexHandle, [FixedDirectedEdgeHandle; 2]) { + new_vertex: FixedVertexHandle, + ) -> [FixedDirectedEdgeHandle; 2] { let edge_handle = self.directed_edge(edge); if edge_handle.is_outer_edge() { - let (new_vertex, [e0, e1]) = - dcel_operations::split_half_edge(self.s_mut(), edge.rev(), new_vertex); - (new_vertex, [e1.rev(), e0.rev()]) + let [e0, e1] = dcel_operations::split_half_edge(self.s_mut(), edge.rev(), new_vertex); + [e1.rev(), e0.rev()] } else if edge_handle.rev().is_outer_edge() { dcel_operations::split_half_edge(self.s_mut(), edge, new_vertex) } else { @@ -412,7 +457,7 @@ pub trait TriangulationExt: Triangulation { &self, start: FixedVertexHandle, position: Point2<::Scalar>, - ) -> VertexHandle { + ) -> VertexHandle<'_, Self::Vertex, Self::DirectedEdge, Self::UndirectedEdge, Self::Face> { let start_position = self.vertex(start).position(); if start_position == position { diff --git a/src/delaunay_triangulation.rs b/src/delaunay_triangulation.rs index 7b4f634..81dc9b3 100644 --- a/src/delaunay_triangulation.rs +++ b/src/delaunay_triangulation.rs @@ -1,10 +1,12 @@ use super::delaunay_core::Dcel; use crate::{ - delaunay_core::bulk_load, handles::VertexHandle, HasPosition, HintGenerator, InsertionError, - LastUsedVertexHintGenerator, NaturalNeighbor, Point2, Triangulation, TriangulationExt, + handles::VertexHandle, ConstrainedDelaunayTriangulation, HasPosition, HintGenerator, + InsertionError, LastUsedVertexHintGenerator, NaturalNeighbor, Point2, Triangulation, + TriangulationExt, }; - +use alloc::vec; use alloc::vec::Vec; + use num_traits::Float; #[cfg(feature = "serde")] @@ -292,7 +294,7 @@ where pub fn nearest_neighbor( &self, position: Point2<::Scalar>, - ) -> Option> { + ) -> Option> { if self.num_vertices() == 0 { return None; } @@ -312,15 +314,19 @@ where /// /// # Duplicate handling /// - /// If two vertices have the same position, only one of them will be included in the final - /// triangulation. It is undefined which of them is discarded. + /// For every set of vertices with identical positions, only the vertex with the lowest index + /// is kept. /// - /// For example, if the input vertices are [v0, v1, v2, v1] (where v1 is duplicated), the - /// resulting triangulation will be either - /// [v0, v1, v2] or [v0, v2, v1] + /// For example, if the input vertices are `[v0, v1, v2, v1]` (where v1 is duplicated), the + /// resulting triangulation will be `[v0, v1, v2]`. /// /// Consider checking the triangulation's size after calling this method to ensure that no - /// duplicates were present. + /// duplicates were present. Removing duplicates requires additional work and slightly + /// increases the run time. + /// + /// # Run Time + /// + /// `O(n*log(n))` for `n` input vertices. Slightly (5% - 10%) slower than [Triangulation::bulk_load]. /// /// # Example /// ``` @@ -345,15 +351,23 @@ where /// # Ok(()) } /// ``` pub fn bulk_load_stable(elements: Vec) -> Result { - let mut result: Self = crate::delaunay_core::bulk_load_stable::< - _, - _, - DelaunayTriangulation<_, _, _, _, _>, - _, - >(bulk_load, elements)?; - *result.hint_generator_mut() = L::initialize_from_triangulation(&result); + let cdt = ConstrainedDelaunayTriangulation::bulk_load_cdt(elements, vec![])?; + let result = Self::from_cdt(cdt); + Ok(result) } + + fn from_cdt( + cdt: ConstrainedDelaunayTriangulation, + ) -> DelaunayTriangulation { + let (dcel, hint_generator, num_constraints) = cdt.into_parts(); + let dcel = dcel.map_undirected_edges(|cdt_edge| cdt_edge.deconstruct().1); + assert_eq!(num_constraints, 0); + DelaunayTriangulation { + dcel, + hint_generator, + } + } } impl Default for DelaunayTriangulation @@ -383,7 +397,7 @@ where { /// Allows using natural neighbor interpolation on this triangulation. Refer to the documentation /// of [NaturalNeighbor] for more information. - pub fn natural_neighbor(&self) -> NaturalNeighbor { + pub fn natural_neighbor(&self) -> NaturalNeighbor<'_, Self> { NaturalNeighbor::new(self) } } diff --git a/src/flood_fill_iterator.rs b/src/flood_fill_iterator.rs index d26ffae..1b5291c 100644 --- a/src/flood_fill_iterator.rs +++ b/src/flood_fill_iterator.rs @@ -1,5 +1,9 @@ use alloc::{collections::VecDeque, vec::Vec}; + +#[cfg(not(feature = "std"))] use hashbrown::HashSet; +#[cfg(feature = "std")] +use std::collections::HashSet; use num_traits::{one, zero, Float}; use smallvec::smallvec; diff --git a/src/intersection_iterator.rs b/src/intersection_iterator.rs index 65dd37f..4224e87 100644 --- a/src/intersection_iterator.rs +++ b/src/intersection_iterator.rs @@ -197,7 +197,7 @@ where delaunay: &T, from: FixedVertexHandle, to: FixedVertexHandle, - ) -> LineIntersectionIterator + ) -> LineIntersectionIterator<'_, V, DE, UE, F> where T: Triangulation, { diff --git a/src/triangulation.rs b/src/triangulation.rs index 384e4d9..0d86ab2 100644 --- a/src/triangulation.rs +++ b/src/triangulation.rs @@ -9,6 +9,7 @@ use crate::flood_fill_iterator::RectangleMetric; use crate::flood_fill_iterator::VerticesInShapeIterator; use crate::iterators::*; use crate::Barycentric; +use crate::DelaunayTriangulation; use crate::HintGenerator; use crate::{delaunay_core::Dcel, handles::*}; use crate::{HasPosition, InsertionError, Point2, TriangulationExt}; @@ -169,9 +170,13 @@ pub trait Triangulation: Default { /// more efficient very quickly. #[doc = include_str!("../images/bulk_load_vs_incremental_graph.svg")] fn bulk_load(elements: Vec) -> Result { - let mut result: Self = crate::delaunay_core::bulk_load(elements)?; - *result.hint_generator_mut() = Self::HintGenerator::initialize_from_triangulation(&result); - Ok(result) + let with_incorrect_hint_generator: DelaunayTriangulation<_, _, _, _> = + crate::delaunay_core::bulk_load(elements)?; + let hint_generator = + Self::HintGenerator::initialize_from_triangulation(&with_incorrect_hint_generator); + let (dcel, _, _) = with_incorrect_hint_generator.into_parts(); + + Ok(Self::from_parts(dcel, hint_generator, 0)) } /// Converts a fixed vertex handle to a reference vertex handle. @@ -180,7 +185,7 @@ pub trait Triangulation: Default { fn vertex( &self, handle: FixedVertexHandle, - ) -> VertexHandle { + ) -> VertexHandle<'_, Self::Vertex, Self::DirectedEdge, Self::UndirectedEdge, Self::Face> { self.s().vertex(handle) } @@ -195,8 +200,14 @@ pub trait Triangulation: Default { fn face( &self, handle: FixedFaceHandle, - ) -> FaceHandle - { + ) -> FaceHandle< + '_, + InnerOuter, + Self::Vertex, + Self::DirectedEdge, + Self::UndirectedEdge, + Self::Face, + > { self.s().face(handle) } @@ -204,6 +215,7 @@ pub trait Triangulation: Default { fn outer_face( &self, ) -> FaceHandle< + '_, PossiblyOuterTag, Self::Vertex, Self::DirectedEdge, @@ -219,7 +231,7 @@ pub trait Triangulation: Default { fn directed_edge( &self, handle: FixedDirectedEdgeHandle, - ) -> DirectedEdgeHandle + ) -> DirectedEdgeHandle<'_, Self::Vertex, Self::DirectedEdge, Self::UndirectedEdge, Self::Face> { DirectedEdgeHandle::new(self.s(), handle) } @@ -230,7 +242,7 @@ pub trait Triangulation: Default { fn undirected_edge( &self, handle: FixedUndirectedEdgeHandle, - ) -> UndirectedEdgeHandle + ) -> UndirectedEdgeHandle<'_, Self::Vertex, Self::DirectedEdge, Self::UndirectedEdge, Self::Face> { UndirectedEdgeHandle::new(self.s(), handle) } @@ -286,7 +298,7 @@ pub trait Triangulation: Default { /// The iterator type is [DirectedEdgeHandle]. fn directed_edges( &self, - ) -> DirectedEdgeIterator + ) -> DirectedEdgeIterator<'_, Self::Vertex, Self::DirectedEdge, Self::UndirectedEdge, Self::Face> { self.s().directed_edges() } @@ -296,8 +308,13 @@ pub trait Triangulation: Default { /// The iterator type is [UndirectedEdgeHandle] fn undirected_edges( &self, - ) -> UndirectedEdgeIterator - { + ) -> UndirectedEdgeIterator< + '_, + Self::Vertex, + Self::DirectedEdge, + Self::UndirectedEdge, + Self::Face, + > { self.s().undirected_edges() } @@ -311,7 +328,8 @@ pub trait Triangulation: Default { /// The iterator type is [VertexHandle] fn vertices( &self, - ) -> VertexIterator { + ) -> VertexIterator<'_, Self::Vertex, Self::DirectedEdge, Self::UndirectedEdge, Self::Face> + { self.s().vertices() } @@ -330,7 +348,7 @@ pub trait Triangulation: Default { fn get_vertex( &self, handle: FixedVertexHandle, - ) -> Option> + ) -> Option> { self.s().get_vertex(handle) } @@ -343,7 +361,7 @@ pub trait Triangulation: Default { /// *See also [inner_faces()](Triangulation::inner_faces())* fn all_faces( &self, - ) -> FaceIterator { + ) -> FaceIterator<'_, Self::Vertex, Self::DirectedEdge, Self::UndirectedEdge, Self::Face> { self.s().faces() } @@ -352,7 +370,8 @@ pub trait Triangulation: Default { /// The iterator type is [FaceHandle](FaceHandle). fn inner_faces( &self, - ) -> InnerFaceIterator { + ) -> InnerFaceIterator<'_, Self::Vertex, Self::DirectedEdge, Self::UndirectedEdge, Self::Face> + { self.s().inner_faces() } @@ -361,7 +380,7 @@ pub trait Triangulation: Default { /// The iterator type is [VoronoiFace] fn voronoi_faces( &self, - ) -> VoronoiFaceIterator + ) -> VoronoiFaceIterator<'_, Self::Vertex, Self::DirectedEdge, Self::UndirectedEdge, Self::Face> { VoronoiFaceIterator::new(self.s()) } @@ -372,6 +391,7 @@ pub trait Triangulation: Default { fn directed_voronoi_edges( &self, ) -> DirectedVoronoiEdgeIterator< + '_, Self::Vertex, Self::DirectedEdge, Self::UndirectedEdge, @@ -386,6 +406,7 @@ pub trait Triangulation: Default { fn undirected_voronoi_edges( &self, ) -> UndirectedVoronoiEdgeIterator< + '_, Self::Vertex, Self::DirectedEdge, Self::UndirectedEdge, @@ -418,7 +439,7 @@ pub trait Triangulation: Default { fn locate_vertex( &self, point: Point2<::Scalar>, - ) -> Option> + ) -> Option> { match self.locate(point) { PositionInTriangulation::OnVertex(vertex) => Some(self.vertex(vertex)), @@ -437,7 +458,7 @@ pub trait Triangulation: Default { from: FixedVertexHandle, to: FixedVertexHandle, ) -> Option< - DirectedEdgeHandle, + DirectedEdgeHandle<'_, Self::Vertex, Self::DirectedEdge, Self::UndirectedEdge, Self::Face>, > { self.s().get_edge_from_neighbors(from, to) } @@ -604,7 +625,7 @@ pub trait Triangulation: Default { /// *See also [convex_hull_size](Triangulation::convex_hull_size)* fn convex_hull( &self, - ) -> HullIterator { + ) -> HullIterator<'_, Self::Vertex, Self::DirectedEdge, Self::UndirectedEdge, Self::Face> { { HullIterator::new(self.s()) } @@ -655,7 +676,8 @@ where &self, lower: Point2<::Scalar>, upper: Point2<::Scalar>, - ) -> EdgesInShapeIterator::Scalar>> { + ) -> EdgesInShapeIterator<'_, Self, RectangleMetric<::Scalar>> + { let distance_metric = RectangleMetric::new(lower, upper); let center = lower.add(upper).mul(0.5f32.into()); EdgesInShapeIterator { @@ -686,7 +708,7 @@ where &self, center: Point2<::Scalar>, radius_2: ::Scalar, - ) -> EdgesInShapeIterator::Scalar>> { + ) -> EdgesInShapeIterator<'_, Self, CircleMetric<::Scalar>> { let metric = CircleMetric::new(center, radius_2); EdgesInShapeIterator { inner_iter: FloodFillIterator::new(self, metric, center), @@ -712,7 +734,8 @@ where &self, lower: Point2<::Scalar>, upper: Point2<::Scalar>, - ) -> VerticesInShapeIterator::Scalar>> { + ) -> VerticesInShapeIterator<'_, Self, RectangleMetric<::Scalar>> + { let distance_metric = RectangleMetric::new(lower, upper); let center = lower.add(upper).mul(0.5f32.into()); @@ -741,7 +764,8 @@ where &self, center: Point2<::Scalar>, radius_2: ::Scalar, - ) -> VerticesInShapeIterator::Scalar>> { + ) -> VerticesInShapeIterator<'_, Self, CircleMetric<::Scalar>> + { let distance_metric = CircleMetric::new(center, radius_2); VerticesInShapeIterator::new(FloodFillIterator::new(self, distance_metric, center)) @@ -752,7 +776,7 @@ where /// /// *Note:* In contrast to the other interpolation algorithms, barycentric interpolation also works /// for [crate::ConstrainedDelaunayTriangulation]s. - fn barycentric(&self) -> Barycentric { + fn barycentric(&self) -> Barycentric<'_, Self> { Barycentric::new(self) } }