ceph/src/boost/boost/geometry/strategies/spherical/distance_cross_track.hpp

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