Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

MultiLineString.contains produces incorrect result #933

Closed
thomcom opened this issue Jun 29, 2023 · 6 comments · Fixed by #937
Closed

MultiLineString.contains produces incorrect result #933

thomcom opened this issue Jun 29, 2023 · 6 comments · Fixed by #937
Labels

Comments

@thomcom
Copy link

thomcom commented Jun 29, 2023

I'm working on rapidsai/cuspatial#1063, which are GPU-accelerated multi-geometry binary predicate operations. I discovered a test case that fails in comparison with GeoPandas/Shapely/GEOS results that I thought might interest you:

import shapely
from shapely.geometry import LineString, MultiLineString
b1 = LineString([(0, 0), (1, 1)])
b2 = LineString([(0.5, 0.5), (1, 0.1), (-1, 0.1)])
m = MultiLineString([b1, b2])
print(m.contains(b1))
print(m.contains(b2))
False
True

Expected result is that m.contains(b1) == True. I'm very interested to learn if I'm doing something wrong or if there's a bug here that has since gone unnoticed. I'll be keeping an eye over here! Thanks!

@pramsey
Copy link
Member

pramsey commented Jun 29, 2023

Looks pretty suspect. The 0-dimensional interior interaction at (0.1 0.1) seems to be the culprit... because it's a multilinestring it can have different kinds of interactions (it has both 0-d and 1-d interior interactions), I wonder if the implication is that there is such a thing as an invalid multilinestring.

The relate matrix is

0F1
F00
1F2

It's definitely finding a 0-d interior/interior interaction, and contains would expect a 1-d interaction. @dr-jts ?

@dr-jts
Copy link
Contributor

dr-jts commented Jun 29, 2023

This is failing because the contains IM pattern of T*****FF* does not match the computed IM of 1F10001F2. So the culprit is the evaluation of Ext(m) x Int(b1) = 1 (in other words, the exterior of m is evaluated as containing a point from the interior of b1).

This is clearly wrong. I suspect a robustness issue due to the noding carried out by the relate algorithm. There's a similar issue here: locationtech/jts#396

I'l dig a bit more and confirm if this is what's happening.

In any case, I agree that this is a bug, and in fact m.contains(b1) = true.

@dr-jts dr-jts added the Bug label Jun 29, 2023
@dr-jts
Copy link
Contributor

dr-jts commented Jun 29, 2023

Here's a simpler reproducer (in JTS only - this appears to work in GEOS):

A: LINESTRING (0 0, 1 1)
B: LINESTRING (0.2 0.2, 0.5 0.5)
image

A.contains(B) currently evaluates to false, which is incorrect.

@dr-jts
Copy link
Contributor

dr-jts commented Jun 29, 2023

It's also interesting to note that this works:

A: MULTILINESTRING ((0 0, 1 1), (0.5 0.5, 1 0.25, -1 0.25))
B: LINESTRING (0 0, 1 1)

A.contains(B) = true

The reason is probably that 0.25 can be converted to an exact binary fraction, so does not cause robustness errors when the intersection is computed.

@dr-jts
Copy link
Contributor

dr-jts commented Jun 29, 2023

As I suspected, this error is due to a non-robust computation of the intersection point between LINESTRING (1 0.1, -1 0.1) and LINESTRING (0 0, 1 1). Using the Intersection::intersection implementation produces a value of POINT(0.09999999999999998, 0.1). Since the orientation predicate is evaluated exactly (or at least with a high degree of precision), this is enough "drift" to cause the effective line segments to not lie on the original ones, which produces the above-mentioned "exterior intersection".

One piece of good news is that the CGAlgorithmsDD.intersection implementation does not have this failure mode. So it can probably be used instead (with some slight performance impact). This also solves the failure case reported in locationtech/jts#396.

Note: this change will likely have some very small impacts on overlay test case results, based on JTS testing.

@dr-jts
Copy link
Contributor

dr-jts commented Jun 30, 2023

I confirmed that in GEOS using CGAlgorithmsDD.intersection fixes this issue.

But it does cause a few failures in other unit tests, probably due to the slight difference for some cases.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging a pull request may close this issue.

3 participants