Use faster method for calculating distances.

CMakeLists.txt

@@ -559,6 +559,13 @@ include_directories(SYSTEM ${RAPIDJSON_INCLUDE_DIR})
 
 set(MICROTAR_INCLUDE_DIR "${CMAKE_CURRENT_SOURCE_DIR}/third_party/microtar/src")
 include_directories(SYSTEM ${MICROTAR_INCLUDE_DIR})
+
+set(MBXGEOM_INCLUDE_DIR "${CMAKE_CURRENT_SOURCE_DIR}/third_party/geometry.hpp-0.9.2/include")
+include_directories(SYSTEM ${MBXGEOM_INCLUDE_DIR})
+
+set(CHEAPRULER_INCLUDE_DIR "${CMAKE_CURRENT_SOURCE_DIR}/third_party/cheap-ruler-cpp-2.5.4/include")
+include_directories(SYSTEM ${CHEAPRULER_INCLUDE_DIR})
+
 add_library(MICROTAR OBJECT "${CMAKE_CURRENT_SOURCE_DIR}/third_party/microtar/src/microtar.c")
 set_property(TARGET MICROTAR PROPERTY POSITION_INDEPENDENT_CODE ON)
 

include/engine/routing_algorithms/routing_base.hpp

@@ -415,7 +415,7 @@ template <typename FacadeT> EdgeDistance computeEdgeDistance(const FacadeT &faca
     auto geometry_range = facade.GetUncompressedForwardGeometry(geometry_index.id);
     for (auto current = geometry_range.begin(); current < geometry_range.end() - 1; ++current)
     {
-        total_distance += util::coordinate_calculation::haversineDistance(
+        total_distance += util::coordinate_calculation::fccApproximateDistance(
             facade.GetCoordinateOfNode(*current), facade.GetCoordinateOfNode(*std::next(current)));
     }
 

include/util/coordinate_calculation.hpp

@@ -43,6 +43,9 @@ inline double radToDeg(const double radian)
 //! Takes the squared euclidean distance of the input coordinates. Does not return meters!
 std::uint64_t squaredEuclideanDistance(const Coordinate lhs, const Coordinate rhs);
 
+double fccApproximateDistance(const Coordinate first_coordinate,
+                              const Coordinate second_coordinate);
+
 double haversineDistance(const Coordinate first_coordinate, const Coordinate second_coordinate);
 
 double greatCircleDistance(const Coordinate first_coordinate, const Coordinate second_coordinate);

src/util/coordinate_calculation.cpp

@@ -5,6 +5,8 @@
 
 #include <boost/assert.hpp>
 
+#include <mapbox/cheap_ruler.hpp>
+
 #include <algorithm>
 #include <iterator>
 #include <limits>
@@ -18,6 +20,30 @@ namespace util
 namespace coordinate_calculation
 {
 
+namespace
+{
+mapbox::cheap_ruler::CheapRuler cheap_ruler_cache[] = {
+    mapbox::cheap_ruler::CheapRuler(-90, mapbox::cheap_ruler::CheapRuler::Meters),
+    mapbox::cheap_ruler::CheapRuler(-80, mapbox::cheap_ruler::CheapRuler::Meters),
+    mapbox::cheap_ruler::CheapRuler(-70, mapbox::cheap_ruler::CheapRuler::Meters),
+    mapbox::cheap_ruler::CheapRuler(-60, mapbox::cheap_ruler::CheapRuler::Meters),
+    mapbox::cheap_ruler::CheapRuler(-50, mapbox::cheap_ruler::CheapRuler::Meters),
+    mapbox::cheap_ruler::CheapRuler(-40, mapbox::cheap_ruler::CheapRuler::Meters),
+    mapbox::cheap_ruler::CheapRuler(-30, mapbox::cheap_ruler::CheapRuler::Meters),
+    mapbox::cheap_ruler::CheapRuler(-20, mapbox::cheap_ruler::CheapRuler::Meters),
+    mapbox::cheap_ruler::CheapRuler(-10, mapbox::cheap_ruler::CheapRuler::Meters),
+    mapbox::cheap_ruler::CheapRuler(0, mapbox::cheap_ruler::CheapRuler::Meters),
+    mapbox::cheap_ruler::CheapRuler(10, mapbox::cheap_ruler::CheapRuler::Meters),
+    mapbox::cheap_ruler::CheapRuler(20, mapbox::cheap_ruler::CheapRuler::Meters),
+    mapbox::cheap_ruler::CheapRuler(30, mapbox::cheap_ruler::CheapRuler::Meters),
+    mapbox::cheap_ruler::CheapRuler(40, mapbox::cheap_ruler::CheapRuler::Meters),
+    mapbox::cheap_ruler::CheapRuler(50, mapbox::cheap_ruler::CheapRuler::Meters),
+    mapbox::cheap_ruler::CheapRuler(60, mapbox::cheap_ruler::CheapRuler::Meters),
+    mapbox::cheap_ruler::CheapRuler(70, mapbox::cheap_ruler::CheapRuler::Meters),
+    mapbox::cheap_ruler::CheapRuler(80, mapbox::cheap_ruler::CheapRuler::Meters),
+    mapbox::cheap_ruler::CheapRuler(90, mapbox::cheap_ruler::CheapRuler::Meters)};
+}
+
 // Does not project the coordinates!
 std::uint64_t squaredEuclideanDistance(const Coordinate lhs, const Coordinate rhs)
 {
@@ -32,6 +58,19 @@ std::uint64_t squaredEuclideanDistance(const Coordinate lhs, const Coordinate rh
     return result;
 }
 
+// Uses method described here:
+// https://www.gpo.gov/fdsys/pkg/CFR-2005-title47-vol4/pdf/CFR-2005-title47-vol4-sec73-208.pdf
+// should be within 0.1% or so of Vincenty method (assuming 19 buckets are enough)
+// Should be more faster and more precise than Haversine
+double fccApproximateDistance(const Coordinate coordinate_1, const Coordinate coordinate_2)
+{
+    const auto lon1 = static_cast<double>(util::toFloating(coordinate_1.lon));
+    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_cache[std::lround(lat1 / 10) + 9].distance({lon1, lat1}, {lon2, lat2});
+}
+
 double haversineDistance(const Coordinate coordinate_1, const Coordinate coordinate_2)
 {
     auto lon1 = static_cast<int>(coordinate_1.lon);

third_party/cheap-ruler-cpp-2.5.4/.clang-format

@@ -0,0 +1,18 @@
+Standard: Cpp11
+IndentWidth: 4
+AccessModifierOffset: -4
+UseTab: Never
+BinPackParameters: false
+AllowShortIfStatementsOnASingleLine: false
+AllowShortLoopsOnASingleLine: false
+AllowShortBlocksOnASingleLine: false
+AllowShortFunctionsOnASingleLine: false
+ConstructorInitializerAllOnOneLineOrOnePerLine: true
+AlwaysBreakTemplateDeclarations: true
+NamespaceIndentation: None
+PointerBindsToType: true
+SpacesInParentheses: false
+BreakBeforeBraces: Attach
+ColumnLimit: 100
+Cpp11BracedListStyle: false
+SpacesBeforeTrailingComments: 1

third_party/cheap-ruler-cpp-2.5.4/.gitignore

@@ -0,0 +1,2 @@
+build/
+mason_packages/

third_party/cheap-ruler-cpp-2.5.4/.travis.yml

@@ -0,0 +1,25 @@
+language: generic
+
+matrix:
+  include:
+    - os: linux
+      env: CXX=g++-4.9
+      sudo: required
+      dist: trusty
+      addons:
+        apt:
+          sources: [ 'ubuntu-toolchain-r-test' ]
+          packages: [ 'g++-4.9', 'cmake', 'cmake-data' ]
+    - os: linux
+      env: CXX=g++-5
+      sudo: required
+      dist: trusty
+      addons:
+        apt:
+          sources: [ 'ubuntu-toolchain-r-test' ]
+          packages: [ 'g++-5', 'cmake', 'cmake-data' ]
+
+cache: apt
+
+script:
+ - cmake . && make && ./cheap_ruler

third_party/cheap-ruler-cpp-2.5.4/CMakeLists.txt

@@ -0,0 +1,27 @@
+cmake_minimum_required(VERSION 3.0.0 FATAL_ERROR)
+project(cheap_ruler LANGUAGES CXX C)
+
+include(cmake/build.cmake)
+include(cmake/mason.cmake)
+
+set(CMAKE_OSX_DEPLOYMENT_TARGET 10.10)
+set(CMAKE_CONFIGURATION_TYPES Debug Release)
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++14 -Wall -Wextra -Wpedantic -Wshadow")
+
+mason_use(geometry VERSION 0.9.2 HEADER_ONLY)
+mason_use(gtest VERSION 1.8.0)
+mason_use(variant VERSION 1.1.4 HEADER_ONLY)
+
+add_definitions(-D_GLIBCXX_USE_CXX11_ABI=0)
+
+add_executable(cheap_ruler
+    ${PROJECT_SOURCE_DIR}/test/cheap_ruler.cpp
+)
+
+target_include_directories(cheap_ruler
+    PUBLIC ${PROJECT_SOURCE_DIR}/include
+)
+
+target_add_mason_package(cheap_ruler PRIVATE geometry)
+target_add_mason_package(cheap_ruler PRIVATE gtest)
+target_add_mason_package(cheap_ruler PRIVATE variant)

third_party/cheap-ruler-cpp-2.5.4/LICENSE

@@ -0,0 +1,15 @@
+ISC License
+
+Copyright (c) 2017, Mapbox
+
+Permission to use, copy, modify, and/or distribute this software for any purpose
+with or without fee is hereby granted, provided that the above copyright notice
+and this permission notice appear in all copies.
+
+THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH
+REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
+FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT,
+INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS
+OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
+TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
+THIS SOFTWARE.

third_party/cheap-ruler-cpp-2.5.4/README.md

@@ -0,0 +1,199 @@
+# cheap-ruler-cpp
+
+Port to C++ of [Cheap Ruler](https://github.com/mapbox/cheap-ruler), a collection of very fast approximations to common geodesic measurements.
+
+[![Build Status](https://travis-ci.org/mapbox/cheap-ruler-cpp.svg?branch=master)](https://travis-ci.org/mapbox/cheap-ruler-cpp)
+
+# Usage
+
+```cpp
+#include <mapbox/cheap_ruler.hpp>
+
+namespace cr = mapbox::cheap_ruler;
+```
+
+All `point`, `line_string`, `polygon`, and `box` references are [mapbox::geometry](https://github.com/mapbox/geometry.hpp) data structures.
+
+## Create a ruler object
+
+#### `CheapRuler(double latitude, Unit unit)`
+
+Creates a ruler object that will approximate measurements around the given latitude with an optional distance unit. Once created, the ruler object has access to the [methods](#methods) below.
+
+```cpp
+auto ruler = cr::CheapRuler(32.8351);
+auto milesRuler = cr::CheapRuler(32.8351, cr::CheapRuler::Miles);
+```
+
+Possible units:
+
+* `cheap_ruler::CheapRuler::Unit`
+* `cheap_ruler::CheapRuler::Kilometers`
+* `cheap_ruler::CheapRuler::Miles`
+* `cheap_ruler::CheapRuler::NauticalMiles`
+* `cheap_ruler::CheapRuler::Meters`
+* `cheap_ruler::CheapRuler::Yards`
+* `cheap_ruler::CheapRuler::Feet`
+* `cheap_ruler::CheapRuler::Inches`
+
+#### `CheapRuler::fromTile(uint32_t y, uint32_t z)`
+
+Creates a ruler object from tile coordinates (`y` and `z` integers). Convenient in tile-reduce scripts.
+
+```cpp
+auto ruler = cr::CheapRuler::fromTile(11041, 15);
+```
+
+## Methods
+
+#### `distance(point a, point b)`
+
+Given two points of the form [x = longitude, y = latitude], returns the distance (`double`).
+
+```cpp
+cr::point point_a{-96.9148, 32.8351};
+cr::point point_b{-96.9146, 32.8386};
+auto distance = ruler.distance(point_a, point_b);
+std::clog << distance; // 0.388595
+```
+
+#### `bearing(point a, point b)`
+
+Returns the bearing (`double`) between two points in angles.
+
+```cpp
+cr::point point_a{-96.9148, 32.8351};
+cr::point point_b{-96.9146, 32.8386};
+auto bearing = ruler.bearing(point_a, point_b);
+std::clog << bearing; // 2.76206
+```
+
+#### `destination(point origin, double distance, double bearing)`
+
+Returns a new point (`point`) given distance and bearing from the starting point.
+
+```cpp
+cr::point point_a{-96.9148, 32.8351};
+auto dest = ruler.destination(point_a, 1.0, -175);
+std::clog << dest.x << ", " << dest.y; // -96.9148, 32.8261
+```
+
+#### `offset(point origin, double dx, double dy)`
+
+Returns a new point (`point`) given easting and northing offsets from the starting point.
+
+```cpp
+cr::point point_a{-96.9148, 32.8351};
+auto os = ruler.offset(point_a, 10.0, -5.0);
+std::clog << os.x << ", " << os.y; // -96.808, 32.79
+```
+
+#### `lineDistance(const line_string& points)`
+
+Given a line (an array of points), returns the total line distance (`double`).
+
+```cpp
+cr::line_string line_a{{ -96.9, 32.8 }, { -96.8, 32.8 }, { -96.2, 32.3 }};
+auto line_distance = ruler.lineDistance(line_a);
+std::clog << line_distance; // 88.2962
+```
+
+#### `area(polygon poly)`
+
+Given a polygon (an array of rings, where each ring is an array of points), returns the area (`double`).
+
+```cpp
+cr::linear_ring ring{{ -96.9, 32.8 }, { -96.8, 32.8 }, { -96.2, 32.3 }, { -96.9, 32.8 }};
+auto area = ruler.area(cr::polygon{ ring });
+std::clog << area; //
+```
+
+#### `along(const line_string& line, double distance)`
+
+Returns the point (`point`) at a specified distance along the line.
+
+```cpp
+cr::linear_ring ring{{ -96.9, 32.8 }, { -96.8, 32.8 }, { -96.2, 32.3 }, { -96.9, 32.8 }};
+auto area = ruler.area(cr::polygon{ ring });
+std::clog << area; // 259.581
+```
+
+#### `pointOnLine(const line_string& line, point p)`
+
+Returns a tuple of the form `std::pair<point, unsigned>` where point is closest point on the line from the given point, index is the start index of the segment with the closest point, and t is a parameter from 0 to 1 that indicates where the closest point is on that segment.
+
+```cpp
+cr::line_string line{{ -96.9, 32.8 }, { -96.8, 32.8 }, { -96.2, 32.3 }};
+cr::point point{-96.9, 32.79};
+auto pol = ruler.pointOnLine(line, point);
+auto point = std::get<0>(pol);
+std::clog << point.x << ", " << point.y; // -96.9, 32.8 (point)
+std::clog << std::get<1>(pol); // 0 (index)
+std::clog << std::get<2>(pol); // 0. (t)
+```
+
+#### `lineSlice(point start, point stop, const line_string& line)`
+
+Returns a part of the given line (`line_string`) between the start and the stop points (or their closest points on the line).
+
+```cpp
+cr::line_string line{{ -96.9, 32.8 }, { -96.8, 32.8 }, { -96.2, 32.3 }};
+cr::point start_point{-96.9, 32.8};
+cr::point stop_point{-96.8, 32.8};
+auto slice = ruler.lineSlice(start_point, stop_point, line);
+std::clog << slice[0].x << ", " << slice[0].y; // -96.9, 32.8
+std::clog << slice[1].x << ", " << slice[1].y; // -96.8, 32.8
+```
+
+#### `lineSliceAlong(double start, double stop, const line_string& line)`
+
+Returns a part of the given line (`line_string`) between the start and the stop points indicated by distance along the line.
+
+```cpp
+cr::line_string line{{ -96.9, 32.8 }, { -96.8, 32.8 }, { -96.2, 32.3 }};
+auto slice = ruler.lineSliceAlong(0.1, 1.2, line);
+```
+
+#### `bufferPoint(point p, double buffer)`
+
+Given a point, returns a bounding box object ([w, s, e, n]) created from the given point buffered by a given distance.
+
+```cpp
+cr::point point{-96.9, 32.8};
+auto box = ruler.bufferPoint(point, 0.1);
+```
+
+#### `bufferBBox(box bbox, double buffer)`
+
+Given a bounding box, returns the box buffered by a given distance.
+
+```cpp
+cr::box bbox({ 30, 38 }, { 40, 39 });
+auto bbox2 = ruler.bufferBBox(bbox, 1);
+```
+
+#### `insideBBox(point p, box bbox)`
+
+Returns true (`bool`) if the given point is inside in the given bounding box, otherwise false.
+
+```cpp
+cr::box bbox({ 30, 38 }, { 40, 39 });
+auto inside = ruler.insideBBox({ 35, 38.5 }, bbox);
+std::clog << inside; // true
+```
+
+# Develop
+
+```shell
+# create targets
+cmake .
+
+# build
+make
+
+# test
+./cheap_ruler
+
+# or just do it all in one!
+cmake . && make && ./cheap_ruler
+```

third_party/cheap-ruler-cpp-2.5.4/cmake/build.cmake

@@ -0,0 +1,11 @@
+# Generate source groups so the files are properly sorted in IDEs like Xcode.
+function(create_source_groups target)
+    get_target_property(sources ${target} SOURCES)
+    foreach(file ${sources})
+        get_filename_component(file "${file}" ABSOLUTE)
+        string(REGEX REPLACE "^${CMAKE_SOURCE_DIR}/" "" group "${file}")
+        get_filename_component(group "${group}" DIRECTORY)
+        string(REPLACE "/" "\\" group "${group}")
+        source_group("${group}" FILES "${file}")
+    endforeach()
+endfunction()