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