From d7b0685261266e7233c9d050c79a8975e25839a5 Mon Sep 17 00:00:00 2001 From: Hanno Braun Date: Sat, 9 Apr 2022 09:56:21 +0200 Subject: [PATCH 1/5] Add dependency on `robust` to `fj-kernel` --- Cargo.lock | 1 + fj-kernel/Cargo.toml | 1 + 2 files changed, 2 insertions(+) diff --git a/Cargo.lock b/Cargo.lock index 75495a9bb..60b9240a6 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -717,6 +717,7 @@ dependencies = [ "parking_lot 0.12.0", "parry2d-f64", "parry3d-f64", + "robust", "slotmap", "spade", "thiserror", diff --git a/fj-kernel/Cargo.toml b/fj-kernel/Cargo.toml index 3db1e76ff..c52e0074c 100644 --- a/fj-kernel/Cargo.toml +++ b/fj-kernel/Cargo.toml @@ -20,6 +20,7 @@ nalgebra = "0.30.0" parking_lot = "0.12.0" parry2d-f64 = "0.8.0" parry3d-f64 = "0.8.0" +robust = "0.2.3" slotmap = "1.0.6" spade = "2.0.0" thiserror = "1.0.30" From 8cf7ea1520e7e22771c0d50fcc9e4e348a1d9340 Mon Sep 17 00:00:00 2001 From: Hanno Braun Date: Fri, 8 Apr 2022 18:30:35 +0200 Subject: [PATCH 2/5] Add `HorizontalRayToTheRight` This struct is specifically made for the point-in-polygon test. It it simpler than the currently used ray from Parry, allowing it to be more robust, and most importantly, provide exactly the feedback that the point-in-polygon test needs in a robust manner. --- fj-kernel/src/algorithms/triangulation/mod.rs | 1 + fj-kernel/src/algorithms/triangulation/ray.rs | 163 ++++++++++++++++++ 2 files changed, 164 insertions(+) create mode 100644 fj-kernel/src/algorithms/triangulation/ray.rs diff --git a/fj-kernel/src/algorithms/triangulation/mod.rs b/fj-kernel/src/algorithms/triangulation/mod.rs index 611380552..3666c0c19 100644 --- a/fj-kernel/src/algorithms/triangulation/mod.rs +++ b/fj-kernel/src/algorithms/triangulation/mod.rs @@ -1,4 +1,5 @@ mod polygon; +mod ray; use fj_debug::DebugInfo; use fj_math::{Scalar, Triangle}; diff --git a/fj-kernel/src/algorithms/triangulation/ray.rs b/fj-kernel/src/algorithms/triangulation/ray.rs new file mode 100644 index 000000000..23adc6598 --- /dev/null +++ b/fj-kernel/src/algorithms/triangulation/ray.rs @@ -0,0 +1,163 @@ +#![cfg(test)] + +use fj_math::{Point, Segment}; + +pub struct HorizontalRayToTheRight { + pub origin: Point<2>, +} + +impl HorizontalRayToTheRight { + pub fn hits_segment(&self, segment: impl Into>) -> Option { + let [a, b] = segment.into().points(); + let [lower, upper] = if a.v <= b.v { [a, b] } else { [b, a] }; + let right = if a.u > b.u { a } else { b }; + + if self.origin.v > upper.v { + // ray is above segment + return None; + } + if self.origin.v < lower.v { + // ray is below segment + return None; + } + + if self.origin.v == lower.v && lower.v == upper.v { + // ray and segment are parallel and at same height + + if self.origin.u > right.u { + return None; + } + + return Some(Hit::Parallel); + } + + let pa = robust::Coord { + x: lower.u, + y: lower.v, + }; + let pb = robust::Coord { + x: upper.u, + y: upper.v, + }; + let pc = robust::Coord { + x: self.origin.u, + y: self.origin.v, + }; + + if robust::orient2d(pa, pb, pc) >= 0. { + // ray starts on the line or left of it + + if self.origin.v == upper.v { + return Some(Hit::UpperVertex); + } + if self.origin.v == lower.v { + return Some(Hit::LowerVertex); + } + + return Some(Hit::Segment); + } + + None + } +} + +impl

From

for HorizontalRayToTheRight +where + P: Into>, +{ + fn from(point: P) -> Self { + Self { + origin: point.into(), + } + } +} + +#[derive(Clone, Copy, Debug, Eq, PartialEq)] +pub enum Hit { + Segment, + + LowerVertex, + UpperVertex, + + Parallel, +} + +#[cfg(test)] +mod tests { + use super::{Hit, HorizontalRayToTheRight}; + + #[test] + fn hits_segment_right() { + let ray = HorizontalRayToTheRight::from([0., 2.]); + + let below = [[1., 0.], [1., 1.]]; + let above = [[1., 3.], [1., 4.]]; + let same_level = [[1., 1.], [1., 3.]]; + + assert!(ray.hits_segment(below).is_none()); + assert!(ray.hits_segment(above).is_none()); + assert!(matches!(ray.hits_segment(same_level), Some(Hit::Segment))); + } + + #[test] + fn hits_segment_left() { + let ray = HorizontalRayToTheRight::from([1., 2.]); + + let same_level = [[0., 1.], [0., 3.]]; + assert!(ray.hits_segment(same_level).is_none()); + } + + #[test] + fn hits_segment_overlapping() { + let ray = HorizontalRayToTheRight::from([1., 1.]); + + let no_hit = [[0., 0.], [2., 3.]]; + + let hit_segment = [[0., 0.], [3., 2.]]; + let hit_upper = [[0., 0.], [2., 1.]]; + let hit_lower = [[0., 2.], [2., 1.]]; + + assert!(ray.hits_segment(no_hit).is_none()); + assert!(matches!(ray.hits_segment(hit_segment), Some(Hit::Segment))); + assert!(matches!( + ray.hits_segment(hit_upper), + Some(Hit::UpperVertex), + )); + assert!(matches!( + ray.hits_segment(hit_lower), + Some(Hit::LowerVertex), + )); + } + + #[test] + fn hits_segment_on_segment() { + let ray = HorizontalRayToTheRight::from([1., 1.]); + + let hit_segment = [[0., 0.], [2., 2.]]; + let hit_upper = [[0., 0.], [1., 1.]]; + let hit_lower = [[1., 1.], [2., 2.]]; + + assert!(matches!(ray.hits_segment(hit_segment), Some(Hit::Segment))); + assert!(matches!( + ray.hits_segment(hit_upper), + Some(Hit::UpperVertex), + )); + assert!(matches!( + ray.hits_segment(hit_lower), + Some(Hit::LowerVertex), + )); + } + + #[test] + fn hits_segment_parallel() { + let ray = HorizontalRayToTheRight::from([2., 0.]); + + let left = [[0., 0.], [1., 0.]]; + let overlapping = [[1., 0.], [3., 0.]]; + let right = [[3., 0.], [4., 0.]]; + + assert!(ray.hits_segment(left).is_none()); + assert!(matches!(ray.hits_segment(overlapping), Some(Hit::Parallel))); + assert!(matches!(ray.hits_segment(right), Some(Hit::Parallel))); + } +} From 2b56f016ceb244910e786235beb3191c0c84b153 Mon Sep 17 00:00:00 2001 From: Hanno Braun Date: Sat, 9 Apr 2022 12:16:23 +0200 Subject: [PATCH 3/5] Fix most edge cases in point-in-polygon test This leaves one known edge case unfixed, and that's passing through a vertex at the seam of the polygon in the presence of interior polygon chains. I'm going to fix that in a follow-up commit, as that's going to require some additional cleanup. --- .../src/algorithms/triangulation/polygon.rs | 204 +++++++++++++----- fj-kernel/src/algorithms/triangulation/ray.rs | 2 - 2 files changed, 152 insertions(+), 54 deletions(-) diff --git a/fj-kernel/src/algorithms/triangulation/polygon.rs b/fj-kernel/src/algorithms/triangulation/polygon.rs index 9d77df150..6ab8c5e59 100644 --- a/fj-kernel/src/algorithms/triangulation/polygon.rs +++ b/fj-kernel/src/algorithms/triangulation/polygon.rs @@ -1,15 +1,13 @@ -use std::collections::BTreeSet; - use fj_debug::{DebugInfo, TriangleEdgeCheck}; -use fj_math::{Point, PolyChain, Scalar, Segment}; -use parry2d_f64::query::{Ray as Ray2, RayCast as _}; +use fj_math::{Point, PolyChain, Segment}; use crate::geometry::Surface; +use super::ray::{Hit, HorizontalRayToTheRight}; + pub struct Polygon { surface: Surface, edges: Vec>, - max: Point<2>, } impl Polygon { @@ -17,7 +15,6 @@ impl Polygon { Self { surface, edges: Vec::new(), - max: Point::origin(), } } @@ -38,21 +35,18 @@ impl Polygon { fn with_approx(mut self, approx: impl Into>) -> Self { for segment in approx.into().segments() { - let segment = segment.points().map(|point| { - if point > self.max { - self.max = point; - } - - point - }); - - let edge = Segment::from(segment); - self.edges.push(edge); + self.edges.push(segment); } self } + #[cfg(test)] + pub fn invert_winding(mut self) -> Self { + self.edges.reverse(); + self + } + pub fn contains_triangle( &self, &[a, b, c]: &[Point<2>; 3], @@ -92,52 +86,158 @@ impl Polygon { point: impl Into>, debug_info: &mut DebugInfo, ) -> bool { - let point = point.into(); - let outside = self.max * 2.; - - let dir = outside - point; - let ray = Ray2 { - origin: point.to_na(), - dir: dir.to_na(), + let ray = HorizontalRayToTheRight { + origin: point.into(), }; - let mut check = - TriangleEdgeCheck::new(self.surface.point_surface_to_model(&point)); + let mut check = TriangleEdgeCheck::new( + self.surface.point_surface_to_model(&ray.origin), + ); - // We need to keep track of where our ray hits the edges. Otherwise, if - // the ray hits a vertex, we might count that hit twice, as every vertex - // is attached to two edges. - let mut hits = BTreeSet::new(); + let mut num_hits = 0; + + // We need to properly detect the ray passing the boundary at the "seam" + // of the polygon, i.e. the vertex between the last and the first + // segment. The logic in the loop properly takes care of that, as long + // as we initialize the `previous_hit` variable with the result of the + // last segment. + let mut previous_hit = self + .edges + .last() + .copied() + .and_then(|edge| ray.hits_segment(edge)); - // Use ray-casting to determine if `center` is within the face-polygon. for &edge in &self.edges { - // Please note that we if we get to this point, then the point is - // not on a polygon edge, due to the check above. We don't need to - // handle any edge cases that would arise from that case. - - let intersection = edge - .to_parry() - .cast_local_ray(&ray, f64::INFINITY, true) - .map(Scalar::from_f64); - - if let Some(t) = intersection { - // Due to slight inaccuracies, we might get different values for - // the same intersections. Let's round `t` before using it. - let eps = 1_000_000.0; - let t = (t * eps).round() / eps; - - if hits.insert(t) { - let edge = - Segment::from_points(edge.points().map(|point| { - self.surface.point_surface_to_model(&point) - })); - check.hits.push(edge); + let hit = ray.hits_segment(edge); + + let count_hit = match (hit, previous_hit) { + (Some(Hit::Segment), _) => { + // We're hitting a segment right-on. Clear case. + true + } + (Some(Hit::UpperVertex), Some(Hit::LowerVertex)) + | (Some(Hit::LowerVertex), Some(Hit::UpperVertex)) => { + // If we're hitting a vertex, only count it if we've hit the + // other kind of vertex right before. + // + // That means, we're passing through the polygon boundary + // at where two edges touch. Depending on the order in which + // edges are checked, we're seeing this as a hit to one + // edge's lower/upper vertex, then the other edge's opposite + // vertex. + // + // If we're seeing two of the same vertices in a row, we're + // not actually passing through the polygon boundary. Then + // we're just touching a vertex without passing through + // anything. + true + } + (Some(Hit::Parallel), _) => { + // A parallel edge must be completely ignored. Its presence + // won't change anything, so we can treat it as if it + // wasn't there, and its neighbors were connected to each + // other. + continue; + } + _ => { + // Any other case is not a valid hit. + false } + }; + + if count_hit { + num_hits += 1; + + let edge = + Segment::from_points(edge.points().map(|point| { + self.surface.point_surface_to_model(&point) + })); + check.hits.push(edge); } + + previous_hit = hit; } debug_info.triangle_edge_checks.push(check); - hits.len() % 2 == 1 + num_hits % 2 == 1 + } +} + +#[cfg(test)] +mod tests { + use fj_debug::DebugInfo; + use fj_math::{Point, PolyChain}; + + use crate::geometry::Surface; + + use super::Polygon; + + #[test] + fn contains_point_ray_hits_vertex_while_passing_outside() { + let a = [0., 0.]; + let b = [2., 1.]; + let c = [0., 2.]; + + let polygon = Polygon::new(Surface::x_y_plane()) + .with_exterior(PolyChain::from([a, b, c]).close()); + + assert_contains_point(polygon, [1., 1.]); + } + + #[test] + fn contains_point_ray_hits_vertex_at_polygon_seam() { + let a = [2., 1.]; + let b = [0., 2.]; + let c = [0., 0.]; + + let polygon = Polygon::new(Surface::x_y_plane()) + .with_exterior(PolyChain::from([a, b, c]).close()); + + assert_contains_point(polygon, [1., 1.]); + } + + #[test] + fn contains_point_ray_hits_vertex_while_staying_inside() { + let a = [0., 0.]; + let b = [2., 1.]; + let c = [3., 0.]; + let d = [3., 4.]; + + let polygon = Polygon::new(Surface::x_y_plane()) + .with_exterior(PolyChain::from([a, b, c, d]).close()); + + assert_contains_point(polygon, [1., 1.]); + } + + #[test] + fn contains_ray_hits_parallel_edge() { + // Ray passes polygon boundary at a vertex. + let a = [0., 0.]; + let b = [2., 1.]; + let c = [3., 1.]; + let d = [0., 2.]; + let polygon = Polygon::new(Surface::x_y_plane()) + .with_exterior(PolyChain::from([a, b, c, d]).close()); + assert_contains_point(polygon, [1., 1.]); + + // Ray hits a vertex, but doesn't pass polygon boundary there. + let a = [0., 0.]; + let b = [2., 1.]; + let c = [3., 1.]; + let d = [4., 0.]; + let e = [4., 5.]; + let polygon = Polygon::new(Surface::x_y_plane()) + .with_exterior(PolyChain::from([a, b, c, d, e]).close()); + assert_contains_point(polygon, [1., 1.]); + } + + fn assert_contains_point(polygon: Polygon, point: impl Into>) { + let point = point.into(); + + assert!(polygon.contains_point(point, &mut DebugInfo::new())); + assert!(polygon + .invert_winding() + .contains_point(point, &mut DebugInfo::new())); } } diff --git a/fj-kernel/src/algorithms/triangulation/ray.rs b/fj-kernel/src/algorithms/triangulation/ray.rs index 23adc6598..e46cbb05c 100644 --- a/fj-kernel/src/algorithms/triangulation/ray.rs +++ b/fj-kernel/src/algorithms/triangulation/ray.rs @@ -1,5 +1,3 @@ -#![cfg(test)] - use fj_math::{Point, Segment}; pub struct HorizontalRayToTheRight { From b2e03dc921024b04c05f2d108642f736ef06bbe0 Mon Sep 17 00:00:00 2001 From: Hanno Braun Date: Sat, 9 Apr 2022 15:18:07 +0200 Subject: [PATCH 4/5] Inline method calls --- .../src/algorithms/triangulation/polygon.rs | 20 +++++++++---------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/fj-kernel/src/algorithms/triangulation/polygon.rs b/fj-kernel/src/algorithms/triangulation/polygon.rs index 6ab8c5e59..0505a646d 100644 --- a/fj-kernel/src/algorithms/triangulation/polygon.rs +++ b/fj-kernel/src/algorithms/triangulation/polygon.rs @@ -18,8 +18,12 @@ impl Polygon { } } - pub fn with_exterior(self, exterior: impl Into>) -> Self { - self.with_approx(exterior) + pub fn with_exterior(mut self, exterior: impl Into>) -> Self { + for segment in exterior.into().segments() { + self.edges.push(segment); + } + + self } pub fn with_interiors( @@ -27,15 +31,9 @@ impl Polygon { interiors: impl IntoIterator>>, ) -> Self { for interior in interiors { - self = self.with_approx(interior); - } - - self - } - - fn with_approx(mut self, approx: impl Into>) -> Self { - for segment in approx.into().segments() { - self.edges.push(segment); + for segment in interior.into().segments() { + self.edges.push(segment); + } } self From 0ad96f6685d61842ca336246ca5d8798d281e428 Mon Sep 17 00:00:00 2001 From: Hanno Braun Date: Sat, 9 Apr 2022 15:35:05 +0200 Subject: [PATCH 5/5] Fix last known point-in-polygon edge case The ray leaving the polygon through the seam, was not detected in the presence of multiple polygonal chains. Now it is. --- .../src/algorithms/triangulation/polygon.rs | 164 ++++++++++-------- 1 file changed, 89 insertions(+), 75 deletions(-) diff --git a/fj-kernel/src/algorithms/triangulation/polygon.rs b/fj-kernel/src/algorithms/triangulation/polygon.rs index 0505a646d..057e8b481 100644 --- a/fj-kernel/src/algorithms/triangulation/polygon.rs +++ b/fj-kernel/src/algorithms/triangulation/polygon.rs @@ -7,22 +7,21 @@ use super::ray::{Hit, HorizontalRayToTheRight}; pub struct Polygon { surface: Surface, - edges: Vec>, + exterior: PolyChain<2>, + interiors: Vec>, } impl Polygon { pub fn new(surface: Surface) -> Self { Self { surface, - edges: Vec::new(), + exterior: PolyChain::new(), + interiors: Vec::new(), } } pub fn with_exterior(mut self, exterior: impl Into>) -> Self { - for segment in exterior.into().segments() { - self.edges.push(segment); - } - + self.exterior = exterior.into(); self } @@ -30,18 +29,18 @@ impl Polygon { mut self, interiors: impl IntoIterator>>, ) -> Self { - for interior in interiors { - for segment in interior.into().segments() { - self.edges.push(segment); - } - } - + self.interiors.extend(interiors.into_iter().map(Into::into)); self } #[cfg(test)] pub fn invert_winding(mut self) -> Self { - self.edges.reverse(); + self.exterior = self.exterior.reverse(); + + for interior in &mut self.interiors { + *interior = interior.clone().reverse(); + } + self } @@ -76,7 +75,14 @@ impl Polygon { } pub fn contains_edge(&self, edge: Segment<2>) -> bool { - self.edges.contains(&edge) || self.edges.contains(&edge.reverse()) + let mut contains = false; + + for chain in Some(&self.exterior).into_iter().chain(&self.interiors) { + contains |= chain.segments().contains(&edge); + contains |= chain.segments().contains(&edge.reverse()); + } + + contains } pub fn contains_point( @@ -94,66 +100,69 @@ impl Polygon { let mut num_hits = 0; - // We need to properly detect the ray passing the boundary at the "seam" - // of the polygon, i.e. the vertex between the last and the first - // segment. The logic in the loop properly takes care of that, as long - // as we initialize the `previous_hit` variable with the result of the - // last segment. - let mut previous_hit = self - .edges - .last() - .copied() - .and_then(|edge| ray.hits_segment(edge)); - - for &edge in &self.edges { - let hit = ray.hits_segment(edge); - - let count_hit = match (hit, previous_hit) { - (Some(Hit::Segment), _) => { - // We're hitting a segment right-on. Clear case. - true - } - (Some(Hit::UpperVertex), Some(Hit::LowerVertex)) - | (Some(Hit::LowerVertex), Some(Hit::UpperVertex)) => { - // If we're hitting a vertex, only count it if we've hit the - // other kind of vertex right before. - // - // That means, we're passing through the polygon boundary - // at where two edges touch. Depending on the order in which - // edges are checked, we're seeing this as a hit to one - // edge's lower/upper vertex, then the other edge's opposite - // vertex. - // - // If we're seeing two of the same vertices in a row, we're - // not actually passing through the polygon boundary. Then - // we're just touching a vertex without passing through - // anything. - true - } - (Some(Hit::Parallel), _) => { - // A parallel edge must be completely ignored. Its presence - // won't change anything, so we can treat it as if it - // wasn't there, and its neighbors were connected to each - // other. - continue; + for chain in Some(&self.exterior).into_iter().chain(&self.interiors) { + let edges = chain.segments(); + + // We need to properly detect the ray passing the boundary at the + // "seam" of the polygon, i.e. the vertex between the last and the + // first segment. The logic in the loop properly takes care of that, + // as long as we initialize the `previous_hit` variable with the + // result of the last segment. + let mut previous_hit = edges + .last() + .copied() + .and_then(|edge| ray.hits_segment(edge)); + + for edge in edges { + let hit = ray.hits_segment(edge); + + let count_hit = match (hit, previous_hit) { + (Some(Hit::Segment), _) => { + // We're hitting a segment right-on. Clear case. + true + } + (Some(Hit::UpperVertex), Some(Hit::LowerVertex)) + | (Some(Hit::LowerVertex), Some(Hit::UpperVertex)) => { + // If we're hitting a vertex, only count it if we've hit the + // other kind of vertex right before. + // + // That means, we're passing through the polygon boundary + // at where two edges touch. Depending on the order in which + // edges are checked, we're seeing this as a hit to one + // edge's lower/upper vertex, then the other edge's opposite + // vertex. + // + // If we're seeing two of the same vertices in a row, we're + // not actually passing through the polygon boundary. Then + // we're just touching a vertex without passing through + // anything. + true + } + (Some(Hit::Parallel), _) => { + // A parallel edge must be completely ignored. Its presence + // won't change anything, so we can treat it as if it + // wasn't there, and its neighbors were connected to each + // other. + continue; + } + _ => { + // Any other case is not a valid hit. + false + } + }; + + if count_hit { + num_hits += 1; + + let edge = + Segment::from_points(edge.points().map(|point| { + self.surface.point_surface_to_model(&point) + })); + check.hits.push(edge); } - _ => { - // Any other case is not a valid hit. - false - } - }; - - if count_hit { - num_hits += 1; - let edge = - Segment::from_points(edge.points().map(|point| { - self.surface.point_surface_to_model(&point) - })); - check.hits.push(edge); + previous_hit = hit; } - - previous_hit = hit; } debug_info.triangle_edge_checks.push(check); @@ -185,14 +194,19 @@ mod tests { #[test] fn contains_point_ray_hits_vertex_at_polygon_seam() { - let a = [2., 1.]; - let b = [0., 2.]; + let a = [4., 2.]; + let b = [0., 4.]; let c = [0., 0.]; + let d = [1., 1.]; + let e = [2., 1.]; + let f = [1., 3.]; + let polygon = Polygon::new(Surface::x_y_plane()) - .with_exterior(PolyChain::from([a, b, c]).close()); + .with_exterior(PolyChain::from([a, b, c]).close()) + .with_interiors([PolyChain::from([d, e, f]).close()]); - assert_contains_point(polygon, [1., 1.]); + assert_contains_point(polygon, [1., 2.]); } #[test]