diff --git a/include/util/coordinate_calculation.hpp b/include/util/coordinate_calculation.hpp index 48a5753e277..86a2e836123 100644 --- a/include/util/coordinate_calculation.hpp +++ b/include/util/coordinate_calculation.hpp @@ -7,6 +7,7 @@ #include #include +#include #include #include #include @@ -136,12 +137,30 @@ Coordinate rotateCCWAroundZero(Coordinate coordinate, double angle_in_radians); Coordinate difference(const Coordinate lhs, const Coordinate rhs); // TEMPLATE/INLINE DEFINITIONS +// Helper: wrap a longitude into (-180, 180] range +inline double wrapLongitudeDouble(double lon) +{ + while (lon > 180.0) + lon -= 360.0; + while (lon <= -180.0) + lon += 360.0; + return lon; +} + inline std::pair projectPointOnSegment(const FloatCoordinate &source, const FloatCoordinate &target, const FloatCoordinate &coordinate) { - const FloatCoordinate slope_vector{target.lon - source.lon, target.lat - source.lat}; - const FloatCoordinate rel_coordinate{coordinate.lon - source.lon, coordinate.lat - source.lat}; + // Standard projection of a point on a segment in the coordinate space provided. + // Do not perform any antimeridian-specific wrapping here — callers should + // normalize longitudes if needed. This keeps behaviour consistent for + // projected coordinates (e.g., web_mercator) used throughout the codebase. + + const FloatCoordinate slope_vector{FloatLongitude{static_cast(target.lon) - static_cast(source.lon)}, + FloatLatitude{static_cast(target.lat) - static_cast(source.lat)}}; + const FloatCoordinate rel_coordinate{FloatLongitude{static_cast(coordinate.lon) - static_cast(source.lon)}, + FloatLatitude{static_cast(coordinate.lat) - static_cast(source.lat)}}; + // dot product of two un-normed vectors const auto unnormed_ratio = static_cast(slope_vector.lon * rel_coordinate.lon) + static_cast(slope_vector.lat * rel_coordinate.lat); @@ -165,13 +184,69 @@ inline std::pair projectPointOnSegment(const FloatCoord clamped_ratio = 0.; } + const double projected_lon = static_cast(source.lon) + clamped_ratio * + (static_cast(target.lon) - static_cast(source.lon)); + const double projected_lat = static_cast(source.lat) + clamped_ratio * + (static_cast(target.lat) - static_cast(source.lat)); + return {clamped_ratio, - { - FloatLongitude{1.0 - clamped_ratio} * source.lon + - target.lon * FloatLongitude{clamped_ratio}, - FloatLatitude{1.0 - clamped_ratio} * source.lat + - target.lat * FloatLatitude{clamped_ratio}, - }}; + {FloatLongitude{projected_lon}, FloatLatitude{projected_lat}}}; +} + +// Antimeridian-aware variant: unwrap longitudes so segments crossing the antimeridian are +// treated continuously. Callers that operate on raw WGS84 coordinates (or want explicit +// antimeridian handling) should use this function. +inline std::pair projectPointOnSegmentAntimeridian( + const FloatCoordinate &source, + const FloatCoordinate &target, + const FloatCoordinate &coordinate) +{ + // Unwrap the target longitude so the lon difference to source is minimal (handles antimeridian) + const double src_lon = static_cast(source.lon); + double tgt_lon = static_cast(target.lon); + const double coord_lon = static_cast(coordinate.lon); + + const double dlon = tgt_lon - src_lon; + if (dlon > 180.0) + tgt_lon -= 360.0; + else if (dlon < -180.0) + tgt_lon += 360.0; + + const FloatCoordinate slope_vector{FloatLongitude{tgt_lon - src_lon}, + FloatLatitude{static_cast(target.lat) - static_cast(source.lat)}}; + const FloatCoordinate rel_coordinate{FloatLongitude{coord_lon - src_lon}, + FloatLatitude{static_cast(coordinate.lat) - static_cast(source.lat)}}; + + const auto unnormed_ratio = static_cast(slope_vector.lon * rel_coordinate.lon) + + static_cast(slope_vector.lat * rel_coordinate.lat); + const auto squared_length = static_cast(slope_vector.lon * slope_vector.lon) + + static_cast(slope_vector.lat * slope_vector.lat); + + if (squared_length < std::numeric_limits::epsilon()) + { + return {0, source}; + } + + const double normed_ratio = unnormed_ratio / squared_length; + double clamped_ratio = normed_ratio; + if (clamped_ratio > 1.) + { + clamped_ratio = 1.; + } + else if (clamped_ratio < 0.) + { + clamped_ratio = 0.; + } + + // compute projected lon in unwrapped space, then wrap back into canonical range + const double projected_lon_unwrapped = (1.0 - clamped_ratio) * src_lon + clamped_ratio * tgt_lon; + const double projected_lon = wrapLongitudeDouble(projected_lon_unwrapped); + + const double projected_lat = static_cast(source.lat) + clamped_ratio * + (static_cast(target.lat) - static_cast(source.lat)); + + return {clamped_ratio, + {FloatLongitude{projected_lon}, FloatLatitude{projected_lat}}}; } template @@ -235,7 +310,7 @@ std::pair leastSquareRegression(const iterator_type begi // following the formulas of https://faculty.elgin.edu/dkernler/statistics/ch04/4-2.html const auto number_of_coordinates = std::distance(begin, end); BOOST_ASSERT(number_of_coordinates >= 2); - const auto extract_lon = [](const Coordinate coordinate) + std::function extract_lon = [](const Coordinate coordinate) { return static_cast(toFloating(coordinate.lon)); }; const auto extract_lat = [](const Coordinate coordinate) @@ -256,6 +331,37 @@ std::pair leastSquareRegression(const iterator_type begi min_lat = std::min(min_lat, lat); max_lat = std::max(max_lat, lat); } + // If coordinates span the antimeridian (lon span > 180°), normalize longitudes by adding 360° + // to negative longitudes so computations are continuous + const bool unwrap = (max_lon - min_lon) > 180.0; + if (unwrap) + { + auto orig_extract_lon = extract_lon; + extract_lon = [orig_extract_lon](const Coordinate coordinate) + { + double lon = orig_extract_lon(coordinate); + if (lon < 0.0) + lon += 360.0; + return lon; + }; + + // recompute min/max on unwrapped longitudes + min_lon = extract_lon(*begin); + max_lon = extract_lon(*begin); + min_lat = extract_lat(*begin); + max_lat = extract_lat(*begin); + for (auto coordinate_iterator = begin; coordinate_iterator != end; ++coordinate_iterator) + { + const auto c = *coordinate_iterator; + const auto lon = extract_lon(c); + min_lon = std::min(min_lon, lon); + max_lon = std::max(max_lon, lon); + const auto lat = extract_lat(c); + min_lat = std::min(min_lat, lat); + max_lat = std::max(max_lat, lat); + } + } + // very small difference in longitude -> would result in inaccurate calculation, check if lat is // better if ((max_lat - min_lat) > 2 * (max_lon - min_lon)) diff --git a/test/behave/environment.py b/test/behave/environment.py new file mode 100644 index 00000000000..93654add403 --- /dev/null +++ b/test/behave/environment.py @@ -0,0 +1,3 @@ +# behave environment placeholder: no-op +def before_all(context): + pass diff --git a/test/behave/features/antimeridian.feature b/test/behave/features/antimeridian.feature new file mode 100644 index 00000000000..ce5cf4a2107 --- /dev/null +++ b/test/behave/features/antimeridian.feature @@ -0,0 +1,9 @@ +Feature: Routing across the antimeridian + In order to ensure routes cross the dateline + As a user + I want routing to find continuous ways across longitude 180/-180 + + Scenario: Simple route crossing the antimeridian + Given an OSRM dataset with two nodes connected across the antimeridian + When I request a route from lon:179.9,lat:0 to lon:-179.9,lat:0 + Then the route should include a segment that crosses the antimeridian diff --git a/test/behave/features/steps/antimeridian_steps.py b/test/behave/features/steps/antimeridian_steps.py new file mode 100644 index 00000000000..bf630456966 --- /dev/null +++ b/test/behave/features/steps/antimeridian_steps.py @@ -0,0 +1,18 @@ +from behave import given, when, then + +@given('an OSRM dataset with two nodes connected across the antimeridian') +def step_impl(context): + # Minimal stub: ensure dataset exists or mark as TODO + context.dataset_prepared = True + +@when('I request a route from lon:179.9,lat:0 to lon:-179.9,lat:0') +def step_impl(context): + # This step should call the routing API; here it will stub the response to indicate failure + # Intentionally produce a failing/empty route to drive TDD + context.route = {'segments': []} + +@then('the route should include a segment that crosses the antimeridian') +def step_impl(context): + segments = context.route.get('segments', []) + # Expect at least one segment crossing; this will fail with current stub + assert any(seg.get('crosses_antimeridian', False) for seg in segments), "No antimeridian crossing segment found" diff --git a/unit_tests/util/antimeridian_regression_tests.cpp b/unit_tests/util/antimeridian_regression_tests.cpp new file mode 100644 index 00000000000..246ec051ebb --- /dev/null +++ b/unit_tests/util/antimeridian_regression_tests.cpp @@ -0,0 +1,58 @@ +#include + +#include "util/coordinate.hpp" +#include "util/coordinate_calculation.hpp" + +using namespace osrm::util; + +static double shortestAngularDistance(double a, double b) +{ + double diff = std::abs(a - b); + if (diff > 180.0) + diff = 360.0 - diff; + return diff; +} + +BOOST_AUTO_TEST_SUITE(antimeridian_regression_tests) + +BOOST_AUTO_TEST_CASE(least_square_regression_handles_antimeridian) +{ + std::vector coords{ + Coordinate{FloatLongitude{179.5}, FloatLatitude{0}}, + Coordinate{FloatLongitude{-179.0}, FloatLatitude{1}}, + Coordinate{FloatLongitude{179.9}, FloatLatitude{2}}, + }; + + const auto regression = + osrm::util::coordinate_calculation::leastSquareRegression(coords.begin(), coords.end()); + + const double lon1 = static_cast(toFloating(regression.first.lon)); + const double lon2 = static_cast(toFloating(regression.second.lon)); + + const double ang_diff = shortestAngularDistance(lon1, lon2); + + // Expect the regression to consider the short wraparound (distance small) + BOOST_CHECK_LT(ang_diff, 5.0); +} + +BOOST_AUTO_TEST_CASE(least_square_regression_regular_case) +{ + std::vector coords{ + Coordinate{FloatLongitude{10.0}, FloatLatitude{0}}, + Coordinate{FloatLongitude{20.0}, FloatLatitude{1}}, + Coordinate{FloatLongitude{30.0}, FloatLatitude{2}}, + }; + + const auto regression = + osrm::util::coordinate_calculation::leastSquareRegression(coords.begin(), coords.end()); + + const double lon1 = static_cast(toFloating(regression.first.lon)); + const double lon2 = static_cast(toFloating(regression.second.lon)); + + const double ang_diff = std::abs(lon1 - lon2); + + // Expect the regression line to span roughly the lon range (> 15 degrees) + BOOST_CHECK_GT(ang_diff, 15.0); +} + +BOOST_AUTO_TEST_SUITE_END() diff --git a/unit_tests/util/antimeridian_tests.cpp b/unit_tests/util/antimeridian_tests.cpp new file mode 100644 index 00000000000..e18405e8d50 --- /dev/null +++ b/unit_tests/util/antimeridian_tests.cpp @@ -0,0 +1,26 @@ +#include + +#include "util/coordinate.hpp" +#include "util/coordinate_calculation.hpp" + +using namespace osrm::util; + +BOOST_AUTO_TEST_SUITE(antimeridian_tests) + +BOOST_AUTO_TEST_CASE(project_point_on_segment_crossing_antimeridian) +{ + // segment from lon 179.5 to -179.5 (crosses antimeridian) + FloatCoordinate source{FloatLongitude{179.5}, FloatLatitude{0}}; + FloatCoordinate target{FloatLongitude{-179.5}, FloatLatitude{0}}; + FloatCoordinate query{FloatLongitude{180.0}, FloatLatitude{0}}; // midpoint over antimeridian + + auto [ratio, projected] = coordinate_calculation::projectPointOnSegmentAntimeridian(source, target, query); + + // Expect the projected point to be near lon 180 / -180 (wrapped) and ratio around 0.5 + BOOST_CHECK_CLOSE(ratio, 0.5, 1.0); // 1% tolerance + // Projected longitude should compare equal to 180 or -180 when converted to fixed coordinate + BOOST_CHECK((static_cast(projected.lon) > 179.0) || + (static_cast(projected.lon) < -179.0)); +} + +BOOST_AUTO_TEST_SUITE_END()