Port of JTS [RelateNG](locationtech/jts#1052), rewrite of boolean pre…

…dicates and relate matrix calculations. * No new functions, just rework of existing functionality * "Prepared" mode now available for all predicates and relate matrix * CAPI functions GEOSPreparedRelate and GEOSPreparedRelatePattern expose new functionality * CAPI implementations of GEOSPreparedTouches, etc, that were previously defaulting into non-prepared implementations now default into the RelateNG prepared implementation * Prepared implementations for Intersects, Covers, still use the older implementations * https://lin-ear-th-inking.blogspot.com/2024/05/jts-topological-relationships-next.html * https://lin-ear-th-inking.blogspot.com/2024/05/relateng-performance.html
libgeos · Aug 13, 2024 · d963c0f · d963c0f
1 parent 85b06ab
commit d963c0f
Show file tree

Hide file tree

Showing 93 changed files with 12,202 additions and 200 deletions.
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
@@ -273,7 +273,7 @@ jobs:
       working-directory: ./build
       run: |
         ctest --output-on-failure \
-          --overwrite MemoryCheckCommandOptions="--leak-check=full --error-exitcode=100" \
+          --overwrite MemoryCheckCommandOptions="--leak-check=full --error-exitcode=100 --num-callers=100" \
           -R ^all-unit -C Valgrind -T memcheck
 
     - name: 'Upload Valgrind Log'

diff --git a/benchmarks/BenchmarkUtils.h b/benchmarks/BenchmarkUtils.h
@@ -58,6 +58,12 @@ createLine(const geom::CoordinateXY& base, double size, std::size_t npts) {
 }
 
 
+std::vector<std::unique_ptr<geom::Geometry>>
+createPolygons(const geom::Envelope& env, std::size_t nItems, double size, std::size_t npts) {
+    return createGeometriesOnGrid(env, nItems, [size, npts](const geom::CoordinateXY& base) {
+        return createSineStar(base, size, npts);
+    });
+}
 std::vector<std::unique_ptr<geom::Geometry>>
 createLines(const geom::Envelope& env, std::size_t nItems, double size, std::size_t npts) {
     return createGeometriesOnGrid(env, nItems, [size, npts](const geom::CoordinateXY& base) {

diff --git a/benchmarks/geom/CMakeLists.txt b/benchmarks/geom/CMakeLists.txt
@@ -32,8 +32,8 @@ IF(benchmark_FOUND)
             benchmark::benchmark geos_cxx_flags)
 endif()
 
-add_executable(perf_prepared_polygon_intersects
-    PreparedPolygonIntersectsPerfTest.cpp)
-target_include_directories(perf_prepared_polygon_intersects PUBLIC
+add_executable(perf_topo_predicate
+    TopologyPredicatePerfTest.cpp)
+target_include_directories(perf_topo_predicate PUBLIC
     $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/benchmarks>)
-target_link_libraries(perf_prepared_polygon_intersects PRIVATE geos)
+target_link_libraries(perf_topo_predicate PRIVATE geos)
diff --git a/benchmarks/geom/PreparedPolygonIntersectsPerfTest.cpp b/benchmarks/geom/PreparedPolygonIntersectsPerfTest.cpp
diff --git a/benchmarks/geom/TopologyPredicatePerfTest.cpp b/benchmarks/geom/TopologyPredicatePerfTest.cpp
@@ -0,0 +1,266 @@
+/******************************************************
+ *   Performance tests for topological predicates
+ * 
+ * Usage: perf_topo_predicate [ intersects | contains ]
+******************************************************/
+
+#include <geos/geom/util/SineStarFactory.h>
+#include <geos/geom/GeometryFactory.h>
+#include <geos/profiler.h>
+#include <geos/geom/IntersectionMatrix.h>
+#include <geos/geom/prep/PreparedGeometryFactory.h>
+#include <geos/operation/relateng/RelateNG.h>
+#include <geos/operation/relateng/RelatePredicate.h>
+#include <geos/io/WKBWriter.h>
+#include <BenchmarkUtils.h>
+
+#include <iomanip>
+
+using namespace geos::geom;
+
+std::size_t MAX_ITER = 10;
+std::size_t NUM_LINES = 10000;
+std::size_t NUM_LINES_PTS = 100;
+
+#define INTERSECTS 0
+#define CONTAINS 1
+#define COVERS 2
+#define TOUCHES 3
+
+std::string predicateName{"intersects"};
+int predicateOp = INTERSECTS;
+
+int testRelateOpIntersects(const Geometry& g, const std::vector<std::unique_ptr<Geometry>>& geoms) {
+    int count = 0;
+    for (const auto& geom : geoms) {
+        auto im = g.relate(geom.get());
+        count += im->isIntersects();
+    }
+    return count;
+}
+
+int testRelateOpContains(const Geometry& g, const std::vector<std::unique_ptr<Geometry>>& geoms) {
+    int count = 0;
+    for (const auto& geom : geoms) {
+        auto im = g.relate(geom.get());
+        count += im->isContains();
+    }
+    return count;
+}
+
+int testRelateOpCovers(const Geometry& g, const std::vector<std::unique_ptr<Geometry>>& geoms) {
+    int count = 0;
+    for (const auto& geom : geoms) {
+        auto im = g.relate(geom.get());
+        count += im->isCovers();
+    }
+    return count;
+}
+
+int testRelateOpTouches(const Geometry& g, const std::vector<std::unique_ptr<Geometry>>& geoms) {
+    int count = 0;
+    for (const auto& geom : geoms) {
+        auto im = g.relate(geom.get());
+        count += im->isTouches(g.getDimension(), geom.get()->getDimension());
+    }
+    return count;
+}
+
+int testGeometryIntersects(const Geometry& g, const std::vector<std::unique_ptr<Geometry>>& geoms) {
+    int count = 0;
+    for (const auto& geom : geoms) {
+        count += g.intersects(geom.get());
+    }
+    return count;
+}
+
+int testGeometryContains(const Geometry& g, const std::vector<std::unique_ptr<Geometry>>& geoms) {
+    int count = 0;
+    for (const auto& geom : geoms) {
+        count += g.contains(geom.get());
+    }
+    return count;
+}
+
+int testGeometryCovers(const Geometry& g, const std::vector<std::unique_ptr<Geometry>>& geoms) {
+    int count = 0;
+    for (const auto& geom : geoms) {
+        count += g.covers(geom.get());
+    }
+    return count;
+}
+
+int testGeometryTouches(const Geometry& g, const std::vector<std::unique_ptr<Geometry>>& geoms) {
+    int count = 0;
+    for (const auto& geom : geoms) {
+        count += g.touches(geom.get());
+    }
+    return count;
+}
+
+int testPrepGeomIntersects(const Geometry& g, const std::vector<std::unique_ptr<Geometry>>& geoms) {
+    int count = 0;
+    auto prep = prep::PreparedGeometryFactory::prepare(&g);
+    for (const auto& geom : geoms) {
+        count += prep->intersects(geom.get());
+    }
+    return count;
+}
+
+int testPrepGeomContains(const Geometry& g, const std::vector<std::unique_ptr<Geometry>>& geoms) {
+    int count = 0;
+    auto prep = prep::PreparedGeometryFactory::prepare(&g);
+    for (const auto& geom : geoms) {
+        count += prep->contains(geom.get());
+    }
+    return count;
+}
+
+int testPrepGeomCovers(const Geometry& g, const std::vector<std::unique_ptr<Geometry>>& geoms) {
+    int count = 0;
+    auto prep = prep::PreparedGeometryFactory::prepare(&g);
+    for (const auto& geom : geoms) {
+        count += prep->covers(geom.get());
+    }
+    return count;
+}
+
+int testRelateNGPreparedIntersects(const Geometry& g, const std::vector<std::unique_ptr<Geometry>>& geoms) {
+    int count = 0;
+    auto prep = geos::operation::relateng::RelateNG::prepare(&g);
+    for (const auto& line : geoms) {
+        count += prep->evaluate(line.get(), *geos::operation::relateng::RelatePredicate::intersects());
+    }
+    return count;
+}
+
+int testRelateNGPreparedContains(const Geometry& g, const std::vector<std::unique_ptr<Geometry>>& geoms) {
+    int count = 0;
+    auto prep = geos::operation::relateng::RelateNG::prepare(&g);
+    for (const auto& line : geoms) {
+        count += prep->evaluate(line.get(), *geos::operation::relateng::RelatePredicate::contains());
+    }
+    return count;
+}
+
+int testRelateNGPreparedCovers(const Geometry& g, const std::vector<std::unique_ptr<Geometry>>& geoms) {
+    int count = 0;
+    auto prep = geos::operation::relateng::RelateNG::prepare(&g);
+    for (const auto& line : geoms) {
+        count += prep->evaluate(line.get(), *geos::operation::relateng::RelatePredicate::covers());
+    }
+    return count;
+}
+
+int testRelateNGPreparedTouches(const Geometry& g, const std::vector<std::unique_ptr<Geometry>>& geoms) {
+    int count = 0;
+    auto prep = geos::operation::relateng::RelateNG::prepare(&g);
+    for (const auto& line : geoms) {
+        count += prep->evaluate(line.get(), *geos::operation::relateng::RelatePredicate::touches());
+    }
+    return count;
+}
+
+template<typename F>
+double test(const Geometry& g, const std::vector<std::unique_ptr<Geometry>>& geoms, const std::string& method, F&& fun, double baseTime)
+{
+    geos::util::Profile sw("TopologuyPredicatePerf");
+    sw.start();
+
+    int count = 0;
+    for (std::size_t i = 0; i < MAX_ITER; i++) {
+        count += fun(g, geoms);
+    }
+
+    sw.stop();
+    double tot = sw.getTot();
+    double timesFaster = baseTime == 0 ? 1 : baseTime / tot;
+    std::cout << std::fixed << std::setprecision(0);
+    std::cout << g.getNumPoints() << "," 
+        << MAX_ITER * geoms.size() << "," 
+        << count << "," << geoms[0]->getGeometryType() << "," 
+        << geoms[0]->getNumPoints() << "," 
+        << method << " - " << predicateName << "," 
+        << tot << ",";
+    std::cout << std::fixed << std::setprecision(1);
+    std::cout << timesFaster 
+        << std::endl;
+    return tot;
+}
+
+void test(int dim, std::size_t npts) {
+
+    auto target = geos::benchmark::createSineStar({0, 0}, 100, npts);
+    std::vector<std::unique_ptr<Geometry>> geoms;
+    switch (dim) {
+    case 0:
+        geoms = geos::benchmark::createPoints(*target->getEnvelopeInternal(), NUM_LINES);
+        break;
+    case 1:
+        geoms = geos::benchmark::createLines(*target->getEnvelopeInternal(), NUM_LINES, 1.0, NUM_LINES_PTS);
+        break;
+    case 2:
+        geoms = geos::benchmark::createPolygons(*target->getEnvelopeInternal(), NUM_LINES, 1.0, NUM_LINES_PTS);
+        break;
+    }
+    double baseTime;
+    switch (predicateOp) {
+    case INTERSECTS:
+        baseTime = test(*target, geoms, "RelateOp", testRelateOpIntersects, 0);
+        test(*target, geoms, "Geometry", testGeometryIntersects, baseTime);
+        test(*target, geoms, "PreparedGeom", testPrepGeomIntersects, baseTime);
+        test(*target, geoms, "RelateNGPrepared", testRelateNGPreparedIntersects, baseTime);
+        break;
+    case CONTAINS:
+        baseTime = test(*target, geoms, "RelateOp", testRelateOpIntersects, 0);
+        test(*target, geoms, "Geometry", testGeometryIntersects, baseTime);
+        test(*target, geoms, "PreparedGeom", testPrepGeomContains, baseTime);
+        test(*target, geoms, "RelateNGPrepared", testRelateNGPreparedContains, baseTime);
+        break;
+    case COVERS:
+        baseTime = test(*target, geoms, "RelateOp", testRelateOpCovers, 0);
+        test(*target, geoms, "Geometry", testGeometryCovers, baseTime);
+        test(*target, geoms, "PreparedGeom", testPrepGeomCovers, baseTime);
+        test(*target, geoms, "RelateNGPrepared", testRelateNGPreparedCovers, baseTime);
+        break;
+    case TOUCHES:
+        baseTime = test(*target, geoms, "RelateOp", testRelateOpTouches, 0);
+        test(*target, geoms, "Geometry", testGeometryTouches, baseTime);
+        test(*target, geoms, "RelateNGPrepared", testRelateNGPreparedTouches, baseTime);
+        break;
+    }
+}
+
+void testAll(int dim)
+{
+    test(dim, 5);
+    test(dim, 10);
+    test(dim, 500);
+    test(dim, 1000);
+    test(dim, 2000);
+    test(dim, 4000);
+    test(dim, 8000);
+    test(dim, 16000);
+}
+
+int main(int argc, char** argv) {
+    predicateName = "intersects";
+    if (argc >= 2) {
+        predicateName = argv[1];
+    }
+    predicateOp = INTERSECTS;
+    if (predicateName == "contains") {
+        predicateOp = CONTAINS;
+    }
+    else if (predicateName == "covers") {
+        predicateOp = COVERS;
+    }
+    else if (predicateName == "touches") {
+        predicateOp = TOUCHES;
+    }
+
+    std::cout << "target_points,num_tests,num_hits,test_type,pts_in_test,method,time,factor" << std::endl;
+    testAll(0);
+    testAll(1);
+    testAll(2);
+}