ceph/src/boost/boost/geometry/strategies/spherical/distance_haversine.hpp

   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
  75     inline haversine()
  76         : m_radius(1.0)
  77     {}
  78
  79     template <typename RadiusOrSphere>
  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
  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
  97     template <typename T1, typename T2>
  98     inline radius_type meridian(T1 lat1, T2 lat2) const
  99     {
 100         return m_radius * (lat1 - lat2);
 101     }
 102
 103     inline radius_type radius() const
 104     {
 105         return m_radius;
 106     }
 107
 108
 109 private :
 110     template <typename R, typename T1, typename T2>
 111     static inline R calculate(T1 const& lon1, T1 const& lat1,
 112                               T2 const& lon2, T2 const& lat2)
 113     {
 114         return math::hav(lat2 - lat1)
 115                 + cos(lat1) * cos(lat2) * math::hav(lon2 - lon1);
 116     }
 117
 118     radius_type m_radius;
 119 };
 120
 121
 122
 123 } // namespace comparable
 124
 125 /*!
 126 \brief Distance calculation for spherical coordinates
 127 on a perfect sphere using haversine
 128 \ingroup strategies
 129 \tparam RadiusTypeOrSphere \tparam_radius_or_sphere
 130 \tparam CalculationType \tparam_calculation
 131 \author Adapted from: http://williams.best.vwh.net/avform.htm
 132 \see http://en.wikipedia.org/wiki/Great-circle_distance
 133 \note (from Wiki:) The great circle distance d between two
 134 points with coordinates {lat1,lon1} and {lat2,lon2} is given by:
 135     d=acos(sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(lon1-lon2))
 136 A mathematically equivalent formula, which is less subject
 137     to rounding error for short distances is:
 138     d=2*asin(sqrt((sin((lat1-lat2) / 2))^2
 139     + cos(lat1)*cos(lat2)*(sin((lon1-lon2) / 2))^2))
 140 \qbk{
 141 [heading See also]
 142 [link geometry.reference.algorithms.distance.distance_3_with_strategy distance (with strategy)]
 143 }
 144 */
 145 template
 146 <
 147     typename RadiusTypeOrSphere = double,
 148     typename CalculationType = void
 149 >
 150 class haversine
 151 {
 152     typedef comparable::haversine<RadiusTypeOrSphere, CalculationType> comparable_type;
 153
 154 public :
 155     template <typename Point1, typename Point2>
 156     struct calculation_type
 157         : services::return_type<comparable_type, Point1, Point2>
 158     {};
 159
 160     typedef typename strategy_detail::get_radius
 161         <
 162             RadiusTypeOrSphere
 163         >::type radius_type;
 164
 165     /*!
 166     \brief Default constructor, radius set to 1.0 for the unit sphere
 167     */
 168     inline haversine()
 169         : m_radius(1.0)
 170     {}
 171
 172     /*!
 173     \brief Constructor
 174     \param radius_or_sphere radius of the sphere or sphere model
 175     */
 176     template <typename RadiusOrSphere>
 177     explicit inline haversine(RadiusOrSphere const& radius_or_sphere)
 178         : m_radius(strategy_detail::get_radius
 179                     <
 180                         RadiusOrSphere
 181                     >::apply(radius_or_sphere))
 182     {}
 183
 184     /*!
 185     \brief applies the distance calculation
 186     \return the calculated distance (including multiplying with radius)
 187     \param p1 first point
 188     \param p2 second point
 189     */
 190     template <typename Point1, typename Point2>
 191     inline typename calculation_type<Point1, Point2>::type
 192     apply(Point1 const& p1, Point2 const& p2) const
 193     {
 194         typedef typename calculation_type<Point1, Point2>::type calculation_type;
 195         calculation_type const a = comparable_type::apply(p1, p2);
 196         calculation_type const c = calculation_type(2.0) * asin(math::sqrt(a));
 197         return calculation_type(m_radius) * c;
 198     }
 199
 200     /*!
 201     \brief meridian distance calculation
 202     \return the calculated distance (including multiplying with radius)
 203     \param p1 first point
 204     \param p2 second point
 205     */
 206
 207     template <typename T1, typename T2>
 208     inline radius_type meridian(T1 lat1, T2 lat2) const
 209     {
 210         return m_radius * (lat1 - lat2);
 211     }
 212
 213     /*!
 214     \brief access to radius value
 215     \return the radius
 216     */
 217     inline radius_type radius() const
 218     {
 219         return m_radius;
 220     }
 221
 222 private :
 223     radius_type m_radius;
 224 };
 225
 226
 227 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
 228 namespace services
 229 {
 230
 231 template <typename RadiusType, typename CalculationType>
 232 struct tag<haversine<RadiusType, CalculationType> >
 233 {
 234     typedef strategy_tag_distance_point_point type;
 235 };
 236
 237
 238 template <typename RadiusType, typename CalculationType, typename P1, typename P2>
 239 struct return_type<haversine<RadiusType, CalculationType>, P1, P2>
 240     : haversine<RadiusType, CalculationType>::template calculation_type<P1, P2>
 241 {};
 242
 243
 244 template <typename RadiusType, typename CalculationType>
 245 struct comparable_type<haversine<RadiusType, CalculationType> >
 246 {
 247     typedef comparable::haversine<RadiusType, CalculationType> type;
 248 };
 249
 250
 251 template <typename RadiusType, typename CalculationType>
 252 struct get_comparable<haversine<RadiusType, CalculationType> >
 253 {
 254 private :
 255     typedef haversine<RadiusType, CalculationType> this_type;
 256     typedef comparable::haversine<RadiusType, CalculationType> comparable_type;
 257 public :
 258     static inline comparable_type apply(this_type const& input)
 259     {
 260         return comparable_type(input.radius());
 261     }
 262 };
 263
 264 template <typename RadiusType, typename CalculationType, typename P1, typename P2>
 265 struct result_from_distance<haversine<RadiusType, CalculationType>, P1, P2>
 266 {
 267 private :
 268     typedef haversine<RadiusType, CalculationType> this_type;
 269     typedef typename return_type<this_type, P1, P2>::type return_type;
 270 public :
 271     template <typename T>
 272     static inline return_type apply(this_type const& , T const& value)
 273     {
 274         return return_type(value);
 275     }
 276 };
 277
 278
 279 // Specializations for comparable::haversine
 280 template <typename RadiusType, typename CalculationType>
 281 struct tag<comparable::haversine<RadiusType, CalculationType> >
 282 {
 283     typedef strategy_tag_distance_point_point type;
 284 };
 285
 286
 287 template <typename RadiusType, typename CalculationType, typename P1, typename P2>
 288 struct return_type<comparable::haversine<RadiusType, CalculationType>, P1, P2>
 289     : comparable::haversine<RadiusType, CalculationType>::template calculation_type<P1, P2>
 290 {};
 291
 292
 293 template <typename RadiusType, typename CalculationType>
 294 struct comparable_type<comparable::haversine<RadiusType, CalculationType> >
 295 {
 296     typedef comparable::haversine<RadiusType, CalculationType> type;
 297 };
 298
 299
 300 template <typename RadiusType, typename CalculationType>
 301 struct get_comparable<comparable::haversine<RadiusType, CalculationType> >
 302 {
 303 private :
 304     typedef comparable::haversine<RadiusType, CalculationType> this_type;
 305 public :
 306     static inline this_type apply(this_type const& input)
 307     {
 308         return input;
 309     }
 310 };
 311
 312
 313 template <typename RadiusType, typename CalculationType, typename P1, typename P2>
 314 struct result_from_distance<comparable::haversine<RadiusType, CalculationType>, P1, P2>
 315 {
 316 private :
 317     typedef comparable::haversine<RadiusType, CalculationType> strategy_type;
 318     typedef typename return_type<strategy_type, P1, P2>::type return_type;
 319 public :
 320     template <typename T>
 321     static inline return_type apply(strategy_type const& strategy, T const& distance)
 322     {
 323         return_type const s = sin((distance / strategy.radius()) / return_type(2));
 324         return s * s;
 325     }
 326 };
 327
 328
 329 // Register it as the default for point-types
 330 // in a spherical equatorial coordinate system
 331 template <typename Point1, typename Point2>
 332 struct default_strategy
 333     <
 334         point_tag, point_tag, Point1, Point2,
 335         spherical_equatorial_tag, spherical_equatorial_tag
 336     >
 337 {
 338     typedef strategy::distance::haversine<typename select_coordinate_type<Point1, Point2>::type> type;
 339 };
 340
 341 // Note: spherical polar coordinate system requires "get_as_radian_equatorial"
 342
 343
 344 } // namespace services
 345 #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
 346
 347
 348 }} // namespace strategy::distance
 349
 350
 351 }} // namespace boost::geometry
 352
 353
 354 #endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_HAVERSINE_HPP