From 421115200bd815019c73675fabb6e6b75bf75b9f Mon Sep 17 00:00:00 2001 From: Michael Krasnyk Date: Wed, 30 Aug 2017 16:33:23 +0200 Subject: [PATCH] Port osmium point-in-polygon function --- include/extractor/location_dependent_data.hpp | 10 +- src/extractor/location_dependent_data.cpp | 112 +++++++++++++++++- .../location_dependent_data_tests.cpp | 68 +++++++++-- 3 files changed, 175 insertions(+), 15 deletions(-) diff --git a/include/extractor/location_dependent_data.hpp b/include/extractor/location_dependent_data.hpp index b1638b093..c76244e46 100644 --- a/include/extractor/location_dependent_data.hpp +++ b/include/extractor/location_dependent_data.hpp @@ -3,6 +3,7 @@ #include #include +#include #include #include @@ -14,11 +15,14 @@ namespace osrm { namespace extractor { + struct LocationDependentData { - using point_t = boost::geometry::model:: - point>; + using point_t = boost::geometry::model::d2:: + point_xy>; + using segment_t = std::pair; using polygon_t = boost::geometry::model::polygon; + using polygon_bands_t = std::vector>; using box_t = boost::geometry::model::box; using polygon_position_t = std::size_t; @@ -42,7 +46,7 @@ struct LocationDependentData void loadLocationDependentData(const boost::filesystem::path &file_path); rtree_t rtree; - std::vector> polygons; + std::vector> polygons; std::vector properties; }; } diff --git a/src/extractor/location_dependent_data.cpp b/src/extractor/location_dependent_data.cpp index c4355d8e4..c11884225 100644 --- a/src/extractor/location_dependent_data.cpp +++ b/src/extractor/location_dependent_data.cpp @@ -9,6 +9,7 @@ #include #include +#include #include #include @@ -109,7 +110,62 @@ void LocationDependentData::loadLocationDependentData(const boost::filesystem::p } auto envelop = boost::geometry::return_envelope(polygon); bounding_boxes.emplace_back(envelop, polygons.size()); - polygons.emplace_back(std::make_pair(polygon, properties_index)); + + // here is a part of ExtractPolygon::ExtractPolygon code + // TODO: remove copy overhead + constexpr const int32_t segments_per_band = 10; + constexpr const int32_t max_bands = 10000; + const auto y_min = envelop.min_corner().y(); + const auto y_max = envelop.max_corner().y(); + + std::vector segments; + auto add_ring = [&segments](const auto &ring) { + auto it = ring.begin(); + const auto end = ring.end(); + + BOOST_ASSERT(it != end); + auto last_it = it++; + while (it != end) + { + segments.emplace_back(*last_it, *it); + last_it = it++; + } + }; + + add_ring(polygon.outer()); + for (const auto &ring : polygon.inners()) + add_ring(ring); + + int32_t num_bands = static_cast(segments.size()) / segments_per_band; + if (num_bands < 1) + { + num_bands = 1; + } + else if (num_bands > max_bands) + { + num_bands = max_bands; + } + + polygon_bands_t bands(num_bands); + const auto dy = (y_max - y_min) / num_bands; + + for (const auto &segment : segments) + { + using coord_t = boost::geometry::traits::coordinate_type::type; + const std::pair mm = + std::minmax(segment.first.y(), segment.second.y()); + const auto band_min = std::min(num_bands - 1, (mm.first - y_min) / dy); + const auto band_max = std::min( + num_bands, ((mm.second - y_min) / dy) + 1); // TODO: use integer coordinates + + for (auto band = band_min; band < band_max; ++band) + { + bands[band].push_back(segment); + } + } + // EOC + + polygons.emplace_back(std::make_pair(bands, properties_index)); }; for (rapidjson::SizeType ifeature = 0; ifeature < features_array.Size(); ifeature++) @@ -156,7 +212,59 @@ LocationDependentData::properties_t LocationDependentData::operator()(const poin // Search the R-tree and collect a Lua table of tags that correspond to the location rtree.query(boost::geometry::index::intersects(point) && boost::geometry::index::satisfies([this, &point](const rtree_t::value_type &v) { - return boost::geometry::within(point, polygons[v.second].first); + + // Simple point-in-polygon algorithm adapted from + // https://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html + + const auto &envelop = v.first; + const auto &bands = polygons[v.second].first; + + const auto y_min = envelop.min_corner().y(); + const auto y_max = envelop.max_corner().y(); + auto dy = (y_max - y_min) / bands.size(); + + std::size_t band = (point.y() - y_min) / dy; + if (band >= bands.size()) + { + band = bands.size() - 1; + } + + bool inside = false; + + for (const auto &segment : bands[band]) + { + const auto point_x = point.x(), point_y = point.y(); + const auto from_x = segment.first.x(), from_y = segment.first.y(); + const auto to_x = segment.second.x(), to_y = segment.second.y(); + + if (to_y == from_y) + { // handle horizontal segments: check if on boundary or skip + if ((to_y == point_y) && + (from_x == point_x || (to_x > point_x) != (from_x > point_x))) + return true; + continue; + } + + if ((to_y > point_y) != (from_y > point_y)) + { + const auto ax = to_x - from_x; + const auto ay = to_y - from_y; + const auto tx = point_x - from_x; + const auto ty = point_y - from_y; + + const auto cross_product = tx * ay - ax * ty; + + if (cross_product == 0) + return true; + + if ((ay > 0) == (cross_product > 0)) + { + inside = !inside; + } + } + } + + return inside; }), boost::make_function_output_iterator(std::ref(merger))); diff --git a/unit_tests/extractor/location_dependent_data_tests.cpp b/unit_tests/extractor/location_dependent_data_tests.cpp index 32972c508..f1d21dcb0 100644 --- a/unit_tests/extractor/location_dependent_data_tests.cpp +++ b/unit_tests/extractor/location_dependent_data_tests.cpp @@ -51,10 +51,10 @@ BOOST_AUTO_TEST_CASE(polygon_tests) LocationDependentData data(fixture.temporary_file); BOOST_CHECK(data(point_t(0, 0)).empty()); - BOOST_CHECK(data(point_t(1, 1)).empty()); - BOOST_CHECK(data(point_t(0, 1)).empty()); - BOOST_CHECK(data(point_t(0.5, -0.5)).empty()); - BOOST_CHECK(data(point_t(0, -3)).empty()); + BOOST_CHECK(!data(point_t(1, 1)).empty()); + BOOST_CHECK(!data(point_t(0, 1)).empty()); + BOOST_CHECK(!data(point_t(0.5, -0.5)).empty()); + BOOST_CHECK(!data(point_t(0, -3)).empty()); BOOST_CHECK(!data(point_t(-0.75, 0.75)).empty()); BOOST_CHECK(!data(point_t(-2, 6)).empty()); BOOST_CHECK_EQUAL(boost::get(data(point_t(2, 0))["answer"]), 42.); @@ -110,13 +110,61 @@ BOOST_AUTO_TEST_CASE(polygon_merging_tests) ]})json"); LocationDependentData data(fixture.temporary_file); - BOOST_CHECK_EQUAL(boost::get(data(point_t(0, 0))["answer"]), "a"); + BOOST_CHECK_EQUAL(boost::get(data(point_t(-3, 3))["answer"]), "a"); + BOOST_CHECK_EQUAL(boost::get(data(point_t(-3, 1))["answer"]), "a"); + BOOST_CHECK_EQUAL(boost::get(data(point_t(-3, -3))["answer"]), "a"); + BOOST_CHECK_EQUAL(boost::get(data(point_t(0, 3))["answer"]), "a"); BOOST_CHECK_EQUAL(boost::get(data(point_t(1, 0))["answer"]), "a"); - BOOST_CHECK_EQUAL(boost::get(data(point_t(2, 0))["answer"]), "a"); - BOOST_CHECK_EQUAL(boost::get(data(point_t(3, 0))["answer"]), "b"); - BOOST_CHECK_EQUAL(boost::get(data(point_t(4, 0))["answer"]), "b"); - BOOST_CHECK_EQUAL(boost::get(data(point_t(6, 0))["answer"]), "b"); - BOOST_CHECK_EQUAL(boost::get(data(point_t(7, 0))["answer"]), "c"); + BOOST_CHECK_EQUAL(boost::get(data(point_t(2, -3))["answer"]), "a"); + BOOST_CHECK_EQUAL(boost::get(data(point_t(3, 0))["answer"]), "a"); + BOOST_CHECK_EQUAL(boost::get(data(point_t(4, 3))["answer"]), "b"); + BOOST_CHECK_EQUAL(boost::get(data(point_t(6, 1))["answer"]), "b"); + BOOST_CHECK_EQUAL(boost::get(data(point_t(7, 0))["answer"]), "b"); + BOOST_CHECK_EQUAL(boost::get(data(point_t(8, 3))["answer"]), "c"); + BOOST_CHECK_EQUAL(boost::get(data(point_t(8, -1))["answer"]), "c"); + BOOST_CHECK_EQUAL(boost::get(data(point_t(8, -3))["answer"]), "c"); +} + +BOOST_AUTO_TEST_CASE(staircase_polygon) +{ + LocationDataFixture fixture(R"json({ +"type": "FeatureCollection", +"features": [ +{ + "type": "Feature", + "properties": { "answer": "a" }, + "geometry": { "type": "Polygon", "coordinates": [ [ [0, 0], [3, 0], [3, 3], [2, 3], [2, 2], [1, 2], [1, 1], [0, 1], [0, 0] ] ] } +} +]})json"); + + LocationDependentData data(fixture.temporary_file); + + // all corners + BOOST_CHECK(!data(point_t(0, 0)).empty()); + BOOST_CHECK(!data(point_t(0, 1)).empty()); + BOOST_CHECK(!data(point_t(1, 1)).empty()); + BOOST_CHECK(!data(point_t(1, 2)).empty()); + BOOST_CHECK(!data(point_t(2, 2)).empty()); + BOOST_CHECK(!data(point_t(2, 3)).empty()); + BOOST_CHECK(!data(point_t(3, 3)).empty()); + BOOST_CHECK(!data(point_t(3, 0)).empty()); + + // at x = 1 + BOOST_CHECK(data(point_t(1, -0.5)).empty()); + BOOST_CHECK(!data(point_t(1, 0)).empty()); + BOOST_CHECK(!data(point_t(1, 0.5)).empty()); + BOOST_CHECK(!data(point_t(1, 1.5)).empty()); + BOOST_CHECK(data(point_t(1, 2.5)).empty()); + BOOST_CHECK(data(point_t(3.5, 2)).empty()); + + // at y = 2 + BOOST_CHECK(data(point_t(0.5, 2)).empty()); + BOOST_CHECK(!data(point_t(1, 2)).empty()); + BOOST_CHECK(!data(point_t(1.5, 2)).empty()); + BOOST_CHECK(!data(point_t(2, 2)).empty()); + BOOST_CHECK(!data(point_t(2.5, 2)).empty()); + BOOST_CHECK(!data(point_t(3, 2)).empty()); + BOOST_CHECK(data(point_t(3.5, 2)).empty()); } BOOST_AUTO_TEST_SUITE_END()