From 6107054967fd11a25a9db3542a57f82fe2ec13af Mon Sep 17 00:00:00 2001 From: Stefan Altmayer Date: Sat, 3 Aug 2024 11:54:41 +0200 Subject: [PATCH 1/3] Adds `ConstrainedDelaunayTriangulation::remove_constraint_edge` --- src/cdt.rs | 46 ++++++++++++++++++++++++++ src/delaunay_core/bulk_load.rs | 2 +- src/delaunay_core/triangulation_ext.rs | 26 ++++++++++----- 3 files changed, 64 insertions(+), 10 deletions(-) diff --git a/src/cdt.rs b/src/cdt.rs index 0304a19..e9c92eb 100644 --- a/src/cdt.rs +++ b/src/cdt.rs @@ -44,6 +44,11 @@ impl CdtEdge { self.0 = true; } + fn unmake_constraint_edge(&mut self) { + assert!(self.is_constraint_edge()); + self.0 = false; + } + /// Returns the wrapped undirected edge data type. pub fn data(&self) -> &UE { &self.1 @@ -739,6 +744,23 @@ where } } + /// Removes a constraint edge. + /// + /// Does nothing and returns `false` if the given edge is not a constraint edge. + /// Otherwise, the edge is unmarked and the Delaunay property is restored in its vicinity. + pub fn remove_constraint_edge(&mut self, edge: FixedUndirectedEdgeHandle) -> bool { + if self.is_constraint_edge(edge) { + self.dcel + .undirected_edge_data_mut(edge) + .unmake_constraint_edge(); + self.num_constraints -= 1; + self.legalize_edge(edge.as_directed(), true); + true + } else { + false + } + } + /// Attempts to add a constraint edge. Leaves the triangulation unchanged if the new edge would /// intersect an existing constraint edge. /// @@ -1813,6 +1835,30 @@ mod test { Ok(()) } + #[test] + fn test_remove_constraint_edge() -> Result<(), InsertionError> { + let mut cdt = get_cdt_for_try_add_constraint()?; + for edge in cdt.fixed_undirected_edges() { + cdt.remove_constraint_edge(edge); + } + assert_eq!(cdt.num_constraints, 0); + cdt.sanity_check(); + + let added_edges = cdt.try_add_constraint( + FixedVertexHandle::from_index(0), + FixedVertexHandle::from_index(1), + ); + assert_eq!(added_edges.len(), 1); + + assert!(cdt.remove_constraint_edge(added_edges.first().unwrap().as_undirected())); + assert_eq!(cdt.num_constraints, 0); + cdt.sanity_check(); + + Ok(()) + } + Ok(()) + } + fn check_returned_edges( cdt: &mut ConstrainedDelaunayTriangulation>, edges: &[FixedDirectedEdgeHandle], diff --git a/src/delaunay_core/bulk_load.rs b/src/delaunay_core/bulk_load.rs index 3643370..db21b23 100644 --- a/src/delaunay_core/bulk_load.rs +++ b/src/delaunay_core/bulk_load.rs @@ -489,7 +489,7 @@ where let cw_walk_start = result.directed_edge(edge).next().rev().fix(); // Check if the edge that was just connected requires legalization - result.legalize_edge(edge); + result.legalize_edge(edge, false); // At this stage the new vertex was successfully inserted. However, insertions like this will end // up in a strongly *star shaped* triangulation instead of a nice nearly-convex blob of faces. diff --git a/src/delaunay_core/triangulation_ext.rs b/src/delaunay_core/triangulation_ext.rs index 34d93c0..8119283 100644 --- a/src/delaunay_core/triangulation_ext.rs +++ b/src/delaunay_core/triangulation_ext.rs @@ -242,7 +242,7 @@ pub trait TriangulationExt: Triangulation { let ccw_walk_start = self.directed_edge(convex_hull_edge).prev().rev().fix(); let cw_walk_start = self.directed_edge(convex_hull_edge).next().rev().fix(); - self.legalize_edge(convex_hull_edge); + self.legalize_edge(convex_hull_edge, false); let mut current_edge = ccw_walk_start; loop { @@ -254,7 +254,7 @@ pub trait TriangulationExt: Triangulation { self.s_mut(), current_edge, ); - self.legalize_edge(current_edge); + self.legalize_edge(current_edge, false); current_edge = new_edge; } else { break; @@ -271,7 +271,7 @@ pub trait TriangulationExt: Triangulation { self.s_mut(), current_edge, ); - self.legalize_edge(next_fix); + self.legalize_edge(next_fix, false); current_edge = new_edge; } else { break; @@ -317,7 +317,7 @@ pub trait TriangulationExt: Triangulation { .collect(); for edge_to_legalize in edges { - self.legalize_edge(edge_to_legalize); + self.legalize_edge(edge_to_legalize, false); } } @@ -337,7 +337,10 @@ pub trait TriangulationExt: Triangulation { /// four sided polygon. /// /// Returns `true` if at least one edge was flipped. This will always include the initial edge. - fn legalize_edge(&mut self, edge: FixedDirectedEdgeHandle) -> bool { + /// + /// Some additional checks are performed if `fully_legalize` is `true`. This is more costly + /// and often not required. + fn legalize_edge(&mut self, edge: FixedDirectedEdgeHandle, fully_legalize: bool) -> bool { let mut edges: SmallVec<[FixedDirectedEdgeHandle; 8]> = Default::default(); edges.push(edge); @@ -379,12 +382,17 @@ pub trait TriangulationExt: Triangulation { result |= should_flip; if should_flip { - let e1 = edge.rev().next().fix(); - let e2 = edge.rev().prev().fix(); + edges.push(edge.rev().next().fix()); + edges.push(edge.rev().prev().fix()); + + if fully_legalize { + // Check all adjacent edges of the input edge. This is more costly but + // can fix Delaunay violations that lie on the left side of this edge. + edges.push(edge.next().fix()); + edges.push(edge.prev().fix()); + } dcel_operations::flip_cw(self.s_mut(), e.as_undirected()); - edges.push(e1); - edges.push(e2); } } } From 94d1172dcc4b6bf9559e3688f81bcdf2397f870b Mon Sep 17 00:00:00 2001 From: Stefan Altmayer Date: Sat, 3 Aug 2024 12:37:31 +0200 Subject: [PATCH 2/3] Bugfix: CDTs could fail to split a constraint edge if the split vertex would be in an unexpected position due to imprecise calculations. This fix adds a slow fallback routine that should be more robust. --- src/cdt.rs | 219 +++++++++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 197 insertions(+), 22 deletions(-) diff --git a/src/cdt.rs b/src/cdt.rs index e9c92eb..fcb142f 100644 --- a/src/cdt.rs +++ b/src/cdt.rs @@ -1,14 +1,17 @@ use alloc::vec; use alloc::vec::Vec; -use num_traits::{zero, Float}; +use num_traits::{Float, NumCast}; #[cfg(feature = "serde")] use serde::{Deserialize, Serialize}; use crate::cdt::GroupEnd::Existing; use crate::delaunay_core::dcel_operations::flip_cw; use crate::delaunay_core::{bulk_load_cdt, bulk_load_stable}; -use crate::{delaunay_core::Dcel, intersection_iterator::LineIntersectionIterator, SpadeNum}; +use crate::{ + delaunay_core::Dcel, intersection_iterator::LineIntersectionIterator, PositionInTriangulation, + SpadeNum, +}; use crate::{handles::*, intersection_iterator::Intersection}; use crate::{ DelaunayTriangulation, HasPosition, HintGenerator, InsertionError, LastUsedVertexHintGenerator, @@ -819,12 +822,75 @@ where // - First, identify potential constraint edge intersections (conflicts). This must be done // beforehand in case that the caller chooses to `ConflictResolution::Cancel` the // operation - no mutation should have happened at this stage. - let initial_conflict_regions = + let (initial_conflict_regions, all_regions_intact) = self.get_conflict_resolutions(from, to, &mut conflict_resolver); // - Second, apply the conflict resolutions, e.g. by inserting new split points and by // rotating non-constraint edges that intersect the new constraint edge (see function // `resolve_conflict_region`). - self.resolve_conflict_groups(to, initial_conflict_regions) + + if all_regions_intact { + self.resolve_conflict_groups(to, initial_conflict_regions) + } else { + self.add_splitting_constraint_edge_fallback(initial_conflict_regions) + } + } + + /// Fallback routine to add splitting constraints that cannot be added via the fast path. + /// + /// This routine simply adds all split vertices first and then adds any missing constraints. + /// This avoids edge cases that can arise when the split point for a constraint + /// intersection lies "too far" off the conflict edge due to imprecise calculations. + fn add_splitting_constraint_edge_fallback( + &mut self, + initial_conflict_regions: Vec>, + ) -> Vec { + let mut vertices_to_connect = Vec::new(); + let mut temporarily_removed = Vec::new(); + + // Phase 1: Add all pending split vertices directly. + for region in initial_conflict_regions { + match region.group_end { + Existing(v) => vertices_to_connect.push(v), + GroupEnd::NewVertex(new_vertex, edge) => { + let new_handle = self + .insert(new_vertex) + .expect("Failed to insert vertex as expected. This is a bug in spade."); + + let [old_from, old_to] = self.directed_edge(edge).vertices().map(|v| v.fix()); + // The conflict edge can prevent the forced insertion to the split vertex. + // It will be removed temporarily + self.remove_constraint_edge(edge.as_undirected()); + + // Re-add the temporarily removed edge later as if it was split by the new + // vertex + temporarily_removed.push([old_from, new_handle]); + temporarily_removed.push([new_handle, old_to]); + vertices_to_connect.push(new_handle); + } + } + } + + let mut result = Vec::new(); + let mut last_vertex = None; + + // Phase 2: Add all constraint edges + for vertex in vertices_to_connect { + if let Some(last_vertex) = last_vertex { + let new_edges = self.try_add_constraint(last_vertex, vertex); + // try_add_constraint should always succeed as any conflicting edge should have been + // removed temporarily + assert_ne!(new_edges, Vec::new()); + result.extend(new_edges); + } + + last_vertex = Some(vertex); + } + + for [from, to] in temporarily_removed { + self.try_add_constraint(from, to); + } + + result } fn get_conflict_resolutions( @@ -832,31 +898,43 @@ where from: FixedVertexHandle, to: FixedVertexHandle, conflict_resolver: &mut R, - ) -> Vec> + ) -> (Vec>, bool) where R: FnMut(DirectedEdgeHandle, F>) -> ConflictResolution, { + let mut all_regions_intact = true; let mut conflict_groups = Vec::new(); let mut current_group = Vec::new(); for intersection in LineIntersectionIterator::new_from_handles(self, from, to) { match intersection { Intersection::EdgeIntersection(edge) => { - if edge.is_constraint_edge() { + if !edge.is_constraint_edge() { + current_group.push(edge.fix()); + } else { match conflict_resolver(edge) { ConflictResolution::Cancel => { - return Vec::new(); + return (Vec::new(), true); } ConflictResolution::Split(new_vertex) => { - let group_end = GroupEnd::NewVertex(new_vertex, edge.fix()); + let position = new_vertex.position(); + let (group_end, is_valid) = + self.verify_split_position(edge, position); + + // A region is considered to be intact if the split vertex lies + // within the region and not outside or on its border. + all_regions_intact &= is_valid; + let conflict_edges = core::mem::take(&mut current_group); + + let group_end = group_end + .unwrap_or(GroupEnd::NewVertex(new_vertex, edge.fix())); + conflict_groups.push(InitialConflictRegion { conflict_edges, group_end, }); } } - } else { - current_group.push(edge.fix()); } } Intersection::VertexIntersection(v) => { @@ -876,7 +954,39 @@ where } } - conflict_groups + (conflict_groups, all_regions_intact) + } + + fn verify_split_position( + &self, + conflict_edge: DirectedEdgeHandle, F>, + split_position: Point2<::Scalar>, + ) -> (Option>, bool) { + // Not every split vertex may lead to a conflict region that will properly contain the + // split vertex. This can happen as not all split positions can be represented precisely. + // + // Instead, these vertices will be handled by a slower fallback routine. + // + // A split position is considered to be valid if it lies either *on* the edge that was split + // or *within any of the edges neighboring faces*. + match self.locate_with_hint(split_position, conflict_edge.from().fix()) { + PositionInTriangulation::OnEdge(real_edge) => { + let is_valid = real_edge.as_undirected() == conflict_edge.fix().as_undirected(); + (None, is_valid) + } + PositionInTriangulation::OnFace(face) => { + let face = face.adjust_inner_outer(); + let is_valid = + face == conflict_edge.face().fix() || face == conflict_edge.rev().face().fix(); + (None, is_valid) + } + PositionInTriangulation::OutsideOfConvexHull(_) => { + let is_valid = conflict_edge.is_part_of_convex_hull(); + (None, is_valid) + } + PositionInTriangulation::OnVertex(v) => (Some(Existing(v)), false), + PositionInTriangulation::NoTriangulation => unreachable!(), + } } fn resolve_conflict_groups( @@ -1029,7 +1139,7 @@ where self.try_add_constraint_inner(from, to, |edge| { let [p0, p1] = edge.positions(); - let line_intersection = line_intersection(p0, p1, from_pos, to_pos); + let line_intersection = get_edge_intersections(p0, p1, from_pos, to_pos); let new_vertex = vertex_constructor(line_intersection); assert_eq!(new_vertex.position(), line_intersection); ConflictResolution::Split(new_vertex) @@ -1056,10 +1166,17 @@ enum ConflictResolution { Split(V), } -fn line_intersection(p1: Point2, p2: Point2, p3: Point2, p4: Point2) -> Point2 -where - S: SpadeNum + Float, -{ +pub fn get_edge_intersections( + p1: Point2, + p2: Point2, + p3: Point2, + p4: Point2, +) -> Point2 { + let p1 = p1.to_f64(); + let p2 = p2.to_f64(); + let p3 = p3.to_f64(); + let p4 = p4.to_f64(); + let a1 = p2.y - p1.y; let b1 = p1.x - p2.x; let c1 = a1 * p1.x + b1 * p1.y; @@ -1070,13 +1187,19 @@ where let determinant = a1 * b2 - a2 * b1; - if determinant == zero() { - panic!("Lines are parallel"); + let x: f64; + let y: f64; + if determinant == 0.0 { + x = f64::infinity(); + y = f64::infinity(); + } else { + x = (b2 * c1 - b1 * c2) / determinant; + y = (a1 * c2 - a2 * c1) / determinant; } - let x = (b2 * c1 - b1 * c2) / determinant; - let y = (a1 * c2 - a2 * c1) / determinant; - Point2::new(x, y) + [x, y] + .map(|s| ::from(s).unwrap_or_else(|| (s as f32).into())) + .into() } #[cfg(test)] @@ -1086,7 +1209,7 @@ mod test { use rand::distributions::{Distribution, Uniform}; use rand::{Rng, SeedableRng}; - use crate::delaunay_core::FixedDirectedEdgeHandle; + use crate::delaunay_core::{FixedDirectedEdgeHandle, TriangulationExt}; use crate::handles::FixedVertexHandle; use crate::test_utilities::*; use crate::{DelaunayTriangulation, InsertionError, Point2, Triangulation}; @@ -1856,6 +1979,58 @@ mod test { Ok(()) } + + #[test] + fn edge_intersection_precision_test() -> Result<(), InsertionError> { + let edges = [ + [ + Point2::new(17.064112, -17.96008), + Point2::new(16.249594, -17.145563), + ], + [ + Point2::new(-25.290726, -24.435482), + Point2::new(-5.6608872, -24.435482), + ], + [ + Point2::new(17.878626, -18.774595), + Point2::new(15.435078, -16.331045), + ], + ]; + let mut cdt: ConstrainedDelaunayTriangulation> = + ConstrainedDelaunayTriangulation::new(); + + for edge in edges.iter() { + let point_a = cdt.insert(edge[0])?; + let point_b = cdt.insert(edge[1])?; + + // The intersection calculation of the last edge is susceptible to floating point + // inaccuracies. Spade has a fallback routine that is more costly but should handle + // these more robustly. This test is set up to trigger this routine. + cdt.add_constraint_and_split(point_a, point_b, |v| v); + cdt.cdt_sanity_check(); + } + + assert_eq!(cdt.num_vertices(), 7); + + // Gather all constraint edges as [from, to] index tuples + let mut constraint_edges = cdt + .undirected_edges() + .filter(|e| e.is_constraint_edge()) + .map(|e| e.vertices().map(|v| v.index())) + .collect::>(); + + // Normalize to make comparison order-independent + for edge_pair in &mut constraint_edges { + edge_pair.sort(); + } + constraint_edges.sort(); + + // Manually checked for correctness... + assert_eq!( + constraint_edges, + vec![[0, 6], [1, 6], [2, 3], [4, 6], [5, 6]] + ); + Ok(()) } From a9ffb06035fb0e906dc6548c8e103b2e3d8d03d0 Mon Sep 17 00:00:00 2001 From: Stefan Altmayer Date: Sat, 3 Aug 2024 12:41:25 +0200 Subject: [PATCH 3/3] Updates CHANGELOG for 2.11 release. --- CHANGELOG.md | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index b82f1c7..3bf6e15 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,16 @@ 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.11.0] - 2024-08-03 + +### Added + +- Added `ConstrainedDelaunayTriangulation::remove_constraint_edge(edge)` + +### Fix + +- Fixes potential crash when using `ConstrainedDelaunayTriangulation::add_constraint_and_split`. See #109 + ## [2.10.0] - 2024-07-25 ### Added @@ -447,6 +457,8 @@ A lot has changed for the 1.0. release, only larger changes are shown. Initial commit +[2.11.0]: https://github.com/Stoeoef/spade/compare/v2.10.0...v2.11.0 + [2.10.0]: https://github.com/Stoeoef/spade/compare/v2.9.0...v2.10.0 [2.9.0]: https://github.com/Stoeoef/spade/compare/v2.8.0...v2.9.0