1 // Boost.Geometry (aka GGL, Generic Geometry Library)
3 // Copyright (c) 2015 Barend Gehrels, Amsterdam, the Netherlands.
5 // This file was modified by Oracle on 2015, 2019.
6 // Modifications copyright (c) 2015-2019, Oracle and/or its affiliates.
8 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
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)
14 #ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP
15 #define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP
20 #include <boost/core/ignore_unused.hpp>
21 #include <boost/type_traits/is_integral.hpp>
22 #include <boost/type_traits/make_unsigned.hpp>
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>
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>
38 namespace boost { namespace geometry
41 namespace strategy { namespace side
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
52 struct multiplicable_integral
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);
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
69 multiplicable_integral(int sign, unsigned_type ms, unsigned_type ls)
70 : m_sign(sign), m_ms(ms), m_ls(ls)
73 explicit multiplicable_integral(T const& val)
75 unsigned_type val_u = val > 0 ?
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;
86 friend multiplicable_integral operator*(multiplicable_integral const& a,
87 multiplicable_integral const& b)
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,
100 friend bool operator<(multiplicable_integral const& a,
101 multiplicable_integral const& b)
103 if ( a.m_sign == b.m_sign )
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);
111 return a.m_sign < b.m_sign;
115 friend bool operator>(multiplicable_integral const& a,
116 multiplicable_integral const& b)
121 #ifdef BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG
122 template <typename CmpVal>
123 void check_value(CmpVal const& cmp_val) const
125 unsigned_type b = base; // a workaround for MinGW - undefined reference base
126 CmpVal val = CmpVal(m_sign) * (CmpVal(m_ms) * CmpVal(b) + CmpVal(m_ls));
127 BOOST_GEOMETRY_ASSERT(cmp_val == val);
129 #endif // BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG
132 } // namespace detail
134 // Calculates the side of the intersection-point (if any) of
135 // of segment a//b w.r.t. segment c
136 // This is calculated without (re)calculating the IP itself again and fully
137 // based on integer mathematics; there are no divisions
138 // It can be used for either integer (rescaled) points, and also for FP
139 class side_of_intersection
142 template <typename T, typename U>
144 int sign_of_product(T const& a, U const& b)
146 return a == 0 || b == 0 ? 0
152 template <typename T>
154 int sign_of_compare(T const& a, T const& b, T const& c, T const& d)
156 // Both a*b and c*d are positive
157 // We have to judge if a*b > c*d
159 using side::detail::multiplicable_integral;
160 multiplicable_integral<T> ab = multiplicable_integral<T>(a)
161 * multiplicable_integral<T>(b);
162 multiplicable_integral<T> cd = multiplicable_integral<T>(c)
163 * multiplicable_integral<T>(d);
165 int result = ab > cd ? 1
170 #ifdef BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG
171 using namespace boost::multiprecision;
172 cpp_int const lab = cpp_int(a) * cpp_int(b);
173 cpp_int const lcd = cpp_int(c) * cpp_int(d);
178 int result2 = lab > lcd ? 1
182 BOOST_GEOMETRY_ASSERT(result == result2);
188 template <typename T>
190 int sign_of_addition_of_two_products(T const& a, T const& b, T const& c, T const& d)
192 // sign of a*b+c*d, 1 if positive, -1 if negative, else 0
193 int const ab = sign_of_product(a, b);
194 int const cd = sign_of_product(c, d);
206 // Both positive or both negative
210 // One is positive, one is negative, both are non zero
211 // If ab is positive, we have to judge if a*b > -c*d (then 1 because sum is positive)
212 // If ab is negative, we have to judge if c*d > -a*b (idem)
214 ? sign_of_compare(a, b, -c, d)
215 : sign_of_compare(c, d, -a, b);
221 // Calculates the side of the intersection-point (if any) of
222 // of segment a//b w.r.t. segment c
223 // This is calculated without (re)calculating the IP itself again and fully
224 // based on integer mathematics
225 template <typename T, typename Segment, typename Point>
226 static inline T side_value(Segment const& a, Segment const& b,
227 Segment const& c, Point const& fallback_point)
229 // The first point of the three segments is reused several times
230 T const ax = get<0, 0>(a);
231 T const ay = get<0, 1>(a);
232 T const bx = get<0, 0>(b);
233 T const by = get<0, 1>(b);
234 T const cx = get<0, 0>(c);
235 T const cy = get<0, 1>(c);
237 T const dx_a = get<1, 0>(a) - ax;
238 T const dy_a = get<1, 1>(a) - ay;
240 T const dx_b = get<1, 0>(b) - bx;
241 T const dy_b = get<1, 1>(b) - by;
243 T const dx_c = get<1, 0>(c) - cx;
244 T const dy_c = get<1, 1>(c) - cy;
246 // Cramer's rule: d (see cart_intersect.hpp)
247 T const d = geometry::detail::determinant<T>
256 // There is no IP of a//b, they are collinear or parallel
257 // Assuming they intersect (this method should be called for
258 // segments known to intersect), they are collinear and overlap.
259 // They have one or two intersection points - we don't know and
260 // have to rely on the fallback intersection point
263 geometry::detail::assign_point_from_index<0>(c, c1);
264 geometry::detail::assign_point_from_index<1>(c, c2);
265 return side_by_triangle<>::apply(c1, c2, fallback_point);
268 // Cramer's rule: da (see cart_intersect.hpp)
269 T const da = geometry::detail::determinant<T>
275 // IP is at (ax + (da/d) * dx_a, ay + (da/d) * dy_a)
276 // Side of IP is w.r.t. c is: determinant(dx_c, dy_c, ipx-cx, ipy-cy)
277 // We replace ipx by expression above and multiply each term by d
279 #ifdef BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG
280 T const result1 = geometry::detail::determinant<T>
283 d * (ax - cx) + dx_a * da, d * (ay - cy) + dy_a * da
286 // Note: result / (d * d)
287 // is identical to the side_value of side_by_triangle
288 // Therefore, the sign is always the same as that result, and the
289 // resulting side (left,right,collinear) is the same
291 // The first row we divide again by d because of determinant multiply rule
292 T const result2 = d * geometry::detail::determinant<T>
295 d * (ax - cx) + dx_a * da, d * (ay - cy) + dy_a * da
298 T const result3 = d * (dx_c * (d * (ay - cy) + dy_a * da)
299 - dy_c * (d * (ax - cx) + dx_a * da));
300 // Write out in braces:
301 T const result4 = d * (dx_c * d * (ay - cy) + dx_c * dy_a * da
302 - dy_c * d * (ax - cx) - dy_c * dx_a * da);
303 // Write in terms of d * XX + da * YY
304 T const result5 = d * (d * (dx_c * (ay - cy) - dy_c * (ax - cx))
305 + da * (dx_c * dy_a - dy_c * dx_a));
307 boost::ignore_unused(result1, result2, result3, result4, result5);
311 // We consider the results separately
312 // (in the end we only have to return the side-value 1,0 or -1)
314 // To avoid multiplications we judge the product (easy, avoids *d)
315 // and the sign of p*q+r*s (more elaborate)
316 T const result = sign_of_product
319 sign_of_addition_of_two_products
321 d, dx_c * (ay - cy) - dy_c * (ax - cx),
322 da, dx_c * dy_a - dy_c * dx_a
330 template <typename Segment, typename Point>
331 static inline int apply(Segment const& a, Segment const& b,
333 Point const& fallback_point)
335 typedef typename geometry::coordinate_type<Segment>::type coordinate_type;
336 coordinate_type const s = side_value<coordinate_type>(a, b, c, fallback_point);
337 coordinate_type const zero = coordinate_type();
338 return math::equals(s, zero) ? 0
346 }} // namespace strategy::side
348 }} // namespace boost::geometry
351 #endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP