Chris@16
|
1 // Boost.Geometry (aka GGL, Generic Geometry Library)
|
Chris@16
|
2
|
Chris@16
|
3 // Copyright (c) 2011-2012 Barend Gehrels, Amsterdam, the Netherlands.
|
Chris@16
|
4
|
Chris@16
|
5 // Use, modification and distribution is subject to the Boost Software License,
|
Chris@16
|
6 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
|
Chris@16
|
7 // http://www.boost.org/LICENSE_1_0.txt)
|
Chris@16
|
8
|
Chris@16
|
9 #ifndef BOOST_GEOMETRY_STRATEGIES_SPHERICAL_SSF_HPP
|
Chris@16
|
10 #define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_SSF_HPP
|
Chris@16
|
11
|
Chris@16
|
12 #include <boost/mpl/if.hpp>
|
Chris@16
|
13 #include <boost/type_traits.hpp>
|
Chris@16
|
14
|
Chris@16
|
15 #include <boost/geometry/core/cs.hpp>
|
Chris@16
|
16 #include <boost/geometry/core/access.hpp>
|
Chris@16
|
17 #include <boost/geometry/core/radian_access.hpp>
|
Chris@16
|
18
|
Chris@16
|
19 #include <boost/geometry/util/select_coordinate_type.hpp>
|
Chris@16
|
20 #include <boost/geometry/util/math.hpp>
|
Chris@16
|
21
|
Chris@16
|
22 #include <boost/geometry/strategies/side.hpp>
|
Chris@16
|
23 //#include <boost/geometry/strategies/concepts/side_concept.hpp>
|
Chris@16
|
24
|
Chris@16
|
25
|
Chris@16
|
26 namespace boost { namespace geometry
|
Chris@16
|
27 {
|
Chris@16
|
28
|
Chris@16
|
29
|
Chris@16
|
30 namespace strategy { namespace side
|
Chris@16
|
31 {
|
Chris@16
|
32
|
Chris@16
|
33
|
Chris@16
|
34 /*!
|
Chris@16
|
35 \brief Check at which side of a Great Circle segment a point lies
|
Chris@16
|
36 left of segment (> 0), right of segment (< 0), on segment (0)
|
Chris@16
|
37 \ingroup strategies
|
Chris@16
|
38 \tparam CalculationType \tparam_calculation
|
Chris@16
|
39 */
|
Chris@16
|
40 template <typename CalculationType = void>
|
Chris@16
|
41 class spherical_side_formula
|
Chris@16
|
42 {
|
Chris@16
|
43
|
Chris@16
|
44 public :
|
Chris@16
|
45 template <typename P1, typename P2, typename P>
|
Chris@16
|
46 static inline int apply(P1 const& p1, P2 const& p2, P const& p)
|
Chris@16
|
47 {
|
Chris@16
|
48 typedef typename boost::mpl::if_c
|
Chris@16
|
49 <
|
Chris@16
|
50 boost::is_void<CalculationType>::type::value,
|
Chris@16
|
51
|
Chris@16
|
52 // Select at least a double...
|
Chris@16
|
53 typename select_most_precise
|
Chris@16
|
54 <
|
Chris@16
|
55 typename select_most_precise
|
Chris@16
|
56 <
|
Chris@16
|
57 typename select_most_precise
|
Chris@16
|
58 <
|
Chris@16
|
59 typename coordinate_type<P1>::type,
|
Chris@16
|
60 typename coordinate_type<P2>::type
|
Chris@16
|
61 >::type,
|
Chris@16
|
62 typename coordinate_type<P>::type
|
Chris@16
|
63 >::type,
|
Chris@16
|
64 double
|
Chris@16
|
65 >::type,
|
Chris@16
|
66 CalculationType
|
Chris@16
|
67 >::type coordinate_type;
|
Chris@16
|
68
|
Chris@16
|
69 // Convenient shortcuts
|
Chris@16
|
70 typedef coordinate_type ct;
|
Chris@16
|
71 ct const lambda1 = get_as_radian<0>(p1);
|
Chris@16
|
72 ct const delta1 = get_as_radian<1>(p1);
|
Chris@16
|
73 ct const lambda2 = get_as_radian<0>(p2);
|
Chris@16
|
74 ct const delta2 = get_as_radian<1>(p2);
|
Chris@16
|
75 ct const lambda = get_as_radian<0>(p);
|
Chris@16
|
76 ct const delta = get_as_radian<1>(p);
|
Chris@16
|
77
|
Chris@16
|
78 // Create temporary points (vectors) on unit a sphere
|
Chris@16
|
79 ct const cos_delta1 = cos(delta1);
|
Chris@16
|
80 ct const c1x = cos_delta1 * cos(lambda1);
|
Chris@16
|
81 ct const c1y = cos_delta1 * sin(lambda1);
|
Chris@16
|
82 ct const c1z = sin(delta1);
|
Chris@16
|
83
|
Chris@16
|
84 ct const cos_delta2 = cos(delta2);
|
Chris@16
|
85 ct const c2x = cos_delta2 * cos(lambda2);
|
Chris@16
|
86 ct const c2y = cos_delta2 * sin(lambda2);
|
Chris@16
|
87 ct const c2z = sin(delta2);
|
Chris@16
|
88
|
Chris@16
|
89 // (Third point is converted directly)
|
Chris@16
|
90 ct const cos_delta = cos(delta);
|
Chris@16
|
91
|
Chris@16
|
92 // Apply the "Spherical Side Formula" as presented on my blog
|
Chris@16
|
93 ct const dist
|
Chris@16
|
94 = (c1y * c2z - c1z * c2y) * cos_delta * cos(lambda)
|
Chris@16
|
95 + (c1z * c2x - c1x * c2z) * cos_delta * sin(lambda)
|
Chris@16
|
96 + (c1x * c2y - c1y * c2x) * sin(delta);
|
Chris@16
|
97
|
Chris@16
|
98 ct zero = ct();
|
Chris@16
|
99 return dist > zero ? 1
|
Chris@16
|
100 : dist < zero ? -1
|
Chris@16
|
101 : 0;
|
Chris@16
|
102 }
|
Chris@16
|
103 };
|
Chris@16
|
104
|
Chris@16
|
105
|
Chris@16
|
106 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
|
Chris@16
|
107 namespace services
|
Chris@16
|
108 {
|
Chris@16
|
109
|
Chris@16
|
110 /*template <typename CalculationType>
|
Chris@16
|
111 struct default_strategy<spherical_polar_tag, CalculationType>
|
Chris@16
|
112 {
|
Chris@16
|
113 typedef spherical_side_formula<CalculationType> type;
|
Chris@16
|
114 };*/
|
Chris@16
|
115
|
Chris@16
|
116 template <typename CalculationType>
|
Chris@16
|
117 struct default_strategy<spherical_equatorial_tag, CalculationType>
|
Chris@16
|
118 {
|
Chris@16
|
119 typedef spherical_side_formula<CalculationType> type;
|
Chris@16
|
120 };
|
Chris@16
|
121
|
Chris@16
|
122 template <typename CalculationType>
|
Chris@16
|
123 struct default_strategy<geographic_tag, CalculationType>
|
Chris@16
|
124 {
|
Chris@16
|
125 typedef spherical_side_formula<CalculationType> type;
|
Chris@16
|
126 };
|
Chris@16
|
127
|
Chris@16
|
128 }
|
Chris@16
|
129 #endif
|
Chris@16
|
130
|
Chris@16
|
131 }} // namespace strategy::side
|
Chris@16
|
132
|
Chris@16
|
133 }} // namespace boost::geometry
|
Chris@16
|
134
|
Chris@16
|
135
|
Chris@16
|
136 #endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_SSF_HPP
|