]>
Commit | Line | Data |
---|---|---|
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. | |
6 | // Modifications copyright (c) 2015, Oracle and/or its affiliates. | |
7 | ||
8 | // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle | |
9 | ||
10 | // Use, modification and distribution is subject to the Boost Software License, | |
11 | // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at | |
12 | // http://www.boost.org/LICENSE_1_0.txt) | |
13 | ||
14 | #ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP | |
15 | #define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP | |
16 | ||
17 | ||
18 | #include <limits> | |
19 | ||
20 | #include <boost/core/ignore_unused.hpp> | |
21 | #include <boost/type_traits/is_integral.hpp> | |
22 | #include <boost/type_traits/make_unsigned.hpp> | |
23 | ||
24 | #include <boost/geometry/arithmetic/determinant.hpp> | |
25 | #include <boost/geometry/core/access.hpp> | |
26 | #include <boost/geometry/core/assert.hpp> | |
27 | #include <boost/geometry/core/coordinate_type.hpp> | |
28 | #include <boost/geometry/algorithms/detail/assign_indexed_point.hpp> | |
29 | #include <boost/geometry/strategies/cartesian/side_by_triangle.hpp> | |
30 | #include <boost/geometry/util/math.hpp> | |
31 | ||
32 | #ifdef BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG | |
33 | #include <boost/math/common_factor_ct.hpp> | |
34 | #include <boost/math/common_factor_rt.hpp> | |
35 | #include <boost/multiprecision/cpp_int.hpp> | |
36 | #endif | |
37 | ||
38 | namespace boost { namespace geometry | |
39 | { | |
40 | ||
41 | namespace strategy { namespace side | |
42 | { | |
43 | ||
44 | namespace detail | |
45 | { | |
46 | ||
47 | // A tool for multiplication of integers avoiding overflow | |
48 | // It's a temporary workaround until we can use Multiprecision | |
49 | // The algorithm is based on Karatsuba algorithm | |
50 | // see: http://en.wikipedia.org/wiki/Karatsuba_algorithm | |
51 | template <typename T> | |
52 | struct multiplicable_integral | |
53 | { | |
54 | // Currently this tool can't be used with non-integral coordinate types. | |
55 | // Also side_of_intersection strategy sign_of_product() and sign_of_compare() | |
56 | // functions would have to be modified to properly support floating-point | |
57 | // types (comparisons and multiplication). | |
58 | BOOST_STATIC_ASSERT(boost::is_integral<T>::value); | |
59 | ||
60 | static const std::size_t bits = CHAR_BIT * sizeof(T); | |
61 | static const std::size_t half_bits = bits / 2; | |
62 | typedef typename boost::make_unsigned<T>::type unsigned_type; | |
63 | static const unsigned_type base = unsigned_type(1) << half_bits; // 2^half_bits | |
64 | ||
65 | int m_sign; | |
66 | unsigned_type m_ms; | |
67 | unsigned_type m_ls; | |
68 | ||
69 | multiplicable_integral(int sign, unsigned_type ms, unsigned_type ls) | |
70 | : m_sign(sign), m_ms(ms), m_ls(ls) | |
71 | {} | |
72 | ||
73 | explicit multiplicable_integral(T const& val) | |
74 | { | |
75 | unsigned_type val_u = val > 0 ? | |
76 | unsigned_type(val) | |
77 | : val == (std::numeric_limits<T>::min)() ? | |
78 | unsigned_type((std::numeric_limits<T>::max)()) + 1 | |
79 | : unsigned_type(-val); | |
80 | // MMLL -> S 00MM 00LL | |
81 | m_sign = math::sign(val); | |
82 | m_ms = val_u >> half_bits; // val_u / base | |
83 | m_ls = val_u - m_ms * base; | |
84 | } | |
85 | ||
86 | friend multiplicable_integral operator*(multiplicable_integral const& a, | |
87 | multiplicable_integral const& b) | |
88 | { | |
89 | // (S 00MM 00LL) * (S 00MM 00LL) -> (S Z2MM 00LL) | |
90 | unsigned_type z2 = a.m_ms * b.m_ms; | |
91 | unsigned_type z0 = a.m_ls * b.m_ls; | |
92 | unsigned_type z1 = (a.m_ms + a.m_ls) * (b.m_ms + b.m_ls) - z2 - z0; | |
93 | // z0 may be >= base so it must be normalized to allow comparison | |
94 | unsigned_type z0_ms = z0 >> half_bits; // z0 / base | |
95 | return multiplicable_integral(a.m_sign * b.m_sign, | |
96 | z2 * base + z1 + z0_ms, | |
97 | z0 - base * z0_ms); | |
98 | } | |
99 | ||
100 | friend bool operator<(multiplicable_integral const& a, | |
101 | multiplicable_integral const& b) | |
102 | { | |
103 | if ( a.m_sign == b.m_sign ) | |
104 | { | |
105 | bool u_less = a.m_ms < b.m_ms | |
106 | || (a.m_ms == b.m_ms && a.m_ls < b.m_ls); | |
107 | return a.m_sign > 0 ? u_less : (! u_less); | |
108 | } | |
109 | else | |
110 | { | |
111 | return a.m_sign < b.m_sign; | |
112 | } | |
113 | } | |
114 | ||
115 | friend bool operator>(multiplicable_integral const& a, | |
116 | multiplicable_integral const& b) | |
117 | { | |
118 | return b < a; | |
119 | } | |
120 | ||
121 | template <typename CmpVal> | |
122 | void check_value(CmpVal const& cmp_val) const | |
123 | { | |
124 | unsigned_type b = base; // a workaround for MinGW - undefined reference base | |
125 | CmpVal val = CmpVal(m_sign) * (CmpVal(m_ms) * CmpVal(b) + CmpVal(m_ls)); | |
126 | BOOST_GEOMETRY_ASSERT(cmp_val == val); | |
127 | } | |
128 | }; | |
129 | ||
130 | } // namespace detail | |
131 | ||
132 | // Calculates the side of the intersection-point (if any) of | |
133 | // of segment a//b w.r.t. segment c | |
134 | // This is calculated without (re)calculating the IP itself again and fully | |
135 | // based on integer mathematics; there are no divisions | |
136 | // It can be used for either integer (rescaled) points, and also for FP | |
137 | class side_of_intersection | |
138 | { | |
139 | private : | |
140 | template <typename T, typename U> | |
141 | static inline | |
142 | int sign_of_product(T const& a, U const& b) | |
143 | { | |
144 | return a == 0 || b == 0 ? 0 | |
145 | : a > 0 && b > 0 ? 1 | |
146 | : a < 0 && b < 0 ? 1 | |
147 | : -1; | |
148 | } | |
149 | ||
150 | template <typename T> | |
151 | static inline | |
152 | int sign_of_compare(T const& a, T const& b, T const& c, T const& d) | |
153 | { | |
154 | // Both a*b and c*d are positive | |
155 | // We have to judge if a*b > c*d | |
156 | ||
157 | using side::detail::multiplicable_integral; | |
158 | multiplicable_integral<T> ab = multiplicable_integral<T>(a) | |
159 | * multiplicable_integral<T>(b); | |
160 | multiplicable_integral<T> cd = multiplicable_integral<T>(c) | |
161 | * multiplicable_integral<T>(d); | |
162 | ||
163 | int result = ab > cd ? 1 | |
164 | : ab < cd ? -1 | |
165 | : 0 | |
166 | ; | |
167 | ||
168 | #ifdef BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG | |
169 | using namespace boost::multiprecision; | |
170 | cpp_int const lab = cpp_int(a) * cpp_int(b); | |
171 | cpp_int const lcd = cpp_int(c) * cpp_int(d); | |
172 | ||
173 | ab.check_value(lab); | |
174 | cd.check_value(lcd); | |
175 | ||
176 | int result2 = lab > lcd ? 1 | |
177 | : lab < lcd ? -1 | |
178 | : 0 | |
179 | ; | |
180 | BOOST_GEOMETRY_ASSERT(result == result2); | |
181 | #endif | |
182 | ||
183 | return result; | |
184 | } | |
185 | ||
186 | template <typename T> | |
187 | static inline | |
188 | int sign_of_addition_of_two_products(T const& a, T const& b, T const& c, T const& d) | |
189 | { | |
190 | // sign of a*b+c*d, 1 if positive, -1 if negative, else 0 | |
191 | int const ab = sign_of_product(a, b); | |
192 | int const cd = sign_of_product(c, d); | |
193 | if (ab == 0) | |
194 | { | |
195 | return cd; | |
196 | } | |
197 | if (cd == 0) | |
198 | { | |
199 | return ab; | |
200 | } | |
201 | ||
202 | if (ab == cd) | |
203 | { | |
204 | // Both positive or both negative | |
205 | return ab; | |
206 | } | |
207 | ||
208 | // One is positive, one is negative, both are non zero | |
209 | // If ab is positive, we have to judge if a*b > -c*d (then 1 because sum is positive) | |
210 | // If ab is negative, we have to judge if c*d > -a*b (idem) | |
211 | return ab == 1 | |
212 | ? sign_of_compare(a, b, -c, d) | |
213 | : sign_of_compare(c, d, -a, b); | |
214 | } | |
215 | ||
216 | ||
217 | public : | |
218 | ||
219 | // Calculates the side of the intersection-point (if any) of | |
220 | // of segment a//b w.r.t. segment c | |
221 | // This is calculated without (re)calculating the IP itself again and fully | |
222 | // based on integer mathematics | |
223 | template <typename T, typename Segment, typename Point> | |
224 | static inline T side_value(Segment const& a, Segment const& b, | |
225 | Segment const& c, Point const& fallback_point) | |
226 | { | |
227 | // The first point of the three segments is reused several times | |
228 | T const ax = get<0, 0>(a); | |
229 | T const ay = get<0, 1>(a); | |
230 | T const bx = get<0, 0>(b); | |
231 | T const by = get<0, 1>(b); | |
232 | T const cx = get<0, 0>(c); | |
233 | T const cy = get<0, 1>(c); | |
234 | ||
235 | T const dx_a = get<1, 0>(a) - ax; | |
236 | T const dy_a = get<1, 1>(a) - ay; | |
237 | ||
238 | T const dx_b = get<1, 0>(b) - bx; | |
239 | T const dy_b = get<1, 1>(b) - by; | |
240 | ||
241 | T const dx_c = get<1, 0>(c) - cx; | |
242 | T const dy_c = get<1, 1>(c) - cy; | |
243 | ||
244 | // Cramer's rule: d (see cart_intersect.hpp) | |
245 | T const d = geometry::detail::determinant<T> | |
246 | ( | |
247 | dx_a, dy_a, | |
248 | dx_b, dy_b | |
249 | ); | |
250 | ||
251 | T const zero = T(); | |
252 | if (d == zero) | |
253 | { | |
254 | // There is no IP of a//b, they are collinear or parallel | |
255 | // Assuming they intersect (this method should be called for | |
256 | // segments known to intersect), they are collinear and overlap. | |
257 | // They have one or two intersection points - we don't know and | |
258 | // have to rely on the fallback intersection point | |
259 | ||
260 | Point c1, c2; | |
261 | geometry::detail::assign_point_from_index<0>(c, c1); | |
262 | geometry::detail::assign_point_from_index<1>(c, c2); | |
263 | return side_by_triangle<>::apply(c1, c2, fallback_point); | |
264 | } | |
265 | ||
266 | // Cramer's rule: da (see cart_intersect.hpp) | |
267 | T const da = geometry::detail::determinant<T> | |
268 | ( | |
269 | dx_b, dy_b, | |
270 | ax - bx, ay - by | |
271 | ); | |
272 | ||
273 | // IP is at (ax + (da/d) * dx_a, ay + (da/d) * dy_a) | |
274 | // Side of IP is w.r.t. c is: determinant(dx_c, dy_c, ipx-cx, ipy-cy) | |
275 | // We replace ipx by expression above and multiply each term by d | |
276 | ||
277 | #ifdef BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG | |
278 | T const result1 = geometry::detail::determinant<T> | |
279 | ( | |
280 | dx_c * d, dy_c * d, | |
281 | d * (ax - cx) + dx_a * da, d * (ay - cy) + dy_a * da | |
282 | ); | |
283 | ||
284 | // Note: result / (d * d) | |
285 | // is identical to the side_value of side_by_triangle | |
286 | // Therefore, the sign is always the same as that result, and the | |
287 | // resulting side (left,right,collinear) is the same | |
288 | ||
289 | // The first row we divide again by d because of determinant multiply rule | |
290 | T const result2 = d * geometry::detail::determinant<T> | |
291 | ( | |
292 | dx_c, dy_c, | |
293 | d * (ax - cx) + dx_a * da, d * (ay - cy) + dy_a * da | |
294 | ); | |
295 | // Write out: | |
296 | T const result3 = d * (dx_c * (d * (ay - cy) + dy_a * da) | |
297 | - dy_c * (d * (ax - cx) + dx_a * da)); | |
298 | // Write out in braces: | |
299 | T const result4 = d * (dx_c * d * (ay - cy) + dx_c * dy_a * da | |
300 | - dy_c * d * (ax - cx) - dy_c * dx_a * da); | |
301 | // Write in terms of d * XX + da * YY | |
302 | T const result5 = d * (d * (dx_c * (ay - cy) - dy_c * (ax - cx)) | |
303 | + da * (dx_c * dy_a - dy_c * dx_a)); | |
304 | ||
305 | boost::ignore_unused(result1, result2, result3, result4, result5); | |
306 | //return result; | |
307 | #endif | |
308 | ||
309 | // We consider the results separately | |
310 | // (in the end we only have to return the side-value 1,0 or -1) | |
311 | ||
312 | // To avoid multiplications we judge the product (easy, avoids *d) | |
313 | // and the sign of p*q+r*s (more elaborate) | |
314 | T const result = sign_of_product | |
315 | ( | |
316 | d, | |
317 | sign_of_addition_of_two_products | |
318 | ( | |
319 | d, dx_c * (ay - cy) - dy_c * (ax - cx), | |
320 | da, dx_c * dy_a - dy_c * dx_a | |
321 | ) | |
322 | ); | |
323 | return result; | |
324 | ||
325 | ||
326 | } | |
327 | ||
328 | template <typename Segment, typename Point> | |
329 | static inline int apply(Segment const& a, Segment const& b, | |
330 | Segment const& c, | |
331 | Point const& fallback_point) | |
332 | { | |
333 | typedef typename geometry::coordinate_type<Segment>::type coordinate_type; | |
334 | coordinate_type const s = side_value<coordinate_type>(a, b, c, fallback_point); | |
335 | coordinate_type const zero = coordinate_type(); | |
336 | return math::equals(s, zero) ? 0 | |
337 | : s > zero ? 1 | |
338 | : -1; | |
339 | } | |
340 | ||
341 | }; | |
342 | ||
343 | ||
344 | }} // namespace strategy::side | |
345 | ||
346 | }} // namespace boost::geometry | |
347 | ||
348 | ||
349 | #endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP |