Prevent merging of circular-shaped roads

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                                            |
+            | waypoints | route                                      | turns                       |
-            | i,h       | Kurfürstendamm,Rathenauplatz,Hubertusallee,Hubertusallee | depart,turn slight left,turn slight right,arrive |
+            | i,h       | Kurfürstendamm,Hubertusallee,Hubertusallee | depart,turn straight,arrive |

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<Coordinate> &polygon);
+
 } // ns coordinate_calculation
 } // ns util
 } // ns osrm

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(

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<Coordinate> &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

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<Coordinate>{{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<Coordinate>{{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()