From a548155f65a45ac0bebbc99232d3c81f891e61f8 Mon Sep 17 00:00:00 2001 From: Daniel Patterson Date: Mon, 28 Mar 2016 13:29:59 -0700 Subject: [PATCH] Use correct perpendicular distance measure when simplifying line geometries. --- src/engine/douglas_peucker.cpp | 59 +++------------------------ unit_tests/engine/douglas_peucker.cpp | 39 +++++++++++++++--- 2 files changed, 40 insertions(+), 58 deletions(-) diff --git a/src/engine/douglas_peucker.cpp b/src/engine/douglas_peucker.cpp index 3333768cb..e7e2f9024 100644 --- a/src/engine/douglas_peucker.cpp +++ b/src/engine/douglas_peucker.cpp @@ -16,53 +16,6 @@ namespace osrm namespace engine { -namespace -{ -// FIXME This algorithm is a very naive approximation that leads to -// problems like (almost) co-linear points not being simplified. -// Switch to real-point-segment distance of projected coordinates -struct CoordinatePairCalculator -{ - CoordinatePairCalculator() = delete; - CoordinatePairCalculator(const util::Coordinate coordinate_a, - const util::Coordinate coordinate_b) - { - using namespace util::coordinate_calculation; - // initialize distance calculator with two fixed coordinates a, b - first_lon = static_cast(toFloating(coordinate_a.lon)) * DEGREE_TO_RAD; - first_lat = static_cast(toFloating(coordinate_a.lat)) * DEGREE_TO_RAD; - second_lon = static_cast(toFloating(coordinate_b.lon)) * DEGREE_TO_RAD; - second_lat = static_cast(toFloating(coordinate_b.lat)) * DEGREE_TO_RAD; - } - - int operator()(const util::Coordinate other) const - { - using namespace util::coordinate_calculation; - // set third coordinate c - const float float_lon1 = static_cast(toFloating(other.lon)) * DEGREE_TO_RAD; - const float float_lat1 = static_cast(toFloating(other.lat)) * DEGREE_TO_RAD; - - // compute distance (a,c) - const float x_value_1 = (first_lon - float_lon1) * cos((float_lat1 + first_lat) / 2.f); - const float y_value_1 = first_lat - float_lat1; - const float dist1 = std::hypot(x_value_1, y_value_1) * EARTH_RADIUS; - - // compute distance (b,c) - const float x_value_2 = (second_lon - float_lon1) * cos((float_lat1 + second_lat) / 2.f); - const float y_value_2 = second_lat - float_lat1; - const float dist2 = std::hypot(x_value_2, y_value_2) * EARTH_RADIUS; - - // return the minimum - return static_cast(std::min(dist1, dist2)); - } - - float first_lat; - float first_lon; - float second_lat; - float second_lon; -}; -} - std::vector douglasPeucker(std::vector::const_iterator begin, std::vector::const_iterator end, const unsigned zoom_level) @@ -98,25 +51,25 @@ std::vector douglasPeucker(std::vector::cons BOOST_ASSERT_MSG(pair.second < size, "right border outside of geometry"); BOOST_ASSERT_MSG(pair.first <= pair.second, "left border on the wrong side"); - int max_int_distance = 0; + double max_distance = 0; auto farthest_entry_index = pair.second; - const CoordinatePairCalculator dist_calc(begin[pair.first], begin[pair.second]); // sweep over range to find the maximum for (auto idx = pair.first + 1; idx != pair.second; ++idx) { - const int distance = dist_calc(begin[idx]); + using namespace util::coordinate_calculation; + const auto distance = perpendicularDistance(begin[pair.first], begin[pair.second], begin[idx]); // found new feasible maximum? - if (distance > max_int_distance && + if (distance > max_distance && distance > detail::DOUGLAS_PEUCKER_THRESHOLDS[zoom_level]) { farthest_entry_index = idx; - max_int_distance = distance; + max_distance = distance; } } // check if maximum violates a zoom level dependent threshold - if (max_int_distance > detail::DOUGLAS_PEUCKER_THRESHOLDS[zoom_level]) + if (max_distance > detail::DOUGLAS_PEUCKER_THRESHOLDS[zoom_level]) { // mark idx as necessary is_necessary[farthest_entry_index] = true; diff --git a/unit_tests/engine/douglas_peucker.cpp b/unit_tests/engine/douglas_peucker.cpp index a1290eba5..2b4c0cf23 100644 --- a/unit_tests/engine/douglas_peucker.cpp +++ b/unit_tests/engine/douglas_peucker.cpp @@ -23,14 +23,11 @@ BOOST_AUTO_TEST_CASE(removed_middle_test) */ std::vector coordinates = { util::Coordinate(util::FloatLongitude(5), util::FloatLatitude(5)), - util::Coordinate(util::FloatLongitude(6), util::FloatLatitude(6)), + util::Coordinate(util::FloatLongitude(7.4999999999999725), util::FloatLatitude(7.50718628974349)), util::Coordinate(util::FloatLongitude(10), util::FloatLatitude(10)), util::Coordinate(util::FloatLongitude(15), util::FloatLatitude(5))}; - // FIXME this test fails for everything below z4 because the DP algorithms - // only used a naive distance measurement - // for (unsigned z = 0; z < detail::DOUGLAS_PEUCKER_THRESHOLDS_SIZE; z++) - for (unsigned z = 0; z < 2; z++) + for (unsigned z = 0; z < detail::DOUGLAS_PEUCKER_THRESHOLDS_SIZE; z++) { auto result = douglasPeucker(coordinates, z); BOOST_CHECK_EQUAL(result.size(), 3); @@ -40,6 +37,38 @@ BOOST_AUTO_TEST_CASE(removed_middle_test) } } +BOOST_AUTO_TEST_CASE(removed_middle_test_zoom_sensitive) +{ + /* + x + / \ + x \ + / \ + x x + */ + std::vector coordinates = { + util::Coordinate(util::FloatLongitude(5), util::FloatLatitude(5)), + util::Coordinate(util::FloatLongitude(6), util::FloatLatitude(6)), + util::Coordinate(util::FloatLongitude(10), util::FloatLatitude(10)), + util::Coordinate(util::FloatLongitude(15), util::FloatLatitude(5))}; + + // Coordinate 6,6 should start getting included at Z12 and higher + // Note that 5,5->6,6->10,10 is *not* a straight line on the surface + // of the earth + for (unsigned z = 0; z < 11; z++) + { + auto result = douglasPeucker(coordinates, z); + BOOST_CHECK_EQUAL(result.size(), 3); + } + // From 12 to max zoom, we should get all coordinates + for (unsigned z = 12; z < detail::DOUGLAS_PEUCKER_THRESHOLDS_SIZE; z++) + { + auto result = douglasPeucker(coordinates, z); + BOOST_CHECK_EQUAL(result.size(), 4); + } +} + + BOOST_AUTO_TEST_CASE(remove_second_node_test) { for (unsigned z = 0; z < detail::DOUGLAS_PEUCKER_THRESHOLDS_SIZE; z++)