GeometryOps.within(poly1, poly2) inconsistent with LibGEOS.within(poly1, poly2)

[BonelessPi]

For certain polygons, within() seems to incorrectly return true. I honestly don't know why this is happening. Maybe it relates to the polygons sharing vertices? I used Point2i to try to avoid floating point issues but maybe that's responsible? Both polygons are simply connected. In my example, I added the within() checks for the point (4,2) to test if the winding order of my points was causing an issue. I'm rather new to the julia geometry packages, so apologies if this is user error. I suspect this is a bug as I can't understand how a certain point could be within rect, not within poly, yet rect is within poly.

Example:

using GeometryBasics, GeometryOps, CairoMakie

rect = Polygon{2, Int64}(Point{2, Int64}[[0, 0], [8, 0], [8, 4], [0, 4]], Vector{Point{2, Int64}}[])
poly = Polygon{2, Int64}(Point{2, Int64}[[8, 4], [0, 4], [0, 3], [7, 3], [7, 1], [0, 2], [0, 0], [8, 0]], Vector{Point{2, Int64}}[])

within(rect, poly) # returns true, should be false
within(poly, rect) # returns true
within(Point2i(4, 2), rect) # returns true
within(Point2i(4, 2), poly) # returns false

f,ax,_ = plot(rect);
plot!(ax,poly);
display(f)

I was able to reproduce this issue in a temporary environment with the following package versions:
[13f3f980] CairoMakie v0.15.9
[5c1252a2] GeometryBasics v0.5.10
[3251bfac] GeometryOps v0.1.38

[BonelessPi]

I read a bit into the tests to see how this package expects to be called. I have tweaked my example to use LibGEOS.Polygon instead of GeometryBasics, changed to a CW winding order, and made sure the list of points given to the polygon is closed.

Based on this example, LibGEOS.within() gave the correct result (as did GO.within(GO.GEOS(),...)). It seems that there is an issue with the native julia implementation in GeometryOps. Maybe this should be a new test case once the issue is fixed?

import GeometryOps as GO
import LibGEOS as LG
using CairoMakie

rect = LG.Polygon([Vector{Float64}[[0, 0], [0, 4], [8, 4], [8, 0], [0, 0]]])
poly = LG.Polygon([Vector{Float64}[[0, 0], [0, 2], [7, 1], [7, 3], [0, 3], [0, 4], [8, 4], [8, 0], [0, 0]]])
@show rect
@show poly
println()
@show GO.within(rect, poly) # returns true, should be false
@show GO.within(poly, rect) # returns true
@show GO.within(GO.GEOS(), rect, poly) # returns false
@show GO.within(GO.GEOS(), poly, rect) # returns true
@show LG.within(rect, poly) # returns false
@show LG.within(poly, rect) # returns true
println()
@show GO.within(LG.Point([4.0, 2.0]), rect) # returns true
@show GO.within(LG.Point([4.0, 2.0]), poly) # returns false

f,ax,_ = plot(rect);
plot!(ax,poly);
display(f)

[asinghvi17]

Thanks for the report and the perfect MWE! Looks like this was a bug in the intersection finding code, should be fixed by the linked PR. I want to put that through some more testing and see if it fixes some other issues that we have had for a while as well.