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