comparison DEPENDENCIES/generic/include/boost/geometry/strategies/spherical/ssf.hpp @ 101:c530137014c0

Update Boost headers (1.58.0)
author Chris Cannam
date Mon, 07 Sep 2015 11:12:49 +0100
parents 2665513ce2d3
children
comparison
equal deleted inserted replaced
100:793467b5e61c 101:c530137014c0
14 14
15 #include <boost/geometry/core/cs.hpp> 15 #include <boost/geometry/core/cs.hpp>
16 #include <boost/geometry/core/access.hpp> 16 #include <boost/geometry/core/access.hpp>
17 #include <boost/geometry/core/radian_access.hpp> 17 #include <boost/geometry/core/radian_access.hpp>
18 18
19 #include <boost/geometry/util/select_coordinate_type.hpp>
20 #include <boost/geometry/util/math.hpp> 19 #include <boost/geometry/util/math.hpp>
20 #include <boost/geometry/util/promote_floating_point.hpp>
21 #include <boost/geometry/util/select_calculation_type.hpp>
21 22
22 #include <boost/geometry/strategies/side.hpp> 23 #include <boost/geometry/strategies/side.hpp>
23 //#include <boost/geometry/strategies/concepts/side_concept.hpp> 24 //#include <boost/geometry/strategies/concepts/side_concept.hpp>
24 25
25 26
28 29
29 30
30 namespace strategy { namespace side 31 namespace strategy { namespace side
31 { 32 {
32 33
34 #ifndef DOXYGEN_NO_DETAIL
35 namespace detail
36 {
37
38 template <typename T>
39 int spherical_side_formula(T const& lambda1, T const& delta1,
40 T const& lambda2, T const& delta2,
41 T const& lambda, T const& delta)
42 {
43 // Create temporary points (vectors) on unit a sphere
44 T const cos_delta1 = cos(delta1);
45 T const c1x = cos_delta1 * cos(lambda1);
46 T const c1y = cos_delta1 * sin(lambda1);
47 T const c1z = sin(delta1);
48
49 T const cos_delta2 = cos(delta2);
50 T const c2x = cos_delta2 * cos(lambda2);
51 T const c2y = cos_delta2 * sin(lambda2);
52 T const c2z = sin(delta2);
53
54 // (Third point is converted directly)
55 T const cos_delta = cos(delta);
56
57 // Apply the "Spherical Side Formula" as presented on my blog
58 T const dist
59 = (c1y * c2z - c1z * c2y) * cos_delta * cos(lambda)
60 + (c1z * c2x - c1x * c2z) * cos_delta * sin(lambda)
61 + (c1x * c2y - c1y * c2x) * sin(delta);
62
63 T zero = T();
64 return dist > zero ? 1
65 : dist < zero ? -1
66 : 0;
67 }
68
69 }
70 #endif // DOXYGEN_NO_DETAIL
33 71
34 /*! 72 /*!
35 \brief Check at which side of a Great Circle segment a point lies 73 \brief Check at which side of a Great Circle segment a point lies
36 left of segment (> 0), right of segment (< 0), on segment (0) 74 left of segment (> 0), right of segment (< 0), on segment (0)
37 \ingroup strategies 75 \ingroup strategies
43 81
44 public : 82 public :
45 template <typename P1, typename P2, typename P> 83 template <typename P1, typename P2, typename P>
46 static inline int apply(P1 const& p1, P2 const& p2, P const& p) 84 static inline int apply(P1 const& p1, P2 const& p2, P const& p)
47 { 85 {
48 typedef typename boost::mpl::if_c 86 typedef typename promote_floating_point
49 < 87 <
50 boost::is_void<CalculationType>::type::value, 88 typename select_calculation_type_alt
89 <
90 CalculationType,
91 P1, P2, P
92 >::type
93 >::type calculation_type;
51 94
52 // Select at least a double... 95 calculation_type const lambda1 = get_as_radian<0>(p1);
53 typename select_most_precise 96 calculation_type const delta1 = get_as_radian<1>(p1);
54 < 97 calculation_type const lambda2 = get_as_radian<0>(p2);
55 typename select_most_precise 98 calculation_type const delta2 = get_as_radian<1>(p2);
56 < 99 calculation_type const lambda = get_as_radian<0>(p);
57 typename select_most_precise 100 calculation_type const delta = get_as_radian<1>(p);
58 <
59 typename coordinate_type<P1>::type,
60 typename coordinate_type<P2>::type
61 >::type,
62 typename coordinate_type<P>::type
63 >::type,
64 double
65 >::type,
66 CalculationType
67 >::type coordinate_type;
68 101
69 // Convenient shortcuts 102 return detail::spherical_side_formula(lambda1, delta1,
70 typedef coordinate_type ct; 103 lambda2, delta2,
71 ct const lambda1 = get_as_radian<0>(p1); 104 lambda, delta);
72 ct const delta1 = get_as_radian<1>(p1);
73 ct const lambda2 = get_as_radian<0>(p2);
74 ct const delta2 = get_as_radian<1>(p2);
75 ct const lambda = get_as_radian<0>(p);
76 ct const delta = get_as_radian<1>(p);
77
78 // Create temporary points (vectors) on unit a sphere
79 ct const cos_delta1 = cos(delta1);
80 ct const c1x = cos_delta1 * cos(lambda1);
81 ct const c1y = cos_delta1 * sin(lambda1);
82 ct const c1z = sin(delta1);
83
84 ct const cos_delta2 = cos(delta2);
85 ct const c2x = cos_delta2 * cos(lambda2);
86 ct const c2y = cos_delta2 * sin(lambda2);
87 ct const c2z = sin(delta2);
88
89 // (Third point is converted directly)
90 ct const cos_delta = cos(delta);
91
92 // Apply the "Spherical Side Formula" as presented on my blog
93 ct const dist
94 = (c1y * c2z - c1z * c2y) * cos_delta * cos(lambda)
95 + (c1z * c2x - c1x * c2z) * cos_delta * sin(lambda)
96 + (c1x * c2y - c1y * c2x) * sin(delta);
97
98 ct zero = ct();
99 return dist > zero ? 1
100 : dist < zero ? -1
101 : 0;
102 } 105 }
103 }; 106 };
104 107
105 108
106 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS 109 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS