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