From afcb6aa70bbd514199a5e24b0d8628d157bcb2c7 Mon Sep 17 00:00:00 2001 From: Stefan Altmayer Date: Sat, 3 Aug 2024 12:37:31 +0200 Subject: [PATCH] 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(()) }