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

Move projection function into own header and inline

Browse Source
This commit is contained in:
Patrick Niklaus
2016-04-09 02:18:47 +02:00
parent 9a617f5d41
commit 5052c4ae3a
11 changed files with 293 additions and 286 deletions
Download Patch File Download Diff File Expand all files Collapse all files
include/util/coordinate_calculation.hpp
+3 -32
View File
@@ -3,7 +3,6 @@
#include "util/coordinate.hpp"
#include <boost/math/constants/constants.hpp>
#include <boost/optional.hpp>
#include <utility>
@@ -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<double>();
// ^ 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<double>());
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
include/util/static_rtree.hpp
+39 -34
View File
@@ -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<uint64_t>(0, m_element_count),
[&input_data_vector, &input_wrapper_vector,
&coordinate_list](const tbb::blocked_range<uint64_t> &range) {
&coordinate_list](const tbb::blocked_range<uint64_t> &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<std::uint32_t>(0, search_tree_size),
[this, &search_tree_size](const tbb::blocked_range<std::uint32_t> &range) {
for (std::uint32_t i = range.begin(), end = range.end(); i != end; ++i)
{
TreeNode &current_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<std::uint32_t>(0, search_tree_size),
[this, &search_tree_size](const tbb::blocked_range<std::uint32_t> &range)
{
for (std::uint32_t i = range.begin(), end = range.end(); i != end;
++i)
{
TreeNode &current_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<EdgeDataT> results;
std::queue<TreeNode> traversal_queue;
@@ -391,8 +394,12 @@ class StaticRTree
std::vector<EdgeDataT> 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<EdgeDataT> 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 &current_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.));
include/util/viewport.hpp
+10 -10
View File
@@ -2,7 +2,7 @@
#define UTIL_VIEWPORT_HPP
#include "util/coordinate.hpp"
#include "util/coordinate_calculation.hpp"
#include "util/web_mercator.hpp"
#include <boost/assert.hpp>
@@ -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<unsigned>(detail::MIN_ZOOM, zoom);
else
return detail::MIN_ZOOM;
return detail::MIN_ZOOM;
}
}
}
}
include/util/web_mercator.hpp
+127
View File
@@ -0,0 +1,127 @@
#ifndef OSRM_WEB_MERCATOR_HPP
#define OSRM_WEB_MERCATOR_HPP
#include "util/coordinate.hpp"
#include <boost/math/constants/constants.hpp>
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<double>();
// ^ 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<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;
}
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<double>(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<double>(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<double>(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<unsigned>(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<double>(clamp(util::FloatLongitude(minx))) * DEGREE_TO_PX;
miny = latToY(clamp(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;
}
}
}
}
#endif