From b0514129a158a9b01ae35a15ba6ddc661708cc08 Mon Sep 17 00:00:00 2001 From: "daniel.korn" Date: Sat, 31 Oct 2026 21:08:26 +0100 Subject: [PATCH] perf: skip intersects() for rings with disjoint envelopes in correct() When make_valid repairs a multipolygon, result_combine merges each corrected ring into the result by testing it against every ring added before it with boost::geometry::intersects(). That is quadratic in the number of rings: Lake Huron (r1205151, 14233 inner rings) spends 88 of its 103 s in it. Keep the envelope of each ring alongside and skip intersects() when the envelopes are disjoint, which rules out an intersection. The rings are merged in the same order, so the output is unchanged. Adds a unit test for correct(). Lake Huron 103 s -> 25 s (-76%); Germany 43.6 s -> 39.3 s wall (-10%); planet relations phase 310 s -> 263 s (-15%). Co-Authored-By: Claude Opus 5.5 (1M context) --- Makefile | 5 +++ include/geometry/correct.hpp | 43 ++++++++++++++------- test/correct.test.cpp | 72 ++++++++++++++++++++++++++++++++++++ 3 files changed, 107 insertions(+), 13 deletions(-) create mode 100644 test/correct.test.cpp diff --git a/Makefile b/Makefile index 57227d4c..467393ea 100644 --- a/Makefile +++ b/Makefile @@ -156,6 +156,7 @@ test: \ test_append_vector \ test_attribute_store \ test_config_validator \ + test_correct \ test_deque_map \ test_geom \ test_helpers \ @@ -193,6 +194,10 @@ src/config_schema.h: resources/config-schema.json src/config_validator.o: src/config_schema.h +test_correct: \ + test/correct.test.o + $(CXX) $(CXXFLAGS) -o test.correct $^ $(INC) $(LIB) $(LDFLAGS) && ./test.correct + test_deque_map: \ test/deque_map.test.o $(CXX) $(CXXFLAGS) -o test.deque_map $^ $(INC) $(LIB) $(LDFLAGS) && ./test.deque_map diff --git a/include/geometry/correct.hpp b/include/geometry/correct.hpp index fba96c69..a385d86c 100644 --- a/include/geometry/correct.hpp +++ b/include/geometry/correct.hpp @@ -23,13 +23,18 @@ namespace geometry { namespace impl { -template -static inline void result_combine(C &result, T &&new_element) +// boxes[i] is the envelope of result[i]. Polygons whose envelopes are disjoint +// cannot intersect, so the box test skips most of the costly intersects() calls +// (otherwise quadratic in the number of rings, e.g. a lake with 14000 islands). +template +static inline void result_combine(C &result, std::vector &boxes, T &&new_element) { result.push_back(std::forward(new_element)); + boxes.push_back(boost::geometry::return_envelope(result.back())); for(std::size_t i = 0; i < result.size() - 1; ) { - if(!boost::geometry::intersects(result[i], result.back())) { + if(boost::geometry::disjoint(boxes[i], boxes.back()) || + !boost::geometry::intersects(result[i], result.back())) { ++i; continue; } @@ -43,8 +48,10 @@ static inline void result_combine(C &result, T &&new_element) } result.back() = std::move(union_result[0]); + boxes.back() = boost::geometry::return_envelope(result.back()); result.erase(result.begin() + i); - } + boxes.erase(boxes.begin() + i); + } } struct pseudo_vertice_key @@ -302,13 +309,15 @@ template< > struct combine_non_zero_winding { - inline void operator()(multi_polygon_t &combined_outers, multi_polygon_t &combined_inners, polygon_t &poly) + using box_t = boost::geometry::model::box; + + inline void operator()(multi_polygon_t &combined_outers, std::vector &outer_boxes, multi_polygon_t &combined_inners, std::vector &inner_boxes, polygon_t &poly) { if(boost::geometry::area(poly) > 0) - result_combine(combined_outers, std::move(poly)); + result_combine(combined_outers, outer_boxes, std::move(poly)); else { std::reverse(poly.outer().begin(), poly.outer().end()); - result_combine(combined_inners, std::move(poly)); + result_combine(combined_inners, inner_boxes, std::move(poly)); } } }; @@ -320,7 +329,9 @@ template< > struct combine_odd_even { - inline void operator()(multi_polygon_t &combined_outers, multi_polygon_t &combined_inners, polygon_t &poly) + using box_t = boost::geometry::model::box; + + inline void operator()(multi_polygon_t &combined_outers, std::vector &, multi_polygon_t &combined_inners, std::vector &, polygon_t &poly) { if(boost::geometry::area(poly) < 0) std::reverse(poly.outer().begin(), poly.outer().end()); @@ -346,11 +357,12 @@ static inline void correct(polygon_t const &input, multi_polygon_t &output, doub // Calculate all outers and combine them if possible multi_polygon_t combined_outers; multi_polygon_t combined_inners; + std::vector> outer_boxes, inner_boxes; for(auto &ring: outer_rings) { polygon_t poly; poly.outer() = std::move(ring); - combine(combined_outers, combined_inners, poly); + combine(combined_outers, outer_boxes, combined_inners, inner_boxes, poly); } // Calculate all inners and combine them if possible @@ -362,8 +374,8 @@ static inline void correct(polygon_t const &input, multi_polygon_t &output, doub correct(poly, new_inners, remove_spike_min_area, combine); for(auto &poly: new_inners) { - result_combine(combined_inners, std::move(poly)); - } + result_combine(combined_inners, inner_boxes, std::move(poly)); + } } // Cut out all inners from all the outers @@ -379,13 +391,18 @@ template< > static inline void correct(multi_polygon_t const &input, multi_polygon_t &output, double remove_spike_min_area, combine_function_t combine) { + using box_t = boost::geometry::model::box; + std::vector output_boxes; + for(auto const &polygon: output) + output_boxes.push_back(boost::geometry::return_envelope(polygon)); + for(auto const &polygon: input) { multi_polygon_t new_polygons; correct(polygon, new_polygons, remove_spike_min_area, combine); - for(auto &new_polygon: new_polygons) - result_combine(output, std::move(new_polygon)); + for(auto &new_polygon: new_polygons) + result_combine(output, output_boxes, std::move(new_polygon)); } } diff --git a/test/correct.test.cpp b/test/correct.test.cpp new file mode 100644 index 00000000..7178d173 --- /dev/null +++ b/test/correct.test.cpp @@ -0,0 +1,72 @@ +#include +#include +#include +#include "external/minunit.h" +#include "geometry/correct.hpp" + +typedef boost::geometry::model::d2::point_xy point_t; +typedef boost::geometry::model::polygon polygon_t; +typedef boost::geometry::model::multi_polygon multi_polygon_t; + +static bool near(double a, double b) { return std::abs(a - b) < 1e-9; } + +MU_TEST(test_correct_merges_overlapping_holes) { + // Two overlapping holes are merged into one; a disjoint hole is kept as is + polygon_t input; + boost::geometry::read_wkt( + "POLYGON((0 0,0 10,10 10,10 0,0 0)," + "(1 1,3 1,3 3,1 3,1 1),(2 2,4 2,4 4,2 4,2 2),(6 6,7 6,7 7,6 7,6 6))", input); + multi_polygon_t output; + geometry::correct(input, output); + + mu_check(output.size() == 1); + mu_check(output[0].inners().size() == 2); + mu_check(near(boost::geometry::area(output), 100 - 7 - 1)); + mu_check(boost::geometry::is_valid(output)); +} + +MU_TEST(test_correct_merges_overlapping_outers) { + multi_polygon_t input; + boost::geometry::read_wkt( + "MULTIPOLYGON(((0 0,0 2,2 2,2 0,0 0)),((1 1,1 3,3 3,3 1,1 1)),((5 5,5 6,6 6,6 5,5 5)))", input); + multi_polygon_t output; + geometry::correct(input, output); + + mu_check(output.size() == 2); + mu_check(near(boost::geometry::area(output), 4 + 4 - 1 + 1)); + mu_check(boost::geometry::is_valid(output)); +} + +MU_TEST(test_correct_many_disjoint_holes) { + // 400 disjoint holes, e.g. a lake with many islands: all are kept + std::string wkt = "POLYGON((0 0,0 100,100 100,100 0,0 0)"; + for (int x = 0; x < 20; x++) + for (int y = 0; y < 20; y++) { + double x0 = 1 + x * 5, y0 = 1 + y * 5; + wkt += ",(" + std::to_string(x0) + " " + std::to_string(y0) + "," + std::to_string(x0 + 2) + " " + std::to_string(y0) + "," + + std::to_string(x0 + 2) + " " + std::to_string(y0 + 2) + "," + std::to_string(x0) + " " + std::to_string(y0 + 2) + "," + + std::to_string(x0) + " " + std::to_string(y0) + ")"; + } + wkt += ")"; + polygon_t input; + boost::geometry::read_wkt(wkt, input); + multi_polygon_t output; + geometry::correct(input, output); + + mu_check(output.size() == 1); + mu_check(output[0].inners().size() == 400); + mu_check(near(boost::geometry::area(output), 100 * 100 - 400 * 4)); + mu_check(boost::geometry::is_valid(output)); +} + +MU_TEST_SUITE(test_suite_correct) { + MU_RUN_TEST(test_correct_merges_overlapping_holes); + MU_RUN_TEST(test_correct_merges_overlapping_outers); + MU_RUN_TEST(test_correct_many_disjoint_holes); +} + +int main() { + MU_RUN_SUITE(test_suite_correct); + MU_REPORT(); + return MU_EXIT_CODE; +}