1 // Boost.Geometry
2
3 // Copyright (c) 2016, Oracle and/or its affiliates.
4 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
5 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
6
7 // Use, modification and distribution is subject to the Boost Software License,
8 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
9 // http://www.boost.org/LICENSE_1_0.txt)
10
11 #ifndef BOOST_GEOMETRY_FORMULAS_SPHERICAL_HPP
12 #define BOOST_GEOMETRY_FORMULAS_SPHERICAL_HPP
13
14 #include <boost/geometry/core/coordinate_system.hpp>
15 #include <boost/geometry/core/coordinate_type.hpp>
16 #include <boost/geometry/core/access.hpp>
17 #include <boost/geometry/core/radian_access.hpp>
18
19 //#include <boost/geometry/arithmetic/arithmetic.hpp>
20 #include <boost/geometry/arithmetic/cross_product.hpp>
21 #include <boost/geometry/arithmetic/dot_product.hpp>
22
23 #include <boost/geometry/util/math.hpp>
24 #include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>
25 #include <boost/geometry/util/select_coordinate_type.hpp>
26
27 namespace boost { namespace geometry {
28
29 namespace formula {
30
31 template <typename T>
32 struct result_spherical
33 {
result_sphericalboost::geometry::formula::result_spherical34 result_spherical()
35 : azimuth(0)
36 , reverse_azimuth(0)
37 {}
38
39 T azimuth;
40 T reverse_azimuth;
41 };
42
43 template <typename T>
sph_to_cart3d(T const & lon,T const & lat,T & x,T & y,T & z)44 static inline void sph_to_cart3d(T const& lon, T const& lat, T & x, T & y, T & z)
45 {
46 T const cos_lat = cos(lat);
47 x = cos_lat * cos(lon);
48 y = cos_lat * sin(lon);
49 z = sin(lat);
50 }
51
52 template <typename Point3d, typename PointSph>
sph_to_cart3d(PointSph const & point_sph)53 static inline Point3d sph_to_cart3d(PointSph const& point_sph)
54 {
55 typedef typename coordinate_type<Point3d>::type calc_t;
56
57 calc_t const lon = get_as_radian<0>(point_sph);
58 calc_t const lat = get_as_radian<1>(point_sph);
59 calc_t x, y, z;
60 sph_to_cart3d(lon, lat, x, y, z);
61
62 Point3d res;
63 set<0>(res, x);
64 set<1>(res, y);
65 set<2>(res, z);
66
67 return res;
68 }
69
70 template <typename T>
cart3d_to_sph(T const & x,T const & y,T const & z,T & lon,T & lat)71 static inline void cart3d_to_sph(T const& x, T const& y, T const& z, T & lon, T & lat)
72 {
73 lon = atan2(y, x);
74 lat = asin(z);
75 }
76
77 template <typename PointSph, typename Point3d>
cart3d_to_sph(Point3d const & point_3d)78 static inline PointSph cart3d_to_sph(Point3d const& point_3d)
79 {
80 typedef typename coordinate_type<PointSph>::type coord_t;
81 typedef typename coordinate_type<Point3d>::type calc_t;
82
83 calc_t const x = get<0>(point_3d);
84 calc_t const y = get<1>(point_3d);
85 calc_t const z = get<2>(point_3d);
86 calc_t lonr, latr;
87 cart3d_to_sph(x, y, z, lonr, latr);
88
89 PointSph res;
90 set_from_radian<0>(res, lonr);
91 set_from_radian<1>(res, latr);
92
93 coord_t lon = get<0>(res);
94 coord_t lat = get<1>(res);
95
96 math::normalize_spheroidal_coordinates
97 <
98 typename coordinate_system<PointSph>::type::units,
99 coord_t
100 >(lon, lat);
101
102 set<0>(res, lon);
103 set<1>(res, lat);
104
105 return res;
106 }
107
108 // -1 right
109 // 1 left
110 // 0 on
111 template <typename Point3d1, typename Point3d2>
sph_side_value(Point3d1 const & norm,Point3d2 const & pt)112 static inline int sph_side_value(Point3d1 const& norm, Point3d2 const& pt)
113 {
114 typedef typename select_coordinate_type<Point3d1, Point3d2>::type calc_t;
115 calc_t c0 = 0;
116 calc_t d = dot_product(norm, pt);
117 return math::equals(d, c0) ? 0
118 : d > c0 ? 1
119 : -1; // d < 0
120 }
121
122 template <typename CT, bool ReverseAzimuth, typename T1, typename T2>
spherical_azimuth(T1 const & lon1,T1 const & lat1,T2 const & lon2,T2 const & lat2)123 static inline result_spherical<CT> spherical_azimuth(T1 const& lon1,
124 T1 const& lat1,
125 T2 const& lon2,
126 T2 const& lat2)
127 {
128 typedef result_spherical<CT> result_type;
129 result_type result;
130
131 // http://williams.best.vwh.net/avform.htm#Crs
132 // https://en.wikipedia.org/wiki/Great-circle_navigation
133 CT dlon = lon2 - lon1;
134
135 // An optimization which should kick in often for Boxes
136 //if ( math::equals(dlon, ReturnType(0)) )
137 //if ( get<0>(p1) == get<0>(p2) )
138 //{
139 // return - sin(get_as_radian<1>(p1)) * cos_p2lat);
140 //}
141
142 CT const cos_dlon = cos(dlon);
143 CT const sin_dlon = sin(dlon);
144 CT const cos_lat1 = cos(lat1);
145 CT const cos_lat2 = cos(lat2);
146 CT const sin_lat1 = sin(lat1);
147 CT const sin_lat2 = sin(lat2);
148
149 {
150 // "An alternative formula, not requiring the pre-computation of d"
151 // In the formula below dlon is used as "d"
152 CT const y = sin_dlon * cos_lat2;
153 CT const x = cos_lat1 * sin_lat2 - sin_lat1 * cos_lat2 * cos_dlon;
154 result.azimuth = atan2(y, x);
155 }
156
157 if (ReverseAzimuth)
158 {
159 CT const y = sin_dlon * cos_lat1;
160 CT const x = sin_lat2 * cos_lat1 * cos_dlon - cos_lat2 * sin_lat1;
161 result.reverse_azimuth = atan2(y, x);
162 }
163
164 return result;
165 }
166
167 template <typename ReturnType, typename T1, typename T2>
spherical_azimuth(T1 const & lon1,T1 const & lat1,T2 const & lon2,T2 const & lat2)168 inline ReturnType spherical_azimuth(T1 const& lon1, T1 const& lat1,
169 T2 const& lon2, T2 const& lat2)
170 {
171 return spherical_azimuth<ReturnType, false>(lon1, lat1, lon2, lat2).azimuth;
172 }
173
174 template <typename T>
spherical_azimuth(T const & lon1,T const & lat1,T const & lon2,T const & lat2)175 inline T spherical_azimuth(T const& lon1, T const& lat1, T const& lon2, T const& lat2)
176 {
177 return spherical_azimuth<T, false>(lon1, lat1, lon2, lat2).azimuth;
178 }
179
180 template <typename T>
azimuth_side_value(T const & azi_a1_p,T const & azi_a1_a2)181 inline int azimuth_side_value(T const& azi_a1_p, T const& azi_a1_a2)
182 {
183 T const pi = math::pi<T>();
184 T const two_pi = math::two_pi<T>();
185
186 // instead of the formula from XTD
187 //calc_t a_diff = asin(sin(azi_a1_p - azi_a1_a2));
188
189 T a_diff = azi_a1_p - azi_a1_a2;
190 // normalize, angle in [-pi, pi]
191 while (a_diff > pi)
192 a_diff -= two_pi;
193 while (a_diff < -pi)
194 a_diff += two_pi;
195
196 // NOTE: in general it shouldn't be required to support the pi/-pi case
197 // because in non-cartesian systems it makes sense to check the side
198 // only "between" the endpoints.
199 // However currently the winding strategy calls the side strategy
200 // for vertical segments to check if the point is "between the endpoints.
201 // This could be avoided since the side strategy is not required for that
202 // because meridian is the shortest path. So a difference of
203 // longitudes would be sufficient (of course normalized to [-pi, pi]).
204
205 // NOTE: with the above said, the pi/-pi check is temporary
206 // however in case if this was required
207 // the geodesics on ellipsoid aren't "symmetrical"
208 // therefore instead of comparing a_diff to pi and -pi
209 // one should probably use inverse azimuths and compare
210 // the difference to 0 as well
211
212 // positive azimuth is on the right side
213 return math::equals(a_diff, 0)
214 || math::equals(a_diff, pi)
215 || math::equals(a_diff, -pi) ? 0
216 : a_diff > 0 ? -1 // right
217 : 1; // left
218 }
219
220 } // namespace formula
221
222 }} // namespace boost::geometry
223
224 #endif // BOOST_GEOMETRY_FORMULAS_SPHERICAL_HPP
225