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 // Use, modification and distribution is subject to the Boost Software License,
   6 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
   7 // http://www.boost.org/LICENSE_1_0.txt)
   8
   9 #ifndef BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_HAVERSINE_HPP
  10 #define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_HAVERSINE_HPP
  11
  12
  13 #include <boost/geometry/core/cs.hpp>
  14 #include <boost/geometry/core/access.hpp>
  15 #include <boost/geometry/core/radian_access.hpp>
  16
  17 #include <boost/geometry/util/math.hpp>
  18 #include <boost/geometry/util/select_calculation_type.hpp>
  19 #include <boost/geometry/util/promote_floating_point.hpp>
  20
  21 #include <boost/geometry/strategies/distance.hpp>
  22
  23
  24
  25 namespace boost { namespace geometry
  26 {
  27
  28
  29 namespace strategy { namespace distance
  30 {
  31
  32
  33 namespace comparable
  34 {
  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 RadiusType,
  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
  62     typedef RadiusType radius_type;
  63
  64     explicit inline haversine(RadiusType const& r = 1.0)
  65         : m_radius(r)
  66     {}
  67
  68     template <typename Point1, typename Point2>
  69     static inline typename calculation_type<Point1, Point2>::type
  70     apply(Point1 const& p1, Point2 const& p2)
  71     {
  72         return calculate<typename calculation_type<Point1, Point2>::type>(
  73                    get_as_radian<0>(p1), get_as_radian<1>(p1),
  74                    get_as_radian<0>(p2), get_as_radian<1>(p2)
  75                );
  76     }
  77
  78     inline RadiusType radius() const
  79     {
  80         return m_radius;
  81     }
  82
  83
  84 private :
  85     template <typename R, typename T1, typename T2>
  86     static inline R calculate(T1 const& lon1, T1 const& lat1,
  87                               T2 const& lon2, T2 const& lat2)
  88     {
  89         return math::hav(lat2 - lat1)
  90                 + cos(lat1) * cos(lat2) * math::hav(lon2 - lon1);
  91     }
  92
  93     RadiusType m_radius;
  94 };
  95
  96
  97
  98 } // namespace comparable
  99
 100 /*!
 101 \brief Distance calculation for spherical coordinates
 102 on a perfect sphere using haversine
 103 \ingroup strategies
 104 \tparam RadiusType \tparam_radius
 105 \tparam CalculationType \tparam_calculation
 106 \author Adapted from: http://williams.best.vwh.net/avform.htm
 107 \see http://en.wikipedia.org/wiki/Great-circle_distance
 108 \note (from Wiki:) The great circle distance d between two
 109 points with coordinates {lat1,lon1} and {lat2,lon2} is given by:
 110     d=acos(sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(lon1-lon2))
 111 A mathematically equivalent formula, which is less subject
 112     to rounding error for short distances is:
 113     d=2*asin(sqrt((sin((lat1-lat2) / 2))^2
 114     + cos(lat1)*cos(lat2)*(sin((lon1-lon2) / 2))^2))
 115
 116
 117 \qbk{
 118 [heading See also]
 119 [link geometry.reference.algorithms.distance.distance_3_with_strategy distance (with strategy)]
 120 }
 121
 122 */
 123 template
 124 <
 125     typename RadiusType,
 126     typename CalculationType = void
 127 >
 128 class haversine
 129 {
 130     typedef comparable::haversine<RadiusType, CalculationType> comparable_type;
 131
 132 public :
 133     template <typename Point1, typename Point2>
 134     struct calculation_type
 135         : services::return_type<comparable_type, Point1, Point2>
 136     {};
 137
 138     typedef RadiusType radius_type;
 139
 140     /*!
 141     \brief Constructor
 142     \param radius radius of the sphere, defaults to 1.0 for the unit sphere
 143     */
 144     inline haversine(RadiusType const& radius = 1.0)
 145         : m_radius(radius)
 146     {}
 147
 148     /*!
 149     \brief applies the distance calculation
 150     \return the calculated distance (including multiplying with radius)
 151     \param p1 first point
 152     \param p2 second point
 153     */
 154     template <typename Point1, typename Point2>
 155     inline typename calculation_type<Point1, Point2>::type
 156     apply(Point1 const& p1, Point2 const& p2) const
 157     {
 158         typedef typename calculation_type<Point1, Point2>::type calculation_type;
 159         calculation_type const a = comparable_type::apply(p1, p2);
 160         calculation_type const c = calculation_type(2.0) * asin(math::sqrt(a));
 161         return calculation_type(m_radius) * c;
 162     }
 163
 164     /*!
 165     \brief access to radius value
 166     \return the radius
 167     */
 168     inline RadiusType radius() const
 169     {
 170         return m_radius;
 171     }
 172
 173 private :
 174     RadiusType m_radius;
 175 };
 176
 177
 178 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
 179 namespace services
 180 {
 181
 182 template <typename RadiusType, typename CalculationType>
 183 struct tag<haversine<RadiusType, CalculationType> >
 184 {
 185     typedef strategy_tag_distance_point_point type;
 186 };
 187
 188
 189 template <typename RadiusType, typename CalculationType, typename P1, typename P2>
 190 struct return_type<haversine<RadiusType, CalculationType>, P1, P2>
 191     : haversine<RadiusType, CalculationType>::template calculation_type<P1, P2>
 192 {};
 193
 194
 195 template <typename RadiusType, typename CalculationType>
 196 struct comparable_type<haversine<RadiusType, CalculationType> >
 197 {
 198     typedef comparable::haversine<RadiusType, CalculationType> type;
 199 };
 200
 201
 202 template <typename RadiusType, typename CalculationType>
 203 struct get_comparable<haversine<RadiusType, CalculationType> >
 204 {
 205 private :
 206     typedef haversine<RadiusType, CalculationType> this_type;
 207     typedef comparable::haversine<RadiusType, CalculationType> comparable_type;
 208 public :
 209     static inline comparable_type apply(this_type const& input)
 210     {
 211         return comparable_type(input.radius());
 212     }
 213 };
 214
 215 template <typename RadiusType, typename CalculationType, typename P1, typename P2>
 216 struct result_from_distance<haversine<RadiusType, CalculationType>, P1, P2>
 217 {
 218 private :
 219     typedef haversine<RadiusType, CalculationType> this_type;
 220     typedef typename return_type<this_type, P1, P2>::type return_type;
 221 public :
 222     template <typename T>
 223     static inline return_type apply(this_type const& , T const& value)
 224     {
 225         return return_type(value);
 226     }
 227 };
 228
 229
 230 // Specializations for comparable::haversine
 231 template <typename RadiusType, typename CalculationType>
 232 struct tag<comparable::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<comparable::haversine<RadiusType, CalculationType>, P1, P2>
 240     : comparable::haversine<RadiusType, CalculationType>::template calculation_type<P1, P2>
 241 {};
 242
 243
 244 template <typename RadiusType, typename CalculationType>
 245 struct comparable_type<comparable::haversine<RadiusType, CalculationType> >
 246 {
 247     typedef comparable::haversine<RadiusType, CalculationType> type;
 248 };
 249
 250
 251 template <typename RadiusType, typename CalculationType>
 252 struct get_comparable<comparable::haversine<RadiusType, CalculationType> >
 253 {
 254 private :
 255     typedef comparable::haversine<RadiusType, CalculationType> this_type;
 256 public :
 257     static inline this_type apply(this_type const& input)
 258     {
 259         return input;
 260     }
 261 };
 262
 263
 264 template <typename RadiusType, typename CalculationType, typename P1, typename P2>
 265 struct result_from_distance<comparable::haversine<RadiusType, CalculationType>, P1, P2>
 266 {
 267 private :
 268     typedef comparable::haversine<RadiusType, CalculationType> strategy_type;
 269     typedef typename return_type<strategy_type, P1, P2>::type return_type;
 270 public :
 271     template <typename T>
 272     static inline return_type apply(strategy_type const& strategy, T const& distance)
 273     {
 274         return_type const s = sin((distance / strategy.radius()) / return_type(2));
 275         return s * s;
 276     }
 277 };
 278
 279
 280 // Register it as the default for point-types
 281 // in a spherical equatorial coordinate system
 282 template <typename Point1, typename Point2>
 283 struct default_strategy
 284     <
 285         point_tag, point_tag, Point1, Point2,
 286         spherical_equatorial_tag, spherical_equatorial_tag
 287     >
 288 {
 289     typedef strategy::distance::haversine<typename select_coordinate_type<Point1, Point2>::type> type;
 290 };
 291
 292 // Note: spherical polar coordinate system requires "get_as_radian_equatorial"
 293
 294
 295 } // namespace services
 296 #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
 297
 298
 299 }} // namespace strategy::distance
 300
 301
 302 }} // namespace boost::geometry
 303
 304
 305 #endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_HAVERSINE_HPP