distance_cross_track.hpp 27 KB

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