Merge pull request #17 from bast/radovan/towards-v0.3.4

Working towards v0.3.4

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "polygons"
-version = "0.3.3"
+version = "0.3.4"
 authors = ["Radovan Bast <[email protected]>"]
 description = "Fast points-in-polygon test and distances to polygons."
 license = "GPL-3.0-only"
@@ -16,8 +16,9 @@ crate-type = ["cdylib", "rlib"]
 default = ["rayon", "pyo3"]
 
 [dependencies]
-rayon = { version = "1.8", optional = true }
-pyo3 = { version = "0.20", features = ["extension-module"], optional = true }
+rayon = { version = "1.10", optional = true }
+pyo3 = { version = "0.22", features = ["extension-module"], optional = true }
+float-cmp = "0.10"
 
 [dev-dependencies]
 rand = "0.8"

README.md

@@ -47,10 +47,10 @@ author(s):
 @misc{polygons,
   author    = {Bast, Radovan},
   title     = {Polygons: Fast points-in-polygon test and distances to polygons},
-  month     = {1},
+  month     = {10},
   year      = {2024},
   publisher = {Zenodo},
-  version   = {v0.3.3},
+  version   = {v0.3.4},
   doi       = {10.5281/zenodo.3825616},
   url       = {https://doi.org/10.5281/zenodo.3825616}
 }

src/intersections.rs

@@ -1,5 +1,9 @@
 use crate::tree::{Edge, Node};
+use float_cmp::approx_eq;
 
+// we count each intersection twice since this makes it possible to deal with the reference point
+// has the same y-coordinate as an edge point and then we can avoid double-counting the
+// intersection
 pub fn num_intersections(node: &Node, n: i32, p: (f64, f64)) -> i32 {
     if skip_box_intersection(p, node) {
         return n;
@@ -17,7 +21,14 @@ pub fn num_intersections(node: &Node, n: i32, p: (f64, f64)) -> i32 {
     if !node.edges.is_empty() {
         for edge in &node.edges {
             if crosses(p, edge) {
-                n_ += 1;
+                // if y-coordinate of reference point is equal to y-coordinate of edge point
+                if (approx_eq!(f64, p.1, edge.p1.y, ulps = 2) && edge.p1.in_between)
+                    || (approx_eq!(f64, p.1, edge.p2.y, ulps = 2) && edge.p2.in_between)
+                {
+                    n_ += 1;
+                } else {
+                    n_ += 2;
+                }
             }
         }
         return n_;

src/lib.rs

@@ -13,4 +13,3 @@ mod intersections;
 #[cfg(feature = "pyo3")]
 mod python;
 mod tree;
-

src/python.rs

@@ -1,5 +1,4 @@
 use pyo3::prelude::*;
-use pyo3::wrap_pyfunction;
 
 use crate::tree;
 use crate::tree::build_search_tree;
@@ -22,7 +21,7 @@ fn distances_nearest_edges(tree: Tree, points: Vec<(f64, f64)>) -> Vec<f64> {
 }
 
 #[pymodule]
-fn polygons(_py: Python, m: &PyModule) -> PyResult<()> {
+fn polygons(m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add("__version__", env!("CARGO_PKG_VERSION"))?;
 
     m.add_function(wrap_pyfunction!(build_search_tree, m)?)?;

src/tree.rs

@@ -15,10 +15,11 @@ pub struct Point {
     pub y: f64,
     pub h: f64,
     pub index: usize,
+    pub in_between: bool,
 }
 
 // edge connects two points
-#[derive(Clone)]
+#[derive(Debug, Clone)]
 pub struct Edge {
     pub p1: Point,
     pub p2: Point,
@@ -113,7 +114,10 @@ pub fn points_are_inside(tree: &Tree, points: &[(f64, f64)]) -> Vec<bool> {
     #[cfg(not(feature = "rayon"))]
     let iter = points.iter();
 
-    iter.map(|p| (intersections::num_intersections(&tree[0], 0, *p) % 2) != 0)
+    // the division by 2 is because we count each edge intersection twice
+    // and the reason for that is that it makes it easier to deal with the case
+    // where the point has the same y coordinate as an edge point
+    iter.map(|p| ((intersections::num_intersections(&tree[0], 0, *p) / 2) % 2) != 0)
         .collect()
 }
 
@@ -245,39 +249,82 @@ fn group_edges(num_edges_children: usize, input: Vec<Edge>) -> Tree {
     parents
 }
 
+fn is_in_between(y1: &f64, y2: &f64, y3: &f64) -> bool {
+    (y1 < y2 && y2 < y3) || (y1 > y2 && y2 > y3)
+}
+
 fn points_to_edges(points: &[(f64, f64, f64)], offset: usize) -> Vec<Edge> {
     let mut edges: Vec<Edge> = points
-        .windows(2)
+        .windows(4)
         .enumerate()
         .map(|(i, t)| Edge {
             p1: Point {
-                x: t[0].0,
-                y: t[0].1,
-                h: t[0].2,
-                index: offset + i,
-            },
-            p2: Point {
                 x: t[1].0,
                 y: t[1].1,
                 h: t[1].2,
                 index: offset + i + 1,
+                in_between: is_in_between(&t[0].1, &t[1].1, &t[2].1),
+            },
+            p2: Point {
+                x: t[2].0,
+                y: t[2].1,
+                h: t[2].2,
+                index: offset + i + 2,
+                in_between: is_in_between(&t[1].1, &t[2].1, &t[3].1),
             },
         })
         .collect();
 
     let n = points.len() - 1;
+
+    edges.push(Edge {
+        p1: Point {
+            x: points[n - 1].0,
+            y: points[n - 1].1,
+            h: points[n - 1].2,
+            index: offset + n - 1,
+            in_between: is_in_between(&points[n - 2].1, &points[n - 1].1, &points[n].1),
+        },
+        p2: Point {
+            x: points[n].0,
+            y: points[n].1,
+            h: points[n].2,
+            index: offset + n,
+            in_between: is_in_between(&points[n - 1].1, &points[n].1, &points[0].1),
+        },
+    });
+
     edges.push(Edge {
         p1: Point {
             x: points[n].0,
             y: points[n].1,
             h: points[n].2,
             index: offset + n,
+            in_between: is_in_between(&points[n - 1].1, &points[n].1, &points[0].1),
         },
         p2: Point {
             x: points[0].0,
             y: points[0].1,
             h: points[0].2,
             index: offset,
+            in_between: is_in_between(&points[n].1, &points[0].1, &points[1].1),
+        },
+    });
+
+    edges.push(Edge {
+        p1: Point {
+            x: points[0].0,
+            y: points[0].1,
+            h: points[0].2,
+            index: offset,
+            in_between: is_in_between(&points[n].1, &points[0].1, &points[1].1),
+        },
+        p2: Point {
+            x: points[1].0,
+            y: points[1].1,
+            h: points[1].2,
+            index: offset + 1,
+            in_between: is_in_between(&points[0].1, &points[1].1, &points[2].1),
         },
     });
 

tests/case-2/boundary.txt

@@ -0,0 +1,21 @@
+20
+0.1673524526250124 0.9178821207993946
+0.011481376675199044 0.5040866625032467
+0.1861643712112249 0.7022171017595733
+0.0 0.13416003308604163
+0.08458961373099898 0.0
+0.1395451240976833 0.282147861601127
+0.34764509629213264 0.3991547290646283
+0.16261961056238453 0.1930772723649014
+0.48383690178366934 0.35669979866141643
+0.3100479670261293 0.09031279796896015
+0.38512891495090623 0.030257905954670757
+0.860185843678497 0.19414501989373728
+1.0 0.10829172533525562
+0.4665000797700417 0.8011822602635931
+0.7467270273354121 0.979553598458782
+0.32280537665797454 1.0
+0.2969309912531665 0.8329022529024206
+0.2868483350347063 0.5808802574078954
+0.17518032707146317 0.5689305796970081
+0.23128431456254747 0.8669948775639155