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