diff options
Diffstat (limited to 'boost/geometry/strategies/cartesian/side_of_intersection.hpp')
| -rw-r--r-- | boost/geometry/strategies/cartesian/side_of_intersection.hpp | 248 |
1 files changed, 239 insertions, 9 deletions
diff --git a/boost/geometry/strategies/cartesian/side_of_intersection.hpp b/boost/geometry/strategies/cartesian/side_of_intersection.hpp index 39487676c1..db57644ad5 100644 --- a/boost/geometry/strategies/cartesian/side_of_intersection.hpp +++ b/boost/geometry/strategies/cartesian/side_of_intersection.hpp @@ -2,6 +2,11 @@ // Copyright (c) 2015 Barend Gehrels, Amsterdam, the Netherlands. +// This file was modified by Oracle on 2015. +// Modifications copyright (c) 2015, 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) @@ -10,11 +15,25 @@ #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 { @@ -22,6 +41,94 @@ 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; + } + + 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); + } +}; + +} // 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 @@ -29,15 +136,93 @@ namespace strategy { namespace side // 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> + template <typename T, typename Segment, typename Point> static inline T side_value(Segment const& a, Segment const& b, - Segment const& c) + 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); @@ -67,9 +252,15 @@ public : if (d == zero) { // There is no IP of a//b, they are collinear or parallel - // We don't have to divide but we can already conclude the side-value - // is meaningless and the resulting determinant will be 0 - return zero; + // 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) @@ -82,7 +273,9 @@ public : // 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 - T const result = geometry::detail::determinant<T> + +#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 @@ -93,15 +286,52 @@ public : // 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> - static inline int apply(Segment const& a, Segment const& b, Segment const& c) + 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); + 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 |