1 // Boost.Geometry (aka GGL, Generic Geometry Library)
2 
3 // Copyright (c) 2015 Barend Gehrels, Amsterdam, the Netherlands.
4 
5 // This file was modified by Oracle on 2015, 2017.
6 // Modifications copyright (c) 2015-2017 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_DIRECITON_CODE_HPP
16 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_DIRECITON_CODE_HPP
17 
18 
19 #include <boost/geometry/core/access.hpp>
20 #include <boost/geometry/util/math.hpp>
21 #include <boost/geometry/util/select_coordinate_type.hpp>
22 #include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>
23 
24 #include <boost/mpl/assert.hpp>
25 
26 
27 namespace boost { namespace geometry
28 {
29 
30 
31 #ifndef DOXYGEN_NO_DETAIL
32 namespace detail
33 {
34 
35 
36 // TODO: remove
37 template <std::size_t Index, typename Point1, typename Point2>
sign_of_difference(Point1 const & point1,Point2 const & point2)38 inline int sign_of_difference(Point1 const& point1, Point2 const& point2)
39 {
40     return
41         math::equals(geometry::get<Index>(point1), geometry::get<Index>(point2))
42         ?
43         0
44         :
45         (geometry::get<Index>(point1) > geometry::get<Index>(point2) ? 1 : -1);
46 }
47 
48 
49 template <typename Point, typename CSTag = typename cs_tag<Point>::type>
50 struct direction_code_impl
51 {
52     BOOST_MPL_ASSERT_MSG((false), NOT_IMPLEMENTED_FOR_THIS_CS, (CSTag));
53 };
54 
55 template <typename Point>
56 struct direction_code_impl<Point, cartesian_tag>
57 {
58     template <typename Point1, typename Point2>
applyboost::geometry::detail::direction_code_impl59     static inline int apply(Point1 const& segment_a, Point1 const& segment_b,
60                             Point2 const& p)
61     {
62         typedef typename geometry::select_coordinate_type
63             <
64                 Point1, Point2
65             >::type calc_t;
66 
67         if ( (math::equals(geometry::get<0>(segment_b), geometry::get<0>(segment_a))
68            && math::equals(geometry::get<1>(segment_b), geometry::get<1>(segment_a)))
69           || (math::equals(geometry::get<0>(segment_b), geometry::get<0>(p))
70            && math::equals(geometry::get<1>(segment_b), geometry::get<1>(p))) )
71         {
72             return 0;
73         }
74 
75         calc_t x1 = geometry::get<0>(segment_b) - geometry::get<0>(segment_a);
76         calc_t y1 = geometry::get<1>(segment_b) - geometry::get<1>(segment_a);
77         calc_t x2 = geometry::get<0>(segment_b) - geometry::get<0>(p);
78         calc_t y2 = geometry::get<1>(segment_b) - geometry::get<1>(p);
79 
80         calc_t ax = (std::min)(math::abs(x1), math::abs(x2));
81         calc_t ay = (std::min)(math::abs(y1), math::abs(y2));
82 
83         int s1 = 0, s2 = 0;
84         if (ax >= ay)
85         {
86             s1 = x1 > 0 ? 1 : -1;
87             s2 = x2 > 0 ? 1 : -1;
88         }
89         else
90         {
91             s1 = y1 > 0 ? 1 : -1;
92             s2 = y2 > 0 ? 1 : -1;
93         }
94 
95         return s1 == s2 ? -1 : 1;
96     }
97 };
98 
99 template <typename Point>
100 struct direction_code_impl<Point, spherical_equatorial_tag>
101 {
102     template <typename Point1, typename Point2>
applyboost::geometry::detail::direction_code_impl103     static inline int apply(Point1 const& segment_a, Point1 const& segment_b,
104                             Point2 const& p)
105     {
106         typedef typename coordinate_type<Point1>::type coord1_t;
107         typedef typename coordinate_type<Point2>::type coord2_t;
108         typedef typename coordinate_system<Point1>::type::units units_t;
109         typedef typename coordinate_system<Point2>::type::units units2_t;
110         BOOST_MPL_ASSERT_MSG((boost::is_same<units_t, units2_t>::value),
111                              NOT_IMPLEMENTED_FOR_DIFFERENT_UNITS,
112                              (units_t, units2_t));
113 
114         typedef typename geometry::select_coordinate_type <Point1, Point2>::type calc_t;
115         typedef math::detail::constants_on_spheroid<coord1_t, units_t> constants1;
116         typedef math::detail::constants_on_spheroid<coord2_t, units_t> constants2;
117         typedef math::detail::constants_on_spheroid<calc_t, units_t> constants;
118 
119         coord1_t const a0 = geometry::get<0>(segment_a);
120         coord1_t const a1 = geometry::get<1>(segment_a);
121         coord1_t const b0 = geometry::get<0>(segment_b);
122         coord1_t const b1 = geometry::get<1>(segment_b);
123         coord2_t const p0 = geometry::get<0>(p);
124         coord2_t const p1 = geometry::get<1>(p);
125         coord1_t const pi_half1 = constants1::max_latitude();
126         coord2_t const pi_half2 = constants2::max_latitude();
127         calc_t const pi = constants::half_period();
128         calc_t const pi_half = constants::max_latitude();
129         calc_t const c0 = 0;
130 
131         if ( (math::equals(b0, a0) && math::equals(b1, a1))
132           || (math::equals(b0, p0) && math::equals(b1, p1)) )
133         {
134             return 0;
135         }
136 
137         bool const is_a_pole = math::equals(pi_half1, math::abs(a1));
138         bool const is_b_pole = math::equals(pi_half1, math::abs(b1));
139         bool const is_p_pole = math::equals(pi_half2, math::abs(p1));
140 
141         if ( is_b_pole && ((is_a_pole && math::sign(b1) == math::sign(a1))
142                         || (is_p_pole && math::sign(b1) == math::sign(p1))) )
143         {
144             return 0;
145         }
146 
147         // NOTE: as opposed to the implementation for cartesian CS
148         // here point b is the origin
149 
150         calc_t const dlon1 = math::longitude_distance_signed<units_t>(b0, a0);
151         calc_t const dlon2 = math::longitude_distance_signed<units_t>(b0, p0);
152 
153         bool is_antilon1 = false, is_antilon2 = false;
154         calc_t const dlat1 = latitude_distance_signed(b1, a1, dlon1, pi, is_antilon1);
155         calc_t const dlat2 = latitude_distance_signed(b1, p1, dlon2, pi, is_antilon2);
156 
157         calc_t mx = is_a_pole || is_b_pole || is_p_pole ?
158                     c0 :
159                     (std::min)(is_antilon1 ? c0 : math::abs(dlon1),
160                                is_antilon2 ? c0 : math::abs(dlon2));
161         calc_t my = (std::min)(math::abs(dlat1),
162                                math::abs(dlat2));
163 
164         int s1 = 0, s2 = 0;
165         if (mx >= my)
166         {
167             s1 = dlon1 > 0 ? 1 : -1;
168             s2 = dlon2 > 0 ? 1 : -1;
169         }
170         else
171         {
172             s1 = dlat1 > 0 ? 1 : -1;
173             s2 = dlat2 > 0 ? 1 : -1;
174         }
175 
176         return s1 == s2 ? -1 : 1;
177     }
178 
179     template <typename T>
latitude_distance_signedboost::geometry::detail::direction_code_impl180     static inline T latitude_distance_signed(T const& lat1, T const& lat2, T const& lon_ds, T const& pi, bool & is_antilon)
181     {
182         T const c0 = 0;
183 
184         T res = lat2 - lat1;
185 
186         is_antilon = math::equals(math::abs(lon_ds), pi);
187         if (is_antilon)
188         {
189             res = lat2 + lat1;
190             if (res >= c0)
191                 res = pi - res;
192             else
193                 res = -pi - res;
194         }
195 
196         return res;
197     }
198 };
199 
200 template <typename Point>
201 struct direction_code_impl<Point, spherical_polar_tag>
202 {
203     template <typename Point1, typename Point2>
applyboost::geometry::detail::direction_code_impl204     static inline int apply(Point1 segment_a, Point1 segment_b,
205                             Point2 p)
206     {
207         typedef math::detail::constants_on_spheroid
208             <
209                 typename coordinate_type<Point1>::type,
210                 typename coordinate_system<Point1>::type::units
211             > constants1;
212         typedef math::detail::constants_on_spheroid
213             <
214                 typename coordinate_type<Point2>::type,
215                 typename coordinate_system<Point2>::type::units
216             > constants2;
217 
218         geometry::set<1>(segment_a,
219             constants1::max_latitude() - geometry::get<1>(segment_a));
220         geometry::set<1>(segment_b,
221             constants1::max_latitude() - geometry::get<1>(segment_b));
222         geometry::set<1>(p,
223             constants2::max_latitude() - geometry::get<1>(p));
224 
225         return direction_code_impl
226                 <
227                     Point, spherical_equatorial_tag
228                 >::apply(segment_a, segment_b, p);
229     }
230 };
231 
232 template <typename Point>
233 struct direction_code_impl<Point, geographic_tag>
234     : direction_code_impl<Point, spherical_equatorial_tag>
235 {};
236 
237 // Gives sense of direction for point p, collinear w.r.t. segment (a,b)
238 // Returns -1 if p goes backward w.r.t (a,b), so goes from b in direction of a
239 // Returns 1 if p goes forward, so extends (a,b)
240 // Returns 0 if p is equal with b, or if (a,b) is degenerate
241 // Note that it does not do any collinearity test, that should be done before
242 template <typename Point1, typename Point2>
direction_code(Point1 const & segment_a,Point1 const & segment_b,Point2 const & p)243 inline int direction_code(Point1 const& segment_a, Point1 const& segment_b,
244                           Point2 const& p)
245 {
246     return direction_code_impl<Point1>::apply(segment_a, segment_b, p);
247 }
248 
249 
250 } // namespace detail
251 #endif //DOXYGEN_NO_DETAIL
252 
253 
254 }} // namespace boost::geometry
255 
256 #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_DIRECITON_CODE_HPP
257