Fix PreparedLineStringDistance for lines within envelope and polygons (…

…libgeos#959)

src/geom/prep/PreparedLineStringDistance.cpp

@@ -32,12 +32,19 @@ PreparedLineStringDistance::distance(const geom::Geometry* g) const
         return DoubleInfinity;
     }
 
-    // TODO: test if this shortcut be any useful
-    //if ( prepLine.intersects(g) ) return 0.0;
-
-    /* Not intersecting, compute distance from facets */
+    /* Compute potential distance from facets */
     operation::distance::IndexedFacetDistance *idf = prepLine.getIndexedFacetDistance();
-    return idf->distance(g);
+    double dist = idf->distance(g);
+    if (dist == 0.0)
+        return 0.0;
+
+    // If any point from prepLine is contained by g, the distance is zero
+    // Do this last because this PIP test is not indexed.
+    if ( g->getDimension() == 2 
+            && prepLine.isAnyTargetComponentInTest(g)) {
+        return 0.0;
+    }
+    return dist;
 }
 
 bool
@@ -49,7 +56,15 @@ PreparedLineStringDistance::isWithinDistance(const geom::Geometry* g, double d)
     }
 
     operation::distance::IndexedFacetDistance *idf = prepLine.getIndexedFacetDistance();
-    return idf->isWithinDistance(g, d);
+    if (idf->isWithinDistance(g, d))
+        return true;
+
+    // If any point from prepLine is contained by g, the distance is zero
+    // Do this last because this PIP test is not indexed.
+    if ( g->getDimension() == 2 ) {
+        return prepLine.isAnyTargetComponentInTest(g);
+    }
+    return false;
 }
 
 

src/operation/distance/IndexedFacetDistance.cpp

@@ -59,17 +59,26 @@ IndexedFacetDistance::distance(const Geometry* g) const
 bool
 IndexedFacetDistance::isWithinDistance(const Geometry* g, double maxDistance) const
 {
+    std::cout << "geom dist = " << distance(g) << std::endl;
+
     // short-circuit check
     double envDist = baseGeometry.getEnvelopeInternal()->distance(*g->getEnvelopeInternal());
     if (envDist > maxDistance) {
         return false;
     }
-
-    auto env2 = g->getEnvelope();
-    if (distance(env2.get()) > maxDistance) {
-        return false;
+//*
+    //-- heuristic: for atomic indexed geom, test distance to envelope of test geom
+    if (baseGeometry.getNumGeometries() == 1
+        && ! g->getEnvelopeInternal()->contains(baseGeometry.getEnvelopeInternal()))
+    {
+        auto env2 = g->getEnvelope();
+std::cout << "env dist = " << distance(env2.get()) << std::endl;
+        if (distance(env2.get()) > maxDistance) {
+std::cout << "env dist > maxdistance of " << maxDistance << std::endl;
+            return false;
+        }
     }
-
+//*/
     auto tree2 = FacetSequenceTreeBuilder::build(g);
     return cachedTree->isWithinDistance<FacetDistance>(*tree2, maxDistance);
 }

tests/unit/capi/GEOSPreparedDistanceTest.cpp

@@ -178,5 +178,44 @@ void object::test<9>
     );
 }
 
+// Prepared line within envelope of test line
+template<>
+template<>
+void object::test<12>
+()
+{
+    checkDistance(
+        "LINESTRING (1 4, 4 7)",
+        "LINESTRING (1 1, 5 5, 5 9)",
+        1
+    );
+}
+
+// Prepared line within polygon
+template<>
+template<>
+void object::test<13>
+()
+{
+    checkDistance(
+        "LINESTRING (30 30, 70 70)",
+        "POLYGON ((0 100, 100 100, 100 0, 0 0, 0 100))",
+        0
+    );
+}
+
+// Prepared multiline with one element within polygon
+template<>
+template<>
+void object::test<14>
+()
+{
+    checkDistance(
+        "MULTILINESTRING ((30 30, 70 70), (170 200, 200 170))",
+        "POLYGON ((0 100, 100 100, 100 0, 0 0, 0 100))",
+        0
+    );
+}
+
 } // namespace tut
 

tests/unit/capi/GEOSPreparedDistanceWithinTest.cpp

@@ -1,5 +1,5 @@
 //
-// Test Suite for C-API GEOSPreparedDistance
+// Test Suite for C-API GEOSPreparedDistanceWithin
 
 #include <tut/tut.hpp>
 // geos
@@ -213,5 +213,76 @@ void object::test<11>
     );
 }
 
+// Prepared line within envelope of test line
+// see https://github.com/libgeos/geos/issues/958
+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
+    );
+}
+
+// Prepared line within test geometry
+// see https://github.com/libgeos/geos/issues/960
+template<>
+template<>
+void object::test<13>
+()
+{
+    checkDistanceWithin(
+        "LINESTRING (30 30, 70 70)",
+        "POLYGON ((0 100, 100 100, 100 0, 0 0, 0 100))",
+        1,
+        1
+    );
+}
+
+// Prepared multiline with one element within Polygon
+template<>
+template<>
+void object::test<14>
+()
+{
+    checkDistanceWithin(
+        "MULTILINESTRING ((30 30, 70 70), (170 200, 200 170))",
+        "POLYGON ((0 100, 100 100, 100 0, 0 0, 0 100))",
+        1,
+        1
+    );
+}
+
+// Prepared multiline with one element within MultiPolygon.
+template<>
+template<>
+void object::test<15>
+()
+{
+    checkDistanceWithin(
+        "MULTILINESTRING ((1 6, 1 1), (15 16, 15 14))",
+        "MULTIPOLYGON (((10 20, 20 20, 20 10, 10 10, 10 20)), ((30 20, 40 20, 40 10, 30 10, 30 20)))",
+        1,
+        1
+    );
+}
+
+// Indexed multiline with one element within line envelope.
+template<>
+template<>
+void object::test<16>()
+{
+    checkDistanceWithin(
+        "MULTILINESTRING ((1 6, 1 1), (11 14, 11 11))",
+        "LINESTRING (10 10, 10 20, 30 20)",
+        2,
+        1
+    );
+}
+
 } // namespace tut
 

tests/unit/operation/distance/IndexedFacetDistanceTest.cpp

@@ -12,6 +12,7 @@
 // tut
 #include <tut/tut.hpp>
 #include <tut/tut_macros.hpp>
+#include <utility.h>
 // geos
 #include <geos/profiler.h>
 #include <geos/constants.h>
@@ -33,6 +34,7 @@
 #endif
 
 using geos::operation::distance::IndexedFacetDistance;
+using geos::geom::Coordinate;
 
 namespace tut {
 //
@@ -59,19 +61,35 @@ struct test_facetdistanceop_data {
     {}
 
     void
-    checkDistanceNearestPoints(std::string wkt1, std::string wkt2, double distance,
-                               geos::geom::Coordinate& p1, geos::geom::Coordinate& p2)
+    checkDistanceNearestPoints(std::string wkt1, std::string wkt2, double distanceExpected,
+                               const geos::geom::Coordinate& p1, 
+                               const geos::geom::Coordinate& p2)
     {
         using geos::operation::distance::IndexedFacetDistance;
         GeomPtr g1(_wktreader.read(wkt1));
         GeomPtr g2(_wktreader.read(wkt2));
         auto pts = IndexedFacetDistance::nearestPoints(g1.get(), g2.get());
-        ensure(fabs((*pts)[0].distance((*pts)[1])-distance) < 1e08);
-        ensure(fabs((*pts)[0].x - p1.x) < 1e-08);
-        ensure(fabs((*pts)[0].y - p1.y) < 1e-08);
-        ensure(fabs((*pts)[1].x - p2.x) < 1e-08);
-        ensure(fabs((*pts)[1].y - p2.y) < 1e-08);
-        return;
+        const Coordinate p1Actual = (*pts)[0];
+        const Coordinate p2Actual = (*pts)[1];
+        double distResult = (*pts)[0].distance((*pts)[1]);
+
+        ensure_equals("Distance", distResult, distanceExpected, 1e-4);
+        ensure_equals_xy(p1Actual, p1, 1e-08);
+        ensure_equals_xy(p2Actual, p2, 1e-08);
+    }
+
+    void
+    checkWithinDistance(std::string wkt1, std::string wkt2, double distance,
+                               bool expected)
+    {
+        using geos::operation::distance::IndexedFacetDistance;
+        std::unique_ptr<geos::geom::Geometry> g1(_wktreader.read(wkt1));
+        std::unique_ptr<geos::geom::Geometry> g2(_wktreader.read(wkt2));
+
+        IndexedFacetDistance ifd(g1.get());
+        bool result = ifd.isWithinDistance(g2.get(), distance);
+std::cout << "result = " << result << "  expected = " << expected << std::endl;
+        ensure_equals(result, expected);
     }
 
     int
@@ -104,8 +122,6 @@ struct test_facetdistanceop_data {
         return ls;
     }
 
-
-
 };
 
 typedef test_group<test_facetdistanceop_data> group;
@@ -408,6 +424,57 @@ void object::test<12>()
 }
 
 
+// Indexed multiline with one element within line envelope.
+template<>
+template<>
+void object::test<13>()
+{
+    checkWithinDistance(
+        "MULTILINESTRING ((1 6, 1 1), (11 14, 11 11))",
+        "LINESTRING (10 10, 10 20, 30 20)",
+        2,
+        true
+    );
+}
+
+// Indexed multiline with one element within multiline envelope and close to inner line.
+template<>
+template<>
+void object::test<14>()
+{
+    checkWithinDistance(
+        "MULTILINESTRING ((1 6, 1 1), (16 16, 19 13))",
+        "MULTILINESTRING ((10 10, 10 20, 30 20), (20 13, 20 16))",
+        2,
+        true
+    );
+}
+
+// Indexed multiline with one element within line envelope.
+template<>
+template<>
+void object::test<15>
+()
+{
+    checkDistanceNearestPoints(
+        "MULTILINESTRING ((1 6, 1 1), (11 15, 13 17))", 
+        "LINESTRING (10 10, 10 20, 30 20)", 
+        1, Coordinate{11, 15}, Coordinate{10, 15}
+    );
+}
+
+// Indexed multiline with one element within line envelope.
+template<>
+template<>
+void object::test<16>
+()
+{
+    checkDistanceNearestPoints(
+        "MULTILINESTRING ((6 6, 7 7), (100 100, 101 101))", 
+        "LINESTRING (0 10, 10 0)", 
+        1.41421, Coordinate{6, 6}, Coordinate{5, 5}
+    );
+}
 
 // TODO: finish the tests by adding:
 // 	LINESTRING - *all*

tests/unit/utility.h

@@ -96,6 +96,23 @@ instanceOf(InstanceType const* instance)
     return dynamic_cast<Type const*>(instance);
 }
 
+inline void
+ensure_equals_xy(geos::geom::Coordinate const& actual,
+                  geos::geom::Coordinate const& expected)
+{
+    ensure_equals("Coordinate X", actual.x, expected.x );
+    ensure_equals("Coordinate Y", actual.y, expected.y );
+}
+
+inline void
+ensure_equals_xy(geos::geom::Coordinate const& actual,
+                  geos::geom::Coordinate const& expected,
+                  double tol)
+{
+    ensure_equals("Coordinate X", actual.x, expected.x, tol );
+    ensure_equals("Coordinate Y", actual.y, expected.y, tol );
+}
+
 inline void
 ensure_equals_xyz(geos::geom::Coordinate const& actual,
                   geos::geom::Coordinate const& expected)