GEOSPreparedDistanceWithin fails when comparing an interior subset of a line to the full line

[brendan-ward]

Given two lines where the subset is at least 3 vertices in from the endpoints of the larger line, GEOSPreparedDistanceWithin returns false when querying the subset line against the larger line, whereas GEOSDistanceWithin returns true.

Example using Shapely 2.0:

line1 = shapely.from_wkt('LINESTRING (0 0, 1 1, 2 2, 3 3, 4 4, 5 5, 6 6, 7 7, 8 8, 9 9)')
line2 = shapely.from_wkt('LINESTRING (2 2, 3 3, 4 4, 5 5, 6 6, 7 7)')

print(shapely.dwithin(line2, line1, 1))  # uses GEOSDistanceWithin_r
# True

shapely.prepare(line2)  # uses GEOSPrepare_r
print(shapely.dwithin(line2, line1, 1))  # uses GEOSPreparedDistanceWithin_r
# False

# longer subset works fine
line3 = shapely.from_wkt('LINESTRING (2 2, 3 3, 4 4, 5 5, 6 6, 7 7, 8 8)')

print(shapely.dwithin(line3, line1, 1))
# True

shapely.prepare(line3)
print(shapely.dwithin(line3, line1, 1))
# True

First observed in shapely #1887 using GEOS main (3.13.0dev), Shapely main, on MacOS.

[pramsey]

Test case for GEOSPreparedDistanceWithinTest.cpp

template<>
template<>
void object::test<12>
()
{
    checkDistanceWithin(
        "LINESTRING (2 2, 3 3, 4 4, 5 5, 6 6, 7 7)",
        "LINESTRING (0 0, 1 1, 2 2, 3 3, 4 4, 5 5, 6 6, 7 7, 8 8, 9 9)",
        1,
        1
    );
}

[dr-jts]

The problem is the incorrect heuristic check in IndexedFacetDistance.isWithinDistance here.

    auto env2 = g->getEnvelope();
    if (distance(env2.get()) > maxDistance) {
        return false;
    }

The check attempts to short-circuit the withinDistance test by checking the distance to the envelope polygon of the tested geometry. However, IndexedFacetDistance only computes distances to linework, not to polygonal geometries. So the test is incorrect in some cases (such as the failing case).

This heuristic is not in the JTS code.

@dbaston: looks like it should be removed here?

[brendan-ward]

Thanks for the fix @dr-jts !