]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // Boost.Geometry (aka GGL, Generic Geometry Library) |
2 | ||
3 | // Copyright (c) 2007-2014 Barend Gehrels, Amsterdam, the Netherlands. | |
4 | ||
5 | // This file was modified by Oracle on 2014. | |
6 | // Modifications copyright (c) 2014, Oracle and/or its affiliates. | |
7 | ||
8 | // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle | |
9 | ||
10 | // Use, modification and distribution is subject to the Boost Software License, | |
11 | // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at | |
12 | // http://www.boost.org/LICENSE_1_0.txt) | |
13 | ||
14 | #ifndef BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_CROSS_TRACK_HPP | |
15 | #define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_CROSS_TRACK_HPP | |
16 | ||
17 | #include <algorithm> | |
18 | ||
19 | #include <boost/config.hpp> | |
20 | #include <boost/concept_check.hpp> | |
21 | #include <boost/mpl/if.hpp> | |
22 | #include <boost/type_traits/is_void.hpp> | |
23 | ||
24 | #include <boost/geometry/core/cs.hpp> | |
25 | #include <boost/geometry/core/access.hpp> | |
26 | #include <boost/geometry/core/radian_access.hpp> | |
27 | #include <boost/geometry/core/tags.hpp> | |
28 | ||
29 | #include <boost/geometry/algorithms/detail/course.hpp> | |
30 | ||
31 | #include <boost/geometry/strategies/distance.hpp> | |
32 | #include <boost/geometry/strategies/concepts/distance_concept.hpp> | |
33 | #include <boost/geometry/strategies/spherical/distance_haversine.hpp> | |
34 | ||
35 | #include <boost/geometry/util/math.hpp> | |
36 | #include <boost/geometry/util/promote_floating_point.hpp> | |
37 | #include <boost/geometry/util/select_calculation_type.hpp> | |
38 | ||
39 | #ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK | |
40 | # include <boost/geometry/io/dsv/write.hpp> | |
41 | #endif | |
42 | ||
43 | ||
44 | namespace boost { namespace geometry | |
45 | { | |
46 | ||
47 | namespace strategy { namespace distance | |
48 | { | |
49 | ||
50 | ||
51 | namespace comparable | |
52 | { | |
53 | ||
54 | /* | |
55 | Given a spherical segment AB and a point D, we are interested in | |
56 | computing the distance of D from AB. This is usually known as the | |
57 | cross track distance. | |
58 | ||
59 | If the projection (along great circles) of the point D lies inside | |
60 | the segment AB, then the distance (cross track error) XTD is given | |
61 | by the formula (see http://williams.best.vwh.net/avform.htm#XTE): | |
62 | ||
63 | XTD = asin( sin(dist_AD) * sin(crs_AD-crs_AB) ) | |
64 | ||
65 | where dist_AD is the great circle distance between the points A and | |
66 | B, and crs_AD, crs_AB is the course (bearing) between the points A, | |
67 | D and A, B, respectively. | |
68 | ||
69 | If the point D does not project inside the arc AB, then the distance | |
70 | of D from AB is the minimum of the two distances dist_AD and dist_BD. | |
71 | ||
72 | Our reference implementation for this procedure is listed below | |
73 | (this was the old Boost.Geometry implementation of the cross track distance), | |
74 | where: | |
75 | * The member variable m_strategy is the underlying haversine strategy. | |
76 | * p stands for the point D. | |
77 | * sp1 stands for the segment endpoint A. | |
78 | * sp2 stands for the segment endpoint B. | |
79 | ||
80 | ================= reference implementation -- start ================= | |
81 | ||
82 | return_type d1 = m_strategy.apply(sp1, p); | |
83 | return_type d3 = m_strategy.apply(sp1, sp2); | |
84 | ||
85 | if (geometry::math::equals(d3, 0.0)) | |
86 | { | |
87 | // "Degenerate" segment, return either d1 or d2 | |
88 | return d1; | |
89 | } | |
90 | ||
91 | return_type d2 = m_strategy.apply(sp2, p); | |
92 | ||
93 | return_type crs_AD = geometry::detail::course<return_type>(sp1, p); | |
94 | return_type crs_AB = geometry::detail::course<return_type>(sp1, sp2); | |
95 | return_type crs_BA = crs_AB - geometry::math::pi<return_type>(); | |
96 | return_type crs_BD = geometry::detail::course<return_type>(sp2, p); | |
97 | return_type d_crs1 = crs_AD - crs_AB; | |
98 | return_type d_crs2 = crs_BD - crs_BA; | |
99 | ||
100 | // d1, d2, d3 are in principle not needed, only the sign matters | |
101 | return_type projection1 = cos( d_crs1 ) * d1 / d3; | |
102 | return_type projection2 = cos( d_crs2 ) * d2 / d3; | |
103 | ||
104 | if (projection1 > 0.0 && projection2 > 0.0) | |
105 | { | |
106 | return_type XTD | |
107 | = radius() * math::abs( asin( sin( d1 / radius() ) * sin( d_crs1 ) )); | |
108 | ||
109 | // Return shortest distance, projected point on segment sp1-sp2 | |
110 | return return_type(XTD); | |
111 | } | |
112 | else | |
113 | { | |
114 | // Return shortest distance, project either on point sp1 or sp2 | |
115 | return return_type( (std::min)( d1 , d2 ) ); | |
116 | } | |
117 | ||
118 | ================= reference implementation -- end ================= | |
119 | ||
120 | ||
121 | Motivation | |
122 | ---------- | |
123 | In what follows we develop a comparable version of the cross track | |
124 | distance strategy, that meets the following goals: | |
125 | * It is more efficient than the original cross track strategy (less | |
126 | operations and less calls to mathematical functions). | |
127 | * Distances using the comparable cross track strategy can not only | |
128 | be compared with other distances using the same strategy, but also with | |
129 | distances computed with the comparable version of the haversine strategy. | |
130 | * It can serve as the basis for the computation of the cross track distance, | |
131 | as it is more efficient to compute its comparable version and | |
132 | transform that to the actual cross track distance, rather than | |
133 | follow/use the reference implementation listed above. | |
134 | ||
135 | Major idea | |
136 | ---------- | |
137 | The idea here is to use the comparable haversine strategy to compute | |
138 | the distances d1, d2 and d3 in the above listing. Once we have done | |
139 | that we need also to make sure that instead of returning XTD (as | |
140 | computed above) that we return a distance CXTD that is compatible | |
141 | with the comparable haversine distance. To achieve this CXTD must satisfy | |
142 | the relation: | |
143 | XTD = 2 * R * asin( sqrt(XTD) ) | |
144 | where R is the sphere's radius. | |
145 | ||
146 | Below we perform the mathematical analysis that show how to compute CXTD. | |
147 | ||
148 | ||
149 | Mathematical analysis | |
150 | --------------------- | |
151 | Below we use the following trigonometric identities: | |
152 | sin(2 * x) = 2 * sin(x) * cos(x) | |
153 | cos(asin(x)) = sqrt(1 - x^2) | |
154 | ||
155 | Observation: | |
156 | The distance d1 needed when the projection of the point D is within the | |
157 | segment must be the true distance. However, comparable::haversine<> | |
158 | returns a comparable distance instead of the one needed. | |
159 | To remedy this, we implicitly compute what is needed. | |
160 | More precisely, we need to compute sin(true_d1): | |
161 | ||
162 | sin(true_d1) = sin(2 * asin(sqrt(d1))) | |
163 | = 2 * sin(asin(sqrt(d1)) * cos(asin(sqrt(d1))) | |
164 | = 2 * sqrt(d1) * sqrt(1-(sqrt(d1))^2) | |
165 | = 2 * sqrt(d1 - d1 * d1) | |
166 | This relation is used below. | |
167 | ||
168 | As we mentioned above the goal is to find CXTD (named "a" below for | |
169 | brevity) such that ("b" below stands for "d1", and "c" for "d_crs1"): | |
170 | ||
171 | 2 * R * asin(sqrt(a)) == R * asin(2 * sqrt(b-b^2) * sin(c)) | |
172 | ||
173 | Analysis: | |
174 | 2 * R * asin(sqrt(a)) == R * asin(2 * sqrt(b-b^2) * sin(c)) | |
175 | <=> 2 * asin(sqrt(a)) == asin(sqrt(b-b^2) * sin(c)) | |
176 | <=> sin(2 * asin(sqrt(a))) == 2 * sqrt(b-b^2) * sin(c) | |
177 | <=> 2 * sin(asin(sqrt(a))) * cos(asin(sqrt(a))) == 2 * sqrt(b-b^2) * sin(c) | |
178 | <=> 2 * sqrt(a) * sqrt(1-a) == 2 * sqrt(b-b^2) * sin(c) | |
179 | <=> sqrt(a) * sqrt(1-a) == sqrt(b-b^2) * sin(c) | |
180 | <=> sqrt(a-a^2) == sqrt(b-b^2) * sin(c) | |
181 | <=> a-a^2 == (b-b^2) * (sin(c))^2 | |
182 | ||
183 | Consider the quadratic equation: x^2-x+p^2 == 0, | |
184 | where p = sqrt(b-b^2) * sin(c); its discriminant is: | |
185 | d = 1 - 4 * p^2 = 1 - 4 * (b-b^2) * (sin(c))^2 | |
186 | ||
187 | The two solutions are: | |
188 | a_1 = (1 - sqrt(d)) / 2 | |
189 | a_2 = (1 + sqrt(d)) / 2 | |
190 | ||
191 | Which one to choose? | |
192 | "a" refers to the distance (on the unit sphere) of D from the | |
193 | supporting great circle Circ(A,B) of the segment AB. | |
194 | The two different values for "a" correspond to the lengths of the two | |
195 | arcs delimited D and the points of intersection of Circ(A,B) and the | |
196 | great circle perperdicular to Circ(A,B) passing through D. | |
197 | Clearly, the value we want is the smallest among these two distances, | |
198 | hence the root we must choose is the smallest root among the two. | |
199 | ||
200 | So the answer is: | |
201 | CXTD = ( 1 - sqrt(1 - 4 * (b-b^2) * (sin(c))^2) ) / 2 | |
202 | ||
203 | Therefore, in order to implement the comparable version of the cross | |
204 | track strategy we need to: | |
205 | (1) Use the comparable version of the haversine strategy instead of | |
206 | the non-comparable one. | |
207 | (2) Instead of return XTD when D projects inside the segment AB, we | |
208 | need to return CXTD, given by the following formula: | |
209 | CXTD = ( 1 - sqrt(1 - 4 * (d1-d1^2) * (sin(d_crs1))^2) ) / 2; | |
210 | ||
211 | ||
212 | Complexity Analysis | |
213 | ------------------- | |
214 | In the analysis that follows we refer to the actual implementation below. | |
215 | In particular, instead of computing CXTD as above, we use the more | |
216 | efficient (operation-wise) computation of CXTD shown here: | |
217 | ||
218 | return_type sin_d_crs1 = sin(d_crs1); | |
219 | return_type d1_x_sin = d1 * sin_d_crs1; | |
220 | return_type d = d1_x_sin * (sin_d_crs1 - d1_x_sin); | |
221 | return d / (0.5 + math::sqrt(0.25 - d)); | |
222 | ||
223 | Notice that instead of computing: | |
224 | 0.5 - 0.5 * sqrt(1 - 4 * d) = 0.5 - sqrt(0.25 - d) | |
225 | we use the following formula instead: | |
226 | d / (0.5 + sqrt(0.25 - d)). | |
227 | This is done for numerical robustness. The expression 0.5 - sqrt(0.25 - x) | |
228 | has large numerical errors for values of x close to 0 (if using doubles | |
229 | the error start to become large even when d is as large as 0.001). | |
230 | To remedy that, we re-write 0.5 - sqrt(0.25 - x) as: | |
231 | 0.5 - sqrt(0.25 - d) | |
232 | = (0.5 - sqrt(0.25 - d) * (0.5 - sqrt(0.25 - d)) / (0.5 + sqrt(0.25 - d)). | |
233 | The numerator is the difference of two squares: | |
234 | (0.5 - sqrt(0.25 - d) * (0.5 - sqrt(0.25 - d)) | |
235 | = 0.5^2 - (sqrt(0.25 - d))^ = 0.25 - (0.25 - d) = d, | |
236 | which gives the expression we use. | |
237 | ||
238 | For the complexity analysis, we distinguish between two cases: | |
239 | (A) The distance is realized between the point D and an | |
240 | endpoint of the segment AB | |
241 | ||
242 | Gains: | |
243 | Since we are using comparable::haversine<> which is called | |
244 | 3 times, we gain: | |
245 | -> 3 calls to sqrt | |
246 | -> 3 calls to asin | |
247 | -> 6 multiplications | |
248 | ||
249 | Loses: None | |
250 | ||
251 | So the net gain is: | |
252 | -> 6 function calls (sqrt/asin) | |
253 | -> 6 arithmetic operations | |
254 | ||
255 | If we use comparable::cross_track<> to compute | |
256 | cross_track<> we need to account for a call to sqrt, a call | |
257 | to asin and 2 multiplications. In this case the net gain is: | |
258 | -> 4 function calls (sqrt/asin) | |
259 | -> 4 arithmetic operations | |
260 | ||
261 | ||
262 | (B) The distance is realized between the point D and an | |
263 | interior point of the segment AB | |
264 | ||
265 | Gains: | |
266 | Since we are using comparable::haversine<> which is called | |
267 | 3 times, we gain: | |
268 | -> 3 calls to sqrt | |
269 | -> 3 calls to asin | |
270 | -> 6 multiplications | |
271 | Also we gain the operations used to compute XTD: | |
272 | -> 2 calls to sin | |
273 | -> 1 call to asin | |
274 | -> 1 call to abs | |
275 | -> 2 multiplications | |
276 | -> 1 division | |
277 | So the total gains are: | |
278 | -> 9 calls to sqrt/sin/asin | |
279 | -> 1 call to abs | |
280 | -> 8 multiplications | |
281 | -> 1 division | |
282 | ||
283 | Loses: | |
284 | To compute a distance compatible with comparable::haversine<> | |
285 | we need to perform a few more operations, namely: | |
286 | -> 1 call to sin | |
287 | -> 1 call to sqrt | |
288 | -> 2 multiplications | |
289 | -> 1 division | |
290 | -> 1 addition | |
291 | -> 2 subtractions | |
292 | ||
293 | So roughly speaking the net gain is: | |
294 | -> 8 fewer function calls and 3 fewer arithmetic operations | |
295 | ||
296 | If we were to implement cross_track directly from the | |
297 | comparable version (much like what haversine<> does using | |
298 | comparable::haversine<>) we need additionally | |
299 | -> 2 function calls (asin/sqrt) | |
300 | -> 2 multiplications | |
301 | ||
302 | So it pays off to re-implement cross_track<> to use | |
303 | comparable::cross_track<>; in this case the net gain would be: | |
304 | -> 6 function calls | |
305 | -> 1 arithmetic operation | |
306 | ||
307 | Summary/Conclusion | |
308 | ------------------ | |
309 | Following the mathematical and complexity analysis above, the | |
310 | comparable cross track strategy (as implemented below) satisfies | |
311 | all the goal mentioned in the beginning: | |
312 | * It is more efficient than its non-comparable counter-part. | |
313 | * Comparable distances using this new strategy can also be compared | |
314 | with comparable distances computed with the comparable haversine | |
315 | strategy. | |
316 | * It turns out to be more efficient to compute the actual cross | |
317 | track distance XTD by first computing CXTD, and then computing | |
318 | XTD by means of the formula: | |
319 | XTD = 2 * R * asin( sqrt(CXTD) ) | |
320 | */ | |
321 | ||
322 | template | |
323 | < | |
324 | typename CalculationType = void, | |
325 | typename Strategy = comparable::haversine<double, CalculationType> | |
326 | > | |
327 | class cross_track | |
328 | { | |
329 | public : | |
330 | template <typename Point, typename PointOfSegment> | |
331 | struct return_type | |
332 | : promote_floating_point | |
333 | < | |
334 | typename select_calculation_type | |
335 | < | |
336 | Point, | |
337 | PointOfSegment, | |
338 | CalculationType | |
339 | >::type | |
340 | > | |
341 | {}; | |
342 | ||
343 | typedef typename Strategy::radius_type radius_type; | |
344 | ||
345 | inline cross_track() | |
346 | {} | |
347 | ||
348 | explicit inline cross_track(typename Strategy::radius_type const& r) | |
349 | : m_strategy(r) | |
350 | {} | |
351 | ||
352 | inline cross_track(Strategy const& s) | |
353 | : m_strategy(s) | |
354 | {} | |
355 | ||
356 | ||
357 | // It might be useful in the future | |
358 | // to overload constructor with strategy info. | |
359 | // crosstrack(...) {} | |
360 | ||
361 | ||
362 | template <typename Point, typename PointOfSegment> | |
363 | inline typename return_type<Point, PointOfSegment>::type | |
364 | apply(Point const& p, PointOfSegment const& sp1, PointOfSegment const& sp2) const | |
365 | { | |
366 | ||
367 | #if !defined(BOOST_MSVC) | |
368 | BOOST_CONCEPT_ASSERT | |
369 | ( | |
370 | (concepts::PointDistanceStrategy<Strategy, Point, PointOfSegment>) | |
371 | ); | |
372 | #endif | |
373 | ||
374 | typedef typename return_type<Point, PointOfSegment>::type return_type; | |
375 | ||
376 | #ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK | |
377 | std::cout << "Course " << dsv(sp1) << " to " << dsv(p) << " " | |
378 | << crs_AD * geometry::math::r2d<return_type>() << std::endl; | |
379 | std::cout << "Course " << dsv(sp1) << " to " << dsv(sp2) << " " | |
380 | << crs_AB * geometry::math::r2d<return_type>() << std::endl; | |
381 | std::cout << "Course " << dsv(sp2) << " to " << dsv(p) << " " | |
382 | << crs_BD * geometry::math::r2d << std::endl; | |
383 | std::cout << "Projection AD-AB " << projection1 << " : " | |
384 | << d_crs1 * geometry::math::r2d<return_type>() << std::endl; | |
385 | std::cout << "Projection BD-BA " << projection2 << " : " | |
386 | << d_crs2 * geometry::math::r2d<return_type>() << std::endl; | |
387 | #endif | |
388 | ||
389 | // http://williams.best.vwh.net/avform.htm#XTE | |
390 | return_type d1 = m_strategy.apply(sp1, p); | |
391 | return_type d3 = m_strategy.apply(sp1, sp2); | |
392 | ||
393 | if (geometry::math::equals(d3, 0.0)) | |
394 | { | |
395 | // "Degenerate" segment, return either d1 or d2 | |
396 | return d1; | |
397 | } | |
398 | ||
399 | return_type d2 = m_strategy.apply(sp2, p); | |
400 | ||
401 | return_type crs_AD = geometry::detail::course<return_type>(sp1, p); | |
402 | return_type crs_AB = geometry::detail::course<return_type>(sp1, sp2); | |
403 | return_type crs_BA = crs_AB - geometry::math::pi<return_type>(); | |
404 | return_type crs_BD = geometry::detail::course<return_type>(sp2, p); | |
405 | return_type d_crs1 = crs_AD - crs_AB; | |
406 | return_type d_crs2 = crs_BD - crs_BA; | |
407 | ||
408 | // d1, d2, d3 are in principle not needed, only the sign matters | |
409 | return_type projection1 = cos( d_crs1 ) * d1 / d3; | |
410 | return_type projection2 = cos( d_crs2 ) * d2 / d3; | |
411 | ||
412 | if (projection1 > 0.0 && projection2 > 0.0) | |
413 | { | |
414 | #ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK | |
415 | return_type XTD = radius() * geometry::math::abs( asin( sin( d1 ) * sin( d_crs1 ) )); | |
416 | ||
417 | std::cout << "Projection ON the segment" << std::endl; | |
418 | std::cout << "XTD: " << XTD | |
419 | << " d1: " << (d1 * radius()) | |
420 | << " d2: " << (d2 * radius()) | |
421 | << std::endl; | |
422 | #endif | |
423 | return_type const half(0.5); | |
424 | return_type const quarter(0.25); | |
425 | ||
426 | return_type sin_d_crs1 = sin(d_crs1); | |
427 | /* | |
428 | This is the straightforward obvious way to continue: | |
429 | ||
430 | return_type discriminant | |
431 | = 1.0 - 4.0 * (d1 - d1 * d1) * sin_d_crs1 * sin_d_crs1; | |
432 | return 0.5 - 0.5 * math::sqrt(discriminant); | |
433 | ||
434 | Below we optimize the number of arithmetic operations | |
435 | and account for numerical robustness: | |
436 | */ | |
437 | return_type d1_x_sin = d1 * sin_d_crs1; | |
438 | return_type d = d1_x_sin * (sin_d_crs1 - d1_x_sin); | |
439 | return d / (half + math::sqrt(quarter - d)); | |
440 | } | |
441 | else | |
442 | { | |
443 | #ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK | |
444 | std::cout << "Projection OUTSIDE the segment" << std::endl; | |
445 | #endif | |
446 | ||
447 | // Return shortest distance, project either on point sp1 or sp2 | |
448 | return return_type( (std::min)( d1 , d2 ) ); | |
449 | } | |
450 | } | |
451 | ||
452 | inline typename Strategy::radius_type radius() const | |
453 | { return m_strategy.radius(); } | |
454 | ||
455 | private : | |
456 | Strategy m_strategy; | |
457 | }; | |
458 | ||
459 | } // namespace comparable | |
460 | ||
461 | ||
462 | /*! | |
463 | \brief Strategy functor for distance point to segment calculation | |
464 | \ingroup strategies | |
465 | \details Class which calculates the distance of a point to a segment, for points on a sphere or globe | |
466 | \see http://williams.best.vwh.net/avform.htm | |
467 | \tparam CalculationType \tparam_calculation | |
468 | \tparam Strategy underlying point-point distance strategy, defaults to haversine | |
469 | ||
470 | \qbk{ | |
471 | [heading See also] | |
472 | [link geometry.reference.algorithms.distance.distance_3_with_strategy distance (with strategy)] | |
473 | } | |
474 | ||
475 | */ | |
476 | template | |
477 | < | |
478 | typename CalculationType = void, | |
479 | typename Strategy = haversine<double, CalculationType> | |
480 | > | |
481 | class cross_track | |
482 | { | |
483 | public : | |
484 | template <typename Point, typename PointOfSegment> | |
485 | struct return_type | |
486 | : promote_floating_point | |
487 | < | |
488 | typename select_calculation_type | |
489 | < | |
490 | Point, | |
491 | PointOfSegment, | |
492 | CalculationType | |
493 | >::type | |
494 | > | |
495 | {}; | |
496 | ||
497 | typedef typename Strategy::radius_type radius_type; | |
498 | ||
499 | inline cross_track() | |
500 | {} | |
501 | ||
502 | explicit inline cross_track(typename Strategy::radius_type const& r) | |
503 | : m_strategy(r) | |
504 | {} | |
505 | ||
506 | inline cross_track(Strategy const& s) | |
507 | : m_strategy(s) | |
508 | {} | |
509 | ||
510 | ||
511 | // It might be useful in the future | |
512 | // to overload constructor with strategy info. | |
513 | // crosstrack(...) {} | |
514 | ||
515 | ||
516 | template <typename Point, typename PointOfSegment> | |
517 | inline typename return_type<Point, PointOfSegment>::type | |
518 | apply(Point const& p, PointOfSegment const& sp1, PointOfSegment const& sp2) const | |
519 | { | |
520 | ||
521 | #if !defined(BOOST_MSVC) | |
522 | BOOST_CONCEPT_ASSERT | |
523 | ( | |
524 | (concepts::PointDistanceStrategy<Strategy, Point, PointOfSegment>) | |
525 | ); | |
526 | #endif | |
527 | typedef typename return_type<Point, PointOfSegment>::type return_type; | |
528 | typedef cross_track<CalculationType, Strategy> this_type; | |
529 | ||
530 | typedef typename services::comparable_type | |
531 | < | |
532 | this_type | |
533 | >::type comparable_type; | |
534 | ||
535 | comparable_type cstrategy | |
536 | = services::get_comparable<this_type>::apply(m_strategy); | |
537 | ||
538 | return_type const a = cstrategy.apply(p, sp1, sp2); | |
539 | return_type const c = return_type(2.0) * asin(math::sqrt(a)); | |
540 | return c * radius(); | |
541 | } | |
542 | ||
543 | inline typename Strategy::radius_type radius() const | |
544 | { return m_strategy.radius(); } | |
545 | ||
546 | private : | |
547 | ||
548 | Strategy m_strategy; | |
549 | }; | |
550 | ||
551 | ||
552 | ||
553 | #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS | |
554 | namespace services | |
555 | { | |
556 | ||
557 | template <typename CalculationType, typename Strategy> | |
558 | struct tag<cross_track<CalculationType, Strategy> > | |
559 | { | |
560 | typedef strategy_tag_distance_point_segment type; | |
561 | }; | |
562 | ||
563 | ||
564 | template <typename CalculationType, typename Strategy, typename P, typename PS> | |
565 | struct return_type<cross_track<CalculationType, Strategy>, P, PS> | |
566 | : cross_track<CalculationType, Strategy>::template return_type<P, PS> | |
567 | {}; | |
568 | ||
569 | ||
570 | template <typename CalculationType, typename Strategy> | |
571 | struct comparable_type<cross_track<CalculationType, Strategy> > | |
572 | { | |
573 | typedef comparable::cross_track | |
574 | < | |
575 | CalculationType, typename comparable_type<Strategy>::type | |
576 | > type; | |
577 | }; | |
578 | ||
579 | ||
580 | template | |
581 | < | |
582 | typename CalculationType, | |
583 | typename Strategy | |
584 | > | |
585 | struct get_comparable<cross_track<CalculationType, Strategy> > | |
586 | { | |
587 | typedef typename comparable_type | |
588 | < | |
589 | cross_track<CalculationType, Strategy> | |
590 | >::type comparable_type; | |
591 | public : | |
592 | static inline comparable_type | |
593 | apply(cross_track<CalculationType, Strategy> const& strategy) | |
594 | { | |
595 | return comparable_type(strategy.radius()); | |
596 | } | |
597 | }; | |
598 | ||
599 | ||
600 | template | |
601 | < | |
602 | typename CalculationType, | |
603 | typename Strategy, | |
604 | typename P, | |
605 | typename PS | |
606 | > | |
607 | struct result_from_distance<cross_track<CalculationType, Strategy>, P, PS> | |
608 | { | |
609 | private : | |
610 | typedef typename cross_track | |
611 | < | |
612 | CalculationType, Strategy | |
613 | >::template return_type<P, PS>::type return_type; | |
614 | public : | |
615 | template <typename T> | |
616 | static inline return_type | |
617 | apply(cross_track<CalculationType, Strategy> const& , T const& distance) | |
618 | { | |
619 | return distance; | |
620 | } | |
621 | }; | |
622 | ||
623 | ||
624 | // Specializations for comparable::cross_track | |
625 | template <typename RadiusType, typename CalculationType> | |
626 | struct tag<comparable::cross_track<RadiusType, CalculationType> > | |
627 | { | |
628 | typedef strategy_tag_distance_point_segment type; | |
629 | }; | |
630 | ||
631 | ||
632 | template | |
633 | < | |
634 | typename RadiusType, | |
635 | typename CalculationType, | |
636 | typename P, | |
637 | typename PS | |
638 | > | |
639 | struct return_type<comparable::cross_track<RadiusType, CalculationType>, P, PS> | |
640 | : comparable::cross_track | |
641 | < | |
642 | RadiusType, CalculationType | |
643 | >::template return_type<P, PS> | |
644 | {}; | |
645 | ||
646 | ||
647 | template <typename RadiusType, typename CalculationType> | |
648 | struct comparable_type<comparable::cross_track<RadiusType, CalculationType> > | |
649 | { | |
650 | typedef comparable::cross_track<RadiusType, CalculationType> type; | |
651 | }; | |
652 | ||
653 | ||
654 | template <typename RadiusType, typename CalculationType> | |
655 | struct get_comparable<comparable::cross_track<RadiusType, CalculationType> > | |
656 | { | |
657 | private : | |
658 | typedef comparable::cross_track<RadiusType, CalculationType> this_type; | |
659 | public : | |
660 | static inline this_type apply(this_type const& input) | |
661 | { | |
662 | return input; | |
663 | } | |
664 | }; | |
665 | ||
666 | ||
667 | template | |
668 | < | |
669 | typename RadiusType, | |
670 | typename CalculationType, | |
671 | typename P, | |
672 | typename PS | |
673 | > | |
674 | struct result_from_distance | |
675 | < | |
676 | comparable::cross_track<RadiusType, CalculationType>, P, PS | |
677 | > | |
678 | { | |
679 | private : | |
680 | typedef comparable::cross_track<RadiusType, CalculationType> strategy_type; | |
681 | typedef typename return_type<strategy_type, P, PS>::type return_type; | |
682 | public : | |
683 | template <typename T> | |
684 | static inline return_type apply(strategy_type const& strategy, | |
685 | T const& distance) | |
686 | { | |
687 | return_type const s | |
688 | = sin( (distance / strategy.radius()) / return_type(2.0) ); | |
689 | return s * s; | |
690 | } | |
691 | }; | |
692 | ||
693 | ||
694 | ||
695 | /* | |
696 | ||
697 | TODO: spherical polar coordinate system requires "get_as_radian_equatorial<>" | |
698 | ||
699 | template <typename Point, typename PointOfSegment, typename Strategy> | |
700 | struct default_strategy | |
701 | < | |
702 | segment_tag, Point, PointOfSegment, | |
703 | spherical_polar_tag, spherical_polar_tag, | |
704 | Strategy | |
705 | > | |
706 | { | |
707 | typedef cross_track | |
708 | < | |
709 | void, | |
710 | typename boost::mpl::if_ | |
711 | < | |
712 | boost::is_void<Strategy>, | |
713 | typename default_strategy | |
714 | < | |
715 | point_tag, Point, PointOfSegment, | |
716 | spherical_polar_tag, spherical_polar_tag | |
717 | >::type, | |
718 | Strategy | |
719 | >::type | |
720 | > type; | |
721 | }; | |
722 | */ | |
723 | ||
724 | template <typename Point, typename PointOfSegment, typename Strategy> | |
725 | struct default_strategy | |
726 | < | |
727 | point_tag, segment_tag, Point, PointOfSegment, | |
728 | spherical_equatorial_tag, spherical_equatorial_tag, | |
729 | Strategy | |
730 | > | |
731 | { | |
732 | typedef cross_track | |
733 | < | |
734 | void, | |
735 | typename boost::mpl::if_ | |
736 | < | |
737 | boost::is_void<Strategy>, | |
738 | typename default_strategy | |
739 | < | |
740 | point_tag, point_tag, Point, PointOfSegment, | |
741 | spherical_equatorial_tag, spherical_equatorial_tag | |
742 | >::type, | |
743 | Strategy | |
744 | >::type | |
745 | > type; | |
746 | }; | |
747 | ||
748 | ||
749 | template <typename PointOfSegment, typename Point, typename Strategy> | |
750 | struct default_strategy | |
751 | < | |
752 | segment_tag, point_tag, PointOfSegment, Point, | |
753 | spherical_equatorial_tag, spherical_equatorial_tag, | |
754 | Strategy | |
755 | > | |
756 | { | |
757 | typedef typename default_strategy | |
758 | < | |
759 | point_tag, segment_tag, Point, PointOfSegment, | |
760 | spherical_equatorial_tag, spherical_equatorial_tag, | |
761 | Strategy | |
762 | >::type type; | |
763 | }; | |
764 | ||
765 | ||
766 | } // namespace services | |
767 | #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS | |
768 | ||
769 | }} // namespace strategy::distance | |
770 | ||
771 | }} // namespace boost::geometry | |
772 | ||
773 | #endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_CROSS_TRACK_HPP |