Port osmium point-in-polygon function

include/extractor/location_dependent_data.hpp

@@ -3,6 +3,7 @@
 
 #include <boost/filesystem/path.hpp>
 #include <boost/geometry.hpp>
+#include <boost/geometry/geometries/point_xy.hpp>
 #include <boost/geometry/index/rtree.hpp>
 
 #include <osmium/osm/way.hpp>
@@ -14,11 +15,14 @@ namespace osrm
 {
 namespace extractor
 {
+
 struct LocationDependentData
 {
-    using point_t = boost::geometry::model::
+    using point_t = boost::geometry::model::d2::
-        point<double, 2, boost::geometry::cs::spherical_equatorial<boost::geometry::degree>>;
+        point_xy<double, boost::geometry::cs::spherical_equatorial<boost::geometry::degree>>;
+    using segment_t = std::pair<point_t, point_t>;
     using polygon_t = boost::geometry::model::polygon<point_t>;
+    using polygon_bands_t = std::vector<std::vector<segment_t>>;
     using box_t = boost::geometry::model::box<point_t>;
 
     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<std::pair<polygon_t, std::size_t>> polygons;
+    std::vector<std::pair<polygon_bands_t, std::size_t>> polygons;
     std::vector<properties_t> properties;
 };
 }

src/extractor/location_dependent_data.cpp

@@ -9,6 +9,7 @@
 
 #include <boost/filesystem.hpp>
 #include <boost/function_output_iterator.hpp>
+#include <boost/geometry/algorithms/equals.hpp>
 
 #include <fstream>
 #include <string>
@@ -109,7 +110,62 @@ void LocationDependentData::loadLocationDependentData(const boost::filesystem::p
         }
         auto envelop = boost::geometry::return_envelope<box_t>(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<segment_t> 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<int32_t>(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<point_t>::type;
+            const std::pair<coord_t, coord_t> mm =
+                std::minmax(segment.first.y(), segment.second.y());
+            const auto band_min = std::min<coord_t>(num_bands - 1, (mm.first - y_min) / dy);
+            const auto band_max = std::min<coord_t>(
+                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)));
 

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(1, 1)).empty());
-    BOOST_CHECK(data(point_t(0, 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.5, -0.5)).empty());
-    BOOST_CHECK(data(point_t(0, -3)).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<double>(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<std::string>(data(point_t(0, 0))["answer"]), "a");
+    BOOST_CHECK_EQUAL(boost::get<std::string>(data(point_t(-3, 3))["answer"]), "a");
+    BOOST_CHECK_EQUAL(boost::get<std::string>(data(point_t(-3, 1))["answer"]), "a");
+    BOOST_CHECK_EQUAL(boost::get<std::string>(data(point_t(-3, -3))["answer"]), "a");
+    BOOST_CHECK_EQUAL(boost::get<std::string>(data(point_t(0, 3))["answer"]), "a");
     BOOST_CHECK_EQUAL(boost::get<std::string>(data(point_t(1, 0))["answer"]), "a");
-    BOOST_CHECK_EQUAL(boost::get<std::string>(data(point_t(2, 0))["answer"]), "a");
+    BOOST_CHECK_EQUAL(boost::get<std::string>(data(point_t(2, -3))["answer"]), "a");
-    BOOST_CHECK_EQUAL(boost::get<std::string>(data(point_t(3, 0))["answer"]), "b");
+    BOOST_CHECK_EQUAL(boost::get<std::string>(data(point_t(3, 0))["answer"]), "a");
-    BOOST_CHECK_EQUAL(boost::get<std::string>(data(point_t(4, 0))["answer"]), "b");
+    BOOST_CHECK_EQUAL(boost::get<std::string>(data(point_t(4, 3))["answer"]), "b");
-    BOOST_CHECK_EQUAL(boost::get<std::string>(data(point_t(6, 0))["answer"]), "b");
+    BOOST_CHECK_EQUAL(boost::get<std::string>(data(point_t(6, 1))["answer"]), "b");
-    BOOST_CHECK_EQUAL(boost::get<std::string>(data(point_t(7, 0))["answer"]), "c");
+    BOOST_CHECK_EQUAL(boost::get<std::string>(data(point_t(7, 0))["answer"]), "b");
+    BOOST_CHECK_EQUAL(boost::get<std::string>(data(point_t(8, 3))["answer"]), "c");
+    BOOST_CHECK_EQUAL(boost::get<std::string>(data(point_t(8, -1))["answer"]), "c");
+    BOOST_CHECK_EQUAL(boost::get<std::string>(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()