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

Use DouglasPeucker with squaredEuclideanDistance

Browse Source
This commit is contained in:
Patrick Niklaus 2016-04-09 00:39:20 +02:00
parent 67834def5f
commit 4886d46d91
3 changed files with 116 additions and 56 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

61
include/engine/douglas_peucker.hpp
View File

@ -12,26 +12,47 @@ namespace engine
{
namespace detail
{
const constexpr int DOUGLAS_PEUCKER_THRESHOLDS[19] = {
512440, // z0
256720, // z1
122560, // z2
56780, // z3
28800, // z4
14400, // z5
7200, // z6
3200, // z7
2400, // z8
1000, // z9
600, // z10
120, // z11
60, // z12
45, // z13
36, // z14
20, // z15
8, // z16
6, // z17
4, // z18
// This is derived from the following formular:
// x = b * (1 + lon/180) => dx = b * dlon/180
// y = b * (1 - lat/180) => dy = b * dlat/180
// dx^2 + dy^2 < min_pixel^2
// => dlon^2 + dlat^2 < min_pixel^2 / b^2 * 180^2
inline std::vector<std::uint64_t> generateThreshold(double min_pixel, unsigned number_of_zoomlevels)
{
std::vector<std::uint64_t> thresholds(number_of_zoomlevels);
for (unsigned zoom = 0; zoom < number_of_zoomlevels; ++zoom)
{
const double shift = (1u << zoom) * 256;
const double b = shift / 2.0;
const double pixel_to_deg = 180. / b;
const std::uint64_t min_deg = min_pixel * pixel_to_deg * COORDINATE_PRECISION;
thresholds[zoom] = min_deg * min_deg;
}
return thresholds;
}
const constexpr std::uint64_t DOUGLAS_PEUCKER_THRESHOLDS[19] = {
49438476562500, // z0
12359619140625, // z1
3089903027344, // z2
772475756836, // z3
193118939209, // z4
48279515076, // z5
12069878769, // z6
3017414761, // z7
754326225, // z8
188567824, // z9
47141956, // z10
11785489, // z11
2944656, // z12
736164, // z13
184041, // z14
45796, // z15
11449, // z16
2809, // z17
676, // z18
};
const constexpr auto DOUGLAS_PEUCKER_THRESHOLDS_SIZE =

31
src/engine/douglas_peucker.cpp
View File

@ -16,6 +16,31 @@ namespace osrm
namespace engine
{
struct FastPerpendicularDistance
{
FastPerpendicularDistance(const util::Coordinate start, const util::Coordinate target)
: projected_start(util::coordinate_calculation::mercator::fromWGS84(start)),
projected_target(util::coordinate_calculation::mercator::fromWGS84(target))
{
}
// Normed to the thresholds table
std::uint64_t operator()(const util::Coordinate coordinate) const
{
auto projected = util::coordinate_calculation::mercator::fromWGS84(coordinate);
util::FloatCoordinate projected_point_on_segment;
std::tie(std::ignore, projected_point_on_segment) =
util::coordinate_calculation::projectPointOnSegment(projected_start, projected_target,
projected);
auto squared_distance = util::coordinate_calculation::squaredEuclideanDistance(projected,
projected_point_on_segment);
return squared_distance;
}
const util::FloatCoordinate projected_start;
const util::FloatCoordinate projected_target;
};
std::vector<util::Coordinate> douglasPeucker(std::vector<util::Coordinate>::const_iterator begin,
std::vector<util::Coordinate>::const_iterator end,
const unsigned zoom_level)
@ -51,14 +76,16 @@ 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.first <= pair.second, "left border on the wrong side");
double max_distance = 0;
std::uint64_t max_distance = 0;
auto farthest_entry_index = pair.second;
FastPerpendicularDistance perpendicular_distance(begin[pair.first], begin[pair.second]);
// sweep over range to find the maximum
for (auto idx = pair.first + 1; idx != pair.second; ++idx)
{
using namespace util::coordinate_calculation;
const auto distance = perpendicularDistance(begin[pair.first], begin[pair.second], begin[idx]);
const auto distance = perpendicular_distance(begin[idx]);
// found new feasible maximum?
if (distance > max_distance &&
distance > detail::DOUGLAS_PEUCKER_THRESHOLDS[zoom_level])

72
unit_tests/engine/douglas_peucker.cpp
View File

@ -1,4 +1,5 @@
#include "engine/douglas_peucker.hpp"
#include "util/coordinate_calculation.hpp"
#include <boost/test/unit_test.hpp>
#include <boost/test/test_case_template.hpp>
@ -12,6 +13,16 @@ BOOST_AUTO_TEST_SUITE(douglas_peucker_simplification)
using namespace osrm;
using namespace osrm::engine;
BOOST_AUTO_TEST_CASE(calibrate_thresholds)
{
auto thresholds = detail::generateThreshold(5, 19);
auto z = 0;
for (auto t : thresholds)
{
BOOST_TEST_MESSAGE(t << ", // z" << z++);
}
}
BOOST_AUTO_TEST_CASE(removed_middle_test)
{
/*
@ -23,9 +34,9 @@ BOOST_AUTO_TEST_CASE(removed_middle_test)
*/
std::vector<util::Coordinate> coordinates = {
util::Coordinate(util::FloatLongitude(5), util::FloatLatitude(5)),
util::Coordinate(util::FloatLongitude(5.995715), util::FloatLatitude(6)),
util::Coordinate(util::FloatLongitude(10), util::FloatLatitude(10)),
util::Coordinate(util::FloatLongitude(15), util::FloatLatitude(5))};
util::Coordinate(util::FloatLongitude(12.5), util::FloatLatitude(12.6096298302)),
util::Coordinate(util::FloatLongitude(20), util::FloatLatitude(20)),
util::Coordinate(util::FloatLongitude(25), util::FloatLatitude(5))};
for (unsigned z = 0; z < detail::DOUGLAS_PEUCKER_THRESHOLDS_SIZE; z++)
{
@ -49,19 +60,19 @@ BOOST_AUTO_TEST_CASE(removed_middle_test_zoom_sensitive)
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))};
util::Coordinate(util::FloatLongitude(20), util::FloatLatitude(20)),
util::Coordinate(util::FloatLongitude(25), util::FloatLatitude(5))};
// Coordinate 6,6 should start getting included at Z12 and higher
// Coordinate 6,6 should start getting included at Z9 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++)
for (unsigned z = 0; z < 9; 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++)
// From 9 to max zoom, we should get all coordinates
for (unsigned z = 9; z < detail::DOUGLAS_PEUCKER_THRESHOLDS_SIZE; z++)
{
auto result = douglasPeucker(coordinates, z);
BOOST_CHECK_EQUAL(result.size(), 4);
@ -71,38 +82,39 @@ BOOST_AUTO_TEST_CASE(removed_middle_test_zoom_sensitive)
BOOST_AUTO_TEST_CASE(remove_second_node_test)
{
for (unsigned z = 0; z < detail::DOUGLAS_PEUCKER_THRESHOLDS_SIZE; z++)
// derived from the degreeToPixel function
const auto delta_pixel_to_delta_degree = [](const int pixel, const unsigned zoom) {
const double shift = (1u << zoom) * 256;
const double b = shift / 2.0;
return pixel * 180. / b;
};
// For zoom level 0 all gets reduced to a line
for (unsigned z = 1; z < detail::DOUGLAS_PEUCKER_THRESHOLDS_SIZE; z++)
{
/*
x--x
x
| \
x-x x
|
x
*/
std::vector<util::Coordinate> input = {
util::Coordinate(util::FixedLongitude(5 * COORDINATE_PRECISION),
util::FixedLatitude(5 * COORDINATE_PRECISION)),
util::Coordinate(util::FixedLongitude(5 * COORDINATE_PRECISION),
util::FixedLatitude(5 * COORDINATE_PRECISION +
detail::DOUGLAS_PEUCKER_THRESHOLDS[z])),
util::Coordinate(util::FixedLongitude(10 * COORDINATE_PRECISION),
util::FixedLatitude(10 * COORDINATE_PRECISION)),
util::Coordinate(util::FixedLongitude(10 * COORDINATE_PRECISION),
util::FixedLatitude(10 + COORDINATE_PRECISION +
detail::DOUGLAS_PEUCKER_THRESHOLDS[z] * 2)),
util::Coordinate(util::FixedLongitude(5 * COORDINATE_PRECISION),
util::FixedLatitude(15 * COORDINATE_PRECISION)),
util::Coordinate(util::FixedLongitude(5 * COORDINATE_PRECISION +
detail::DOUGLAS_PEUCKER_THRESHOLDS[z]),
util::FixedLatitude(15 * COORDINATE_PRECISION))};
BOOST_TEST_MESSAGE("Threshold (" << z << "): " << detail::DOUGLAS_PEUCKER_THRESHOLDS[z]);
util::Coordinate(util::FloatLongitude(5),
util::FloatLatitude(5)),
util::Coordinate(util::FloatLongitude(5 + delta_pixel_to_delta_degree(2, z)),
util::FloatLatitude(5)),
util::Coordinate(util::FloatLongitude(10),
util::FloatLatitude(10)),
util::Coordinate(util::FloatLongitude(5),
util::FloatLatitude(15)),
util::Coordinate(util::FloatLongitude(5),
util::FloatLatitude(15 + delta_pixel_to_delta_degree(2, z)))};
BOOST_TEST_MESSAGE("Delta (" << z << "): " << delta_pixel_to_delta_degree(2, z));
auto result = douglasPeucker(input, z);
BOOST_CHECK_EQUAL(result.size(), 4);
BOOST_CHECK_EQUAL(result.size(), 3);
BOOST_CHECK_EQUAL(input[0], result[0]);
BOOST_CHECK_EQUAL(input[2], result[1]);
BOOST_CHECK_EQUAL(input[3], result[2]);
BOOST_CHECK_EQUAL(input[5], result[3]);
BOOST_CHECK_EQUAL(input[4], result[2]);
}
}