Merge pull request #448 from hannobraun/ray

Fix triangulation issues regarding point-in-polygon test

Cargo.lock

@@ -717,6 +717,7 @@ dependencies = [
  "parking_lot 0.12.0",
  "parry2d-f64",
  "parry3d-f64",
+ "robust",
  "slotmap",
  "spade",
  "thiserror",

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"

fj-kernel/src/algorithms/triangulation/mod.rs

@@ -1,4 +1,5 @@
 mod polygon;
+mod ray;
 
 use fj_debug::DebugInfo;
 use fj_math::{Scalar, Triangle};

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<Segment<2>>,
-    max: Point<2>,
+    exterior: PolyChain<2>,
+    interiors: Vec<PolyChain<2>>,
 }
 
 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<PolyChain<2>>) -> Self {
-        self.with_approx(exterior)
+    pub fn with_exterior(mut self, exterior: impl Into<PolyChain<2>>) -> Self {
+        self.exterior = exterior.into();
+        self
     }
 
     pub fn with_interiors(
         mut self,
         interiors: impl IntoIterator<Item = impl Into<PolyChain<2>>>,
     ) -> 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<PolyChain<2>>) -> Self {
-        for segment in approx.into().segments() {
-            let segment = segment.points().map(|point| {
-                if point > self.max {
-                    self.max = point;
-                }
+    #[cfg(test)]
+    pub fn invert_winding(mut self) -> Self {
+        self.exterior = self.exterior.reverse();
 
-                point
-            });
-
-            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<Point<2>>,
         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;
 
-        // 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.
+        for chain in Some(&self.exterior).into_iter().chain(&self.interiors) {
+            let edges = chain.segments();
 
-            let intersection = edge
-                .to_parry()
-                .cast_local_ray(&ray, f64::INFINITY, true)
-                .map(Scalar::from_f64);
+            // 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));
 
-            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;
+            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 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<Point<2>>) {
+        let point = point.into();
+
+        assert!(polygon.contains_point(point, &mut DebugInfo::new()));
+        assert!(polygon
+            .invert_winding()
+            .contains_point(point, &mut DebugInfo::new()));
     }
 }

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<Segment<2>>) -> Option<Hit> {
+        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<P> From<P> for HorizontalRayToTheRight
+where
+    P: Into<Point<2>>,
+{
+    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)));
+    }
+}