distance_cross_track.hpp 26 KB

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