Fix PreparedLineStringDistance for lines within envelope and polygons

[dr-jts]

This PR fixes PreparedLineStringDistance to handle the cases:

1. the prepared line is contained within the envelope of the test geometry
2. the prepared line is contained within a polygonal geometry

Fixing the first situation requires removing an invalid heuristic check from IndexedFacetDistance.isWithinDistance code:

    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 #958).

This heuristic is not in the JTS code.

Fixing the second situation requires an additional point-in-polygon check in
PreparedLineStringDistance.distance and PreparedLineStringDistance.isWithinDistance to test for the prepared line (or element of) lying within a test polygonal geometry.

Fixes #958.
Fixes #960.

[dr-jts]

The heuristic could be kept if it is only applied to test geometries whose envelope does not contain the indexed geometry envelope. Not sure if this effort is worth it though?

[dbaston]

Is it sufficient to guard the test with

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

[dr-jts]

Is it sufficient to guard the test with...

No, because the indexed geometry baseGeometry might have two elements. So its envelope might not be contained in the envelope of a polygon g, but one of the elements could still be contained in g and hence have distance = 0.

[dbaston]

Do you have a test case that would fail? I'm sorry, I've spent a good 20 minutes going over this and I can't understand the problem.

[dr-jts]

Do you have a test case that would fail?

My bad - this does work, with the addition of the PIP check in PreparedLineStringDistance.isWithinDistance. I will add the test case in GEOSPreparedDistanceWithin.

    checkDistanceWithin(
        "MULTILINESTRING ((1 6, 1 1), (14 16, 16 14))",
        "MULTIPOLYGON (((10 20, 20 20, 20 10, 10 10, 10 20)), ((30 20, 40 20, 40 10, 30 10, 30 20)))",
        1,
        1
    );

It's a bit confusing, because the IndexedFacetDistance.isWithinDistance actually returns false for this case. It would be good to have unit tests for IndexedFacetDistance.isWithinDistance.

[dr-jts]

Do you have a test case that would fail?

Actually while this heuristic works in the context of PreparedLineStringDistance (which has the extra PIP test), I'm not convinced it works in all cases for IndexedFacetDistance.isWithinDistance itself. I'll explore this further and make a IndexedFacetDistance test case.

[dbaston]

I think baseGeom is either:

• inside the envelope of g: skip the heuristic. This check was missing in the original implementation.
• partially inside the envelope of g: computed distance will be zero, which is not greater thanmaxDistance
• outside the envelope of g: if the envelope of g is > maxDistance from baseGeom, than any point on g is >= maxDistance from baseGeom so we can return false and avoid building the whole tree on g.

[dr-jts]

• partially inside the envelope of g: computed distance will be zero, which is not greater thanmaxDistance

If the indexed geometry has multiple elements, the computed distance to g.env may be > 0.

[dr-jts]

I'd like to commit this particular fix, and move the issue of the correctness of IndexedFacetDistance to a new issue. For now I will remove the heuristic (since IMO there is still a question of correctness). It can be reinstated if it can be demonstrated to work later on.

[dbaston]

If the indexed geometry has multiple elements, the computed distance to g.env may be > 0.

You're right, the heuristic can be wrong if baseGeom is a multipart. I would suggest changing it to

    if (baseGeometry.getNumGeometries() == 1 && !g->getEnvelopeInternal()->contains(baseGeometry.getEnvelopeInternal()))
    {
        auto env2 = g->getEnvelope();
        if (distance(env2.get()) > maxDistance) {
            return false;
        }
    }

and adding the following test to IndexedFacetDistance (demonstrating the need for the getNumGeometries() condition)

template<>
template<>
void object::test<13>()
{
    auto g1 = _wktreader.read("MULTILINESTRING ((6 6, 6.1 6.1), (100 100, 101 101))");
    auto g2 = _wktreader.read("LINESTRING (0 10, 10 0)");

    IndexedFacetDistance ifd(g1.get());

    ensure(ifd.isWithinDistance(g2.get(), 2));
}