Open-Harbor/osrm-backend
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 3 Wiki Activity

Use correct perpendicular distance measure when simplifying line geometries.

Browse Source
This commit is contained in:
Daniel Patterson 2016-03-28 13:29:59 -07:00 committed by Patrick Niklaus
parent 0a4015690f
commit a548155f65
2 changed files with 40 additions and 58 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

59
src/engine/douglas_peucker.cpp
View File

@ -16,53 +16,6 @@ namespace osrm
namespace engine 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<double>(toFloating(coordinate_a.lon)) * DEGREE_TO_RAD;
first_lat = static_cast<double>(toFloating(coordinate_a.lat)) * DEGREE_TO_RAD;
second_lon = static_cast<double>(toFloating(coordinate_b.lon)) * DEGREE_TO_RAD;
second_lat = static_cast<double>(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<double>(toFloating(other.lon)) * DEGREE_TO_RAD;
const float float_lat1 = static_cast<double>(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<int>(std::min(dist1, dist2));
}
float first_lat;
float first_lon;
float second_lat;
float second_lon;
};
}
std::vector<util::Coordinate> douglasPeucker(std::vector<util::Coordinate>::const_iterator begin, std::vector<util::Coordinate> douglasPeucker(std::vector<util::Coordinate>::const_iterator begin,
std::vector<util::Coordinate>::const_iterator end, std::vector<util::Coordinate>::const_iterator end,
const unsigned zoom_level) const unsigned zoom_level)
@ -98,25 +51,25 @@ std::vector<util::Coordinate> douglasPeucker(std::vector<util::Coordinate>::cons
BOOST_ASSERT_MSG(pair.second < size, "right border outside of geometry"); BOOST_ASSERT_MSG(pair.second < size, "right border outside of geometry");
BOOST_ASSERT_MSG(pair.first <= pair.second, "left border on the wrong side"); 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; auto farthest_entry_index = pair.second;
const CoordinatePairCalculator dist_calc(begin[pair.first], begin[pair.second]);
// sweep over range to find the maximum // sweep over range to find the maximum
for (auto idx = pair.first + 1; idx != pair.second; ++idx) 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? // found new feasible maximum?
if (distance > max_int_distance && if (distance > max_distance &&
distance > detail::DOUGLAS_PEUCKER_THRESHOLDS[zoom_level]) distance > detail::DOUGLAS_PEUCKER_THRESHOLDS[zoom_level])
{ {
farthest_entry_index = idx; farthest_entry_index = idx;
max_int_distance = distance; max_distance = distance;
} }
} }
// check if maximum violates a zoom level dependent threshold // 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 // mark idx as necessary
is_necessary[farthest_entry_index] = true; is_necessary[farthest_entry_index] = true;

39
unit_tests/engine/douglas_peucker.cpp
View File

@ -23,14 +23,11 @@ BOOST_AUTO_TEST_CASE(removed_middle_test)
*/ */
std::vector<util::Coordinate> coordinates = { std::vector<util::Coordinate> coordinates = {
util::Coordinate(util::FloatLongitude(5), util::FloatLatitude(5)), 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(10), util::FloatLatitude(10)),
util::Coordinate(util::FloatLongitude(15), util::FloatLatitude(5))}; util::Coordinate(util::FloatLongitude(15), util::FloatLatitude(5))};
// FIXME this test fails for everything below z4 because the DP algorithms for (unsigned z = 0; z < detail::DOUGLAS_PEUCKER_THRESHOLDS_SIZE; z++)
// only used a naive distance measurement
// for (unsigned z = 0; z < detail::DOUGLAS_PEUCKER_THRESHOLDS_SIZE; z++)
for (unsigned z = 0; z < 2; z++)
{ {
auto result = douglasPeucker(coordinates, z); auto result = douglasPeucker(coordinates, z);
BOOST_CHECK_EQUAL(result.size(), 3); 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<util::Coordinate> 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) BOOST_AUTO_TEST_CASE(remove_second_node_test)
{ {
for (unsigned z = 0; z < detail::DOUGLAS_PEUCKER_THRESHOLDS_SIZE; z++) for (unsigned z = 0; z < detail::DOUGLAS_PEUCKER_THRESHOLDS_SIZE; z++)