Improve Polygon.contains(LineString) predicate logic. (#1186)

The previous logic for handling the case of a LineString that shares boundary points with a Polygon and has no interior points was wrong, it was just finding the center point of such a LineString and testing that for containment. The correct logic is as follows:

A LineString is contained in a polygon if:
The sum of the points of the LineString that are contained in the Polygon, plus the number of vertices of the linestring that are on the polygon boundary, is equal to the length of the LineString less the number of points in the LineString that are contained in the boundary of the Polygon.

```
polygon.contains(linestring) = contains + point_intersection_count == rhs.sizes - linestring_intersection_count
```

The new logic counts the number of point intersections and linestring intersection that the LineString shares with the polygon. Point intersections occur when the LineString shares an edge point with the Polygon exterior, implying either a LineString segment that approaches from the exterior or the interior. An interior point will be counted, an exterior point will be left in the size of the LineString during comparison, with the final count being less than the size of the LineString.

The hardest part of this implementation was in writing `_pli_features_rebuild_offsets` that takes the results of a `pairwise_linestring_intersection` result and splits them into row-appropriate points in one set and linestrings in another set. This is potentially a good place to move the logic into C++, though I don't think it will be a major profiling issue initially.

Also improves `touches` and `crosses` where the predicate had been written too tightly to the test cases.

---

In addition, a bug with `pairwise_multipoint_equals_count` is fixed when `lhs` contains more than 1 multipoints, but all multipoints are empty. The bug causes the API to raise a cuda invalid configuration error.

Authors:
  - H. Thomson Comer (https://github.com/thomcom)
  - Michael Wang (https://github.com/isVoid)

Approvers:
  - Mark Harris (https://github.com/harrism)
  - Michael Wang (https://github.com/isVoid)

URL: #1186

cpp/include/cuspatial/detail/pairwise_multipoint_equals_count.cuh

@@ -97,12 +97,14 @@ OutputIt pairwise_multipoint_equals_count(MultiPointRangeA lhs,
     rhs.offsets_begin(), rhs.offsets_end(), rhs_point_sorted.begin(), rhs_point_sorted.end()};
 
   detail::zero_data_async(output, output + lhs.size(), stream);
-  auto [tpb, n_blocks] = grid_1d(lhs.num_points());
-  detail::pairwise_multipoint_equals_count_kernel<<<n_blocks, tpb, 0, stream.value()>>>(
-    lhs, rhs_sorted, output);
 
-  CUSPATIAL_CHECK_CUDA(stream.value());
+  if (lhs.num_points() > 0) {
+    auto [tpb, n_blocks] = grid_1d(lhs.num_points());
+    detail::pairwise_multipoint_equals_count_kernel<<<n_blocks, tpb, 0, stream.value()>>>(
+      lhs, rhs_sorted, output);
 
+    CUSPATIAL_CHECK_CUDA(stream.value());
+  }
   return output + lhs.size();
 }
 

python/cuspatial/cuspatial/core/binpreds/feature_contains.py

@@ -24,9 +24,9 @@
     MultiPoint,
     Point,
     Polygon,
-    _false_series,
-    _linestrings_to_center_point,
     _open_polygon_rings,
+    _pli_lines_to_multipoints,
+    _pli_points_to_multipoints,
     _points_and_lines_to_multipoints,
     _zero_series,
 )
@@ -49,7 +49,23 @@ def _preprocess(self, lhs, rhs):
         preprocessor_result = super()._preprocess_multipoint_rhs(lhs, rhs)
         return self._compute_predicate(lhs, rhs, preprocessor_result)
 
-    def _intersection_results_for_contains(self, lhs, rhs):
+    def _intersection_results_for_contains_linestring(self, lhs, rhs):
+        pli = _basic_intersects_pli(lhs, rhs)
+
+        # Convert the pli into points and multipoint intersections.
+        multipoint_points = _pli_points_to_multipoints(pli)
+        multipoint_lines = _pli_lines_to_multipoints(pli)
+
+        # Count the point intersections that are equal to points in the
+        # LineString
+        # Count the linestring intersections that are equal to points in
+        # the LineString
+        return (
+            _basic_equals_count(rhs, multipoint_points),
+            _basic_equals_count(rhs, multipoint_lines),
+        )
+
+    def _intersection_results_for_contains_polygon(self, lhs, rhs):
         pli = _basic_intersects_pli(lhs, rhs)
         pli_features = pli[1]
         if len(pli_features) == 0:
@@ -62,19 +78,15 @@ def _intersection_results_for_contains(self, lhs, rhs):
             pli_features, pli_offsets
         )
 
-        # A point in the rhs can be one of three possible states:
-        # 1. It is in the interior of the lhs
-        # 2. It is in the exterior of the lhs
-        # 3. It is on the boundary of the lhs
-        # This function tests if the point in the rhs is in the boundary
-        # of the lhs
         intersect_equals_count = _basic_equals_count(rhs, multipoints)
         return intersect_equals_count
 
     def _compute_polygon_polygon_contains(self, lhs, rhs, preprocessor_result):
         lines_rhs = _open_polygon_rings(rhs)
         contains = _basic_contains_count(lhs, lines_rhs).reset_index(drop=True)
-        intersects = self._intersection_results_for_contains(lhs, lines_rhs)
+        intersects = self._intersection_results_for_contains_polygon(
+            lhs, lines_rhs
+        )
         # A closed polygon has an extra line segment that is not used in
         # counting the number of points. We need to subtract this from the
         # number of points in the polygon.
@@ -87,47 +99,35 @@ def _compute_polygon_polygon_contains(self, lhs, rhs, preprocessor_result):
         result = contains + intersects >= rhs.sizes - polygon_size_reduction
         return result
 
-    def _test_interior(self, lhs, rhs):
-        # We only need to test linestrings that are length 2.
-        # Divide the linestring in half and test the point for containment
-        # in the polygon.
-        size_two = rhs.sizes == 2
-        if (size_two).any():
-            center_points = _linestrings_to_center_point(rhs[size_two])
-            size_two_results = _false_series(len(lhs))
-            size_two_results.iloc[rhs.index[size_two]] = (
-                _basic_contains_count(lhs, center_points) > 0
-            )
-            return size_two_results
-        else:
-            return _false_series(len(lhs))
-
     def _compute_polygon_linestring_contains(
         self, lhs, rhs, preprocessor_result
     ):
+        # Count the number of points in lhs that are properly contained by
+        # rhs
         contains = _basic_contains_count(lhs, rhs).reset_index(drop=True)
-        intersects = self._intersection_results_for_contains(lhs, rhs)
-
-        # If a linestring has intersection but not containment, we need to
-        # test if the linestring is in the interior of the polygon.
-        final_result = _false_series(len(lhs))
-        intersection_with_no_containment = (contains == 0) & (intersects != 0)
-        interior_tests = self._test_interior(
-            lhs[intersection_with_no_containment].reset_index(drop=True),
-            rhs[intersection_with_no_containment].reset_index(drop=True),
+
+        # Count the number of point intersections (line crossings) between
+        # lhs and rhs.
+        # Also count the number of perfectly overlapping linestring sections.
+        # Each linestring overlap counts as two point overlaps.
+        (
+            point_intersects_count,
+            linestring_intersects_count,
+        ) = self._intersection_results_for_contains_linestring(lhs, rhs)
+
+        # Subtract the length of the linestring intersections from the length
+        # of the rhs linestring, then test that the sum of contained points
+        # is equal to that adjusted rhs length.
+        rhs_sizes_less_line_intersection_size = (
+            rhs.sizes - linestring_intersects_count
         )
-        interior_tests.index = intersection_with_no_containment[
-            intersection_with_no_containment
-        ].index
-        # LineStrings that have intersection but no containment are set
-        # according to the `intersection_with_no_containment` mask.
-        final_result[intersection_with_no_containment] = interior_tests
-        # LineStrings that do not are contained if the sum of intersecting
-        # and containing points is greater than or equal to the number of
-        # points that make up the linestring.
-        final_result[~intersection_with_no_containment] = (
-            contains + intersects >= rhs.sizes
+        rhs_sizes_less_line_intersection_size[
+            rhs_sizes_less_line_intersection_size <= 0
+        ] = 1
+        final_result = contains + point_intersects_count == (
+            rhs_sizes_less_line_intersection_size
         )
+
         return final_result
 
     def _compute_predicate(self, lhs, rhs, preprocessor_result):

python/cuspatial/cuspatial/core/binpreds/feature_crosses.py

@@ -54,7 +54,7 @@ def _compute_predicate(self, lhs, rhs, preprocessor_result):
 
 class LineStringPolygonCrosses(BinPred):
     def _preprocess(self, lhs, rhs):
-        intersects = _basic_intersects_count(rhs, lhs) > 1
+        intersects = _basic_intersects_count(rhs, lhs) > 0
         touches = rhs.touches(lhs)
         contains = rhs.contains(lhs)
         return ~touches & intersects & ~contains

python/cuspatial/cuspatial/core/binpreds/feature_touches.py

@@ -25,6 +25,7 @@
     Point,
     Polygon,
     _false_series,
+    _pli_points_to_multipoints,
     _points_and_lines_to_multipoints,
 )
 
@@ -97,18 +98,19 @@ def _preprocess(self, lhs, rhs):
 
 class LineStringPolygonTouches(BinPred):
     def _preprocess(self, lhs, rhs):
+        intersects = _basic_intersects_count(lhs, rhs) > 0
+        contains = rhs.contains(lhs)
+        contains_any = _basic_contains_properly_any(rhs, lhs)
+
         pli = _basic_intersects_pli(lhs, rhs)
         if len(pli[1]) == 0:
             return _false_series(len(lhs))
-        intersections = _points_and_lines_to_multipoints(pli[1], pli[0])
+        points = _pli_points_to_multipoints(pli)
         # A touch can only occur if the point in the intersection
         # is equal to a point in the linestring: it must
         # terminate in the boundary of the polygon.
-        equals = _basic_equals_count(intersections, lhs) > 0
-        intersects = _basic_intersects_count(lhs, rhs)
-        intersects = (intersects == 1) | (intersects == 2)
-        contains = rhs.contains(lhs)
-        contains_any = _basic_contains_properly_any(rhs, lhs)
+        equals = _basic_equals_count(points, lhs) == points.sizes
+
         return equals & intersects & ~contains & ~contains_any
 
 

python/cuspatial/cuspatial/core/geoseries.py

@@ -191,23 +191,30 @@ def sizes(self):
             full_sizes = self.polygons.ring_offset.take(
                 self.polygons.part_offset.take(self.polygons.geometry_offset)
             )
-            return cudf.Series(full_sizes[1:] - full_sizes[:-1])
+            return cudf.Series(
+                full_sizes[1:] - full_sizes[:-1], index=self.index
+            )
         elif contains_only_linestrings(self):
             # Not supporting multilinestring yet
             full_sizes = self.lines.part_offset.take(
                 self.lines.geometry_offset
             )
-            return cudf.Series(full_sizes[1:] - full_sizes[:-1])
+            return cudf.Series(
+                full_sizes[1:] - full_sizes[:-1], index=self.index
+            )
         elif contains_only_multipoints(self):
-            return (
+            return cudf.Series(
                 self.multipoints.geometry_offset[1:]
-                - self.multipoints.geometry_offset[:-1]
+                - self.multipoints.geometry_offset[:-1],
+                index=self.index,
             )
         elif contains_only_points(self):
-            return cudf.Series(cp.repeat(cp.array(1), len(self)))
+            return cudf.Series(
+                cp.repeat(cp.array(1), len(self)), index=self.index
+            )
         else:
             if len(self) == 0:
-                return cudf.Series([0], dtype="int32")
+                return cudf.Series([0], dtype="int32", index=self.index)
             raise NotImplementedError(
                 "GeoSeries must contain only Points, MultiPoints, Lines, or "
                 "Polygons to return sizes."

python/cuspatial/cuspatial/tests/binpreds/binpred_test_dispatch.py

@@ -326,6 +326,78 @@ def predicate(request):
         LineString([(0.5, 1.25), (0.5, -0.25)]),
         point_polygon,
     ),
+    "linestring-polygon-cross-concave-edge": (
+        """
+    x  x  x
+    |\\ | /|
+    | xx- |
+    |  |  |
+    ---x---
+    """,
+        LineString([(0.5, 0.0), (0.5, 1.0)]),
+        Polygon([(0, 0), (0, 1), (0.3, 0.4), (1, 1), (1, 0)]),
+    ),
+    "linestring-polygon-half-in": (
+        """
+    -----
+    |   |
+    | x |
+    |/ \\|
+    xx-xx
+    """,
+        LineString(
+            [(0.0, 0.0), (0.25, 0.0), (0.5, 0.5), (0.75, 0.0), (1.0, 0.0)]
+        ),
+        point_polygon,
+    ),
+    "linestring-polygon-half-out": (
+        """
+    -----
+    |   |
+    |   |
+    |   |
+    xx-xx
+     \\/
+      x
+    """,
+        LineString(
+            [(0.0, 0.0), (0.25, 0.0), (0.5, -0.5), (0.75, 0.0), (1.0, 0.0)]
+        ),
+        point_polygon,
+    ),
+    "linestring-polygon-two-edges": (
+        """
+    x----
+    |   |
+    |   |
+    |   |
+    x---x
+    """,
+        LineString([(0.0, 1.0), (0.0, 0.0), (1.0, 0.0)]),
+        point_polygon,
+    ),
+    "linestring-polygon-edge-to-interior": (
+        """
+    x----
+    |   |
+    |  -x
+    |-/ |
+    x----
+    """,
+        LineString([(0.0, 1.0), (0.0, 0.0), (1.0, 0.5)]),
+        point_polygon,
+    ),
+    "linestring-polygon-edge-cross-to-exterior": (
+        """
+    x------
+    |     |
+    |    ---x
+    | --- |
+    x------
+    """,
+        LineString([(0.0, 1.0), (0.0, 0.0), (1.5, 0.5)]),
+        point_polygon,
+    ),
     "polygon-polygon-disjoint": (
         """
     Polygon polygon tests use a triangle for the lhs and a square for the rhs.
@@ -519,6 +591,12 @@ def predicate(request):
     "linestring-polygon-edge-interior",
     "linestring-polygon-in",
     "linestring-polygon-crosses",
+    "linestring-polygon-cross-concave-edge",
+    "linestring-polygon-half-in",
+    "linestring-polygon-half-out",
+    "linestring-polygon-two-edges",
+    "linestring-polygon-edge-to-interior",
+    "linestring-polygon-edge-cross-to-exterior",
 ]
 
 polygon_polygon_dispatch_list = [