1 // Boost.Geometry (aka GGL, Generic Geometry Library)
2 // Unit Test
3
4 // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
5 // Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
6 // Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
7
8 // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
9 // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
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
16 #include <geometry_test_common.hpp>
17
18 #include <boost/concept/requires.hpp>
19 #include <boost/concept_check.hpp>
20 #include <boost/core/ignore_unused.hpp>
21
22 #include <boost/geometry/algorithms/assign.hpp>
23 #include <boost/geometry/strategies/spherical/distance_haversine.hpp>
24 #include <boost/geometry/strategies/concepts/distance_concept.hpp>
25
26
27 #include <boost/geometry/geometries/point.hpp>
28
29 #ifdef HAVE_TTMATH
30 # include <boost/geometry/extensions/contrib/ttmath_stub.hpp>
31 #endif
32
33
34
35 double const average_earth_radius = 6372795.0;
36
37
38 template <typename Point, typename LatitudePolicy>
39 struct test_distance
40 {
41 typedef bg::strategy::distance::haversine<double> haversine_type;
42 typedef typename bg::strategy::distance::services::return_type<haversine_type, Point, Point>::type return_type;
43
44 BOOST_CONCEPT_ASSERT
45 (
46 (bg::concepts::PointDistanceStrategy<haversine_type, Point, Point>)
47 );
48
49
testtest_distance50 static void test(double lon1, double lat1, double lon2, double lat2,
51 double radius, double expected, double tolerance)
52 {
53 haversine_type strategy(radius);
54
55 Point p1, p2;
56 bg::assign_values(p1, lon1, LatitudePolicy::apply(lat1));
57 bg::assign_values(p2, lon2, LatitudePolicy::apply(lat2));
58 return_type d = strategy.apply(p1, p2);
59
60 BOOST_CHECK_CLOSE(d, expected, tolerance);
61 }
62 };
63
64 template <typename Point, typename LatitudePolicy>
test_all()65 void test_all()
66 {
67 // earth to unit-sphere -> divide by earth circumference, then it is from 0-1,
68 // then multiply with 2 PI, so effectively just divide by earth radius
69 double e2u = 1.0 / average_earth_radius;
70
71 // ~ Amsterdam/Paris, 467 kilometers
72 double const a_p = 467.2704 * 1000.0;
73 test_distance<Point, LatitudePolicy>::test(4, 52, 2, 48, average_earth_radius, a_p, 1.0);
74 test_distance<Point, LatitudePolicy>::test(2, 48, 4, 52, average_earth_radius, a_p, 1.0);
75 test_distance<Point, LatitudePolicy>::test(4, 52, 2, 48, 1.0, a_p * e2u, 0.001);
76
77 // ~ Amsterdam/Barcelona
78 double const a_b = 1232.9065 * 1000.0;
79 test_distance<Point, LatitudePolicy>::test(4, 52, 2, 41, average_earth_radius, a_b, 1.0);
80 test_distance<Point, LatitudePolicy>::test(2, 41, 4, 52, average_earth_radius, a_b, 1.0);
81 test_distance<Point, LatitudePolicy>::test(4, 52, 2, 41, 1.0, a_b * e2u, 0.001);
82 }
83
84
85 template <typename P1, typename P2, typename CalculationType, typename LatitudePolicy>
test_services()86 void test_services()
87 {
88 namespace bgsd = bg::strategy::distance;
89 namespace services = bg::strategy::distance::services;
90
91 {
92
93 // Compile-check if there is a strategy for this type
94 typedef typename services::default_strategy
95 <
96 bg::point_tag, bg::point_tag, P1, P2
97 >::type haversine_strategy_type;
98
99 boost::ignore_unused<haversine_strategy_type>();
100 }
101
102 P1 p1;
103 bg::assign_values(p1, 4, 52);
104
105 P2 p2;
106 bg::assign_values(p2, 2, 48);
107
108 // ~ Amsterdam/Paris, 467 kilometers
109 double const km = 1000.0;
110 double const a_p = 467.2704 * km;
111 double const expected = a_p;
112
113 double const expected_lower = 460.0 * km;
114 double const expected_higher = 470.0 * km;
115
116 // 1: normal, calculate distance:
117
118 typedef bgsd::haversine<double, CalculationType> strategy_type;
119 typedef typename bgsd::services::return_type<strategy_type, P1, P2>::type return_type;
120
121 strategy_type strategy(average_earth_radius);
122 return_type result = strategy.apply(p1, p2);
123 BOOST_CHECK_CLOSE(result, return_type(expected), 0.001);
124
125 // 2: the strategy should return the same result if we reverse parameters
126 result = strategy.apply(p2, p1);
127 BOOST_CHECK_CLOSE(result, return_type(expected), 0.001);
128
129
130 // 3: "comparable" to construct a "comparable strategy" for P1/P2
131 // a "comparable strategy" is a strategy which does not calculate the exact distance, but
132 // which returns results which can be mutually compared (e.g. avoid sqrt)
133
134 // 3a: "comparable_type"
135 typedef typename services::comparable_type<strategy_type>::type comparable_type;
136
137 // 3b: "get_comparable"
138 comparable_type comparable = bgsd::services::get_comparable<strategy_type>::apply(strategy);
139
140 // Check vice versa:
141 // First the result of the comparable strategy
142 return_type c_result = comparable.apply(p1, p2);
143 // Second the comparable result of the expected distance
144 return_type c_expected = services::result_from_distance<comparable_type, P1, P2>::apply(comparable, expected);
145 // And that one should be equa.
146 BOOST_CHECK_CLOSE(c_result, return_type(c_expected), 0.001);
147
148 // 4: the comparable_type should have a distance_strategy_constructor as well,
149 // knowing how to compare something with a fixed distance
150 return_type c_dist_lower = services::result_from_distance<comparable_type, P1, P2>::apply(comparable, expected_lower);
151 return_type c_dist_higher = services::result_from_distance<comparable_type, P1, P2>::apply(comparable, expected_higher);
152
153 // If this is the case:
154 BOOST_CHECK(c_dist_lower < c_result && c_result < c_dist_higher);
155
156 // Calculate the Haversine by hand here:
157 return_type c_check = return_type(2.0) * asin(sqrt(c_result)) * average_earth_radius;
158 BOOST_CHECK_CLOSE(c_check, expected, 0.001);
159
160 // This should also be the case
161 return_type dist_lower = services::result_from_distance<strategy_type, P1, P2>::apply(strategy, expected_lower);
162 return_type dist_higher = services::result_from_distance<strategy_type, P1, P2>::apply(strategy, expected_higher);
163 BOOST_CHECK(dist_lower < result && result < dist_higher);
164 }
165
166 /****
167 template <typename P, typename Strategy>
168 void time_compare_s(int const n)
169 {
170 boost::timer t;
171 P p1, p2;
172 bg::assign_values(p1, 1, 1);
173 bg::assign_values(p2, 2, 2);
174 Strategy strategy;
175 typename Strategy::return_type s = 0;
176 for (int i = 0; i < n; i++)
177 {
178 for (int j = 0; j < n; j++)
179 {
180 s += strategy.apply(p1, p2);
181 }
182 }
183 std::cout << "s: " << s << " t: " << t.elapsed() << std::endl;
184 }
185
186 template <typename P>
187 void time_compare(int const n)
188 {
189 time_compare_s<P, bg::strategy::distance::haversine<double> >(n);
190 time_compare_s<P, bg::strategy::distance::comparable::haversine<double> >(n);
191 }
192
193 #include <time.h>
194 double time_sqrt(int n)
195 {
196 clock_t start = clock();
197
198 double v = 2.0;
199 double s = 0.0;
200 for (int i = 0; i < n; i++)
201 {
202 for (int j = 0; j < n; j++)
203 {
204 s += sqrt(v);
205 v += 1.0e-10;
206 }
207 }
208 clock_t end = clock();
209 double diff = double(end - start) / CLOCKS_PER_SEC;
210
211 std::cout << "Check: " << s << " Time: " << diff << std::endl;
212 return diff;
213 }
214
215 double time_normal(int n)
216 {
217 clock_t start = clock();
218
219 double v = 2.0;
220 double s = 0.0;
221 for (int i = 0; i < n; i++)
222 {
223 for (int j = 0; j < n; j++)
224 {
225 s += v;
226 v += 1.0e-10;
227 }
228 }
229 clock_t end = clock();
230 double diff = double(end - start) / CLOCKS_PER_SEC;
231
232 std::cout << "Check: " << s << " Time: " << diff << std::endl;
233 return diff;
234 }
235 ***/
236
test_main(int,char * [])237 int test_main(int, char* [])
238 {
239 test_all<bg::model::point<int, 2, bg::cs::spherical_equatorial<bg::degree> >, geographic_policy>();
240 test_all<bg::model::point<float, 2, bg::cs::spherical_equatorial<bg::degree> >, geographic_policy>();
241 test_all<bg::model::point<double, 2, bg::cs::spherical_equatorial<bg::degree> >, geographic_policy>();
242
243 // NYI: haversine for mathematical spherical coordinate systems
244 // test_all<bg::model::point<double, 2, bg::cs::spherical<bg::degree> >, mathematical_policy>();
245
246 //double t1 = time_sqrt(20000);
247 //double t2 = time_normal(20000);
248 //std::cout << "Factor: " << (t1 / t2) << std::endl;
249 //time_compare<bg::model::point<double, 2, bg::cs::spherical<bg::radian> > >(10000);
250
251 #if defined(HAVE_TTMATH)
252 typedef ttmath::Big<1,4> tt;
253 test_all<bg::model::point<tt, 2, bg::cs::spherical_equatorial<bg::degree> >, geographic_policy>();
254 #endif
255
256
257 test_services
258 <
259 bg::model::point<double, 2, bg::cs::spherical_equatorial<bg::degree> >,
260 bg::model::point<double, 2, bg::cs::spherical_equatorial<bg::degree> >,
261 double,
262 geographic_policy
263 >();
264
265 return 0;
266 }
267