distance_haversine.hpp 9.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335
  1. // Boost.Geometry (aka GGL, Generic Geometry Library)
  2. // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
  3. // This file was modified by Oracle on 2017, 2018.
  4. // Modifications copyright (c) 2017-2018, Oracle and/or its affiliates.
  5. // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
  6. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
  7. // Use, modification and distribution is subject to the Boost Software License,
  8. // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  9. // http://www.boost.org/LICENSE_1_0.txt)
  10. #ifndef BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_HAVERSINE_HPP
  11. #define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_HAVERSINE_HPP
  12. #include <boost/geometry/core/access.hpp>
  13. #include <boost/geometry/core/cs.hpp>
  14. #include <boost/geometry/core/radian_access.hpp>
  15. #include <boost/geometry/srs/sphere.hpp>
  16. #include <boost/geometry/strategies/distance.hpp>
  17. #include <boost/geometry/strategies/spherical/get_radius.hpp>
  18. #include <boost/geometry/util/math.hpp>
  19. #include <boost/geometry/util/promote_floating_point.hpp>
  20. #include <boost/geometry/util/select_calculation_type.hpp>
  21. namespace boost { namespace geometry
  22. {
  23. namespace strategy { namespace distance
  24. {
  25. namespace comparable
  26. {
  27. // Comparable haversine.
  28. // To compare distances, we can avoid:
  29. // - multiplication with radius and 2.0
  30. // - applying sqrt
  31. // - applying asin (which is strictly (monotone) increasing)
  32. template
  33. <
  34. typename RadiusTypeOrSphere = double,
  35. typename CalculationType = void
  36. >
  37. class haversine
  38. {
  39. public :
  40. template <typename Point1, typename Point2>
  41. struct calculation_type
  42. : promote_floating_point
  43. <
  44. typename select_calculation_type
  45. <
  46. Point1,
  47. Point2,
  48. CalculationType
  49. >::type
  50. >
  51. {};
  52. typedef typename strategy_detail::get_radius
  53. <
  54. RadiusTypeOrSphere
  55. >::type radius_type;
  56. inline haversine()
  57. : m_radius(1.0)
  58. {}
  59. template <typename RadiusOrSphere>
  60. explicit inline haversine(RadiusOrSphere const& radius_or_sphere)
  61. : m_radius(strategy_detail::get_radius
  62. <
  63. RadiusOrSphere
  64. >::apply(radius_or_sphere))
  65. {}
  66. template <typename Point1, typename Point2>
  67. static inline typename calculation_type<Point1, Point2>::type
  68. apply(Point1 const& p1, Point2 const& p2)
  69. {
  70. return calculate<typename calculation_type<Point1, Point2>::type>(
  71. get_as_radian<0>(p1), get_as_radian<1>(p1),
  72. get_as_radian<0>(p2), get_as_radian<1>(p2)
  73. );
  74. }
  75. inline radius_type radius() const
  76. {
  77. return m_radius;
  78. }
  79. private :
  80. template <typename R, typename T1, typename T2>
  81. static inline R calculate(T1 const& lon1, T1 const& lat1,
  82. T2 const& lon2, T2 const& lat2)
  83. {
  84. return math::hav(lat2 - lat1)
  85. + cos(lat1) * cos(lat2) * math::hav(lon2 - lon1);
  86. }
  87. radius_type m_radius;
  88. };
  89. } // namespace comparable
  90. /*!
  91. \brief Distance calculation for spherical coordinates
  92. on a perfect sphere using haversine
  93. \ingroup strategies
  94. \tparam RadiusTypeOrSphere \tparam_radius_or_sphere
  95. \tparam CalculationType \tparam_calculation
  96. \author Adapted from: http://williams.best.vwh.net/avform.htm
  97. \see http://en.wikipedia.org/wiki/Great-circle_distance
  98. \note (from Wiki:) The great circle distance d between two
  99. points with coordinates {lat1,lon1} and {lat2,lon2} is given by:
  100. d=acos(sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(lon1-lon2))
  101. A mathematically equivalent formula, which is less subject
  102. to rounding error for short distances is:
  103. d=2*asin(sqrt((sin((lat1-lat2) / 2))^2
  104. + cos(lat1)*cos(lat2)*(sin((lon1-lon2) / 2))^2))
  105. \qbk{
  106. [heading See also]
  107. [link geometry.reference.algorithms.distance.distance_3_with_strategy distance (with strategy)]
  108. }
  109. */
  110. template
  111. <
  112. typename RadiusTypeOrSphere = double,
  113. typename CalculationType = void
  114. >
  115. class haversine
  116. {
  117. typedef comparable::haversine<RadiusTypeOrSphere, CalculationType> comparable_type;
  118. public :
  119. template <typename Point1, typename Point2>
  120. struct calculation_type
  121. : services::return_type<comparable_type, Point1, Point2>
  122. {};
  123. typedef typename strategy_detail::get_radius
  124. <
  125. RadiusTypeOrSphere
  126. >::type radius_type;
  127. /*!
  128. \brief Default constructor, radius set to 1.0 for the unit sphere
  129. */
  130. inline haversine()
  131. : m_radius(1.0)
  132. {}
  133. /*!
  134. \brief Constructor
  135. \param radius_or_sphere radius of the sphere or sphere model
  136. */
  137. template <typename RadiusOrSphere>
  138. explicit inline haversine(RadiusOrSphere const& radius_or_sphere)
  139. : m_radius(strategy_detail::get_radius
  140. <
  141. RadiusOrSphere
  142. >::apply(radius_or_sphere))
  143. {}
  144. /*!
  145. \brief applies the distance calculation
  146. \return the calculated distance (including multiplying with radius)
  147. \param p1 first point
  148. \param p2 second point
  149. */
  150. template <typename Point1, typename Point2>
  151. inline typename calculation_type<Point1, Point2>::type
  152. apply(Point1 const& p1, Point2 const& p2) const
  153. {
  154. typedef typename calculation_type<Point1, Point2>::type calculation_type;
  155. calculation_type const a = comparable_type::apply(p1, p2);
  156. calculation_type const c = calculation_type(2.0) * asin(math::sqrt(a));
  157. return calculation_type(m_radius) * c;
  158. }
  159. /*!
  160. \brief access to radius value
  161. \return the radius
  162. */
  163. inline radius_type radius() const
  164. {
  165. return m_radius;
  166. }
  167. private :
  168. radius_type m_radius;
  169. };
  170. #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
  171. namespace services
  172. {
  173. template <typename RadiusType, typename CalculationType>
  174. struct tag<haversine<RadiusType, CalculationType> >
  175. {
  176. typedef strategy_tag_distance_point_point type;
  177. };
  178. template <typename RadiusType, typename CalculationType, typename P1, typename P2>
  179. struct return_type<haversine<RadiusType, CalculationType>, P1, P2>
  180. : haversine<RadiusType, CalculationType>::template calculation_type<P1, P2>
  181. {};
  182. template <typename RadiusType, typename CalculationType>
  183. struct comparable_type<haversine<RadiusType, CalculationType> >
  184. {
  185. typedef comparable::haversine<RadiusType, CalculationType> type;
  186. };
  187. template <typename RadiusType, typename CalculationType>
  188. struct get_comparable<haversine<RadiusType, CalculationType> >
  189. {
  190. private :
  191. typedef haversine<RadiusType, CalculationType> this_type;
  192. typedef comparable::haversine<RadiusType, CalculationType> comparable_type;
  193. public :
  194. static inline comparable_type apply(this_type const& input)
  195. {
  196. return comparable_type(input.radius());
  197. }
  198. };
  199. template <typename RadiusType, typename CalculationType, typename P1, typename P2>
  200. struct result_from_distance<haversine<RadiusType, CalculationType>, P1, P2>
  201. {
  202. private :
  203. typedef haversine<RadiusType, CalculationType> this_type;
  204. typedef typename return_type<this_type, P1, P2>::type return_type;
  205. public :
  206. template <typename T>
  207. static inline return_type apply(this_type const& , T const& value)
  208. {
  209. return return_type(value);
  210. }
  211. };
  212. // Specializations for comparable::haversine
  213. template <typename RadiusType, typename CalculationType>
  214. struct tag<comparable::haversine<RadiusType, CalculationType> >
  215. {
  216. typedef strategy_tag_distance_point_point type;
  217. };
  218. template <typename RadiusType, typename CalculationType, typename P1, typename P2>
  219. struct return_type<comparable::haversine<RadiusType, CalculationType>, P1, P2>
  220. : comparable::haversine<RadiusType, CalculationType>::template calculation_type<P1, P2>
  221. {};
  222. template <typename RadiusType, typename CalculationType>
  223. struct comparable_type<comparable::haversine<RadiusType, CalculationType> >
  224. {
  225. typedef comparable::haversine<RadiusType, CalculationType> type;
  226. };
  227. template <typename RadiusType, typename CalculationType>
  228. struct get_comparable<comparable::haversine<RadiusType, CalculationType> >
  229. {
  230. private :
  231. typedef comparable::haversine<RadiusType, CalculationType> this_type;
  232. public :
  233. static inline this_type apply(this_type const& input)
  234. {
  235. return input;
  236. }
  237. };
  238. template <typename RadiusType, typename CalculationType, typename P1, typename P2>
  239. struct result_from_distance<comparable::haversine<RadiusType, CalculationType>, P1, P2>
  240. {
  241. private :
  242. typedef comparable::haversine<RadiusType, CalculationType> strategy_type;
  243. typedef typename return_type<strategy_type, P1, P2>::type return_type;
  244. public :
  245. template <typename T>
  246. static inline return_type apply(strategy_type const& strategy, T const& distance)
  247. {
  248. return_type const s = sin((distance / strategy.radius()) / return_type(2));
  249. return s * s;
  250. }
  251. };
  252. // Register it as the default for point-types
  253. // in a spherical equatorial coordinate system
  254. template <typename Point1, typename Point2>
  255. struct default_strategy
  256. <
  257. point_tag, point_tag, Point1, Point2,
  258. spherical_equatorial_tag, spherical_equatorial_tag
  259. >
  260. {
  261. typedef strategy::distance::haversine<typename select_coordinate_type<Point1, Point2>::type> type;
  262. };
  263. // Note: spherical polar coordinate system requires "get_as_radian_equatorial"
  264. } // namespace services
  265. #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
  266. }} // namespace strategy::distance
  267. }} // namespace boost::geometry
  268. #endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_HAVERSINE_HPP