strategies/spherical/distance_cross_track.hpp

0001 // Boost.Geometry (aka GGL, Generic Geometry Library)
0002
0003 // Copyright (c) 2007-2014 Barend Gehrels, Amsterdam, the Netherlands.
0004
0005 // This file was modified by Oracle on 2014-2021.
0006 // Modifications copyright (c) 2014-2021, Oracle and/or its affiliates.
0007
0008 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
0009 // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
0010 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0011
0012 // Use, modification and distribution is subject to the Boost Software License,
0013 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0014 // http://www.boost.org/LICENSE_1_0.txt)
0015
0016 #ifndef BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_CROSS_TRACK_HPP
0017 #define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_CROSS_TRACK_HPP
0018
0019 #include <algorithm>
0020 #include <type_traits>
0021
0022 #include <boost/config.hpp>
0023 #include <boost/concept_check.hpp>
0024
0025 #include <boost/geometry/core/cs.hpp>
0026 #include <boost/geometry/core/access.hpp>
0027 #include <boost/geometry/core/coordinate_promotion.hpp>
0028 #include <boost/geometry/core/radian_access.hpp>
0029 #include <boost/geometry/core/tags.hpp>
0030
0031 #include <boost/geometry/formulas/spherical.hpp>
0032
0033 #include <boost/geometry/strategies/distance.hpp>
0034 #include <boost/geometry/strategies/concepts/distance_concept.hpp>
0035 #include <boost/geometry/strategies/spherical/distance_haversine.hpp>
0036 #include <boost/geometry/strategies/spherical/point_in_point.hpp>
0037 #include <boost/geometry/strategies/spherical/intersection.hpp>
0038
0039 #include <boost/geometry/util/math.hpp>
0040 #include <boost/geometry/util/select_calculation_type.hpp>
0041
0042 #ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK
0043 #  include <boost/geometry/io/dsv/write.hpp>
0044 #endif
0045
0046
0047 namespace boost { namespace geometry
0048 {
0049
0050 namespace strategy { namespace distance
0051 {
0052
0053 #ifndef DOXYGEN_NO_DETAIL
0054 namespace detail
0055 {
0056     template <typename CalculationType>
0057     struct compute_cross_track_pair
0058     {
0059         template <typename Point, typename PointOfSegment>
0060         static inline auto apply(Point const& p,
0061                                  PointOfSegment const& sp1,
0062                                  PointOfSegment const& sp2)
0063         {
0064             CalculationType lon1 = geometry::get_as_radian<0>(sp1);
0065             CalculationType lat1 = geometry::get_as_radian<1>(sp1);
0066             CalculationType lon2 = geometry::get_as_radian<0>(sp2);
0067             CalculationType lat2 = geometry::get_as_radian<1>(sp2);
0068             CalculationType lon = geometry::get_as_radian<0>(p);
0069             CalculationType lat = geometry::get_as_radian<1>(p);
0070
0071             CalculationType const crs_AD = geometry::formula::spherical_azimuth
0072                 <
0073                     CalculationType,
0074                     false
0075                 >(lon1, lat1, lon, lat).azimuth;
0076
0077             auto result = geometry::formula::spherical_azimuth
0078                 <
0079                     CalculationType,
0080                     true
0081                 >(lon1, lat1, lon2, lat2);
0082
0083             CalculationType crs_AB = result.azimuth;
0084             CalculationType crs_BA = result.reverse_azimuth -
0085                 geometry::math::pi<CalculationType>();
0086
0087             CalculationType crs_BD = geometry::formula::spherical_azimuth
0088                 <
0089                     CalculationType,
0090                     false
0091                 >(lon2, lat2, lon, lat).azimuth;
0092
0093             CalculationType d_crs1 = crs_AD - crs_AB;
0094             CalculationType d_crs2 = crs_BD - crs_BA;
0095
0096             return std::pair<CalculationType, CalculationType>(d_crs1, d_crs2);
0097         }
0098     };
0099
0100     struct compute_cross_track_distance
0101     {
0102         template <typename CalculationType>
0103         static inline auto apply(CalculationType const& d_crs1,
0104                                  CalculationType const& d1)
0105         {
0106             CalculationType const half(0.5);
0107             CalculationType const quarter(0.25);
0108
0109             CalculationType sin_d_crs1 = sin(d_crs1);
0110             /*
0111               This is the straightforward obvious way to continue:
0112
0113               return_type discriminant
0114                   = 1.0 - 4.0 * (d1 - d1 * d1) * sin_d_crs1 * sin_d_crs1;
0115               return 0.5 - 0.5 * math::sqrt(discriminant);
0116
0117               Below we optimize the number of arithmetic operations
0118               and account for numerical robustness:
0119             */
0120             CalculationType d1_x_sin = d1 * sin_d_crs1;
0121             CalculationType d = d1_x_sin * (sin_d_crs1 - d1_x_sin);
0122             return d / (half + math::sqrt(quarter - d));
0123         }
0124     };
0125
0126 }
0127 #endif // DOXYGEN_NO_DETAIL
0128
0129
0130 namespace comparable
0131 {
0132
0133 /*
0134   Given a spherical segment AB and a point D, we are interested in
0135   computing the distance of D from AB. This is usually known as the
0136   cross track distance.
0137
0138   If the projection (along great circles) of the point D lies inside
0139   the segment AB, then the distance (cross track error) XTD is given
0140   by the formula (see http://williams.best.vwh.net/avform.htm#XTE):
0141
0142   XTD = asin( sin(dist_AD) * sin(crs_AD-crs_AB) )
0143
0144   where dist_AD is the great circle distance between the points A and
0145   B, and crs_AD, crs_AB is the course (bearing) between the points A,
0146   D and A, B, respectively.
0147
0148   If the point D does not project inside the arc AB, then the distance
0149   of D from AB is the minimum of the two distances dist_AD and dist_BD.
0150
0151   Our reference implementation for this procedure is listed below
0152   (this was the old Boost.Geometry implementation of the cross track distance),
0153   where:
0154   * The member variable m_strategy is the underlying haversine strategy.
0155   * p stands for the point D.
0156   * sp1 stands for the segment endpoint A.
0157   * sp2 stands for the segment endpoint B.
0158
0159   ================= reference implementation -- start =================
0160
0161   return_type d1 = m_strategy.apply(sp1, p);
0162   return_type d3 = m_strategy.apply(sp1, sp2);
0163
0164   if (geometry::math::equals(d3, 0.0))
0165   {
0166       // "Degenerate" segment, return either d1 or d2
0167       return d1;
0168   }
0169
0170   return_type d2 = m_strategy.apply(sp2, p);
0171
0172   return_type crs_AD = geometry::detail::course<return_type>(sp1, p);
0173   return_type crs_AB = geometry::detail::course<return_type>(sp1, sp2);
0174   return_type crs_BA = crs_AB - geometry::math::pi<return_type>();
0175   return_type crs_BD = geometry::detail::course<return_type>(sp2, p);
0176   return_type d_crs1 = crs_AD - crs_AB;
0177   return_type d_crs2 = crs_BD - crs_BA;
0178
0179   // d1, d2, d3 are in principle not needed, only the sign matters
0180   return_type projection1 = cos( d_crs1 ) * d1 / d3;
0181   return_type projection2 = cos( d_crs2 ) * d2 / d3;
0182
0183   if (projection1 > 0.0 && projection2 > 0.0)
0184   {
0185       return_type XTD
0186           = radius() * math::abs( asin( sin( d1 / radius() ) * sin( d_crs1 ) ));
0187
0188       // Return shortest distance, projected point on segment sp1-sp2
0189       return return_type(XTD);
0190   }
0191   else
0192   {
0193       // Return shortest distance, project either on point sp1 or sp2
0194       return return_type( (std::min)( d1 , d2 ) );
0195   }
0196
0197   ================= reference implementation -- end =================
0198
0199
0200   Motivation
0201   ----------
0202   In what follows we develop a comparable version of the cross track
0203   distance strategy, that meets the following goals:
0204   * It is more efficient than the original cross track strategy (less
0205     operations and less calls to mathematical functions).
0206   * Distances using the comparable cross track strategy can not only
0207     be compared with other distances using the same strategy, but also with
0208     distances computed with the comparable version of the haversine strategy.
0209   * It can serve as the basis for the computation of the cross track distance,
0210     as it is more efficient to compute its comparable version and
0211     transform that to the actual cross track distance, rather than
0212     follow/use the reference implementation listed above.
0213
0214   Major idea
0215   ----------
0216   The idea here is to use the comparable haversine strategy to compute
0217   the distances d1, d2 and d3 in the above listing. Once we have done
0218   that we need also to make sure that instead of returning XTD (as
0219   computed above) that we return a distance CXTD that is compatible
0220   with the comparable haversine distance. To achieve this CXTD must satisfy
0221   the relation:
0222       XTD = 2 * R * asin( sqrt(XTD) )
0223   where R is the sphere's radius.
0224
0225   Below we perform the mathematical analysis that show how to compute CXTD.
0226
0227
0228   Mathematical analysis
0229   ---------------------
0230   Below we use the following trigonometric identities:
0231       sin(2 * x) = 2 * sin(x) * cos(x)
0232       cos(asin(x)) = sqrt(1 - x^2)
0233
0234   Observation:
0235   The distance d1 needed when the projection of the point D is within the
0236   segment must be the true distance. However, comparable::haversine<>
0237   returns a comparable distance instead of the one needed.
0238   To remedy this, we implicitly compute what is needed.
0239   More precisely, we need to compute sin(true_d1):
0240
0241   sin(true_d1) = sin(2 * asin(sqrt(d1)))
0242                = 2 * sin(asin(sqrt(d1)) * cos(asin(sqrt(d1)))
0243                = 2 * sqrt(d1) * sqrt(1-(sqrt(d1))^2)
0244                = 2 * sqrt(d1 - d1 * d1)
0245   This relation is used below.
0246
0247   As we mentioned above the goal is to find CXTD (named "a" below for
0248   brevity) such that ("b" below stands for "d1", and "c" for "d_crs1"):
0249
0250       2 * R * asin(sqrt(a)) == R * asin(2 * sqrt(b-b^2) * sin(c))
0251
0252   Analysis:
0253       2 * R * asin(sqrt(a)) == R * asin(2 * sqrt(b-b^2) * sin(c))
0254   <=> 2 * asin(sqrt(a)) == asin(sqrt(b-b^2) * sin(c))
0255   <=> sin(2 * asin(sqrt(a))) == 2 * sqrt(b-b^2) * sin(c)
0256   <=> 2 * sin(asin(sqrt(a))) * cos(asin(sqrt(a))) == 2 * sqrt(b-b^2) * sin(c)
0257   <=> 2 * sqrt(a) * sqrt(1-a) == 2 * sqrt(b-b^2) * sin(c)
0258   <=> sqrt(a) * sqrt(1-a) == sqrt(b-b^2) * sin(c)
0259   <=> sqrt(a-a^2) == sqrt(b-b^2) * sin(c)
0260   <=> a-a^2 == (b-b^2) * (sin(c))^2
0261
0262   Consider the quadratic equation: x^2-x+p^2 == 0,
0263   where p = sqrt(b-b^2) * sin(c); its discriminant is:
0264       d = 1 - 4 * p^2 = 1 - 4 * (b-b^2) * (sin(c))^2
0265
0266   The two solutions are:
0267       a_1 = (1 - sqrt(d)) / 2
0268       a_2 = (1 + sqrt(d)) / 2
0269
0270   Which one to choose?
0271   "a" refers to the distance (on the unit sphere) of D from the
0272   supporting great circle Circ(A,B) of the segment AB.
0273   The two different values for "a" correspond to the lengths of the two
0274   arcs delimited D and the points of intersection of Circ(A,B) and the
0275   great circle perperdicular to Circ(A,B) passing through D.
0276   Clearly, the value we want is the smallest among these two distances,
0277   hence the root we must choose is the smallest root among the two.
0278
0279   So the answer is:
0280       CXTD = ( 1 - sqrt(1 - 4 * (b-b^2) * (sin(c))^2) ) / 2
0281
0282   Therefore, in order to implement the comparable version of the cross
0283   track strategy we need to:
0284   (1) Use the comparable version of the haversine strategy instead of
0285       the non-comparable one.
0286   (2) Instead of return XTD when D projects inside the segment AB, we
0287       need to return CXTD, given by the following formula:
0288           CXTD = ( 1 - sqrt(1 - 4 * (d1-d1^2) * (sin(d_crs1))^2) ) / 2;
0289
0290
0291   Complexity Analysis
0292   -------------------
0293   In the analysis that follows we refer to the actual implementation below.
0294   In particular, instead of computing CXTD as above, we use the more
0295   efficient (operation-wise) computation of CXTD shown here:
0296
0297       return_type sin_d_crs1 = sin(d_crs1);
0298       return_type d1_x_sin = d1 * sin_d_crs1;
0299       return_type d = d1_x_sin * (sin_d_crs1 - d1_x_sin);
0300       return d / (0.5 + math::sqrt(0.25 - d));
0301
0302   Notice that instead of computing:
0303       0.5 - 0.5 * sqrt(1 - 4 * d) = 0.5 - sqrt(0.25 - d)
0304   we use the following formula instead:
0305       d / (0.5 + sqrt(0.25 - d)).
0306   This is done for numerical robustness. The expression 0.5 - sqrt(0.25 - x)
0307   has large numerical errors for values of x close to 0 (if using doubles
0308   the error start to become large even when d is as large as 0.001).
0309   To remedy that, we re-write 0.5 - sqrt(0.25 - x) as:
0310       0.5 - sqrt(0.25 - d)
0311       = (0.5 - sqrt(0.25 - d) * (0.5 - sqrt(0.25 - d)) / (0.5 + sqrt(0.25 - d)).
0312   The numerator is the difference of two squares:
0313       (0.5 - sqrt(0.25 - d) * (0.5 - sqrt(0.25 - d))
0314       = 0.5^2 - (sqrt(0.25 - d))^ = 0.25 - (0.25 - d) = d,
0315   which gives the expression we use.
0316
0317   For the complexity analysis, we distinguish between two cases:
0318   (A) The distance is realized between the point D and an
0319       endpoint of the segment AB
0320
0321       Gains:
0322       Since we are using comparable::haversine<> which is called
0323       3 times, we gain:
0324       -> 3 calls to sqrt
0325       -> 3 calls to asin
0326       -> 6 multiplications
0327
0328       Loses: None
0329
0330       So the net gain is:
0331       -> 6 function calls (sqrt/asin)
0332       -> 6 arithmetic operations
0333
0334       If we use comparable::cross_track<> to compute
0335       cross_track<> we need to account for a call to sqrt, a call
0336       to asin and 2 multiplications. In this case the net gain is:
0337       -> 4 function calls (sqrt/asin)
0338       -> 4 arithmetic operations
0339
0340
0341   (B) The distance is realized between the point D and an
0342       interior point of the segment AB
0343
0344       Gains:
0345       Since we are using comparable::haversine<> which is called
0346       3 times, we gain:
0347       -> 3 calls to sqrt
0348       -> 3 calls to asin
0349       -> 6 multiplications
0350       Also we gain the operations used to compute XTD:
0351       -> 2 calls to sin
0352       -> 1 call to asin
0353       -> 1 call to abs
0354       -> 2 multiplications
0355       -> 1 division
0356       So the total gains are:
0357       -> 9 calls to sqrt/sin/asin
0358       -> 1 call to abs
0359       -> 8 multiplications
0360       -> 1 division
0361
0362       Loses:
0363       To compute a distance compatible with comparable::haversine<>
0364       we need to perform a few more operations, namely:
0365       -> 1 call to sin
0366       -> 1 call to sqrt
0367       -> 2 multiplications
0368       -> 1 division
0369       -> 1 addition
0370       -> 2 subtractions
0371
0372       So roughly speaking the net gain is:
0373       -> 8 fewer function calls and 3 fewer arithmetic operations
0374
0375       If we were to implement cross_track directly from the
0376       comparable version (much like what haversine<> does using
0377       comparable::haversine<>) we need additionally
0378       -> 2 function calls (asin/sqrt)
0379       -> 2 multiplications
0380
0381       So it pays off to re-implement cross_track<> to use
0382       comparable::cross_track<>; in this case the net gain would be:
0383       -> 6 function calls
0384       -> 1 arithmetic operation
0385
0386    Summary/Conclusion
0387    ------------------
0388    Following the mathematical and complexity analysis above, the
0389    comparable cross track strategy (as implemented below) satisfies
0390    all the goal mentioned in the beginning:
0391    * It is more efficient than its non-comparable counter-part.
0392    * Comparable distances using this new strategy can also be compared
0393      with comparable distances computed with the comparable haversine
0394      strategy.
0395    * It turns out to be more efficient to compute the actual cross
0396      track distance XTD by first computing CXTD, and then computing
0397      XTD by means of the formula:
0398                 XTD = 2 * R * asin( sqrt(CXTD) )
0399 */
0400
0401 template
0402 <
0403     typename CalculationType = void,
0404     typename Strategy = comparable::haversine<double, CalculationType>
0405 >
0406 class cross_track
0407 {
0408 public:
0409     template <typename Point, typename PointOfSegment>
0410     struct return_type
0411         : promote_floating_point
0412           <
0413               typename select_calculation_type
0414                   <
0415                       Point,
0416                       PointOfSegment,
0417                       CalculationType
0418                   >::type
0419           >
0420     {};
0421
0422     using radius_type = typename Strategy::radius_type;
0423
0424     cross_track() = default;
0425
0426     explicit inline cross_track(typename Strategy::radius_type const& r)
0427         : m_strategy(r)
0428     {}
0429
0430     inline cross_track(Strategy const& s)
0431         : m_strategy(s)
0432     {}
0433
0434     // It might be useful in the future
0435     // to overload constructor with strategy info.
0436     // crosstrack(...) {}
0437
0438
0439     template <typename Point, typename PointOfSegment>
0440     inline typename return_type<Point, PointOfSegment>::type
0441     apply(Point const& p, PointOfSegment const& sp1, PointOfSegment const& sp2) const
0442     {
0443
0444 #if !defined(BOOST_MSVC)
0445         BOOST_CONCEPT_ASSERT
0446             (
0447                 (concepts::PointDistanceStrategy<Strategy, Point, PointOfSegment>)
0448             );
0449 #endif
0450
0451         using return_type = typename return_type<Point, PointOfSegment>::type;
0452
0453         // http://williams.best.vwh.net/avform.htm#XTE
0454         return_type d1 = m_strategy.apply(sp1, p);
0455         return_type d3 = m_strategy.apply(sp1, sp2);
0456
0457         if (geometry::math::equals(d3, 0.0))
0458         {
0459             // "Degenerate" segment, return either d1 or d2
0460             return d1;
0461         }
0462
0463         return_type d2 = m_strategy.apply(sp2, p);
0464
0465         auto d_crs_pair = detail::compute_cross_track_pair<return_type>::apply(
0466             p, sp1, sp2);
0467
0468         // d1, d2, d3 are in principle not needed, only the sign matters
0469         return_type projection1 = cos(d_crs_pair.first) * d1 / d3;
0470         return_type projection2 = cos(d_crs_pair.second) * d2 / d3;
0471
0472 #ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK
0473         std::cout << "Course " << dsv(sp1) << " to " << dsv(p) << " "
0474                   << crs_AD * geometry::math::r2d<return_type>() << std::endl;
0475         std::cout << "Course " << dsv(sp1) << " to " << dsv(sp2) << " "
0476                   << crs_AB * geometry::math::r2d<return_type>() << std::endl;
0477         std::cout << "Course " << dsv(sp2) << " to " << dsv(sp1) << " "
0478                   << crs_BA * geometry::math::r2d<return_type>() << std::endl;
0479         std::cout << "Course " << dsv(sp2) << " to " << dsv(p) << " "
0480                   << crs_BD * geometry::math::r2d<return_type>() << std::endl;
0481         std::cout << "Projection AD-AB " << projection1 << " : "
0482                   << d_crs1 * geometry::math::r2d<return_type>() << std::endl;
0483         std::cout << "Projection BD-BA " << projection2 << " : "
0484                   << d_crs2 * geometry::math::r2d<return_type>() << std::endl;
0485         std::cout << " d1: " << (d1 )
0486                   << " d2: " << (d2 )
0487                   << std::endl;
0488 #endif
0489
0490         if (projection1 > 0.0 && projection2 > 0.0)
0491         {
0492 #ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK
0493             return_type XTD = radius() * geometry::math::abs( asin( sin( d1 ) * sin( d_crs1 ) ));
0494
0495             std::cout << "Projection ON the segment" << std::endl;
0496             std::cout << "XTD: " << XTD
0497                       << " d1: " << (d1 * radius())
0498                       << " d2: " << (d2 * radius())
0499                       << std::endl;
0500 #endif
0501             return detail::compute_cross_track_distance::apply(
0502                 d_crs_pair.first, d1);
0503         }
0504         else
0505         {
0506 #ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK
0507             std::cout << "Projection OUTSIDE the segment" << std::endl;
0508 #endif
0509             // Return shortest distance, project either on point sp1 or sp2
0510             return return_type( (std::min)( d1 , d2 ) );
0511         }
0512     }
0513
0514     template <typename T1, typename T2>
0515     inline radius_type vertical_or_meridian(T1 lat1, T2 lat2) const
0516     {
0517         return m_strategy.radius() * (lat1 - lat2);
0518     }
0519
0520     inline typename Strategy::radius_type radius() const
0521     { return m_strategy.radius(); }
0522
0523 private :
0524     Strategy m_strategy;
0525 };
0526
0527 } // namespace comparable
0528
0529
0530 /*!
0531 \brief Strategy functor for distance point to segment calculation
0532 \ingroup strategies
0533 \details Class which calculates the distance of a point to a segment, for points on a sphere or globe
0534 \see http://williams.best.vwh.net/avform.htm
0535 \tparam CalculationType \tparam_calculation
0536 \tparam Strategy underlying point-point distance strategy, defaults to haversine
0537
0538 \qbk{
0539 [heading See also]
0540 [link geometry.reference.algorithms.distance.distance_3_with_strategy distance (with strategy)]
0541 }
0542
0543 */
0544 template
0545 <
0546     typename CalculationType = void,
0547     typename Strategy = haversine<double, CalculationType>
0548 >
0549 class cross_track
0550 {
0551 public :
0552
0553     template <typename Point, typename PointOfSegment>
0554     struct return_type
0555         : promote_floating_point
0556           <
0557               typename select_calculation_type
0558                   <
0559                       Point,
0560                       PointOfSegment,
0561                       CalculationType
0562                   >::type
0563           >
0564     {};
0565
0566     using radius_type = typename Strategy::radius_type;
0567
0568     inline cross_track()
0569     {}
0570
0571     explicit inline cross_track(typename Strategy::radius_type const& r)
0572         : m_strategy(r)
0573     {}
0574
0575     inline cross_track(Strategy const& s)
0576         : m_strategy(s)
0577     {}
0578
0579     // It might be useful in the future
0580     // to overload constructor with strategy info.
0581     // crosstrack(...) {}
0582
0583
0584     template <typename Point, typename PointOfSegment>
0585     inline auto apply(Point const& p,
0586                       PointOfSegment const& sp1,
0587                       PointOfSegment const& sp2) const
0588     {
0589
0590 #if !defined(BOOST_MSVC)
0591         BOOST_CONCEPT_ASSERT
0592             (
0593                 (concepts::PointDistanceStrategy<Strategy, Point, PointOfSegment>)
0594             );
0595 #endif
0596         using return_type = typename return_type<Point, PointOfSegment>::type;
0597         using this_type = cross_track<CalculationType, Strategy>;
0598
0599         using comparable_type = typename services::comparable_type
0600             <
0601                 this_type
0602             >::type;
0603
0604         comparable_type cstrategy
0605             = services::get_comparable<this_type>::apply(m_strategy);
0606
0607         return_type const a = cstrategy.apply(p, sp1, sp2);
0608         return_type const c = return_type(2.0) * asin(math::sqrt(a));
0609         return c * radius();
0610     }
0611
0612     template <typename T1, typename T2>
0613     inline radius_type vertical_or_meridian(T1 lat1, T2 lat2) const
0614     {
0615         return m_strategy.radius() * (lat1 - lat2);
0616     }
0617
0618     inline typename Strategy::radius_type radius() const
0619     { return m_strategy.radius(); }
0620
0621 private :
0622
0623     Strategy m_strategy;
0624 };
0625
0626
0627
0628 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
0629 namespace services
0630 {
0631
0632 template <typename CalculationType, typename Strategy>
0633 struct tag<cross_track<CalculationType, Strategy> >
0634 {
0635     using type = strategy_tag_distance_point_segment;
0636 };
0637
0638
0639 template <typename CalculationType, typename Strategy, typename P, typename PS>
0640 struct return_type<cross_track<CalculationType, Strategy>, P, PS>
0641     : cross_track<CalculationType, Strategy>::template return_type<P, PS>
0642 {};
0643
0644
0645 template <typename CalculationType, typename Strategy>
0646 struct comparable_type<cross_track<CalculationType, Strategy> >
0647 {
0648     using type = comparable::cross_track
0649         <
0650             CalculationType, typename comparable_type<Strategy>::type
0651         > ;
0652 };
0653
0654
0655 template
0656 <
0657     typename CalculationType,
0658     typename Strategy
0659 >
0660 struct get_comparable<cross_track<CalculationType, Strategy> >
0661 {
0662     using comparable_type = typename comparable_type
0663         <
0664             cross_track<CalculationType, Strategy>
0665         >::type;
0666 public :
0667     static inline comparable_type
0668     apply(cross_track<CalculationType, Strategy> const& strategy)
0669     {
0670         return comparable_type(strategy.radius());
0671     }
0672 };
0673
0674
0675 template
0676 <
0677     typename CalculationType,
0678     typename Strategy,
0679     typename P,
0680     typename PS
0681 >
0682 struct result_from_distance<cross_track<CalculationType, Strategy>, P, PS>
0683 {
0684 private :
0685     using return_type = typename cross_track
0686         <
0687             CalculationType, Strategy
0688         >::template return_type<P, PS>::type;
0689 public :
0690     template <typename T>
0691     static inline return_type
0692     apply(cross_track<CalculationType, Strategy> const& , T const& distance)
0693     {
0694         return distance;
0695     }
0696 };
0697
0698
0699 // Specializations for comparable::cross_track
0700 template <typename RadiusType, typename CalculationType>
0701 struct tag<comparable::cross_track<RadiusType, CalculationType> >
0702 {
0703     using type = strategy_tag_distance_point_segment;
0704 };
0705
0706
0707 template
0708 <
0709     typename RadiusType,
0710     typename CalculationType,
0711     typename P,
0712     typename PS
0713 >
0714 struct return_type<comparable::cross_track<RadiusType, CalculationType>, P, PS>
0715     : comparable::cross_track
0716         <
0717             RadiusType, CalculationType
0718         >::template return_type<P, PS>
0719 {};
0720
0721
0722 template <typename RadiusType, typename CalculationType>
0723 struct comparable_type<comparable::cross_track<RadiusType, CalculationType> >
0724 {
0725     using type = comparable::cross_track<RadiusType, CalculationType>;
0726 };
0727
0728
0729 template <typename RadiusType, typename CalculationType>
0730 struct get_comparable<comparable::cross_track<RadiusType, CalculationType> >
0731 {
0732 private :
0733     using this_type = comparable::cross_track<RadiusType, CalculationType>;
0734 public :
0735     static inline this_type apply(this_type const& input)
0736     {
0737         return input;
0738     }
0739 };
0740
0741
0742 template
0743 <
0744     typename RadiusType,
0745     typename CalculationType,
0746     typename P,
0747     typename PS
0748 >
0749 struct result_from_distance
0750     <
0751         comparable::cross_track<RadiusType, CalculationType>, P, PS
0752     >
0753 {
0754 private :
0755     using strategy_type = comparable::cross_track<RadiusType, CalculationType>;
0756     using return_type = typename return_type<strategy_type, P, PS>::type;
0757 public :
0758     template <typename T>
0759     static inline return_type apply(strategy_type const& strategy,
0760                                     T const& distance)
0761     {
0762         return_type const s
0763             = sin( (distance / strategy.radius()) / return_type(2.0) );
0764         return s * s;
0765     }
0766 };
0767
0768
0769
0770 /*
0771
0772 TODO:  spherical polar coordinate system requires "get_as_radian_equatorial<>"
0773
0774 template <typename Point, typename PointOfSegment, typename Strategy>
0775 struct default_strategy
0776     <
0777         segment_tag, Point, PointOfSegment,
0778         spherical_polar_tag, spherical_polar_tag,
0779         Strategy
0780     >
0781 {
0782     typedef cross_track
0783         <
0784             void,
0785             std::conditional_t
0786                 <
0787                     std::is_void<Strategy>::value,
0788                     typename default_strategy
0789                         <
0790                             point_tag, Point, PointOfSegment,
0791                             spherical_polar_tag, spherical_polar_tag
0792                         >::type,
0793                     Strategy
0794                 >
0795         > type;
0796 };
0797 */
0798
0799 template <typename Point, typename PointOfSegment, typename Strategy>
0800 struct default_strategy
0801     <
0802         point_tag, segment_tag, Point, PointOfSegment,
0803         spherical_equatorial_tag, spherical_equatorial_tag,
0804         Strategy
0805     >
0806 {
0807     using type = cross_track
0808         <
0809             void,
0810             std::conditional_t
0811                 <
0812                     std::is_void<Strategy>::value,
0813                     typename default_strategy
0814                         <
0815                             point_tag, point_tag, Point, PointOfSegment,
0816                             spherical_equatorial_tag, spherical_equatorial_tag
0817                         >::type,
0818                     Strategy
0819                 >
0820         >;
0821 };
0822
0823
0824 template <typename PointOfSegment, typename Point, typename Strategy>
0825 struct default_strategy
0826     <
0827         segment_tag, point_tag, PointOfSegment, Point,
0828         spherical_equatorial_tag, spherical_equatorial_tag,
0829         Strategy
0830     >
0831 {
0832     using type = typename default_strategy
0833         <
0834             point_tag, segment_tag, Point, PointOfSegment,
0835             spherical_equatorial_tag, spherical_equatorial_tag,
0836             Strategy
0837         >::type;
0838 };
0839
0840
0841 } // namespace services
0842 #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
0843
0844 }} // namespace strategy::distance
0845
0846 }} // namespace boost::geometry
0847
0848 #endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_CROSS_TRACK_HPP