From 93a2e667042dcb19dd4ad723cf1b59048528c60d Mon Sep 17 00:00:00 2001 From: Mortada Mehyar Date: Sat, 26 Dec 2015 11:12:10 -0800 Subject: [PATCH] use double precision calculations instead of mixing double and float --- include/engine/geospatial_query.hpp | 26 ++--- .../routing_algorithms/map_matching.hpp | 2 +- include/osrm/coordinate.hpp | 4 +- include/util/coordinate_calculation.hpp | 22 ++--- include/util/rectangle.hpp | 10 +- src/util/coordinate.cpp | 2 +- src/util/coordinate_calculation.cpp | 94 +++++++++---------- unit_tests/util/coordinate.cpp | 8 +- unit_tests/util/static_rtree.cpp | 20 ++-- 9 files changed, 94 insertions(+), 94 deletions(-) diff --git a/include/engine/geospatial_query.hpp b/include/engine/geospatial_query.hpp index 05105554e..95a41da93 100644 --- a/include/engine/geospatial_query.hpp +++ b/include/engine/geospatial_query.hpp @@ -30,7 +30,7 @@ template class GeospatialQuery // Does not filter by small/big component! std::vector NearestPhantomNodesInRange(const FixedPointCoordinate &input_coordinate, - const float max_distance, + const double max_distance, const int bearing = 0, const int bearing_range = 180) { @@ -40,7 +40,7 @@ template class GeospatialQuery { return checkSegmentBearing(data, bearing, bearing_range); }, - [max_distance](const std::size_t, const float min_dist) + [max_distance](const std::size_t, const double min_dist) { return min_dist > max_distance; }); @@ -61,7 +61,7 @@ template class GeospatialQuery { return checkSegmentBearing(data, bearing, bearing_range); }, - [max_results](const std::size_t num_results, const float) + [max_results](const std::size_t num_results, const double) { return num_results >= max_results; }); @@ -99,7 +99,7 @@ template class GeospatialQuery return use_directions; }, - [&has_big_component](const std::size_t num_results, const float) + [&has_big_component](const std::size_t num_results, const double) { return num_results > 0 && has_big_component; }); @@ -132,7 +132,7 @@ template class GeospatialQuery const EdgeData &data) const { FixedPointCoordinate point_on_segment; - float ratio; + double ratio; const auto current_perpendicular_distance = coordinate_calculation::perpendicular_distance( coordinates->at(data.u), coordinates->at(data.v), input_coordinate, point_on_segment, ratio); @@ -140,7 +140,7 @@ template class GeospatialQuery auto transformed = PhantomNodeWithDistance { PhantomNode{data, point_on_segment}, current_perpendicular_distance }; - ratio = std::min(1.f, std::max(0.f, ratio)); + ratio = std::min(1.0, std::max(0.0, ratio)); if (SPECIAL_NODEID != transformed.phantom_node.forward_node_id) { @@ -148,27 +148,27 @@ template class GeospatialQuery } if (SPECIAL_NODEID != transformed.phantom_node.reverse_node_id) { - transformed.phantom_node.reverse_weight *= 1.f - ratio; + transformed.phantom_node.reverse_weight *= 1.0 - ratio; } return transformed; } std::pair checkSegmentBearing(const EdgeData &segment, - const float filter_bearing, - const float filter_bearing_range) + const int filter_bearing, + const int filter_bearing_range) { - const float forward_edge_bearing = + const double forward_edge_bearing = coordinate_calculation::bearing(coordinates->at(segment.u), coordinates->at(segment.v)); - const float backward_edge_bearing = (forward_edge_bearing + 180) > 360 + const double backward_edge_bearing = (forward_edge_bearing + 180) > 360 ? (forward_edge_bearing - 180) : (forward_edge_bearing + 180); const bool forward_bearing_valid = - bearing::CheckInBounds(forward_edge_bearing, filter_bearing, filter_bearing_range) && + bearing::CheckInBounds(std::round(forward_edge_bearing), filter_bearing, filter_bearing_range) && segment.forward_edge_based_node_id != SPECIAL_NODEID; const bool backward_bearing_valid = - bearing::CheckInBounds(backward_edge_bearing, filter_bearing, filter_bearing_range) && + bearing::CheckInBounds(std::round(backward_edge_bearing), filter_bearing, filter_bearing_range) && segment.reverse_edge_based_node_id != SPECIAL_NODEID; return std::make_pair(forward_bearing_valid, backward_bearing_valid); } diff --git a/include/engine/routing_algorithms/map_matching.hpp b/include/engine/routing_algorithms/map_matching.hpp index 9afc2653a..00db48eff 100644 --- a/include/engine/routing_algorithms/map_matching.hpp +++ b/include/engine/routing_algorithms/map_matching.hpp @@ -359,7 +359,7 @@ class MapMatching final : public BasicRoutingInterface(0u, reconstructed_indices.size())) diff --git a/include/osrm/coordinate.hpp b/include/osrm/coordinate.hpp index 6318465e1..46029cbc9 100644 --- a/include/osrm/coordinate.hpp +++ b/include/osrm/coordinate.hpp @@ -34,7 +34,7 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. namespace { -constexpr static const float COORDINATE_PRECISION = 1000000.f; +constexpr static const double COORDINATE_PRECISION = 1000000.0; } struct FixedPointCoordinate @@ -58,7 +58,7 @@ struct FixedPointCoordinate bool is_valid() const; bool operator==(const FixedPointCoordinate &other) const; - float bearing(const FixedPointCoordinate &other) const; + double bearing(const FixedPointCoordinate &other) const; void output(std::ostream &out) const; }; diff --git a/include/util/coordinate_calculation.hpp b/include/util/coordinate_calculation.hpp index 80ec7c289..0ab9fe617 100644 --- a/include/util/coordinate_calculation.hpp +++ b/include/util/coordinate_calculation.hpp @@ -41,41 +41,41 @@ namespace coordinate_calculation double haversine_distance(const FixedPointCoordinate &first_coordinate, const FixedPointCoordinate &second_coordinate); - float great_circle_distance(const FixedPointCoordinate &first_coordinate, + double great_circle_distance(const FixedPointCoordinate &first_coordinate, const FixedPointCoordinate &second_coordinate); - float great_circle_distance(const int lat1, const int lon1, const int lat2, const int lon2); + double great_circle_distance(const int lat1, const int lon1, const int lat2, const int lon2); void lat_or_lon_to_string(const int value, std::string &output); - float perpendicular_distance(const FixedPointCoordinate &segment_source, + double perpendicular_distance(const FixedPointCoordinate &segment_source, const FixedPointCoordinate &segment_target, const FixedPointCoordinate &query_location); - float perpendicular_distance(const FixedPointCoordinate &segment_source, + double perpendicular_distance(const FixedPointCoordinate &segment_source, const FixedPointCoordinate &segment_target, const FixedPointCoordinate &query_location, FixedPointCoordinate &nearest_location, - float &ratio); + double &ratio); - float perpendicular_distance_from_projected_coordinate( + double perpendicular_distance_from_projected_coordinate( const FixedPointCoordinate &segment_source, const FixedPointCoordinate &segment_target, const FixedPointCoordinate &query_location, const std::pair &projected_coordinate); - float perpendicular_distance_from_projected_coordinate( + double perpendicular_distance_from_projected_coordinate( const FixedPointCoordinate &segment_source, const FixedPointCoordinate &segment_target, const FixedPointCoordinate &query_location, const std::pair &projected_coordinate, FixedPointCoordinate &nearest_location, - float &ratio); + double &ratio); - float deg_to_rad(const float degree); - float rad_to_deg(const float radian); + double deg_to_rad(const double degree); + double rad_to_deg(const double radian); - float bearing(const FixedPointCoordinate &first_coordinate, + double bearing(const FixedPointCoordinate &first_coordinate, const FixedPointCoordinate &second_coordinate); } diff --git a/include/util/rectangle.hpp b/include/util/rectangle.hpp index f4dd5280d..9aa99ac2c 100644 --- a/include/util/rectangle.hpp +++ b/include/util/rectangle.hpp @@ -84,7 +84,7 @@ struct RectangleInt2D Contains(lower_left)); } - float GetMinDist(const FixedPointCoordinate &location) const + double GetMinDist(const FixedPointCoordinate &location) const { const bool is_contained = Contains(location); if (is_contained) @@ -117,7 +117,7 @@ struct RectangleInt2D BOOST_ASSERT(d != INVALID); - float min_dist = std::numeric_limits::max(); + double min_dist = std::numeric_limits::max(); switch (d) { case NORTH: @@ -156,14 +156,14 @@ struct RectangleInt2D break; } - BOOST_ASSERT(min_dist < std::numeric_limits::max()); + BOOST_ASSERT(min_dist < std::numeric_limits::max()); return min_dist; } - float GetMinMaxDist(const FixedPointCoordinate &location) const + double GetMinMaxDist(const FixedPointCoordinate &location) const { - float min_max_dist = std::numeric_limits::max(); + double min_max_dist = std::numeric_limits::max(); // Get minmax distance to each of the four sides const FixedPointCoordinate upper_left(max_lat, min_lon); const FixedPointCoordinate upper_right(max_lat, max_lon); diff --git a/src/util/coordinate.cpp b/src/util/coordinate.cpp index 63f81aadb..6124dd8fe 100644 --- a/src/util/coordinate.cpp +++ b/src/util/coordinate.cpp @@ -81,7 +81,7 @@ void FixedPointCoordinate::output(std::ostream &out) const out << "(" << lat / COORDINATE_PRECISION << "," << lon / COORDINATE_PRECISION << ")"; } -float FixedPointCoordinate::bearing(const FixedPointCoordinate &other) const +double FixedPointCoordinate::bearing(const FixedPointCoordinate &other) const { return coordinate_calculation::bearing(other, *this); } diff --git a/src/util/coordinate_calculation.cpp b/src/util/coordinate_calculation.cpp index 54ebaf947..e1975ddeb 100644 --- a/src/util/coordinate_calculation.cpp +++ b/src/util/coordinate_calculation.cpp @@ -40,10 +40,10 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. namespace { -constexpr static const float RAD = 0.017453292519943295769236907684886f; +constexpr static const double RAD = 0.017453292519943295769236907684886; // 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) -constexpr static const float earth_radius = 6372797.560856f; +constexpr static const double earth_radius = 6372797.560856; } namespace coordinate_calculation @@ -84,14 +84,14 @@ double haversine_distance(const FixedPointCoordinate &coordinate_1, coordinate_2.lon); } -float great_circle_distance(const FixedPointCoordinate &coordinate_1, +double great_circle_distance(const FixedPointCoordinate &coordinate_1, const FixedPointCoordinate &coordinate_2) { return great_circle_distance(coordinate_1.lat, coordinate_1.lon, coordinate_2.lat, coordinate_2.lon); } -float great_circle_distance(const int lat1, +double great_circle_distance(const int lat1, const int lon1, const int lat2, const int lon2) @@ -101,32 +101,32 @@ float great_circle_distance(const int lat1, BOOST_ASSERT(lat2 != std::numeric_limits::min()); BOOST_ASSERT(lon2 != std::numeric_limits::min()); - const float float_lat1 = (lat1 / COORDINATE_PRECISION) * RAD; - const float float_lon1 = (lon1 / COORDINATE_PRECISION) * RAD; - const float float_lat2 = (lat2 / COORDINATE_PRECISION) * RAD; - const float float_lon2 = (lon2 / COORDINATE_PRECISION) * RAD; + const double float_lat1 = (lat1 / COORDINATE_PRECISION) * RAD; + const double float_lon1 = (lon1 / COORDINATE_PRECISION) * RAD; + const double float_lat2 = (lat2 / COORDINATE_PRECISION) * RAD; + const double float_lon2 = (lon2 / COORDINATE_PRECISION) * RAD; - const float x_value = (float_lon2 - float_lon1) * std::cos((float_lat1 + float_lat2) / 2.f); - const float y_value = float_lat2 - float_lat1; + 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; } -float perpendicular_distance(const FixedPointCoordinate &source_coordinate, +double perpendicular_distance(const FixedPointCoordinate &source_coordinate, const FixedPointCoordinate &target_coordinate, const FixedPointCoordinate &query_location) { - float ratio; + double ratio; FixedPointCoordinate nearest_location; return perpendicular_distance(source_coordinate, target_coordinate, query_location, nearest_location, ratio); } -float perpendicular_distance(const FixedPointCoordinate &segment_source, +double perpendicular_distance(const FixedPointCoordinate &segment_source, const FixedPointCoordinate &segment_target, const FixedPointCoordinate &query_location, FixedPointCoordinate &nearest_location, - float &ratio) + double &ratio) { return perpendicular_distance_from_projected_coordinate( segment_source, segment_target, query_location, @@ -135,13 +135,13 @@ float perpendicular_distance(const FixedPointCoordinate &segment_source, nearest_location, ratio); } -float perpendicular_distance_from_projected_coordinate( +double perpendicular_distance_from_projected_coordinate( const FixedPointCoordinate &source_coordinate, const FixedPointCoordinate &target_coordinate, const FixedPointCoordinate &query_location, const std::pair &projected_coordinate) { - float ratio; + double ratio; FixedPointCoordinate nearest_location; return perpendicular_distance_from_projected_coordinate(source_coordinate, target_coordinate, @@ -149,13 +149,13 @@ float perpendicular_distance_from_projected_coordinate( nearest_location, ratio); } -float perpendicular_distance_from_projected_coordinate( +double perpendicular_distance_from_projected_coordinate( const FixedPointCoordinate &segment_source, const FixedPointCoordinate &segment_target, const FixedPointCoordinate &query_location, const std::pair &projected_coordinate, FixedPointCoordinate &nearest_location, - float &ratio) + double &ratio) { BOOST_ASSERT(query_location.is_valid()); @@ -171,7 +171,7 @@ float perpendicular_distance_from_projected_coordinate( { const double m = (d - b) / (c - a); // slope // Projection of (x,y) on line joining (a,b) and (c,d) - p = ((x + (m * y)) + (m * m * a - m * b)) / (1.f + m * m); + p = ((x + (m * y)) + (m * m * a - m * b)) / (1.0 + m * m); q = b + m * (p - a); } else @@ -182,36 +182,36 @@ float perpendicular_distance_from_projected_coordinate( nY = (d * p - c * q) / (a * d - b * c); // discretize the result to coordinate precision. it's a hack! - if (std::abs(nY) < (1.f / COORDINATE_PRECISION)) + if (std::abs(nY) < (1.0 / COORDINATE_PRECISION)) { - nY = 0.f; + nY = 0.0; } // compute ratio ratio = - static_cast((p - nY * a) / c); // These values are actually n/m+n and m/m+n , we need + static_cast((p - nY * a) / c); // These values are actually n/m+n and m/m+n , we need // not calculate the explicit values of m an n as we // are just interested in the ratio if (std::isnan(ratio)) { - ratio = (segment_target == query_location ? 1.f : 0.f); + ratio = (segment_target == query_location ? 1.0 : 0.0); } - else if (std::abs(ratio) <= std::numeric_limits::epsilon()) + else if (std::abs(ratio) <= std::numeric_limits::epsilon()) { - ratio = 0.f; + ratio = 0.0; } - else if (std::abs(ratio - 1.f) <= std::numeric_limits::epsilon()) + else if (std::abs(ratio - 1.0) <= std::numeric_limits::epsilon()) { - ratio = 1.f; + ratio = 1.0; } // compute nearest location BOOST_ASSERT(!std::isnan(ratio)); - if (ratio <= 0.f) + if (ratio <= 0.0) { nearest_location = segment_source; } - else if (ratio >= 1.f) + else if (ratio >= 1.0) { nearest_location = segment_target; } @@ -223,9 +223,9 @@ float perpendicular_distance_from_projected_coordinate( } BOOST_ASSERT(nearest_location.is_valid()); - const float approximate_distance = + const double approximate_distance = great_circle_distance(query_location, nearest_location); - BOOST_ASSERT(0.f <= approximate_distance); + BOOST_ASSERT(0.0 <= approximate_distance); return approximate_distance; } @@ -236,36 +236,36 @@ void lat_or_lon_to_string(const int value, std::string &output) output = printInt<11, 6>(buffer, value); } -float deg_to_rad(const float degree) +double deg_to_rad(const double degree) { - return degree * (static_cast(M_PI) / 180.f); + return degree * (static_cast(M_PI) / 180.0); } -float rad_to_deg(const float radian) +double rad_to_deg(const double radian) { - return radian * (180.f * static_cast(M_1_PI)); + return radian * (180.0 * static_cast(M_1_PI)); } -float bearing(const FixedPointCoordinate &first_coordinate, +double bearing(const FixedPointCoordinate &first_coordinate, const FixedPointCoordinate &second_coordinate) { - const float lon_diff = + const double lon_diff = second_coordinate.lon / COORDINATE_PRECISION - first_coordinate.lon / COORDINATE_PRECISION; - const float lon_delta = deg_to_rad(lon_diff); - const float lat1 = deg_to_rad(first_coordinate.lat / COORDINATE_PRECISION); - const float lat2 = deg_to_rad(second_coordinate.lat / COORDINATE_PRECISION); - const float y = std::sin(lon_delta) * std::cos(lat2); - const float x = + const double lon_delta = deg_to_rad(lon_diff); + const double lat1 = deg_to_rad(first_coordinate.lat / COORDINATE_PRECISION); + const double lat2 = deg_to_rad(second_coordinate.lat / COORDINATE_PRECISION); + const double y = std::sin(lon_delta) * std::cos(lat2); + const double x = std::cos(lat1) * std::sin(lat2) - std::sin(lat1) * std::cos(lat2) * std::cos(lon_delta); - float result = rad_to_deg(std::atan2(y, x)); - while (result < 0.f) + double result = rad_to_deg(std::atan2(y, x)); + while (result < 0.0) { - result += 360.f; + result += 360.0; } - while (result >= 360.f) + while (result >= 360.0) { - result -= 360.f; + result -= 360.0; } return result; } diff --git a/unit_tests/util/coordinate.cpp b/unit_tests/util/coordinate.cpp index f153175ab..c4f361921 100644 --- a/unit_tests/util/coordinate.cpp +++ b/unit_tests/util/coordinate.cpp @@ -40,11 +40,11 @@ BOOST_AUTO_TEST_CASE(regression_test_1347) FixedPointCoordinate v(10.001 * COORDINATE_PRECISION, -100.002 * COORDINATE_PRECISION); FixedPointCoordinate q(10.002 * COORDINATE_PRECISION, -100.001 * COORDINATE_PRECISION); - float d1 = coordinate_calculation::perpendicular_distance(u, v, q); + double d1 = coordinate_calculation::perpendicular_distance(u, v, q); - float ratio; + double ratio; FixedPointCoordinate nearest_location; - float d2 = coordinate_calculation::perpendicular_distance(u, v, q, nearest_location, ratio); + double d2 = coordinate_calculation::perpendicular_distance(u, v, q, nearest_location, ratio); - BOOST_CHECK_LE(std::abs(d1 - d2), 0.01f); + BOOST_CHECK_LE(std::abs(d1 - d2), 0.01); } diff --git a/unit_tests/util/static_rtree.cpp b/unit_tests/util/static_rtree.cpp index 33c9eed74..93b81b906 100644 --- a/unit_tests/util/static_rtree.cpp +++ b/unit_tests/util/static_rtree.cpp @@ -88,11 +88,11 @@ template class LinearSearchNN local_edges.begin(), local_edges.begin() + num_results, local_edges.end(), [this, &input_coordinate](const DataT &lhs, const DataT &rhs) { - float current_ratio = 0.; + double current_ratio = 0.; FixedPointCoordinate nearest; - const float lhs_dist = coordinate_calculation::perpendicular_distance( + const double lhs_dist = coordinate_calculation::perpendicular_distance( coords->at(lhs.u), coords->at(lhs.v), input_coordinate, nearest, current_ratio); - const float rhs_dist = coordinate_calculation::perpendicular_distance( + const double rhs_dist = coordinate_calculation::perpendicular_distance( coords->at(rhs.u), coords->at(rhs.v), input_coordinate, nearest, current_ratio); return lhs_dist < rhs_dist; }); @@ -165,7 +165,7 @@ template struct RandomGraphFixture struct GraphFixture { - GraphFixture(const std::vector> &input_coords, + GraphFixture(const std::vector> &input_coords, const std::vector> &input_edges) : coords(std::make_shared>()) { @@ -253,13 +253,13 @@ void sampling_verify_rtree(RTreeT &rtree, LinearSearchNN &lsnn, const auto lsnn_u = result_lsnn.back().u; auto lsnn_v = result_lsnn.back().v; - float current_ratio = 0.; + double current_ratio = 0.; FixedPointCoordinate nearest; - const float rtree_dist = coordinate_calculation::perpendicular_distance( + const double rtree_dist = coordinate_calculation::perpendicular_distance( coords[rtree_u], coords[rtree_v], q, nearest, current_ratio); - const float lsnn_dist = coordinate_calculation::perpendicular_distance( + const double lsnn_dist = coordinate_calculation::perpendicular_distance( coords[lsnn_u], coords[lsnn_v], q, nearest, current_ratio); - BOOST_CHECK_LE(std::abs(rtree_dist - lsnn_dist), std::numeric_limits::epsilon()); + BOOST_CHECK_LE(std::abs(rtree_dist - lsnn_dist), std::numeric_limits::epsilon()); } } @@ -324,7 +324,7 @@ BOOST_FIXTURE_TEST_CASE(construct_multiple_levels_test, TestRandomGraphFixture_M // one BB will be pruned, even if it could contain a nearer match. BOOST_AUTO_TEST_CASE(regression_test) { - using Coord = std::pair; + using Coord = std::pair; using Edge = std::pair; GraphFixture fixture( { @@ -420,7 +420,7 @@ BOOST_AUTO_TEST_CASE(rectangle_test) BOOST_AUTO_TEST_CASE(bearing_tests) { - using Coord = std::pair; + using Coord = std::pair; using Edge = std::pair; GraphFixture fixture( {