From 088d4edc6b761d2642a9b1ab7d85b701d3165fb2 Mon Sep 17 00:00:00 2001 From: Michael Krasnyk Date: Mon, 16 Oct 2017 10:47:14 +0200 Subject: [PATCH] Prevent merging of circular-shaped roads --- .../guidance/merge-segregated-roads.feature | 4 +- include/util/coordinate_calculation.hpp | 2 + .../guidance/mergable_road_detector.cpp | 32 +++++++++++++- src/util/coordinate_calculation.cpp | 42 +++++++++++++++++++ unit_tests/util/coordinate_calculation.cpp | 25 +++++++++++ 5 files changed, 102 insertions(+), 3 deletions(-) diff --git a/features/guidance/merge-segregated-roads.feature b/features/guidance/merge-segregated-roads.feature index 55a15e49e..441dd73b3 100644 --- a/features/guidance/merge-segregated-roads.feature +++ b/features/guidance/merge-segregated-roads.feature @@ -560,5 +560,5 @@ Feature: Merge Segregated Roads | jd | Hubertusallee | yes | When I route I should get - | waypoints | route | turns | - | i,h | Kurfürstendamm,Rathenauplatz,Hubertusallee,Hubertusallee | depart,turn slight left,turn slight right,arrive | + | waypoints | route | turns | + | i,h | Kurfürstendamm,Hubertusallee,Hubertusallee | depart,turn straight,arrive | diff --git a/include/util/coordinate_calculation.hpp b/include/util/coordinate_calculation.hpp index a824ba8df..43f331a9f 100644 --- a/include/util/coordinate_calculation.hpp +++ b/include/util/coordinate_calculation.hpp @@ -378,6 +378,8 @@ bool areParallel(const iterator_type lhs_begin, return std::abs(slope_rhs) < 0.20; // twenty percent incline at the most } +double computeArea(const std::vector &polygon); + } // ns coordinate_calculation } // ns util } // ns osrm diff --git a/src/extractor/guidance/mergable_road_detector.cpp b/src/extractor/guidance/mergable_road_detector.cpp index 1248b2870..05fc729f5 100644 --- a/src/extractor/guidance/mergable_road_detector.cpp +++ b/src/extractor/guidance/mergable_road_detector.cpp @@ -297,7 +297,7 @@ bool MergableRoadDetector::HaveSameDirection(const NodeID intersection_node, // many roads only do short parallel segments. To get a good impression of how `parallel` two // roads are, we look 100 meters down the road (wich can be quite short for very broad roads). - const double constexpr distance_to_extract = 100; + const double constexpr distance_to_extract = 150; std::tie(distance_traversed_to_the_left, coordinates_to_the_left) = getCoordinatesAlongWay(lhs.eid, distance_to_extract); @@ -317,6 +317,36 @@ bool MergableRoadDetector::HaveSameDirection(const NodeID intersection_node, const auto connect_again = (coordinates_to_the_left.back() == coordinates_to_the_right.back()); + // Tuning parameter to detect and don't merge roads close to circular shapes + // if the area to squared circumference ratio is between the lower bound and 1/(4π) + // that correspond to isoperimetric inequality 4πA ≤ L² or lower bound ≤ A/L² ≤ 1/(4π). + // The lower bound must be larger enough to allow merging of square-shaped intersections + // with A/L² = 1/16 or 78.6% + // The condition suppresses roads merging for intersections like + // . . + // . . + // ---- ---- + // . . + // . . + // but will allow roads merging for intersections like + // ------- + // / \  + // ---- ---- + // \ / + // ------- + const auto constexpr CIRCULAR_POLYGON_ISOPERIMETRIC_LOWER_BOUND = 0.85 / (4 * M_PI); + if (connect_again && coordinates_to_the_left.front() == coordinates_to_the_left.back()) + { // if the left and right roads connect again and are closed polygons ... + const auto area = util::coordinate_calculation::computeArea(coordinates_to_the_left); + const auto perimeter = distance_traversed_to_the_left; + const auto area_to_squared_perimeter_ratio = std::abs(area) / (perimeter * perimeter); + + // then don't merge roads if A/L² is greater than the lower bound + BOOST_ASSERT(area_to_squared_perimeter_ratio <= 1. / (4 * M_PI)); + if (area_to_squared_perimeter_ratio >= CIRCULAR_POLYGON_ISOPERIMETRIC_LOWER_BOUND) + return false; + } + // sampling to correctly weight longer segments in regression calculations const auto constexpr SAMPLE_INTERVAL = 5; coordinates_to_the_left = coordinate_extractor.SampleCoordinates( diff --git a/src/util/coordinate_calculation.cpp b/src/util/coordinate_calculation.cpp index 015aec0d3..ae968203d 100644 --- a/src/util/coordinate_calculation.cpp +++ b/src/util/coordinate_calculation.cpp @@ -381,6 +381,48 @@ Coordinate difference(const Coordinate lhs, const Coordinate rhs) return {util::FixedLongitude{lon_diff_int}, util::FixedLatitude{lat_diff_int}}; } +double computeArea(const std::vector &polygon) +{ + using util::coordinate_calculation::haversineDistance; + + if (polygon.empty()) + return 0.; + + BOOST_ASSERT(polygon.front() == polygon.back()); + + // Take the reference point with the smallest latitude. + // ⚠ ref_latitude is the standard parallel for the equirectangular projection + // that is not an area-preserving projection + const auto ref_point = + std::min_element(polygon.begin(), polygon.end(), [](const auto &lhs, const auto &rhs) { + return lhs.lat < rhs.lat; + }); + const auto ref_latitude = ref_point->lat; + + // Compute area of under a curve and a line that is parallel the equator with ref_latitude + // For closed curves it corresponds to the shoelace algorithm for polygon areas + double area = 0.; + auto first = polygon.begin(); + auto previous_base = util::Coordinate{first->lon, ref_latitude}; + auto previous_y = haversineDistance(previous_base, *first); + for (++first; first != polygon.end(); ++first) + { + BOOST_ASSERT(first->lat >= ref_latitude); + + const auto current_base = util::Coordinate{first->lon, ref_latitude}; + const auto current_y = haversineDistance(current_base, *first); + const auto chunk_area = + haversineDistance(previous_base, current_base) * (previous_y + current_y); + + area += (current_base.lon >= previous_base.lon) ? chunk_area : -chunk_area; + + previous_base = current_base; + previous_y = current_y; + } + + return area / 2.; +} + } // ns coordinate_calculation } // ns util } // ns osrm diff --git a/unit_tests/util/coordinate_calculation.cpp b/unit_tests/util/coordinate_calculation.cpp index 5dcf71c3a..add13d2a1 100644 --- a/unit_tests/util/coordinate_calculation.cpp +++ b/unit_tests/util/coordinate_calculation.cpp @@ -408,4 +408,29 @@ BOOST_AUTO_TEST_CASE(regression_test_3516) BOOST_CHECK_EQUAL(nearest_location, v); } +BOOST_AUTO_TEST_CASE(computeArea) +{ + using osrm::util::coordinate_calculation::computeArea; + + // + auto rhombus = std::vector{{FloatLongitude{.00}, FloatLatitude{.00}}, + {FloatLongitude{.01}, FloatLatitude{.01}}, + {FloatLongitude{.02}, FloatLatitude{.00}}, + {FloatLongitude{.01}, FloatLatitude{-.01}}, + {FloatLongitude{.00}, FloatLatitude{.00}}}; + + BOOST_CHECK_CLOSE(2 * 1112.263 * 1112.263, computeArea(rhombus), 1e-3); + + // edge cases + auto self_intersection = std::vector{{FloatLongitude{.00}, FloatLatitude{.00}}, + {FloatLongitude{.00}, FloatLatitude{.02}}, + {FloatLongitude{.01}, FloatLatitude{.01}}, + {FloatLongitude{.02}, FloatLatitude{.00}}, + {FloatLongitude{.02}, FloatLatitude{.02}}, + {FloatLongitude{.01}, FloatLatitude{.01}}, + {FloatLongitude{.00}, FloatLatitude{.00}}}; + BOOST_CHECK(computeArea(self_intersection) < 1e-3); + BOOST_CHECK_CLOSE(0, computeArea({}), 1e-3); +} + BOOST_AUTO_TEST_SUITE_END()