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