1 // Boost.Geometry (aka GGL, Generic Geometry Library)
2
3 // Copyright (c) 2015-2020 Barend Gehrels, Amsterdam, the Netherlands.
4
5 // This file was modified by Oracle on 2015, 2017, 2019.
6 // Modifications copyright (c) 2015-2019 Oracle and/or its affiliates.
7
8 // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
9 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
10
11 // Use, modification and distribution is subject to the Boost Software License,
12 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
13 // http://www.boost.org/LICENSE_1_0.txt)
14
15 #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_DIRECTION_CODE_HPP
16 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_DIRECTION_CODE_HPP
17
18
19 #include <boost/geometry/core/access.hpp>
20 #include <boost/geometry/arithmetic/infinite_line_functions.hpp>
21 #include <boost/geometry/algorithms/detail/make/make.hpp>
22 #include <boost/geometry/util/math.hpp>
23 #include <boost/geometry/util/select_coordinate_type.hpp>
24 #include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>
25
26 #include <boost/mpl/assert.hpp>
27
28
29 namespace boost { namespace geometry
30 {
31
32
33 #ifndef DOXYGEN_NO_DETAIL
34 namespace detail
35 {
36
37 template <typename CSTag>
38 struct direction_code_impl
39 {
40 BOOST_MPL_ASSERT_MSG((false), NOT_IMPLEMENTED_FOR_THIS_CS, (CSTag));
41 };
42
43 template <>
44 struct direction_code_impl<cartesian_tag>
45 {
46 template <typename Point1, typename Point2>
applyboost::geometry::detail::direction_code_impl47 static inline int apply(Point1 const& segment_a, Point1 const& segment_b,
48 Point2 const& point)
49 {
50 typedef typename geometry::select_coordinate_type
51 <
52 Point1, Point2
53 >::type calc_t;
54
55 typedef model::infinite_line<calc_t> line_type;
56
57 // Situation and construction of perpendicular line
58 //
59 // P1 a--------------->b P2
60 // |
61 // |
62 // v
63 //
64 // P1 is located right of the (directional) perpendicular line
65 // and therefore gets a negative side_value, and returns -1.
66 // P2 is to the left of the perpendicular line and returns 1.
67 // If the specified point is located on top of b, it returns 0.
68
69 line_type const line
70 = detail::make::make_perpendicular_line<calc_t>(segment_a,
71 segment_b, segment_b);
72
73 if (arithmetic::is_degenerate(line))
74 {
75 return 0;
76 }
77
78 calc_t const sv = arithmetic::side_value(line, point);
79 return sv == 0 ? 0 : sv > 0 ? 1 : -1;
80 }
81 };
82
83 template <>
84 struct direction_code_impl<spherical_equatorial_tag>
85 {
86 template <typename Point1, typename Point2>
applyboost::geometry::detail::direction_code_impl87 static inline int apply(Point1 const& segment_a, Point1 const& segment_b,
88 Point2 const& p)
89 {
90 typedef typename coordinate_type<Point1>::type coord1_t;
91 typedef typename coordinate_type<Point2>::type coord2_t;
92 typedef typename cs_angular_units<Point1>::type units_t;
93 typedef typename cs_angular_units<Point2>::type units2_t;
94 BOOST_MPL_ASSERT_MSG((boost::is_same<units_t, units2_t>::value),
95 NOT_IMPLEMENTED_FOR_DIFFERENT_UNITS,
96 (units_t, units2_t));
97
98 typedef typename geometry::select_coordinate_type <Point1, Point2>::type calc_t;
99 typedef math::detail::constants_on_spheroid<coord1_t, units_t> constants1;
100 typedef math::detail::constants_on_spheroid<coord2_t, units_t> constants2;
101 static coord1_t const pi_half1 = constants1::max_latitude();
102 static coord2_t const pi_half2 = constants2::max_latitude();
103 static calc_t const c0 = 0;
104
105 coord1_t const a0 = geometry::get<0>(segment_a);
106 coord1_t const a1 = geometry::get<1>(segment_a);
107 coord1_t const b0 = geometry::get<0>(segment_b);
108 coord1_t const b1 = geometry::get<1>(segment_b);
109 coord2_t const p0 = geometry::get<0>(p);
110 coord2_t const p1 = geometry::get<1>(p);
111
112 if ( (math::equals(b0, a0) && math::equals(b1, a1))
113 || (math::equals(b0, p0) && math::equals(b1, p1)) )
114 {
115 return 0;
116 }
117
118 bool const is_a_pole = math::equals(pi_half1, math::abs(a1));
119 bool const is_b_pole = math::equals(pi_half1, math::abs(b1));
120 bool const is_p_pole = math::equals(pi_half2, math::abs(p1));
121
122 if ( is_b_pole && ((is_a_pole && math::sign(b1) == math::sign(a1))
123 || (is_p_pole && math::sign(b1) == math::sign(p1))) )
124 {
125 return 0;
126 }
127
128 // NOTE: as opposed to the implementation for cartesian CS
129 // here point b is the origin
130
131 calc_t const dlon1 = math::longitude_distance_signed<units_t, calc_t>(b0, a0);
132 calc_t const dlon2 = math::longitude_distance_signed<units_t, calc_t>(b0, p0);
133
134 bool is_antilon1 = false, is_antilon2 = false;
135 calc_t const dlat1 = latitude_distance_signed<units_t, calc_t>(b1, a1, dlon1, is_antilon1);
136 calc_t const dlat2 = latitude_distance_signed<units_t, calc_t>(b1, p1, dlon2, is_antilon2);
137
138 calc_t mx = is_a_pole || is_b_pole || is_p_pole ?
139 c0 :
140 (std::min)(is_antilon1 ? c0 : math::abs(dlon1),
141 is_antilon2 ? c0 : math::abs(dlon2));
142 calc_t my = (std::min)(math::abs(dlat1),
143 math::abs(dlat2));
144
145 int s1 = 0, s2 = 0;
146 if (mx >= my)
147 {
148 s1 = dlon1 > 0 ? 1 : -1;
149 s2 = dlon2 > 0 ? 1 : -1;
150 }
151 else
152 {
153 s1 = dlat1 > 0 ? 1 : -1;
154 s2 = dlat2 > 0 ? 1 : -1;
155 }
156
157 return s1 == s2 ? -1 : 1;
158 }
159
160 template <typename Units, typename T>
latitude_distance_signedboost::geometry::detail::direction_code_impl161 static inline T latitude_distance_signed(T const& lat1, T const& lat2, T const& lon_ds, bool & is_antilon)
162 {
163 typedef math::detail::constants_on_spheroid<T, Units> constants;
164 static T const pi = constants::half_period();
165 static T const c0 = 0;
166
167 T res = lat2 - lat1;
168
169 is_antilon = math::equals(math::abs(lon_ds), pi);
170 if (is_antilon)
171 {
172 res = lat2 + lat1;
173 if (res >= c0)
174 res = pi - res;
175 else
176 res = -pi - res;
177 }
178
179 return res;
180 }
181 };
182
183 template <>
184 struct direction_code_impl<spherical_polar_tag>
185 {
186 template <typename Point1, typename Point2>
applyboost::geometry::detail::direction_code_impl187 static inline int apply(Point1 segment_a, Point1 segment_b,
188 Point2 p)
189 {
190 typedef math::detail::constants_on_spheroid
191 <
192 typename coordinate_type<Point1>::type,
193 typename cs_angular_units<Point1>::type
194 > constants1;
195 typedef math::detail::constants_on_spheroid
196 <
197 typename coordinate_type<Point2>::type,
198 typename cs_angular_units<Point2>::type
199 > constants2;
200
201 geometry::set<1>(segment_a,
202 constants1::max_latitude() - geometry::get<1>(segment_a));
203 geometry::set<1>(segment_b,
204 constants1::max_latitude() - geometry::get<1>(segment_b));
205 geometry::set<1>(p,
206 constants2::max_latitude() - geometry::get<1>(p));
207
208 return direction_code_impl
209 <
210 spherical_equatorial_tag
211 >::apply(segment_a, segment_b, p);
212 }
213 };
214
215 // if spherical_tag is passed then pick cs_tag based on Point1 type
216 // with spherical_equatorial_tag as the default
217 template <>
218 struct direction_code_impl<spherical_tag>
219 {
220 template <typename Point1, typename Point2>
applyboost::geometry::detail::direction_code_impl221 static inline int apply(Point1 segment_a, Point1 segment_b,
222 Point2 p)
223 {
224 return direction_code_impl
225 <
226 typename boost::mpl::if_c
227 <
228 boost::is_same
229 <
230 typename geometry::cs_tag<Point1>::type,
231 spherical_polar_tag
232 >::value,
233 spherical_polar_tag,
234 spherical_equatorial_tag
235 >::type
236 >::apply(segment_a, segment_b, p);
237 }
238 };
239
240 template <>
241 struct direction_code_impl<geographic_tag>
242 : direction_code_impl<spherical_equatorial_tag>
243 {};
244
245 // Gives sense of direction for point p, collinear w.r.t. segment (a,b)
246 // Returns -1 if p goes backward w.r.t (a,b), so goes from b in direction of a
247 // Returns 1 if p goes forward, so extends (a,b)
248 // Returns 0 if p is equal with b, or if (a,b) is degenerate
249 // Note that it does not do any collinearity test, that should be done before
250 template <typename CSTag, typename Point1, typename Point2>
direction_code(Point1 const & segment_a,Point1 const & segment_b,Point2 const & p)251 inline int direction_code(Point1 const& segment_a, Point1 const& segment_b,
252 Point2 const& p)
253 {
254 return direction_code_impl<CSTag>::apply(segment_a, segment_b, p);
255 }
256
257
258 } // namespace detail
259 #endif //DOXYGEN_NO_DETAIL
260
261
262 }} // namespace boost::geometry
263
264 #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_DIRECTION_CODE_HPP
265