From 5052c4ae3a789e68f8d6a72bfd9ede833ed36d12 Mon Sep 17 00:00:00 2001 From: Patrick Niklaus Date: Sat, 9 Apr 2016 02:18:47 +0200 Subject: [PATCH] Move projection function into own header and inline --- include/engine/geospatial_query.hpp | 8 +- include/util/coordinate_calculation.hpp | 35 +----- include/util/static_rtree.hpp | 73 +++++------ include/util/viewport.hpp | 20 +-- include/util/web_mercator.hpp | 127 +++++++++++++++++++ src/engine/douglas_peucker.cpp | 7 +- src/engine/plugins/tile.cpp | 24 ++-- src/util/coordinate_calculation.cpp | 134 +++------------------ unit_tests/util/coordinate_calculation.cpp | 66 ---------- unit_tests/util/static_rtree.cpp | 11 +- unit_tests/util/web_mercator.cpp | 74 ++++++++++++ 11 files changed, 293 insertions(+), 286 deletions(-) create mode 100644 include/util/web_mercator.hpp create mode 100644 unit_tests/util/web_mercator.cpp diff --git a/include/engine/geospatial_query.hpp b/include/engine/geospatial_query.hpp index e880139f3..d2e90b887 100644 --- a/include/engine/geospatial_query.hpp +++ b/include/engine/geospatial_query.hpp @@ -6,6 +6,7 @@ #include "engine/phantom_node.hpp" #include "util/bearing.hpp" #include "util/rectangle.hpp" +#include "util/web_mercator.hpp" #include "osrm/coordinate.hpp" @@ -424,10 +425,11 @@ template class GeospatialQuery BOOST_ASSERT(segment.data.reverse_segment_id.id != SPECIAL_SEGMENTID || !segment.data.reverse_segment_id.enabled); - Coordinate wsg84_coordinate = util::coordinate_calculation::mercator::toWGS84( - segment.fixed_projected_coordinate); + Coordinate wsg84_coordinate = + util::web_mercator::toWGS84(segment.fixed_projected_coordinate); - return util::coordinate_calculation::haversineDistance(input_coordinate, wsg84_coordinate) > max_distance; + return util::coordinate_calculation::haversineDistance(input_coordinate, wsg84_coordinate) > + max_distance; } std::pair checkSegmentBearing(const CandidateSegment &segment, diff --git a/include/util/coordinate_calculation.hpp b/include/util/coordinate_calculation.hpp index 14480d778..afa6c6532 100644 --- a/include/util/coordinate_calculation.hpp +++ b/include/util/coordinate_calculation.hpp @@ -3,7 +3,6 @@ #include "util/coordinate.hpp" -#include #include #include @@ -15,24 +14,16 @@ namespace util namespace coordinate_calculation { +namespace detail +{ const constexpr long double DEGREE_TO_RAD = 0.017453292519943295769236907684886; const constexpr long double RAD_TO_DEGREE = 1. / DEGREE_TO_RAD; // earth radius varies between 6,356.750-6,378.135 km (3,949.901-3,963.189mi) // The IUGG value for the equatorial radius is 6378.137 km (3963.19 miles) const constexpr long double EARTH_RADIUS = 6372797.560856; -// radius used by WGS84 -const constexpr double EARTH_RADIUS_WGS84 = 6378137.0; - -namespace detail -{ -// earth circumference devided by 2 -const constexpr double MAXEXTENT = EARTH_RADIUS_WGS84 * boost::math::constants::pi(); -// ^ math functions are not constexpr since they have side-effects (setting errno) :( -const double MAX_LATITUDE = RAD_TO_DEGREE * (2.0 * std::atan(std::exp(180.0 * DEGREE_TO_RAD)) - - boost::math::constants::half_pi()); -const constexpr double MAX_LONGITUDE = 180.0; } + //! Takes the squared euclidean distance of the input coordinates. Does not return meters! std::uint64_t squaredEuclideanDistance(const Coordinate &lhs, const Coordinate &rhs); @@ -76,26 +67,6 @@ double circleRadius(const Coordinate first_coordinate, // returns to Coordinate interpolateLinear(double factor, const Coordinate from, const Coordinate to); -namespace mercator -{ -// This is the global default tile size for all Mapbox Vector Tiles -const constexpr double TILE_SIZE = 256.0; -// Converts projected mercator degrees to PX -const constexpr double DEGREE_TO_PX = detail::MAXEXTENT / 180.0; - -double degreeToPixel(FloatLatitude lat, unsigned zoom); -double degreeToPixel(FloatLongitude lon, unsigned zoom); -FloatLatitude yToLat(const double value); -double latToY(const FloatLatitude latitude); - -FloatCoordinate fromWGS84(const FloatCoordinate &wgs84_coordinate); -FloatCoordinate toWGS84(const FloatCoordinate &mercator_coordinate); - -void xyzToMercator( - const int x, const int y, const int z, double &minx, double &miny, double &maxx, double &maxy); -void xyzToWGS84( - const int x, const int y, const int z, double &minx, double &miny, double &maxx, double &maxy); -} // ns mercator } // ns coordinate_calculation } // ns util } // ns osrm diff --git a/include/util/static_rtree.hpp b/include/util/static_rtree.hpp index a87e103c9..32bd60bd0 100644 --- a/include/util/static_rtree.hpp +++ b/include/util/static_rtree.hpp @@ -5,11 +5,12 @@ #include "util/hilbert_value.hpp" #include "util/rectangle.hpp" #include "util/shared_memory_vector_wrapper.hpp" - #include "util/bearing.hpp" #include "util/exception.hpp" #include "util/integer_range.hpp" #include "util/typedefs.hpp" +#include "util/coordinate_calculation.hpp" +#include "util/web_mercator.hpp" #include "osrm/coordinate.hpp" @@ -129,7 +130,8 @@ class StaticRTree tbb::parallel_for( tbb::blocked_range(0, m_element_count), [&input_data_vector, &input_wrapper_vector, - &coordinate_list](const tbb::blocked_range &range) { + &coordinate_list](const tbb::blocked_range &range) + { for (uint64_t element_counter = range.begin(), end = range.end(); element_counter != end; ++element_counter) { @@ -144,9 +146,9 @@ class StaticRTree Coordinate current_centroid = coordinate_calculation::centroid( coordinate_list[current_element.u], coordinate_list[current_element.v]); - current_centroid.lat = FixedLatitude( - COORDINATE_PRECISION * - coordinate_calculation::mercator::latToY(toFloating(current_centroid.lat))); + current_centroid.lat = + FixedLatitude(COORDINATE_PRECISION * + web_mercator::latToY(toFloating(current_centroid.lat))); current_wrapper.m_hilbert_value = hilbertCode(current_centroid); } @@ -233,20 +235,21 @@ class StaticRTree std::reverse(m_search_tree.begin(), m_search_tree.end()); std::uint32_t search_tree_size = m_search_tree.size(); - tbb::parallel_for( - tbb::blocked_range(0, search_tree_size), - [this, &search_tree_size](const tbb::blocked_range &range) { - for (std::uint32_t i = range.begin(), end = range.end(); i != end; ++i) - { - TreeNode ¤t_tree_node = this->m_search_tree[i]; - for (std::uint32_t j = 0; j < current_tree_node.child_count; ++j) - { - const std::uint32_t old_id = current_tree_node.children[j]; - const std::uint32_t new_id = search_tree_size - old_id - 1; - current_tree_node.children[j] = new_id; - } - } - }); + tbb::parallel_for(tbb::blocked_range(0, search_tree_size), + [this, &search_tree_size](const tbb::blocked_range &range) + { + for (std::uint32_t i = range.begin(), end = range.end(); i != end; + ++i) + { + TreeNode ¤t_tree_node = this->m_search_tree[i]; + for (std::uint32_t j = 0; j < current_tree_node.child_count; ++j) + { + const std::uint32_t old_id = current_tree_node.children[j]; + const std::uint32_t new_id = search_tree_size - old_id - 1; + current_tree_node.children[j] = new_id; + } + } + }); // open tree file boost::filesystem::ofstream tree_node_file(tree_node_filename, std::ios::binary); @@ -325,10 +328,10 @@ class StaticRTree { const Rectangle projected_rectangle{ search_rectangle.min_lon, search_rectangle.max_lon, - toFixed(FloatLatitude{coordinate_calculation::mercator::latToY( - toFloating(FixedLatitude(search_rectangle.min_lat)))}), - toFixed(FloatLatitude{coordinate_calculation::mercator::latToY( - toFloating(FixedLatitude(search_rectangle.max_lat)))})}; + toFixed(FloatLatitude{ + web_mercator::latToY(toFloating(FixedLatitude(search_rectangle.min_lat)))}), + toFixed(FloatLatitude{ + web_mercator::latToY(toFloating(FixedLatitude(search_rectangle.max_lat)))})}; std::vector results; std::queue traversal_queue; @@ -391,8 +394,12 @@ class StaticRTree std::vector Nearest(const Coordinate input_coordinate, const std::size_t max_results) { return Nearest(input_coordinate, - [](const CandidateSegment &) { return std::make_pair(true, true); }, - [max_results](const std::size_t num_results, const CandidateSegment &) { + [](const CandidateSegment &) + { + return std::make_pair(true, true); + }, + [max_results](const std::size_t num_results, const CandidateSegment &) + { return num_results >= max_results; }); } @@ -403,7 +410,7 @@ class StaticRTree Nearest(const Coordinate input_coordinate, const FilterT filter, const TerminationT terminate) { std::vector results; - auto projected_coordinate = coordinate_calculation::mercator::fromWGS84(input_coordinate); + auto projected_coordinate = web_mercator::fromWGS84(input_coordinate); Coordinate fixed_projected_coordinate{projected_coordinate}; // initialize queue with root element @@ -470,10 +477,8 @@ class StaticRTree for (const auto i : irange(0u, current_leaf_node.object_count)) { const auto ¤t_edge = current_leaf_node.objects[i]; - const auto projected_u = - coordinate_calculation::mercator::fromWGS84((*m_coordinate_list)[current_edge.u]); - const auto projected_v = - coordinate_calculation::mercator::fromWGS84((*m_coordinate_list)[current_edge.v]); + const auto projected_u = web_mercator::fromWGS84((*m_coordinate_list)[current_edge.u]); + const auto projected_v = web_mercator::fromWGS84((*m_coordinate_list)[current_edge.v]); FloatCoordinate projected_nearest; std::tie(std::ignore, projected_nearest) = @@ -535,10 +540,10 @@ class StaticRTree BOOST_ASSERT(objects[i].u < coordinate_list.size()); BOOST_ASSERT(objects[i].v < coordinate_list.size()); - Coordinate projected_u{coordinate_calculation::mercator::fromWGS84( - Coordinate{coordinate_list[objects[i].u]})}; - Coordinate projected_v{coordinate_calculation::mercator::fromWGS84( - Coordinate{coordinate_list[objects[i].v]})}; + Coordinate projected_u{ + web_mercator::fromWGS84(Coordinate{coordinate_list[objects[i].u]})}; + Coordinate projected_v{ + web_mercator::fromWGS84(Coordinate{coordinate_list[objects[i].v]})}; BOOST_ASSERT(toFloating(projected_u.lon) <= FloatLongitude(180.)); BOOST_ASSERT(toFloating(projected_u.lon) >= FloatLongitude(-180.)); diff --git a/include/util/viewport.hpp b/include/util/viewport.hpp index fd0ac73f6..d4ed889e0 100644 --- a/include/util/viewport.hpp +++ b/include/util/viewport.hpp @@ -2,7 +2,7 @@ #define UTIL_VIEWPORT_HPP #include "util/coordinate.hpp" -#include "util/coordinate_calculation.hpp" +#include "util/web_mercator.hpp" #include @@ -23,27 +23,27 @@ namespace detail static constexpr unsigned MAX_ZOOM = 18; static constexpr unsigned MIN_ZOOM = 1; // this is an upper bound to current display sizes -static constexpr double VIEWPORT_WIDTH = 8 * coordinate_calculation::mercator::TILE_SIZE; -static constexpr double VIEWPORT_HEIGHT = 5 * coordinate_calculation::mercator::TILE_SIZE; +static constexpr double VIEWPORT_WIDTH = 8 * web_mercator::TILE_SIZE; +static constexpr double VIEWPORT_HEIGHT = 5 * web_mercator::TILE_SIZE; static double INV_LOG_2 = 1. / std::log(2); } unsigned getFittedZoom(util::Coordinate south_west, util::Coordinate north_east) { - const auto min_x = coordinate_calculation::mercator::degreeToPixel(toFloating(south_west.lon), detail::MAX_ZOOM); - const auto max_y = coordinate_calculation::mercator::degreeToPixel(toFloating(south_west.lat), detail::MAX_ZOOM); - const auto max_x = coordinate_calculation::mercator::degreeToPixel(toFloating(north_east.lon), detail::MAX_ZOOM); - const auto min_y = coordinate_calculation::mercator::degreeToPixel(toFloating(north_east.lat), detail::MAX_ZOOM); + const auto min_x = web_mercator::degreeToPixel(toFloating(south_west.lon), detail::MAX_ZOOM); + const auto max_y = web_mercator::degreeToPixel(toFloating(south_west.lat), detail::MAX_ZOOM); + const auto max_x = web_mercator::degreeToPixel(toFloating(north_east.lon), detail::MAX_ZOOM); + const auto min_y = web_mercator::degreeToPixel(toFloating(north_east.lat), detail::MAX_ZOOM); const double width_ratio = (max_x - min_x) / detail::VIEWPORT_WIDTH; const double height_ratio = (max_y - min_y) / detail::VIEWPORT_HEIGHT; - const auto zoom = detail::MAX_ZOOM - std::max(std::log(width_ratio), std::log(height_ratio)) * detail::INV_LOG_2; + const auto zoom = detail::MAX_ZOOM - + std::max(std::log(width_ratio), std::log(height_ratio)) * detail::INV_LOG_2; if (std::isfinite(zoom)) return std::max(detail::MIN_ZOOM, zoom); else - return detail::MIN_ZOOM; + return detail::MIN_ZOOM; } - } } } diff --git a/include/util/web_mercator.hpp b/include/util/web_mercator.hpp new file mode 100644 index 000000000..8c23de143 --- /dev/null +++ b/include/util/web_mercator.hpp @@ -0,0 +1,127 @@ +#ifndef OSRM_WEB_MERCATOR_HPP +#define OSRM_WEB_MERCATOR_HPP + +#include "util/coordinate.hpp" + +#include + +namespace osrm +{ +namespace util +{ +namespace web_mercator +{ +namespace detail +{ +const constexpr long double DEGREE_TO_RAD = 0.017453292519943295769236907684886; +const constexpr long double RAD_TO_DEGREE = 1. / DEGREE_TO_RAD; +// radius used by WGS84 +const constexpr double EARTH_RADIUS_WGS84 = 6378137.0; +// earth circumference devided by 2 +const constexpr double MAXEXTENT = EARTH_RADIUS_WGS84 * boost::math::constants::pi(); +// ^ math functions are not constexpr since they have side-effects (setting errno) :( +const constexpr double MAX_LATITUDE = 85.; +const constexpr double MAX_LONGITUDE = 180.0; +} + +// Converts projected mercator degrees to PX +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 yToLat(const double y) +{ + const auto clamped_y = std::max(-180., std::min(180., y)); + const double normalized_lat = + detail::RAD_TO_DEGREE * 2. * std::atan(std::exp(clamped_y * detail::DEGREE_TO_RAD)); + + return FloatLatitude(normalized_lat - 90.); +} + +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(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; +} + +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; + x = (x - b) / shift * 360.0; + // FIXME needs to be simplified + const double g = (y - b) / -(shift / (2 * M_PI)) / detail::DEGREE_TO_RAD; + static_assert(detail::DEGREE_TO_RAD / (2 * M_PI) - 1 / 360. < 0.0001, ""); + y = static_cast(yToLat(g)); +} + +inline double degreeToPixel(FloatLongitude lon, unsigned zoom) +{ + const double shift = (1u << zoom) * TILE_SIZE; + const double b = shift / 2.0; + const double x = b * (1 + static_cast(lon) / 180.0); + return x; +} + +inline double degreeToPixel(FloatLatitude lat, unsigned zoom) +{ + const double shift = (1u << zoom) * TILE_SIZE; + const double b = shift / 2.0; + const double y = b * (1. - latToY(lat) / 180.); + return y; +} + +inline FloatCoordinate fromWGS84(const FloatCoordinate &wgs84_coordinate) +{ + return {wgs84_coordinate.lon, FloatLatitude{latToY(wgs84_coordinate.lat)}}; +} + +inline FloatCoordinate toWGS84(const FloatCoordinate &mercator_coordinate) +{ + return {mercator_coordinate.lon, yToLat(static_cast(mercator_coordinate.lat))}; +} + +// Converts a WMS tile coordinate (z,x,y) into a wgs bounding box +inline void xyzToWGS84( + const int x, const int y, const int z, double &minx, double &miny, double &maxx, double &maxy) +{ + minx = x * TILE_SIZE; + miny = (y + 1.0) * TILE_SIZE; + maxx = (x + 1.0) * TILE_SIZE; + maxy = y * TILE_SIZE; + // 2^z * TILE_SIZE + const double shift = (1u << static_cast(z)) * TILE_SIZE; + pixelToDegree(shift, minx, miny); + pixelToDegree(shift, maxx, maxy); +} + +// Converts a WMS tile coordinate (z,x,y) into a mercator bounding box +inline void xyzToMercator( + const int x, const int y, const int z, double &minx, double &miny, double &maxx, double &maxy) +{ + xyzToWGS84(x, y, z, minx, miny, maxx, maxy); + + minx = static_cast(clamp(util::FloatLongitude(minx))) * DEGREE_TO_PX; + miny = latToY(clamp(util::FloatLatitude(miny))) * DEGREE_TO_PX; + maxx = static_cast(clamp(util::FloatLongitude(maxx))) * DEGREE_TO_PX; + maxy = latToY(clamp(util::FloatLatitude(maxy))) * DEGREE_TO_PX; +} +} +} +} + +#endif diff --git a/src/engine/douglas_peucker.cpp b/src/engine/douglas_peucker.cpp index 5f10e1db7..24df636b3 100644 --- a/src/engine/douglas_peucker.cpp +++ b/src/engine/douglas_peucker.cpp @@ -1,5 +1,6 @@ #include "engine/douglas_peucker.hpp" #include "util/coordinate_calculation.hpp" +#include "util/web_mercator.hpp" #include "util/coordinate.hpp" #include @@ -19,15 +20,15 @@ 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)) + : projected_start(util::web_mercator::fromWGS84(start)), + projected_target(util::web_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); + auto projected = util::web_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, diff --git a/src/engine/plugins/tile.cpp b/src/engine/plugins/tile.cpp index 7d18d9029..87586ea47 100644 --- a/src/engine/plugins/tile.cpp +++ b/src/engine/plugins/tile.cpp @@ -2,6 +2,7 @@ #include "engine/plugins/tile.hpp" #include "util/coordinate_calculation.hpp" +#include "util/web_mercator.hpp" #include #include @@ -119,7 +120,6 @@ FixedLine coordinatesToTileLine(const util::Coordinate start, const util::Coordinate target, const detail::BBox &tile_bbox) { - using namespace util::coordinate_calculation; FloatLine geo_line; geo_line.emplace_back(static_cast(util::toFloating(start.lon)), static_cast(util::toFloating(start.lat))); @@ -130,15 +130,16 @@ FixedLine coordinatesToTileLine(const util::Coordinate start, for (auto const &pt : geo_line) { - double px_merc = pt.x * mercator::DEGREE_TO_PX; - double py_merc = mercator::latToY(util::FloatLatitude(pt.y)) * mercator::DEGREE_TO_PX; + double px_merc = pt.x * util::web_mercator::DEGREE_TO_PX; + double py_merc = util::web_mercator::latToY(util::FloatLatitude(pt.y)) * + util::web_mercator::DEGREE_TO_PX; // convert lon/lat to tile coordinates const auto px = std::round( - ((px_merc - tile_bbox.minx) * mercator::TILE_SIZE / tile_bbox.width()) * - detail::VECTOR_TILE_EXTENT / util::coordinate_calculation::mercator::TILE_SIZE); + ((px_merc - tile_bbox.minx) * util::web_mercator::TILE_SIZE / tile_bbox.width()) * + detail::VECTOR_TILE_EXTENT / util::web_mercator::TILE_SIZE); const auto py = std::round( - ((tile_bbox.maxy - py_merc) * mercator::TILE_SIZE / tile_bbox.height()) * - detail::VECTOR_TILE_EXTENT / util::coordinate_calculation::mercator::TILE_SIZE); + ((tile_bbox.maxy - py_merc) * util::web_mercator::TILE_SIZE / tile_bbox.height()) * + detail::VECTOR_TILE_EXTENT / util::web_mercator::TILE_SIZE); boost::geometry::append(unclipped_line, point_t(px, py)); } @@ -170,12 +171,11 @@ Status TilePlugin::HandleRequest(const api::TileParameters ¶meters, std::str { BOOST_ASSERT(parameters.IsValid()); - using namespace util::coordinate_calculation; double min_lon, min_lat, max_lon, max_lat; // Convert the z,x,y mercator tile coordinates into WGS84 lon/lat values - mercator::xyzToWGS84(parameters.x, parameters.y, parameters.z, min_lon, min_lat, max_lon, - max_lat); + util::web_mercator::xyzToWGS84(parameters.x, parameters.y, parameters.z, min_lon, min_lat, + max_lon, max_lat); util::Coordinate southwest{util::FloatLongitude(min_lon), util::FloatLatitude(min_lat)}; util::Coordinate northeast{util::FloatLongitude(max_lon), util::FloatLatitude(max_lat)}; @@ -245,8 +245,8 @@ Status TilePlugin::HandleRequest(const api::TileParameters ¶meters, std::str // TODO: extract speed values for compressed and uncompressed geometries // Convert tile coordinates into mercator coordinates - mercator::xyzToMercator(parameters.x, parameters.y, parameters.z, min_lon, min_lat, max_lon, - max_lat); + util::web_mercator::xyzToMercator(parameters.x, parameters.y, parameters.z, min_lon, min_lat, + max_lon, max_lat); const detail::BBox tile_bbox{min_lon, min_lat, max_lon, max_lat}; // Protobuf serialized blocks when objects go out of scope, hence diff --git a/src/util/coordinate_calculation.cpp b/src/util/coordinate_calculation.cpp index 766a03684..8af5ac2a4 100644 --- a/src/util/coordinate_calculation.cpp +++ b/src/util/coordinate_calculation.cpp @@ -1,6 +1,7 @@ #include "util/coordinate.hpp" #include "util/coordinate_calculation.hpp" #include "util/trigonometry_table.hpp" +#include "util/web_mercator.hpp" #include @@ -40,11 +41,11 @@ double haversineDistance(const Coordinate coordinate_1, const Coordinate coordin const double ln1 = lon1 / COORDINATE_PRECISION; const double lt2 = lat2 / COORDINATE_PRECISION; const double ln2 = lon2 / COORDINATE_PRECISION; - const double dlat1 = lt1 * DEGREE_TO_RAD; + const double dlat1 = lt1 * detail::DEGREE_TO_RAD; - const double dlong1 = ln1 * DEGREE_TO_RAD; - const double dlat2 = lt2 * DEGREE_TO_RAD; - const double dlong2 = ln2 * DEGREE_TO_RAD; + const double dlong1 = ln1 * detail::DEGREE_TO_RAD; + const double dlat2 = lt2 * detail::DEGREE_TO_RAD; + const double dlong2 = ln2 * detail::DEGREE_TO_RAD; const double dlong = dlong1 - dlong2; const double dlat = dlat1 - dlat2; @@ -52,7 +53,7 @@ double haversineDistance(const Coordinate coordinate_1, const Coordinate coordin const double aharv = std::pow(std::sin(dlat / 2.0), 2.0) + std::cos(dlat1) * std::cos(dlat2) * std::pow(std::sin(dlong / 2.), 2); const double charv = 2. * std::atan2(std::sqrt(aharv), std::sqrt(1.0 - aharv)); - return EARTH_RADIUS * charv; + return detail::EARTH_RADIUS * charv; } double greatCircleDistance(const Coordinate coordinate_1, const Coordinate coordinate_2) @@ -66,14 +67,14 @@ double greatCircleDistance(const Coordinate coordinate_1, const Coordinate coord BOOST_ASSERT(lat2 != std::numeric_limits::min()); BOOST_ASSERT(lon2 != std::numeric_limits::min()); - const double float_lat1 = (lat1 / COORDINATE_PRECISION) * DEGREE_TO_RAD; - const double float_lon1 = (lon1 / COORDINATE_PRECISION) * DEGREE_TO_RAD; - const double float_lat2 = (lat2 / COORDINATE_PRECISION) * DEGREE_TO_RAD; - const double float_lon2 = (lon2 / COORDINATE_PRECISION) * DEGREE_TO_RAD; + const double float_lat1 = (lat1 / COORDINATE_PRECISION) * detail::DEGREE_TO_RAD; + const double float_lon1 = (lon1 / COORDINATE_PRECISION) * detail::DEGREE_TO_RAD; + const double float_lat2 = (lat2 / COORDINATE_PRECISION) * detail::DEGREE_TO_RAD; + const double float_lon2 = (lon2 / COORDINATE_PRECISION) * detail::DEGREE_TO_RAD; const double x_value = (float_lon2 - float_lon1) * std::cos((float_lat1 + float_lat2) / 2.0); const double y_value = float_lat2 - float_lat1; - return std::hypot(x_value, y_value) * EARTH_RADIUS; + return std::hypot(x_value, y_value) * detail::EARTH_RADIUS; } std::pair projectPointOnSegment(const FloatCoordinate &source, @@ -125,10 +126,10 @@ double perpendicularDistance(const Coordinate segment_source, BOOST_ASSERT(query_location.IsValid()); FloatCoordinate projected_nearest; - std::tie(ratio, projected_nearest) = projectPointOnSegment(mercator::fromWGS84(segment_source), - mercator::fromWGS84(segment_target), - mercator::fromWGS84(query_location)); - nearest_location = mercator::toWGS84(projected_nearest); + std::tie(ratio, projected_nearest) = projectPointOnSegment( + web_mercator::fromWGS84(segment_source), web_mercator::fromWGS84(segment_target), + web_mercator::fromWGS84(query_location)); + nearest_location = web_mercator::toWGS84(projected_nearest); const double approximate_distance = greatCircleDistance(query_location, nearest_location); BOOST_ASSERT(0.0 <= approximate_distance); @@ -198,10 +199,10 @@ double computeAngle(const Coordinate first, const Coordinate second, const Coord const double v1x = static_cast(toFloating(first.lon - second.lon)); const double v1y = - mercator::latToY(toFloating(first.lat)) - mercator::latToY(toFloating(second.lat)); + web_mercator::latToY(toFloating(first.lat)) - web_mercator::latToY(toFloating(second.lat)); const double v2x = static_cast(toFloating(third.lon - second.lon)); const double v2y = - mercator::latToY(toFloating(third.lat)) - mercator::latToY(toFloating(second.lat)); + web_mercator::latToY(toFloating(third.lat)) - web_mercator::latToY(toFloating(second.lat)); double angle = (atan2_lookup(v2y, v2x) - atan2_lookup(v1y, v1x)) * 180. / pi(); @@ -302,107 +303,6 @@ Coordinate interpolateLinear(double factor, const Coordinate from, const Coordin return {std::move(interpolated_lon), std::move(interpolated_lat)}; } -namespace mercator -{ -FloatLatitude yToLat(const double y) -{ - const auto clamped_y = std::max(-180., std::min(180., y)); - const double normalized_lat = - RAD_TO_DEGREE * 2. * std::atan(std::exp(clamped_y * DEGREE_TO_RAD)); - - return FloatLatitude(normalized_lat - 90.); -} - -double latToY(const FloatLatitude latitude) -{ - // apparently this is the (faster) version of the canonical log(tan()) version - const double f = std::sin(DEGREE_TO_RAD * static_cast(latitude)); - const double y = 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; -} - -FloatLatitude clamp(const FloatLatitude lat) -{ - return std::max(std::min(lat, FloatLatitude(detail::MAX_LATITUDE)), - FloatLatitude(-detail::MAX_LATITUDE)); -} - -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; - x = (x - b) / shift * 360.0; - // FIXME needs to be simplified - const double g = (y - b) / -(shift / (2 * M_PI)) / DEGREE_TO_RAD; - static_assert(DEGREE_TO_RAD / (2 * M_PI) - 1 / 360. < 0.0001, ""); - y = static_cast(util::coordinate_calculation::mercator::yToLat(g)); -} - -double degreeToPixel(FloatLongitude lon, unsigned zoom) -{ - const double shift = (1u << zoom) * TILE_SIZE; - const double b = shift / 2.0; - const double x = b * (1 + static_cast(lon) / 180.0); - return x; -} - -double degreeToPixel(FloatLatitude lat, unsigned zoom) -{ - const double shift = (1u << zoom) * TILE_SIZE; - const double b = shift / 2.0; - const double y = b * (1. - latToY(lat) / 180.); - return y; -} - -FloatCoordinate fromWGS84(const FloatCoordinate &wgs84_coordinate) -{ - return {wgs84_coordinate.lon, - FloatLatitude{coordinate_calculation::mercator::latToY(wgs84_coordinate.lat)}}; -} - -FloatCoordinate toWGS84(const FloatCoordinate &mercator_coordinate) -{ - return {mercator_coordinate.lon, - coordinate_calculation::mercator::yToLat(static_cast(mercator_coordinate.lat))}; -} - -// Converts a WMS tile coordinate (z,x,y) into a wgs bounding box -void xyzToWGS84( - const int x, const int y, const int z, double &minx, double &miny, double &maxx, double &maxy) -{ - using util::coordinate_calculation::mercator::TILE_SIZE; - - minx = x * TILE_SIZE; - miny = (y + 1.0) * TILE_SIZE; - maxx = (x + 1.0) * TILE_SIZE; - maxy = y * TILE_SIZE; - // 2^z * TILE_SIZE - const double shift = (1u << static_cast(z)) * TILE_SIZE; - pixelToDegree(shift, minx, miny); - pixelToDegree(shift, maxx, maxy); -} - -// Converts a WMS tile coordinate (z,x,y) into a mercator bounding box -void xyzToMercator( - const int x, const int y, const int z, double &minx, double &miny, double &maxx, double &maxy) -{ - using namespace util::coordinate_calculation::mercator; - - xyzToWGS84(x, y, z, minx, miny, maxx, maxy); - - minx = static_cast(clamp(util::FloatLongitude(minx))) * DEGREE_TO_PX; - miny = latToY(clamp(util::FloatLatitude(miny))) * DEGREE_TO_PX; - maxx = static_cast(clamp(util::FloatLongitude(maxx))) * DEGREE_TO_PX; - maxy = latToY(clamp(util::FloatLatitude(maxy))) * DEGREE_TO_PX; -} - -} // ns mercato // ns mercatorr } // ns coordinate_calculation } // ns util } // ns osrm diff --git a/unit_tests/util/coordinate_calculation.cpp b/unit_tests/util/coordinate_calculation.cpp index d159441be..c14280f5f 100644 --- a/unit_tests/util/coordinate_calculation.cpp +++ b/unit_tests/util/coordinate_calculation.cpp @@ -27,72 +27,6 @@ BOOST_AUTO_TEST_CASE(regression_test_1347) BOOST_CHECK_LE(std::abs(d1 - d2), 0.01); } -BOOST_AUTO_TEST_CASE(lon_to_pixel) -{ - using namespace coordinate_calculation; - BOOST_CHECK_CLOSE(7.416042 * mercator::DEGREE_TO_PX, 825550.019142, 0.1); - BOOST_CHECK_CLOSE(7.415892 * mercator::DEGREE_TO_PX, 825533.321218, 0.1); - BOOST_CHECK_CLOSE(7.416016 * mercator::DEGREE_TO_PX, 825547.124835, 0.1); - BOOST_CHECK_CLOSE(7.41577 * mercator::DEGREE_TO_PX, 825519.74024, 0.1); - BOOST_CHECK_CLOSE(7.415808 * mercator::DEGREE_TO_PX, 825523.970381, 0.1); -} - -BOOST_AUTO_TEST_CASE(lat_to_pixel) -{ - using namespace coordinate_calculation; - BOOST_CHECK_CLOSE(mercator::latToY(util::FloatLatitude(43.733947)) * mercator::DEGREE_TO_PX, - 5424361.75863, 0.1); - BOOST_CHECK_CLOSE(mercator::latToY(util::FloatLatitude(43.733799)) * mercator::DEGREE_TO_PX, - 5424338.95731, 0.1); - BOOST_CHECK_CLOSE(mercator::latToY(util::FloatLatitude(43.733922)) * mercator::DEGREE_TO_PX, - 5424357.90705, 0.1); - BOOST_CHECK_CLOSE(mercator::latToY(util::FloatLatitude(43.733697)) * mercator::DEGREE_TO_PX, - 5424323.24293, 0.1); - BOOST_CHECK_CLOSE(mercator::latToY(util::FloatLatitude(43.733729)) * mercator::DEGREE_TO_PX, - 5424328.17293, 0.1); -} - -BOOST_AUTO_TEST_CASE(xyz_to_wgs84) -{ - using namespace coordinate_calculation; - - double minx_1; - double miny_1; - double maxx_1; - double maxy_1; - mercator::xyzToWGS84(2, 2, 1, minx_1, miny_1, maxx_1, maxy_1); - BOOST_CHECK_CLOSE(minx_1, 180, 0.0001); - BOOST_CHECK_CLOSE(miny_1, -85.0511, 0.0001); - BOOST_CHECK_CLOSE(maxx_1, 360, 0.0001); - BOOST_CHECK_CLOSE(maxy_1, -85.0511, 0.0001); - - double minx_2; - double miny_2; - double maxx_2; - double maxy_2; - mercator::xyzToWGS84(100, 0, 13, minx_2, miny_2, maxx_2, maxy_2); - BOOST_CHECK_CLOSE(minx_2, -175.6054, 0.0001); - BOOST_CHECK_CLOSE(miny_2, 85.0473, 0.0001); - BOOST_CHECK_CLOSE(maxx_2, -175.5615, 0.0001); - BOOST_CHECK_CLOSE(maxy_2, 85.0511, 0.0001); -} - -BOOST_AUTO_TEST_CASE(xyz_to_mercator) -{ - using namespace coordinate_calculation; - - double minx; - double miny; - double maxx; - double maxy; - mercator::xyzToMercator(100, 0, 13, minx, miny, maxx, maxy); - - BOOST_CHECK_CLOSE(minx, -19548311.361764118075, 0.0001); - BOOST_CHECK_CLOSE(miny, 20032616.372979003936, 0.0001); - BOOST_CHECK_CLOSE(maxx, -19543419.391953866929, 0.0001); - BOOST_CHECK_CLOSE(maxy, 20037508.342789277434, 0.0001); -} - BOOST_AUTO_TEST_CASE(regression_point_on_segment) { // ^ diff --git a/unit_tests/util/static_rtree.cpp b/unit_tests/util/static_rtree.cpp index 11f183088..b3c2c9e31 100644 --- a/unit_tests/util/static_rtree.cpp +++ b/unit_tests/util/static_rtree.cpp @@ -62,10 +62,10 @@ template class LinearSearchNN { std::vector local_edges(edges); - auto projected_input = coordinate_calculation::mercator::fromWGS84(input_coordinate); + auto projected_input = web_mercator::fromWGS84(input_coordinate); const auto segment_comparator = [this, &projected_input](const DataT &lhs, const DataT &rhs) { - using coordinate_calculation::mercator::fromWGS84; + using web_mercator::fromWGS84; const auto lhs_result = coordinate_calculation::projectPointOnSegment( fromWGS84(coords->at(lhs.u)), fromWGS84(coords->at(lhs.v)), projected_input); const auto rhs_result = coordinate_calculation::projectPointOnSegment( @@ -328,13 +328,6 @@ BOOST_AUTO_TEST_CASE(regression_test) auto result_rtree = rtree.Nearest(input, 1); auto result_ls = lsnn.Nearest(input, 1); - auto distance_rtree = coordinate_calculation::perpendicularDistance( - fixture.coords->at(result_rtree.front().u), fixture.coords->at(result_rtree.front().v), - input); - - auto distance_lsnn = coordinate_calculation::perpendicularDistance( - fixture.coords->at(result_ls.front().u), fixture.coords->at(result_ls.front().v), input); - BOOST_CHECK(result_rtree.size() == 1); BOOST_CHECK(result_ls.size() == 1); diff --git a/unit_tests/util/web_mercator.cpp b/unit_tests/util/web_mercator.cpp new file mode 100644 index 000000000..e29cc1a8c --- /dev/null +++ b/unit_tests/util/web_mercator.cpp @@ -0,0 +1,74 @@ +#include + +#include "util/web_mercator.hpp" + +#include + +#include + +using namespace osrm; +using namespace osrm::util; + +BOOST_AUTO_TEST_SUITE(web_mercator_tests) + +BOOST_AUTO_TEST_CASE(lon_to_pixel) +{ + BOOST_CHECK_CLOSE(7.416042 * web_mercator::DEGREE_TO_PX, 825550.019142, 0.1); + BOOST_CHECK_CLOSE(7.415892 * web_mercator::DEGREE_TO_PX, 825533.321218, 0.1); + BOOST_CHECK_CLOSE(7.416016 * web_mercator::DEGREE_TO_PX, 825547.124835, 0.1); + BOOST_CHECK_CLOSE(7.41577 * web_mercator::DEGREE_TO_PX, 825519.74024, 0.1); + BOOST_CHECK_CLOSE(7.415808 * web_mercator::DEGREE_TO_PX, 825523.970381, 0.1); +} + +BOOST_AUTO_TEST_CASE(lat_to_pixel) +{ + BOOST_CHECK_CLOSE(web_mercator::latToY(util::FloatLatitude(43.733947)) * web_mercator::DEGREE_TO_PX, + 5424361.75863, 0.1); + BOOST_CHECK_CLOSE(web_mercator::latToY(util::FloatLatitude(43.733799)) * web_mercator::DEGREE_TO_PX, + 5424338.95731, 0.1); + BOOST_CHECK_CLOSE(web_mercator::latToY(util::FloatLatitude(43.733922)) * web_mercator::DEGREE_TO_PX, + 5424357.90705, 0.1); + BOOST_CHECK_CLOSE(web_mercator::latToY(util::FloatLatitude(43.733697)) * web_mercator::DEGREE_TO_PX, + 5424323.24293, 0.1); + BOOST_CHECK_CLOSE(web_mercator::latToY(util::FloatLatitude(43.733729)) * web_mercator::DEGREE_TO_PX, + 5424328.17293, 0.1); +} + +BOOST_AUTO_TEST_CASE(xyz_to_wgs84) +{ + double minx_1; + double miny_1; + double maxx_1; + double maxy_1; + web_mercator::xyzToWGS84(2, 2, 1, minx_1, miny_1, maxx_1, maxy_1); + BOOST_CHECK_CLOSE(minx_1, 180, 0.0001); + BOOST_CHECK_CLOSE(miny_1, -85.0511, 0.0001); + BOOST_CHECK_CLOSE(maxx_1, 360, 0.0001); + BOOST_CHECK_CLOSE(maxy_1, -85.0511, 0.0001); + + double minx_2; + double miny_2; + double maxx_2; + double maxy_2; + web_mercator::xyzToWGS84(100, 0, 13, minx_2, miny_2, maxx_2, maxy_2); + BOOST_CHECK_CLOSE(minx_2, -175.6054, 0.0001); + BOOST_CHECK_CLOSE(miny_2, 85.0473, 0.0001); + BOOST_CHECK_CLOSE(maxx_2, -175.5615, 0.0001); + BOOST_CHECK_CLOSE(maxy_2, 85.0511, 0.0001); +} + +BOOST_AUTO_TEST_CASE(xyz_to_mercator) +{ + double minx; + double miny; + double maxx; + double maxy; + web_mercator::xyzToMercator(100, 0, 13, minx, miny, maxx, maxy); + + BOOST_CHECK_CLOSE(minx, -19548311.361764118075, 0.0001); + BOOST_CHECK_CLOSE(miny, 20032616.372979003936, 0.0001); + BOOST_CHECK_CLOSE(maxx, -19543419.391953866929, 0.0001); + BOOST_CHECK_CLOSE(maxy, 20037508.342789277434, 0.0001); +} + +BOOST_AUTO_TEST_SUITE_END()