From e4a7be5f12c1a2d5f445cd1fcdaab064c4553199 Mon Sep 17 00:00:00 2001 From: Stefan Altmayer Date: Sun, 14 Jan 2024 22:04:22 +0100 Subject: [PATCH] Implements `DelaunayTriangulation::bulk_load_stable` --- CHANGELOG.md | 8 ++ src/cdt.rs | 53 +++++--- src/delaunay_core/bulk_load.rs | 185 ++++++++++++++++++++++++++- src/delaunay_core/dcel.rs | 50 ++++++++ src/delaunay_core/dcel_operations.rs | 13 +- src/delaunay_core/mod.rs | 2 +- src/delaunay_triangulation.rs | 71 +++++++++- src/triangulation.rs | 16 +++ 8 files changed, 376 insertions(+), 22 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index e6e08a5..f2bc928 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,12 @@ 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.6.0] - 2024-01-14 + +### Added + - Implements `DelaunayTriangulation::bulk_load_stable` (See #77) + - Implements `FixedVertexHandle::from_index(usize)` + ## [2.5.1] - 2023-12-27 ### Fix @@ -302,6 +308,8 @@ A lot has changed for the 1.0. release, only larger changes are shown. ## 0.1.0 - 2016-09-23 Initial commit +[2.6.0]: https://github.com/Stoeoef/spade/compare/v2.5.1...v2.6.0 + [2.5.1]: https://github.com/Stoeoef/spade/compare/v2.5.0...v2.5.1 [2.5.0]: https://github.com/Stoeoef/spade/compare/v2.4.1...v2.5.0 diff --git a/src/cdt.rs b/src/cdt.rs index 0412e2c..7980a49 100644 --- a/src/cdt.rs +++ b/src/cdt.rs @@ -144,9 +144,9 @@ pub struct ConstrainedDelaunayTriangulation< F: Default, L: HintGenerator<::Scalar>, { - s: Dcel, F>, + dcel: Dcel, F>, num_constraints: usize, - lookup: L, + hint_generator: L, } impl Default for ConstrainedDelaunayTriangulation @@ -159,9 +159,9 @@ where { fn default() -> Self { ConstrainedDelaunayTriangulation { - s: Default::default(), + dcel: Default::default(), num_constraints: 0, - lookup: Default::default(), + hint_generator: Default::default(), } } } @@ -181,11 +181,11 @@ where type HintGenerator = L; fn s(&self) -> &Dcel, F> { - &self.s + &self.dcel } fn s_mut(&mut self) -> &mut Dcel, F> { - &mut self.s + &mut self.dcel } fn is_defined_legal(&self, edge: FixedUndirectedEdgeHandle) -> bool { @@ -196,7 +196,7 @@ where self.num_constraints += 1; for handle in handles.iter().map(|e| e.as_undirected()) { if !self.is_constraint_edge(handle) { - self.s + self.dcel .undirected_edge_data_mut(handle) .make_constraint_edge(); } @@ -204,11 +204,11 @@ where } fn hint_generator(&self) -> &Self::HintGenerator { - &self.lookup + &self.hint_generator } fn hint_generator_mut(&mut self) -> &mut Self::HintGenerator { - &mut self.lookup + &mut self.hint_generator } fn clear(&mut self) { @@ -217,6 +217,27 @@ where let new_hint_generator = HintGenerator::initialize_from_triangulation(self); *self.hint_generator_mut() = new_hint_generator; } + + fn from_parts( + dcel: Dcel, + hint_generator: Self::HintGenerator, + num_constraints: usize, + ) -> Self { + Self { + dcel, + num_constraints, + hint_generator, + } + } + + fn into_parts( + self, + ) -> ( + Dcel, + Self::HintGenerator, + ) { + todo!() + } } impl From> @@ -234,9 +255,9 @@ where let lookup = value.hint_generator; ConstrainedDelaunayTriangulation { - s, + dcel: s, num_constraints: 0, - lookup, + hint_generator: lookup, } } } @@ -258,7 +279,7 @@ where /// This method will invalidate all vertex, edge and face handles. pub fn remove(&mut self, vertex: FixedVertexHandle) -> V { let num_removed_constraints = self - .s + .dcel .vertex(vertex) .out_edges() .map(|edge| edge.is_constraint_edge()) @@ -275,7 +296,7 @@ where /// Returns `true` if a given edge is a constraint edge. pub fn is_constraint_edge(&self, edge: FixedUndirectedEdgeHandle) -> bool { - self.s.undirected_edge_data(edge).is_constraint_edge() + self.dcel.undirected_edge_data(edge).is_constraint_edge() } /// Checks if two vertices are connected by a constraint edge. @@ -591,7 +612,9 @@ where fn make_constraint_edge(&mut self, edge: FixedUndirectedEdgeHandle) -> bool { if !self.is_constraint_edge(edge) { - self.s.undirected_edge_data_mut(edge).make_constraint_edge(); + self.dcel + .undirected_edge_data_mut(edge) + .make_constraint_edge(); self.num_constraints += 1; true } else { @@ -602,7 +625,7 @@ where #[cfg(test)] pub fn cdt_sanity_check(&self) { let num_undirected_edges = self - .s + .dcel .undirected_edges() .filter(|e| e.is_constraint_edge()) .count(); diff --git a/src/delaunay_core/bulk_load.rs b/src/delaunay_core/bulk_load.rs index 63d4f0a..7fb181e 100644 --- a/src/delaunay_core/bulk_load.rs +++ b/src/delaunay_core/bulk_load.rs @@ -1,8 +1,9 @@ -use core::cmp::{Ordering, Reverse}; - use crate::{HasPosition, InsertionError, Point2, Triangulation, TriangulationExt}; +use core::cmp::{Ordering, Reverse}; -use super::{dcel_operations, FixedDirectedEdgeHandle, FixedUndirectedEdgeHandle}; +use super::{ + dcel_operations, FixedDirectedEdgeHandle, FixedUndirectedEdgeHandle, FixedVertexHandle, +}; use alloc::vec::Vec; @@ -117,7 +118,6 @@ where } // Get new center that is guaranteed to be within the convex hull - // let center_positions = || { result .vertices() @@ -170,6 +170,111 @@ where Ok(result) } +pub struct PointWithIndex { + data: V, + index: usize, +} + +impl HasPosition for PointWithIndex +where + V: HasPosition, +{ + type Scalar = V::Scalar; + fn position(&self) -> Point2 { + self.data.position() + } +} + +pub fn bulk_load_stable(elements: Vec) -> Result +where + V: HasPosition, + T: Triangulation, + T2: Triangulation< + Vertex = PointWithIndex, + DirectedEdge = T::DirectedEdge, + UndirectedEdge = T::UndirectedEdge, + Face = T::Face, + HintGenerator = T::HintGenerator, + >, +{ + let elements = elements + .into_iter() + .enumerate() + .map(|(index, data)| PointWithIndex { index, data }) + .collect::>(); + + let num_original_elements = elements.len(); + + let mut with_indices = bulk_load::, T2>(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 + }); + + // 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; + } + } + + // 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 { + // 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; + } else { + with_indices + .s_mut() + .swap_vertices(FixedVertexHandle::new(old_index), new_index); + } + + if current_index >= with_indices.num_vertices() { + break; + } + } + + let (dcel, hint_generator) = with_indices.into_parts(); + let dcel = dcel.map_vertices(|point_with_index| point_with_index.data); + + Ok(T::from_parts(dcel, hint_generator, 0)) +} + #[inline(never)] // Prevent inlining for better profiling data fn single_bulk_insertion_step( result: &mut TR, @@ -825,6 +930,63 @@ mod test { Ok(()) } + #[test] + fn test_bulk_load_stable() -> Result<(), InsertionError> { + const SIZE: usize = 200; + let mut vertices = random_points_with_seed(SIZE, SEED2); + + vertices.push(Point2::new(4.0, 4.0)); + vertices.push(Point2::new(4.0, -4.0)); + vertices.push(Point2::new(-4.0, 4.0)); + vertices.push(Point2::new(-4.0, -4.0)); + + vertices.push(Point2::new(5.0, 5.0)); + vertices.push(Point2::new(5.0, -5.0)); + vertices.push(Point2::new(-5.0, 5.0)); + vertices.push(Point2::new(-5.0, -5.0)); + + vertices.push(Point2::new(6.0, 6.0)); + vertices.push(Point2::new(6.0, -6.0)); + vertices.push(Point2::new(-6.0, 6.0)); + vertices.push(Point2::new(-6.0, -6.0)); + + let num_vertices = vertices.len(); + + let triangulation = DelaunayTriangulation::<_>::bulk_load_stable(vertices.clone())?; + triangulation.sanity_check(); + assert_eq!(triangulation.num_vertices(), num_vertices); + + for (inserted, original) in triangulation.vertices().zip(vertices) { + assert_eq!(inserted.data(), &original); + } + + triangulation.sanity_check(); + + Ok(()) + } + + #[test] + fn test_bulk_load_stable_with_duplicates() -> Result<(), InsertionError> { + const SIZE: usize = 200; + let mut vertices = random_points_with_seed(SIZE, SEED2); + let original = vertices.clone(); + let duplicates = vertices.iter().copied().take(SIZE / 10).collect::>(); + for (index, d) in duplicates.into_iter().enumerate() { + vertices.insert(index * 2, d); + } + + let triangulation = DelaunayTriangulation::<_>::bulk_load_stable(vertices)?; + triangulation.sanity_check(); + assert_eq!(triangulation.num_vertices(), SIZE); + + for (inserted, original) in triangulation.vertices().zip(original) { + assert_eq!(inserted.data(), &original); + } + + triangulation.sanity_check(); + Ok(()) + } + #[test] fn test_bulk_load() -> Result<(), InsertionError> { const SIZE: usize = 9000; @@ -853,6 +1015,21 @@ mod test { Ok(()) } + #[test] + fn test_same_vertex_bulk_load() -> Result<(), InsertionError> { + const SIZE: usize = 100; + let mut vertices = random_points_with_seed(SIZE, SEED2); + + for i in 0..SIZE - 5 { + vertices.insert(i * 2, Point2::new(0.5, 0.2)); + } + + let triangulation = DelaunayTriangulation::>::bulk_load(vertices)?; + triangulation.sanity_check(); + assert_eq!(triangulation.num_vertices(), SIZE + 1); + Ok(()) + } + #[test] fn test_hull() -> Result<(), InsertionError> { let mut triangulation = DelaunayTriangulation::<_>::new(); diff --git a/src/delaunay_core/dcel.rs b/src/delaunay_core/dcel.rs index e643c35..17850bb 100644 --- a/src/delaunay_core/dcel.rs +++ b/src/delaunay_core/dcel.rs @@ -133,6 +133,24 @@ impl Dcel { self.faces.truncate(1); // Keep outer face } + pub fn map_vertices(self, f: M) -> Dcel + where + M: Fn(V) -> V2, + { + Dcel { + vertices: self + .vertices + .into_iter() + .map(|vertex_data| VertexEntry { + data: f(vertex_data.data), + out_edge: vertex_data.out_edge, + }) + .collect(), + faces: self.faces, + edges: self.edges, + } + } + pub fn map_undirected_edges(self, f: M) -> Dcel where M: Fn(UE) -> UE2, @@ -336,6 +354,38 @@ impl Dcel { FixedFaceIterator::new(self.num_faces()) } + pub fn swap_vertices(&mut self, v0: FixedVertexHandle, v1: FixedVertexHandle) { + let out_0 = self.vertices[v0.index()].out_edge; + let out_1 = self.vertices[v1.index()].out_edge; + self.vertices.swap(v0.index(), v1.index()); + + if let Some(mut out_0) = out_0 { + loop { + let out_entry = self.half_edge_mut(out_0); + if out_entry.origin == v1 { + break; + } + + out_entry.origin = v1; + + out_0 = out_entry.prev.rev(); + } + } + + if let Some(mut out_1) = out_1 { + loop { + let out_entry = self.half_edge_mut(out_1); + if out_entry.origin == v0 { + break; + } + + out_entry.origin = v0; + + out_1 = out_entry.prev.rev(); + } + } + } + #[cfg(any(test, fuzzing))] pub fn sanity_check(&self) { if self.num_vertices() <= 1 { diff --git a/src/delaunay_core/dcel_operations.rs b/src/delaunay_core/dcel_operations.rs index 57e0f0d..3c56542 100644 --- a/src/delaunay_core/dcel_operations.rs +++ b/src/delaunay_core/dcel_operations.rs @@ -311,7 +311,7 @@ where // After: // // new_prev.rev - // | ^ + // | ^ <---- new // | / \ // | / \ // V/ new \<---new_next.rev @@ -1504,6 +1504,17 @@ mod test { assert_eq!(out_edges, reversed); } + #[test] + fn test_swap() { + let mut dcel = default_triangle(); + super::insert_into_triangle(&mut dcel, 3, 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)); + dcel.sanity_check(); + } + } + #[test] fn test_flip() { let mut dcel = default_triangle(); diff --git a/src/delaunay_core/mod.rs b/src/delaunay_core/mod.rs index 800537b..75bb6dc 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; +pub use bulk_load::{bulk_load, bulk_load_stable}; pub use triangulation_ext::{RemovalResult, TriangulationExt}; diff --git a/src/delaunay_triangulation.rs b/src/delaunay_triangulation.rs index c062a6d..f5c233d 100644 --- a/src/delaunay_triangulation.rs +++ b/src/delaunay_triangulation.rs @@ -1,9 +1,10 @@ use super::delaunay_core::Dcel; use crate::{ - handles::VertexHandle, HasPosition, HintGenerator, LastUsedVertexHintGenerator, + handles::VertexHandle, HasPosition, HintGenerator, InsertionError, LastUsedVertexHintGenerator, NaturalNeighbor, Point2, Triangulation, TriangulationExt, }; +use alloc::vec::Vec; use num_traits::Float; #[cfg(feature = "serde")] @@ -303,6 +304,53 @@ where self.hint_generator().notify_vertex_lookup(vertex.fix()); Some(vertex) } + + /// Creates a new delaunay triangulation with an efficient bulk loading strategy. + /// + /// In contrast to [Triangulation::bulk_load], this method will create a triangulation with + /// vertices returned *in the same order* as the input vertices. + /// + /// # 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 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] + /// + /// Consider checking the triangulation's size after calling this method to ensure that no + /// duplicates were present. + /// + /// # Example + /// ``` + /// # use spade::InsertionError; + /// use spade::{DelaunayTriangulation, Point2, Triangulation}; + /// + /// # fn main() -> Result<(), InsertionError> { + /// let vertices = vec![ + /// Point2::new(0.5, 1.0), + /// Point2::new(-0.5, 2.0), + /// Point2::new(0.2, 3.0), + /// Point2::new(0.0, 4.0), + /// Point2::new(-0.2, 5.0) + /// ]; + /// let triangulation = DelaunayTriangulation::>::bulk_load_stable(vertices.clone())?; + /// // This assert will not hold for regular bulk loading! + /// assert_eq!(triangulation.vertices().map(|v| *v.data()).collect::>(), vertices); + /// + /// // This is how you would check for duplicates: + /// let duplicates_found = triangulation.num_vertices() < vertices.len(); + /// assert!(!duplicates_found); + /// # Ok(()) } + /// ``` + pub fn bulk_load_stable(elements: Vec) -> Result { + let result: Self = + crate::delaunay_core::bulk_load_stable::<_, _, DelaunayTriangulation<_, _, _, _, _>>( + elements, + )?; + Ok(result) + } } impl Default for DelaunayTriangulation @@ -366,6 +414,27 @@ where fn hint_generator_mut(&mut self) -> &mut Self::HintGenerator { &mut self.hint_generator } + + fn from_parts( + dcel: Dcel, + hint_generator: Self::HintGenerator, + num_constraints: usize, + ) -> Self { + assert_eq!(num_constraints, 0); + Self { + dcel, + hint_generator, + } + } + + fn into_parts( + self, + ) -> ( + Dcel, + Self::HintGenerator, + ) { + (self.dcel, self.hint_generator) + } } #[cfg(test)] diff --git a/src/triangulation.rs b/src/triangulation.rs index 63218f7..bc8000f 100644 --- a/src/triangulation.rs +++ b/src/triangulation.rs @@ -81,6 +81,22 @@ pub trait Triangulation: Default { #[doc(hidden)] fn hint_generator_mut(&mut self) -> &mut Self::HintGenerator; + #[doc(hidden)] + fn from_parts( + dcel: Dcel, + hint_generator: Self::HintGenerator, + num_constraints: usize, + ) -> Self; + + #[doc(hidden)] + #[allow(clippy::type_complexity)] + fn into_parts( + self, + ) -> ( + Dcel, + Self::HintGenerator, + ); + /// Creates a new triangulation. /// /// A newly created triangulation contains no vertices, no edges and the single