123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351 |
- // Boost.Geometry (aka GGL, Generic Geometry Library)
- // Copyright (c) 2015 Barend Gehrels, Amsterdam, the Netherlands.
- // This file was modified by Oracle on 2015, 2019.
- // Modifications copyright (c) 2015-2019, Oracle and/or its affiliates.
- // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
- // Use, modification and distribution is subject to the Boost Software License,
- // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
- // http://www.boost.org/LICENSE_1_0.txt)
- #ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP
- #define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP
- #include <limits>
- #include <boost/core/ignore_unused.hpp>
- #include <boost/type_traits/is_integral.hpp>
- #include <boost/type_traits/make_unsigned.hpp>
- #include <boost/geometry/arithmetic/determinant.hpp>
- #include <boost/geometry/core/access.hpp>
- #include <boost/geometry/core/assert.hpp>
- #include <boost/geometry/core/coordinate_type.hpp>
- #include <boost/geometry/algorithms/detail/assign_indexed_point.hpp>
- #include <boost/geometry/strategies/cartesian/side_by_triangle.hpp>
- #include <boost/geometry/util/math.hpp>
- #ifdef BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG
- #include <boost/math/common_factor_ct.hpp>
- #include <boost/math/common_factor_rt.hpp>
- #include <boost/multiprecision/cpp_int.hpp>
- #endif
- namespace boost { namespace geometry
- {
- namespace strategy { namespace side
- {
- namespace detail
- {
- // A tool for multiplication of integers avoiding overflow
- // It's a temporary workaround until we can use Multiprecision
- // The algorithm is based on Karatsuba algorithm
- // see: http://en.wikipedia.org/wiki/Karatsuba_algorithm
- template <typename T>
- struct multiplicable_integral
- {
- // Currently this tool can't be used with non-integral coordinate types.
- // Also side_of_intersection strategy sign_of_product() and sign_of_compare()
- // functions would have to be modified to properly support floating-point
- // types (comparisons and multiplication).
- BOOST_STATIC_ASSERT(boost::is_integral<T>::value);
- static const std::size_t bits = CHAR_BIT * sizeof(T);
- static const std::size_t half_bits = bits / 2;
- typedef typename boost::make_unsigned<T>::type unsigned_type;
- static const unsigned_type base = unsigned_type(1) << half_bits; // 2^half_bits
- int m_sign;
- unsigned_type m_ms;
- unsigned_type m_ls;
- multiplicable_integral(int sign, unsigned_type ms, unsigned_type ls)
- : m_sign(sign), m_ms(ms), m_ls(ls)
- {}
- explicit multiplicable_integral(T const& val)
- {
- unsigned_type val_u = val > 0 ?
- unsigned_type(val)
- : val == (std::numeric_limits<T>::min)() ?
- unsigned_type((std::numeric_limits<T>::max)()) + 1
- : unsigned_type(-val);
- // MMLL -> S 00MM 00LL
- m_sign = math::sign(val);
- m_ms = val_u >> half_bits; // val_u / base
- m_ls = val_u - m_ms * base;
- }
-
- friend multiplicable_integral operator*(multiplicable_integral const& a,
- multiplicable_integral const& b)
- {
- // (S 00MM 00LL) * (S 00MM 00LL) -> (S Z2MM 00LL)
- unsigned_type z2 = a.m_ms * b.m_ms;
- unsigned_type z0 = a.m_ls * b.m_ls;
- unsigned_type z1 = (a.m_ms + a.m_ls) * (b.m_ms + b.m_ls) - z2 - z0;
- // z0 may be >= base so it must be normalized to allow comparison
- unsigned_type z0_ms = z0 >> half_bits; // z0 / base
- return multiplicable_integral(a.m_sign * b.m_sign,
- z2 * base + z1 + z0_ms,
- z0 - base * z0_ms);
- }
- friend bool operator<(multiplicable_integral const& a,
- multiplicable_integral const& b)
- {
- if ( a.m_sign == b.m_sign )
- {
- bool u_less = a.m_ms < b.m_ms
- || (a.m_ms == b.m_ms && a.m_ls < b.m_ls);
- return a.m_sign > 0 ? u_less : (! u_less);
- }
- else
- {
- return a.m_sign < b.m_sign;
- }
- }
- friend bool operator>(multiplicable_integral const& a,
- multiplicable_integral const& b)
- {
- return b < a;
- }
- #ifdef BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG
- template <typename CmpVal>
- void check_value(CmpVal const& cmp_val) const
- {
- unsigned_type b = base; // a workaround for MinGW - undefined reference base
- CmpVal val = CmpVal(m_sign) * (CmpVal(m_ms) * CmpVal(b) + CmpVal(m_ls));
- BOOST_GEOMETRY_ASSERT(cmp_val == val);
- }
- #endif // BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG
- };
- } // namespace detail
- // Calculates the side of the intersection-point (if any) of
- // of segment a//b w.r.t. segment c
- // This is calculated without (re)calculating the IP itself again and fully
- // based on integer mathematics; there are no divisions
- // It can be used for either integer (rescaled) points, and also for FP
- class side_of_intersection
- {
- private :
- template <typename T, typename U>
- static inline
- int sign_of_product(T const& a, U const& b)
- {
- return a == 0 || b == 0 ? 0
- : a > 0 && b > 0 ? 1
- : a < 0 && b < 0 ? 1
- : -1;
- }
- template <typename T>
- static inline
- int sign_of_compare(T const& a, T const& b, T const& c, T const& d)
- {
- // Both a*b and c*d are positive
- // We have to judge if a*b > c*d
- using side::detail::multiplicable_integral;
- multiplicable_integral<T> ab = multiplicable_integral<T>(a)
- * multiplicable_integral<T>(b);
- multiplicable_integral<T> cd = multiplicable_integral<T>(c)
- * multiplicable_integral<T>(d);
-
- int result = ab > cd ? 1
- : ab < cd ? -1
- : 0
- ;
- #ifdef BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG
- using namespace boost::multiprecision;
- cpp_int const lab = cpp_int(a) * cpp_int(b);
- cpp_int const lcd = cpp_int(c) * cpp_int(d);
- ab.check_value(lab);
- cd.check_value(lcd);
- int result2 = lab > lcd ? 1
- : lab < lcd ? -1
- : 0
- ;
- BOOST_GEOMETRY_ASSERT(result == result2);
- #endif
- return result;
- }
- template <typename T>
- static inline
- int sign_of_addition_of_two_products(T const& a, T const& b, T const& c, T const& d)
- {
- // sign of a*b+c*d, 1 if positive, -1 if negative, else 0
- int const ab = sign_of_product(a, b);
- int const cd = sign_of_product(c, d);
- if (ab == 0)
- {
- return cd;
- }
- if (cd == 0)
- {
- return ab;
- }
- if (ab == cd)
- {
- // Both positive or both negative
- return ab;
- }
- // One is positive, one is negative, both are non zero
- // If ab is positive, we have to judge if a*b > -c*d (then 1 because sum is positive)
- // If ab is negative, we have to judge if c*d > -a*b (idem)
- return ab == 1
- ? sign_of_compare(a, b, -c, d)
- : sign_of_compare(c, d, -a, b);
- }
- public :
- // Calculates the side of the intersection-point (if any) of
- // of segment a//b w.r.t. segment c
- // This is calculated without (re)calculating the IP itself again and fully
- // based on integer mathematics
- template <typename T, typename Segment, typename Point>
- static inline T side_value(Segment const& a, Segment const& b,
- Segment const& c, Point const& fallback_point)
- {
- // The first point of the three segments is reused several times
- T const ax = get<0, 0>(a);
- T const ay = get<0, 1>(a);
- T const bx = get<0, 0>(b);
- T const by = get<0, 1>(b);
- T const cx = get<0, 0>(c);
- T const cy = get<0, 1>(c);
- T const dx_a = get<1, 0>(a) - ax;
- T const dy_a = get<1, 1>(a) - ay;
- T const dx_b = get<1, 0>(b) - bx;
- T const dy_b = get<1, 1>(b) - by;
- T const dx_c = get<1, 0>(c) - cx;
- T const dy_c = get<1, 1>(c) - cy;
- // Cramer's rule: d (see cart_intersect.hpp)
- T const d = geometry::detail::determinant<T>
- (
- dx_a, dy_a,
- dx_b, dy_b
- );
- T const zero = T();
- if (d == zero)
- {
- // There is no IP of a//b, they are collinear or parallel
- // Assuming they intersect (this method should be called for
- // segments known to intersect), they are collinear and overlap.
- // They have one or two intersection points - we don't know and
- // have to rely on the fallback intersection point
- Point c1, c2;
- geometry::detail::assign_point_from_index<0>(c, c1);
- geometry::detail::assign_point_from_index<1>(c, c2);
- return side_by_triangle<>::apply(c1, c2, fallback_point);
- }
- // Cramer's rule: da (see cart_intersect.hpp)
- T const da = geometry::detail::determinant<T>
- (
- dx_b, dy_b,
- ax - bx, ay - by
- );
- // IP is at (ax + (da/d) * dx_a, ay + (da/d) * dy_a)
- // Side of IP is w.r.t. c is: determinant(dx_c, dy_c, ipx-cx, ipy-cy)
- // We replace ipx by expression above and multiply each term by d
- #ifdef BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG
- T const result1 = geometry::detail::determinant<T>
- (
- dx_c * d, dy_c * d,
- d * (ax - cx) + dx_a * da, d * (ay - cy) + dy_a * da
- );
- // Note: result / (d * d)
- // is identical to the side_value of side_by_triangle
- // Therefore, the sign is always the same as that result, and the
- // resulting side (left,right,collinear) is the same
- // The first row we divide again by d because of determinant multiply rule
- T const result2 = d * geometry::detail::determinant<T>
- (
- dx_c, dy_c,
- d * (ax - cx) + dx_a * da, d * (ay - cy) + dy_a * da
- );
- // Write out:
- T const result3 = d * (dx_c * (d * (ay - cy) + dy_a * da)
- - dy_c * (d * (ax - cx) + dx_a * da));
- // Write out in braces:
- T const result4 = d * (dx_c * d * (ay - cy) + dx_c * dy_a * da
- - dy_c * d * (ax - cx) - dy_c * dx_a * da);
- // Write in terms of d * XX + da * YY
- T const result5 = d * (d * (dx_c * (ay - cy) - dy_c * (ax - cx))
- + da * (dx_c * dy_a - dy_c * dx_a));
- boost::ignore_unused(result1, result2, result3, result4, result5);
- //return result;
- #endif
- // We consider the results separately
- // (in the end we only have to return the side-value 1,0 or -1)
- // To avoid multiplications we judge the product (easy, avoids *d)
- // and the sign of p*q+r*s (more elaborate)
- T const result = sign_of_product
- (
- d,
- sign_of_addition_of_two_products
- (
- d, dx_c * (ay - cy) - dy_c * (ax - cx),
- da, dx_c * dy_a - dy_c * dx_a
- )
- );
- return result;
- }
- template <typename Segment, typename Point>
- static inline int apply(Segment const& a, Segment const& b,
- Segment const& c,
- Point const& fallback_point)
- {
- typedef typename geometry::coordinate_type<Segment>::type coordinate_type;
- coordinate_type const s = side_value<coordinate_type>(a, b, c, fallback_point);
- coordinate_type const zero = coordinate_type();
- return math::equals(s, zero) ? 0
- : s > zero ? 1
- : -1;
- }
- };
- }} // namespace strategy::side
- }} // namespace boost::geometry
- #endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP
|