Server : Apache/2.4.41 (Ubuntu) System : Linux absol.cf 5.4.0-198-generic #218-Ubuntu SMP Fri Sep 27 20:18:53 UTC 2024 x86_64 User : www-data ( 33) PHP Version : 7.4.33 Disable Function : pcntl_alarm,pcntl_fork,pcntl_waitpid,pcntl_wait,pcntl_wifexited,pcntl_wifstopped,pcntl_wifsignaled,pcntl_wifcontinued,pcntl_wexitstatus,pcntl_wtermsig,pcntl_wstopsig,pcntl_signal,pcntl_signal_get_handler,pcntl_signal_dispatch,pcntl_get_last_error,pcntl_strerror,pcntl_sigprocmask,pcntl_sigwaitinfo,pcntl_sigtimedwait,pcntl_exec,pcntl_getpriority,pcntl_setpriority,pcntl_async_signals,pcntl_unshare, Directory : /usr/include/boost/geometry/strategies/cartesian/
Upload File :
Current File : //usr/include/boost/geometry/strategies/cartesian/side_of_intersection.hpp
// 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
Server : Apache/2.4.41 (Ubuntu)