[BUG]: Basic predicate operations for Polygon.contains(LineString) fail in small specific case.

Version

23.08

On which installation method(s) does this occur?

Rapids-Compose

Describe the issue

It turns out that our strategy of using three basic predicates with pip, intersection, and point equality does not work with a particular arrangement of Polygon.contains(Linestring)

That arrangement is visualized here:

In this example, a LineString that shares only border vertices with a containing polygon can't be distinguished as to whether or not the interior segment is inside or outside of the polygon. This effect occurs regardless of the length of the LineString, so long as the LineString segments are all on the border of the polygon.

The normal approach works well if any vertex is added to the LineString that is not on the border. We'll have to think of a deeper solution to this particular test case.

Minimum reproducible example

concave = cuspatial.GeoSeries([
    Polygon([
        (0, 0),
        (0, 1),
        (0.25, 1.0),
        (0.5, 0.5),
        (0.75, 1.0),
        (1, 1),
        (1, 0),
        (0, 0)
    ])
])
convex = cuspatial.GeoSeries([
    Polygon([
        (0, 0),
        (0, 1),
        (0.25, 1.0),
        (0.5, 1.5),
        (0.75, 1.0),
        (1, 1),
        (1, 0),
        (0, 0)
    ])
])
rhs = cuspatial.GeoSeries([LineString([(0, 1), (1, 1)])])
print(concave.contains(rhs))
print(convex.contains(rhs))

Relevant log output

0    False
dtype: bool
0    False
dtype: bool

Environment details

No response

Other/Misc.

No response

The following algorithm describes the approach we defined for correcting Polygon.contains. This whole routine should be written in libcuspatial.

lhs.contains(rhs) where lhs (left-hand side) is a Polygon and rhs (right-hand side) is a Polygon or LineString does not depend on the vertices that define the rhs geometry, but depends on the line segments.

The operation is the same for polygons or linestrings in the rhs. A polygon is simply a closed LineSegment and a LineString is, by definition, unclosed.

The algorithm for computing lhs.contains(rhs) is as follows:

intersection = linestring_intersection_for_contains(lhs, rhs)

linestring_intersection_for_contains should compute two results for the lhs and rhs, for many lhs and rhs in parallel:

boolean: whether or not a proper intersection in the lhs and rhs exists. A proper intersection between the lhs and rhs is a point of intersection that is not a vertex of the lhs or rhs.
proper intersections on the lhs: vertices that are in the lhs but not in the rhs should be used to divide the rhs segment before calling pip on the segments
intersection segments: a set of vertex pairs defining all of the overlapping segments in the intersection between lhs and rhs.

possibles = ~intersection[0]

possibles is the subset of lhs/rhs pairs for which the rhs may be contained in the lhs, because there is no proper intersection between them.

split = add_vertex_to_linestring(rhs[possibles], intersection[1])

add_vertex_to_linestring splits the segment containing each vertex into two segments at that vertex.

difference = pairwise_linestring_difference(split, intersection[2])

pairwise_linestring_difference subtracts each segment in the second argument from the linestring in the first argument. Split and difference could be combined.

centers = segments_to_centroids(difference)

pip_count = point_in_polygon_count(lhs, centers)

return pip_count == centers.sizes

rapidsai / cuspatial