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" 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/polygon.rs b/fj-kernel/src/algorithms/triangulation/polygon.rs index 9d77df150..057e8b481 100644 --- a/fj-kernel/src/algorithms/triangulation/polygon.rs +++ b/fj-kernel/src/algorithms/triangulation/polygon.rs @@ -1,53 +1,44 @@ -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>, + exterior: PolyChain<2>, + interiors: Vec>, } impl Polygon { pub fn new(surface: Surface) -> Self { Self { surface, - edges: Vec::new(), - max: Point::origin(), + exterior: PolyChain::new(), + interiors: Vec::new(), } } - pub fn with_exterior(self, exterior: impl Into>) -> Self { - self.with_approx(exterior) + pub fn with_exterior(mut self, exterior: impl Into>) -> Self { + self.exterior = exterior.into(); + self } pub fn with_interiors( mut self, interiors: impl IntoIterator>>, ) -> Self { - for interior in interiors { - self = self.with_approx(interior); - } - + self.interiors.extend(interiors.into_iter().map(Into::into)); self } - 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 - }); + #[cfg(test)] + pub fn invert_winding(mut self) -> Self { + self.exterior = self.exterior.reverse(); - let edge = Segment::from(segment); - self.edges.push(edge); + for interior in &mut self.interiors { + *interior = interior.clone().reverse(); } self @@ -84,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( @@ -92,52 +90,166 @@ 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)); - - // 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(); - - // 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); + let mut check = TriangleEdgeCheck::new( + self.surface.point_surface_to_model(&ray.origin), + ); + + let mut num_hits = 0; + + 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; - 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); } + + 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 = [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_interiors([PolyChain::from([d, e, f]).close()]); + + assert_contains_point(polygon, [1., 2.]); + } + + #[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 new file mode 100644 index 000000000..e46cbb05c --- /dev/null +++ b/fj-kernel/src/algorithms/triangulation/ray.rs @@ -0,0 +1,161 @@ +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))); + } +}