distance_haversine.hpp 9.5 KB

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