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

use double precision calculations instead of mixing double and float

Browse Source
This commit is contained in:
Mortada Mehyar 2015-12-26 11:12:10 -08:00
parent facbe2c012
commit 93a2e66704
9 changed files with 94 additions and 94 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
include/engine/geospatial_query.hpp
View File

@ -30,7 +30,7 @@ template <typename RTreeT> class GeospatialQuery
// Does not filter by small/big component!
std::vector<PhantomNodeWithDistance>
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 <typename RTreeT> 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 <typename RTreeT> 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 <typename RTreeT> 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 <typename RTreeT> 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 <typename RTreeT> 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 <typename RTreeT> 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<bool, bool> 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);
}

2
include/engine/routing_algorithms/map_matching.hpp
View File

@ -359,7 +359,7 @@ class MapMatching final : public BasicRoutingInterface<DataFacadeT, MapMatching<
continue;
}
matching.length = 0.0f;
matching.length = 0.0;
matching.nodes.resize(reconstructed_indices.size());
matching.indices.resize(reconstructed_indices.size());
for (const auto i : osrm::irange<std::size_t>(0u, reconstructed_indices.size()))

4
include/osrm/coordinate.hpp
View File

@ -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;
};

22
include/util/coordinate_calculation.hpp
View File

@ -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<double, double> &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<double, double> &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);
}

10
include/util/rectangle.hpp
View File

@ -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<float>::max();
double min_dist = std::numeric_limits<double>::max();
switch (d)
{
case NORTH:
@ -156,14 +156,14 @@ struct RectangleInt2D
break;
}
BOOST_ASSERT(min_dist < std::numeric_limits<float>::max());
BOOST_ASSERT(min_dist < std::numeric_limits<double>::max());
return min_dist;
}
float GetMinMaxDist(const FixedPointCoordinate &location) const
double GetMinMaxDist(const FixedPointCoordinate &location) const
{
float min_max_dist = std::numeric_limits<float>::max();
double min_max_dist = std::numeric_limits<double>::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);

2
src/util/coordinate.cpp
View File

@ -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);
}

94
src/util/coordinate_calculation.cpp
View File

@ -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<int>::min());
BOOST_ASSERT(lon2 != std::numeric_limits<int>::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<double, double> &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<double, double> &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<float>((p - nY * a) / c); // These values are actually n/m+n and m/m+n , we need
static_cast<double>((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<float>::epsilon())
else if (std::abs(ratio) <= std::numeric_limits<double>::epsilon())
{
ratio = 0.f;
ratio = 0.0;
}
else if (std::abs(ratio - 1.f) <= std::numeric_limits<float>::epsilon())
else if (std::abs(ratio - 1.0) <= std::numeric_limits<double>::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<float>(M_PI) / 180.f);
return degree * (static_cast<double>(M_PI) / 180.0);
}
float rad_to_deg(const float radian)
double rad_to_deg(const double radian)
{
return radian * (180.f * static_cast<float>(M_1_PI));
return radian * (180.0 * static_cast<double>(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;
}

8
unit_tests/util/coordinate.cpp
View File

@ -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);
}

20
unit_tests/util/static_rtree.cpp
View File

@ -88,11 +88,11 @@ template <typename DataT> 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 <unsigned NUM_NODES, unsigned NUM_EDGES> struct RandomGraphFixture
struct GraphFixture
{
GraphFixture(const std::vector<std::pair<float, float>> &input_coords,
GraphFixture(const std::vector<std::pair<double, double>> &input_coords,
const std::vector<std::pair<unsigned, unsigned>> &input_edges)
: coords(std::make_shared<std::vector<FixedPointCoordinate>>())
{
@ -253,13 +253,13 @@ void sampling_verify_rtree(RTreeT &rtree, LinearSearchNN<TestData> &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<float>::epsilon());
BOOST_CHECK_LE(std::abs(rtree_dist - lsnn_dist), std::numeric_limits<double>::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<float, float>;
using Coord = std::pair<double, double>;
using Edge = std::pair<unsigned, unsigned>;
GraphFixture fixture(
{
@ -420,7 +420,7 @@ BOOST_AUTO_TEST_CASE(rectangle_test)
BOOST_AUTO_TEST_CASE(bearing_tests)
{
using Coord = std::pair<float, float>;
using Coord = std::pair<double, double>;
using Edge = std::pair<unsigned, unsigned>;
GraphFixture fixture(
{