From e8047c0a4c823e2d5005ff430ab6d50e7997f408 Mon Sep 17 00:00:00 2001
From: Mark Harris <mharris@nvidia.com>
Date: Thu, 8 Dec 2022 12:40:30 +1100
Subject: [PATCH] New strict crossings multiply PIP algorithm

---
 .../detail/is_point_in_polygon.cuh            | 63 +++++++++++++++++++
 .../experimental/detail/point_in_polygon.cuh  | 14 ++---
 2 files changed, 70 insertions(+), 7 deletions(-)

diff --git a/cpp/include/cuspatial/experimental/detail/is_point_in_polygon.cuh b/cpp/include/cuspatial/experimental/detail/is_point_in_polygon.cuh
index b9c7d5ac6..2b3ee2a79 100644
--- a/cpp/include/cuspatial/experimental/detail/is_point_in_polygon.cuh
+++ b/cpp/include/cuspatial/experimental/detail/is_point_in_polygon.cuh
@@ -100,5 +100,68 @@ __device__ inline bool is_point_in_polygon(Cart2d const& test_point,
 
   return point_is_within;
 }
+
+// In progress -- failing TwoPolygonsOneRingEach
+template <class Cart2d,
+          class OffsetType,
+          class OffsetIterator,
+          class Cart2dIt,
+          class OffsetItDiffType = typename std::iterator_traits<OffsetIterator>::difference_type,
+          class Cart2dItDiffType = typename std::iterator_traits<Cart2dIt>::difference_type>
+__device__ inline bool is_point_in_polygon2(Cart2d const& test_point,
+                                            OffsetType poly_begin,
+                                            OffsetType poly_end,
+                                            OffsetIterator ring_offsets_first,
+                                            OffsetItDiffType const& num_rings,
+                                            Cart2dIt poly_points_first,
+                                            Cart2dItDiffType const& num_poly_points)
+{
+  using T = iterator_vec_base_type<Cart2dIt>;
+
+  bool inside_right = false;
+  bool inside_left  = false;
+
+  // for each ring
+  for (auto ring_idx = poly_begin; ring_idx < poly_end; ring_idx++) {
+    int32_t ring_idx_next = ring_idx + 1;
+    int32_t ring_begin    = ring_offsets_first[ring_idx];
+    int32_t ring_end =
+      (ring_idx_next < num_rings) ? ring_offsets_first[ring_idx_next] : num_poly_points;
+
+    Cart2d a = poly_points_first[ring_end - 1];
+
+    // for each line segment, starting with the segment between the last and first vertex
+    for (auto point_idx = ring_begin; point_idx < ring_end; point_idx++) {
+      Cart2d const b = poly_points_first[point_idx];
+
+      T run  = b.x - a.x;
+      T rise = b.y - a.y;
+
+      // Points on the line segment are the same, so intersection is impossible.
+      // This is possible because we allow closed or unclosed polygons.
+      T constexpr zero = 0.0;
+      if (float_equal(run, zero) && float_equal(rise, zero)) continue;
+
+      bool in_y = (a.y <= test_point.y && test_point.y <= b.y) ||
+                  (b.y <= test_point.y && test_point.y <= a.y);
+      if (in_y) {
+        // Transform the following inequality to avoid division
+        //  test_point.x < (run / rise) * (test_point.y - a.y) + a.x
+        auto lhs = run * (test_point.y - a.y);
+        auto rhs = rise * (test_point.x - a.x);
+        if (float_equal(lhs, rhs)) {  // point is on an edge
+          return false;
+        } else if ((lhs < rhs) == (rise > 0)) {
+          inside_left = !inside_left;
+        } else if ((lhs > rhs) == (rise > 0)) {
+          inside_right = !inside_right;
+        }
+      }
+      a = b;
+    }
+  }
+
+  return inside_right && inside_left;
+}
 }  // namespace detail
 }  // namespace cuspatial
diff --git a/cpp/include/cuspatial/experimental/detail/point_in_polygon.cuh b/cpp/include/cuspatial/experimental/detail/point_in_polygon.cuh
index 5232cae1f..034245b4a 100644
--- a/cpp/include/cuspatial/experimental/detail/point_in_polygon.cuh
+++ b/cpp/include/cuspatial/experimental/detail/point_in_polygon.cuh
@@ -68,13 +68,13 @@ __global__ void point_in_polygon_kernel(Cart2dItA test_points_first,
     OffsetType poly_end =
       (poly_idx_next < num_polys) ? poly_offsets_first[poly_idx_next] : num_rings;
 
-    bool const point_is_within = is_point_in_polygon(test_point,
-                                                     poly_begin,
-                                                     poly_end,
-                                                     ring_offsets_first,
-                                                     num_rings,
-                                                     poly_points_first,
-                                                     num_poly_points);
+    bool const point_is_within = is_point_in_polygon2(test_point,
+                                                      poly_begin,
+                                                      poly_end,
+                                                      ring_offsets_first,
+                                                      num_rings,
+                                                      poly_points_first,
+                                                      num_poly_points);
 
     hit_mask |= point_is_within << poly_idx;
   }