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

Fix possible division by zero by clamping latitude to 85.05°

Browse Source 
Resolves #3530
This commit is contained in:
Michael Krasnyk 2017-01-09 21:48:40 +01:00 committed by Patrick Niklaus
parent f82740ed93
commit 1d4d3b80b5
3 changed files with 29 additions and 31 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

37
include/util/web_mercator.hpp
View File

@ -20,7 +20,7 @@ const constexpr double EARTH_RADIUS_WGS84 = 6378137.0;
// earth circumference devided by 2
const constexpr double MAXEXTENT = EARTH_RADIUS_WGS84 * boost::math::constants::pi<double>();
// ^ math functions are not constexpr since they have side-effects (setting errno) :(
const constexpr double MAX_LATITUDE = 85.;
const constexpr double EPSG3857_MAX_LATITUDE = 85.051128779806592378; // 90(4*atan(exp(pi))/pi-1)
const constexpr double MAX_LONGITUDE = 180.0;
}
@ -29,6 +29,18 @@ const constexpr double DEGREE_TO_PX = detail::MAXEXTENT / 180.0;
// This is the global default tile size for all Mapbox Vector Tiles
const constexpr double TILE_SIZE = 256.0;
inline FloatLatitude clamp(const FloatLatitude lat)
{
return std::max(std::min(lat, FloatLatitude{detail::EPSG3857_MAX_LATITUDE}),
FloatLatitude{-detail::EPSG3857_MAX_LATITUDE});
}
inline FloatLongitude clamp(const FloatLongitude lon)
{
return std::max(std::min(lon, FloatLongitude{detail::MAX_LONGITUDE}),
FloatLongitude{-detail::MAX_LONGITUDE});
}
inline FloatLatitude yToLat(const double y)
{
const auto clamped_y = std::max(-180., std::min(180., y));
@ -41,10 +53,9 @@ inline FloatLatitude yToLat(const double y)
inline double latToY(const FloatLatitude latitude)
{
// apparently this is the (faster) version of the canonical log(tan()) version
const double f = std::sin(detail::DEGREE_TO_RAD * static_cast<double>(latitude));
const double y = detail::RAD_TO_DEGREE * 0.5 * std::log((1 + f) / (1 - f));
const auto clamped_y = std::max(-180., std::min(180., y));
return clamped_y;
const auto clamped_latitude = clamp(latitude);
const double f = std::sin(detail::DEGREE_TO_RAD * static_cast<double>(clamped_latitude));
return detail::RAD_TO_DEGREE * 0.5 * std::log((1 + f) / (1 - f));
}
template <typename T> constexpr double horner(double, T an) { return an; }
@ -91,18 +102,6 @@ inline double latToYapprox(const FloatLatitude latitude)
-3.23083224835967391884404730e-28);
}
inline FloatLatitude clamp(const FloatLatitude lat)
{
return std::max(std::min(lat, FloatLatitude{detail::MAX_LATITUDE}),
FloatLatitude{-detail::MAX_LATITUDE});
}
inline FloatLongitude clamp(const FloatLongitude lon)
{
return std::max(std::min(lon, FloatLongitude{detail::MAX_LONGITUDE}),
FloatLongitude{-detail::MAX_LONGITUDE});
}
inline void pixelToDegree(const double shift, double &x, double &y)
{
const double b = shift / 2.0;
@ -166,9 +165,9 @@ inline void xyzToMercator(
xyzToWGS84(x, y, z, minx, miny, maxx, maxy);
minx = static_cast<double>(clamp(util::FloatLongitude{minx})) * DEGREE_TO_PX;
miny = latToY(clamp(util::FloatLatitude{miny})) * DEGREE_TO_PX;
miny = latToY(util::FloatLatitude{miny}) * DEGREE_TO_PX;
maxx = static_cast<double>(clamp(util::FloatLongitude{maxx})) * DEGREE_TO_PX;
maxy = latToY(clamp(util::FloatLatitude{maxy})) * DEGREE_TO_PX;
maxy = latToY(util::FloatLatitude{maxy}) * DEGREE_TO_PX;
}
}
}

17
unit_tests/util/coordinate_calculation.cpp
View File

@ -118,8 +118,7 @@ BOOST_AUTO_TEST_CASE(compute_angle)
BOOST_CHECK_EQUAL(angle, 180);
// Tiny changes below our calculation resolution
// This should be equivalent to having two points on the same
// spot.
// This should be equivalent to having two points on the same spot.
first = Coordinate{FloatLongitude{0}, FloatLatitude{0}};
middle = Coordinate{FloatLongitude{1}, FloatLatitude{0}};
end = Coordinate{FloatLongitude{1 + std::numeric_limits<double>::epsilon()}, FloatLatitude{0}};
@ -127,14 +126,12 @@ BOOST_AUTO_TEST_CASE(compute_angle)
BOOST_CHECK_EQUAL(angle, 180);
// Invalid values
/* TODO: Enable this when I figure out how to use BOOST_CHECK_THROW
* and not have the whole test case fail...
first = Coordinate(FloatLongitude{0}, FloatLatitude{0});
middle = Coordinate(FloatLongitude{1}, FloatLatitude{0});
end = Coordinate(FloatLongitude(std::numeric_limits<double>::max()), FloatLatitude{0});
BOOST_CHECK_THROW( coordinate_calculation::computeAngle(first,middle,end),
boost::numeric::positive_overflow);
*/
BOOST_CHECK_THROW(
coordinate_calculation::computeAngle(
Coordinate(FloatLongitude{0}, FloatLatitude{0}),
Coordinate(FloatLongitude{1}, FloatLatitude{0}),
Coordinate(FloatLongitude{std::numeric_limits<double>::max()}, FloatLatitude{0})),
boost::numeric::positive_overflow);
}
// Regression test for bug captured in #1347

6
unit_tests/util/web_mercator.cpp
View File

@ -73,12 +73,14 @@ BOOST_AUTO_TEST_CASE(xyz_to_mercator)
double miny;
double maxx;
double maxy;
// http://tools.geofabrik.de/map/#13/85.0500/-175.5876&type=Geofabrik_Standard&grid=1
web_mercator::xyzToMercator(100, 0, 13, minx, miny, maxx, maxy);
BOOST_CHECK_CLOSE(minx, -19548311.361764118075, 0.0001);
BOOST_CHECK_CLOSE(miny, 19971868.8804085782, 0.0001);
BOOST_CHECK_CLOSE(miny, 20032616.372979045, 0.0001);
BOOST_CHECK_CLOSE(maxx, -19543419.391953866929, 0.0001);
BOOST_CHECK_CLOSE(maxy, 19971868.880408578, 0.0001);
BOOST_CHECK_CLOSE(maxy, 20037508.342789277, 0.0001); // Mercator 6378137*pi, WGS 85.0511
}
BOOST_AUTO_TEST_SUITE_END()