From be24689b0f50d66553e583dce3f68eaf4169c06b Mon Sep 17 00:00:00 2001
From: Patrick Niklaus <patrick.niklaus@student.kit.edu>
Date: Thu, 26 Apr 2018 20:38:17 +0000
Subject: [PATCH] Optimize distance calculation by getting rid of rounding

---
 src/util/coordinate_calculation.cpp | 33 +++++++++++++++++++----------
 1 file changed, 22 insertions(+), 11 deletions(-)

diff --git a/src/util/coordinate_calculation.cpp b/src/util/coordinate_calculation.cpp
index 388d41ef8..a45eb5209 100644
--- a/src/util/coordinate_calculation.cpp
+++ b/src/util/coordinate_calculation.cpp
@@ -27,26 +27,36 @@ class CheapRulerContainer
   public:
     CheapRulerContainer(const int number_of_rulers)
         : cheap_ruler_cache(number_of_rulers, mapbox::cheap_ruler::CheapRuler(0)),
-          step(90.0 / number_of_rulers)
+          step(90.0 * COORDINATE_PRECISION / number_of_rulers)
     {
         for (int n = 0; n < number_of_rulers; n++)
         {
             cheap_ruler_cache[n] = mapbox::cheap_ruler::CheapRuler(
-                step * (n + 0.5), mapbox::cheap_ruler::CheapRuler::Meters);
+                step * (n + 0.5) / COORDINATE_PRECISION, mapbox::cheap_ruler::CheapRuler::Meters);
         }
     };
-    mapbox::cheap_ruler::CheapRuler &getRuler(const double lat)
+
+    mapbox::cheap_ruler::CheapRuler &getRuler(const FixedLatitude lat_1, const FixedLatitude lat_2)
     {
-        return cheap_ruler_cache[std::min((int)std::floor(abs(lat) / step),
-                                          (int)cheap_ruler_cache.size() - 1)];
+        auto lat = (lat_1 + lat_2) / util::FixedLatitude{2};
+        return getRuler(lat);
+    }
+
+    mapbox::cheap_ruler::CheapRuler &getRuler(const FixedLatitude lat)
+    {
+        BOOST_ASSERT(step > 2);
+        // the |lat| > 0  ->  |lat|-1 > -1  ->  (|lat|-1)/step > -1/step > -1/2 >= -1  ->  bin >= 0
+        std::size_t bin = (std::abs(static_cast<int>(lat)) - 1) / step;
+        BOOST_ASSERT(bin < cheap_ruler_cache.size());
+        return cheap_ruler_cache[bin];
     };
 
   private:
     std::vector<mapbox::cheap_ruler::CheapRuler> cheap_ruler_cache;
-    const double step;
+    const int step;
 };
 CheapRulerContainer cheap_ruler_container(1800);
-}
+} // namespace
 
 // Does not project the coordinates!
 std::uint64_t squaredEuclideanDistance(const Coordinate lhs, const Coordinate rhs)
@@ -72,7 +82,8 @@ double fccApproximateDistance(const Coordinate coordinate_1, const Coordinate co
     const auto lat1 = static_cast<double>(util::toFloating(coordinate_1.lat));
     const auto lon2 = static_cast<double>(util::toFloating(coordinate_2.lon));
     const auto lat2 = static_cast<double>(util::toFloating(coordinate_2.lat));
-    return cheap_ruler_container.getRuler(0.5 * (lat1 + lat2)).distance({lon1, lat1}, {lon2, lat2});
+    return cheap_ruler_container.getRuler(coordinate_1.lat, coordinate_2.lat)
+        .distance({lon1, lat1}, {lon2, lat2});
 }
 
 double haversineDistance(const Coordinate coordinate_1, const Coordinate coordinate_2)
@@ -466,6 +477,6 @@ double computeArea(const std::vector<Coordinate> &polygon)
     return area / 2.;
 }
 
-} // ns coordinate_calculation
-} // ns util
-} // ns osrm
+} // namespace coordinate_calculation
+} // namespace util
+} // namespace osrm