boost/geometry
2
0
Fork 0
mirror of https://github.com/boostorg/geometry.git synced 2026-02-13 12:32:09 +00:00
Code Issues Projects Releases Wiki Activity

Replace side_by_triangle with side_robust

Browse Source
This commit is contained in:
Vissarion Fisikopoulos
2021-06-29 13:20:07 +03:00
parent 174c380da2
commit fb9b1d40ff
14 changed files with 267 additions and 67 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
include/boost/geometry/algorithms/detail/equals/collect_vectors.hpp
View File

@@ -40,7 +40,7 @@
#include <boost/geometry/views/detail/normalized_view.hpp>
#include <boost/geometry/strategies/cartesian/side_by_triangle.hpp>
#include <boost/geometry/strategy/cartesian/side_robust.hpp>
#include <boost/geometry/strategies/spherical/ssf.hpp>
#include <boost/geometry/strategies/normalize.hpp>
@@ -62,11 +62,11 @@ struct collected_vector
: nyi::not_implemented_tag
{};
// compatible with side_by_triangle cartesian strategy
// compatible with side_robust cartesian strategy
template <typename T, typename Geometry, typename CT, typename CSTag>
struct collected_vector
<
T, Geometry, strategy::side::side_by_triangle<CT>, CSTag
T, Geometry, strategy::side::side_robust<CT>, CSTag
>
{
typedef T type;

6
include/boost/geometry/strategies/agnostic/point_in_box_by_side.hpp
View File

@@ -25,7 +25,7 @@
#include <boost/geometry/algorithms/assign.hpp>
#include <boost/geometry/strategies/covered_by.hpp>
#include <boost/geometry/strategies/within.hpp>
#include <boost/geometry/strategies/cartesian/side_by_triangle.hpp>
#include <boost/geometry/strategy/cartesian/side_robust.hpp>
#include <boost/geometry/strategies/geographic/side.hpp>
#include <boost/geometry/strategies/spherical/ssf.hpp>
@@ -118,7 +118,7 @@ struct cartesian_point_box_by_side
<
within::detail::decide_within
>(point, box,
strategy::side::side_by_triangle<CalculationType>());
strategy::side::side_robust<CalculationType>());
}
};
@@ -189,7 +189,7 @@ struct cartesian_point_box_by_side
<
within::detail::decide_covered_by
>(point, box,
strategy::side::side_by_triangle<CalculationType>());
strategy::side::side_robust<CalculationType>());
}
};

1
include/boost/geometry/strategies/cartesian/distance_segment_box.hpp
View File

@@ -17,7 +17,6 @@
#include <boost/geometry/strategies/cartesian/distance_pythagoras.hpp>
#include <boost/geometry/strategies/cartesian/distance_pythagoras_point_box.hpp>
#include <boost/geometry/strategies/cartesian/point_in_point.hpp>
#include <boost/geometry/strategies/cartesian/side_by_triangle.hpp>
namespace boost { namespace geometry
{

4
include/boost/geometry/strategies/cartesian/intersection.hpp
View File

@@ -44,7 +44,7 @@
#include <boost/geometry/strategies/cartesian/distance_pythagoras.hpp>
#include <boost/geometry/strategies/cartesian/point_in_point.hpp>
#include <boost/geometry/strategies/cartesian/point_in_poly_winding.hpp>
#include <boost/geometry/strategies/cartesian/side_by_triangle.hpp>
#include <boost/geometry/strategy/cartesian/side_robust.hpp>
#include <boost/geometry/strategies/covered_by.hpp>
#include <boost/geometry/strategies/intersection.hpp>
#include <boost/geometry/strategies/intersection_result.hpp>
@@ -402,7 +402,7 @@ struct cartesian_segments
return Policy::disjoint();
}
typedef side::side_by_triangle<CalculationType> side_strategy_type;
typedef side::side_robust<CalculationType> side_strategy_type;
side_info sides;
sides.set<0>(side_strategy_type::apply(q1, q2, p1),
side_strategy_type::apply(q1, q2, p2));

4
include/boost/geometry/strategies/cartesian/point_in_poly_winding.hpp
View File

@@ -27,7 +27,7 @@
#include <boost/geometry/strategies/cartesian/point_in_box.hpp>
#include <boost/geometry/strategies/cartesian/disjoint_box_box.hpp>
#include <boost/geometry/strategies/cartesian/side_by_triangle.hpp>
#include <boost/geometry/strategy/cartesian/side_robust.hpp>
#include <boost/geometry/strategies/covered_by.hpp>
#include <boost/geometry/strategies/within.hpp>
@@ -116,7 +116,7 @@ public:
else // count == 2 || count == -2
{
// 1 left, -1 right
typedef side::side_by_triangle<CalculationType> side_strategy_type;
typedef side::side_robust<CalculationType> side_strategy_type;
side = side_strategy_type::apply(s1, s2, point);
}

12
include/boost/geometry/strategies/relate/cartesian.hpp
View File

@@ -20,6 +20,7 @@
#include <boost/geometry/strategies/cartesian/point_in_poly_franklin.hpp>
#include <boost/geometry/strategies/cartesian/point_in_poly_winding.hpp>
#include <boost/geometry/strategies/cartesian/disjoint_box_box.hpp>
#include <boost/geometry/strategies/cartesian/side_by_triangle.hpp>
#include <boost/geometry/strategies/envelope/cartesian.hpp>
#include <boost/geometry/strategies/relate/services.hpp>
@@ -147,7 +148,7 @@ public:
static auto side()
{
return strategy::side::side_by_triangle<CalculationType>();
return strategy::side::side_robust<CalculationType>();
}
// within
@@ -372,6 +373,15 @@ struct strategy_converter<strategy::side::side_by_triangle<CalculationType>>
}
};
template <typename CalculationType>
struct strategy_converter<strategy::side::side_robust<CalculationType>>
{
static auto get(strategy::side::side_robust<CalculationType> const&)
{
return strategies::relate::cartesian<CalculationType>();
}
};
} // namespace services

1
include/boost/geometry/strategy/cartesian/side_non_robust.hpp
View File

@@ -62,6 +62,7 @@ public:
* (promoted_type(get<1>(p2)) - promoted_type(get<1>(p)));
auto detright = (promoted_type(get<1>(p1)) - promoted_type(get<1>(p)))
* (promoted_type(get<0>(p2)) - promoted_type(get<0>(p)));
return detleft > detright ? 1 : (detleft < detright ? -1 : 0 );
}

119
include/boost/geometry/strategy/cartesian/side_robust.hpp
View File

@@ -5,6 +5,11 @@
// Contributed and/or modified by Tinko Bartels,
// as part of Google Summer of Code 2019 program.
// This file was modified by Oracle on 2021.
// Modifications copyright (c) 2021, Oracle and/or its affiliates.
// Contributed and/or modified by Vissarion Fysikopoulos, 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)
@@ -12,9 +17,16 @@
#ifndef BOOST_GEOMETRY_STRATEGY_CARTESIAN_SIDE_ROBUST_HPP
#define BOOST_GEOMETRY_STRATEGY_CARTESIAN_SIDE_ROBUST_HPP
#include <iostream>
#include <boost/rational.hpp>
#include <boost/geometry/strategy/cartesian/side_non_robust.hpp>
#include <boost/geometry/util/select_most_precise.hpp>
#include <boost/geometry/util/select_calculation_type.hpp>
#include <boost/geometry/util/precise_math.hpp>
#include <boost/geometry/util/math.hpp>
namespace boost { namespace geometry
{
@@ -22,6 +34,28 @@ namespace boost { namespace geometry
namespace strategy { namespace side
{
struct eps_equals_policy
{
public:
template <typename Policy, typename T1, typename T2>
static auto apply(T1 a, T2 b, Policy policy)
{
return boost::geometry::math::detail::equals_by_policy(a, b, policy);
}
};
struct fp_equals_policy
{
public:
template <typename Policy, typename T1, typename T2>
static auto apply(T1 a, T2 b, Policy)
{
return a == b;
}
};
/*!
\brief Adaptive precision predicate to check at which side of a segment a point lies:
left of segment (>0), right of segment (< 0), on segment (0).
@@ -33,57 +67,106 @@ namespace strategy { namespace side
template
<
typename CalculationType = void,
typename EpsPolicy = eps_equals_policy,
std::size_t Robustness = 3
>
struct side_robust
{
template <typename Policy>
struct eps_policy
{
eps_policy() {}
template <typename Type>
eps_policy(Type const& a, Type const& b, Type const& c, Type const& d)
: policy(a, b, c, d)
{}
Policy policy;
};
public:
//! \brief Computes double the signed area of the CCW triangle p1, p2, p
typedef cartesian_tag cs_tag;
//! \brief Computes the sign of the CCW triangle p1, p2, p
template
<
typename PromotedType,
typename P1,
typename P2,
typename P
typename P,
typename EpsPolicyInternal,
std::enable_if_t<std::is_fundamental<PromotedType>::value, int> = 0
>
static inline PromotedType side_value(P1 const& p1, P2 const& p2,
P const& p)
static inline PromotedType side_value(P1 const& p1,
P2 const& p2,
P const& p,
EpsPolicyInternal& eps_policy)
{
typedef ::boost::geometry::detail::precise_math::vec2d<PromotedType> vec2d;
vec2d pa { get<0>(p1), get<1>(p1) };
vec2d pb { get<0>(p2), get<1>(p2) };
vec2d pc { get<0>(p), get<1>(p) };
using vec2d = ::boost::geometry::detail::precise_math::vec2d<PromotedType>;
vec2d pa;
pa.x = get<0>(p1);
pa.y = get<1>(p1);
vec2d pb;
pb.x = get<0>(p2);
pb.y = get<1>(p2);
vec2d pc;
pc.x = get<0>(p);
pc.y = get<1>(p);
return ::boost::geometry::detail::precise_math::orient2d
<PromotedType, Robustness>(pa, pb, pc);
<PromotedType, Robustness>(pa, pb, pc, eps_policy);
}
template
<
typename PromotedType,
typename P1,
typename P2,
typename P,
typename EpsPolicyInternal,
std::enable_if_t<!std::is_fundamental<PromotedType>::value, int> = 0
>
static inline auto side_value(P1 const& p1, P2 const& p2, P const& p,
EpsPolicyInternal&)
{
return side_non_robust<>::apply(p1, p2, p);
}
#ifndef DOXYGEN_SHOULD_SKIP_THIS
template
<
<
typename P1,
typename P2,
typename P
>
static inline int apply(P1 const& p1, P2 const& p2, P const& p)
{
typedef typename select_calculation_type_alt
using coordinate_type = typename select_calculation_type_alt
<
CalculationType,
P1,
P2,
P
>::type coordinate_type;
typedef typename select_most_precise
>::type;
using promoted_type = typename select_most_precise
<
coordinate_type,
double
>::type promoted_type;
>::type;
promoted_type sv = side_value<promoted_type>(p1, p2, p);
return sv > 0 ? 1
: sv < 0 ? -1
: 0;
eps_policy
<
boost::geometry::math::detail::equals_factor_policy<promoted_type>
> epsp;
promoted_type sv = side_value<promoted_type>(p1, p2, p, epsp);
promoted_type const zero = promoted_type();
return EpsPolicy::apply(sv, zero, epsp.policy) ? 0
: sv > zero ? 1
: -1;
}
#endif
};

55
include/boost/geometry/util/precise_math.hpp
View File

@@ -5,6 +5,10 @@
// Contributed and/or modified by Tinko Bartels,
// as part of Google Summer of Code 2019 program.
// This file was modified by Oracle on 2021.
// Modifications copyright (c) 2021, Oracle and/or its affiliates.
// Contributed and/or modified by Vissarion Fisikopoulos, 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)
@@ -157,7 +161,9 @@ inline int fast_expansion_sum_zeroelim(
int n = InSize2)
{
std::array<RealNumber, 2> Qh;
int i_e = 0, i_f = 0, i_h = 0;
int i_e = 0;
int i_f = 0;
int i_h = 0;
if (std::abs(f[0]) > std::abs(e[0]))
{
Qh[0] = e[i_e++];
@@ -282,14 +288,16 @@ inline RealNumber orient2dtail(vec2d<RealNumber> const& p1,
std::array<RealNumber, 2>& t4,
std::array<RealNumber, 2>& t5_01,
std::array<RealNumber, 2>& t6_01,
RealNumber const& magnitude
)
RealNumber const& magnitude)
{
t5_01[1] = two_product_tail(t1[0], t2[0], t5_01[0]);
t6_01[1] = two_product_tail(t3[0], t4[0], t6_01[0]);
std::array<RealNumber, 4> tA_03 = two_two_expansion_diff(t5_01, t6_01);
RealNumber det = std::accumulate(tA_03.begin(), tA_03.end(), static_cast<RealNumber>(0));
if(Robustness == 1) return det;
if (Robustness == 1)
{
return det;
}
// see p.39, mind the different definition of epsilon for error bound
RealNumber B_relative_bound =
(1 + 3 * std::numeric_limits<RealNumber>::epsilon())
@@ -347,29 +355,51 @@ inline RealNumber orient2dtail(vec2d<RealNumber> const& p1,
return(D[D_nz - 1]);
}
// see page 38, Figure 21 for the calculations, notation follows the notation in the figure.
// see page 38, Figure 21 for the calculations, notation follows the notation
// in the figure.
template
<
typename RealNumber,
std::size_t Robustness = 3
std::size_t Robustness = 3,
typename EpsPolicy
>
inline RealNumber orient2d(vec2d<RealNumber> const& p1,
vec2d<RealNumber> const& p2,
vec2d<RealNumber> const& p3)
vec2d<RealNumber> const& p3,
EpsPolicy& eps_policy)
{
if(Robustness == 0)
{
return (p1.x - p3.x) * (p2.y - p3.y) - (p1.y - p3.y) * (p2.x - p3.x);
}
auto const x = p3.x;
auto const y = p3.y;
auto const sx1 = p1.x;
auto const sy1 = p1.y;
auto const sx2 = p2.x;
auto const sy2 = p2.y;
auto const dx = sx2 - sx1;
auto const dy = sy2 - sy1;
auto const dpx = x - sx1;
auto const dpy = y - sy1;
eps_policy = EpsPolicy(dx, dy, dpx, dpy);
std::array<RealNumber, 2> t1, t2, t3, t4;
t1[0] = p1.x - p3.x;
t2[0] = p2.y - p3.y;
t3[0] = p1.y - p3.y;
t4[0] = p2.x - p3.x;
std::array<RealNumber, 2> t5_01, t6_01;
t5_01[0] = t1[0] * t2[0];
t6_01[0] = t3[0] * t4[0];
RealNumber det = t5_01[0] - t6_01[0];
if (Robustness == 0)
{
return det;
}
RealNumber const magnitude = std::abs(t5_01[0]) + std::abs(t6_01[0]);
// see p.39, mind the different definition of epsilon for error bound
@@ -388,7 +418,8 @@ inline RealNumber orient2d(vec2d<RealNumber> const& p1,
//obvious
return det;
}
return orient2dtail<RealNumber, Robustness>(p1, p2, p3, t1, t2, t3, t4, t5_01, t6_01, magnitude);
return orient2dtail<RealNumber, Robustness>(p1, p2, p3, t1, t2, t3, t4,
t5_01, t6_01, magnitude);
}
// This method adaptively computes increasingly precise approximations of the following

40
include/boost/geometry/util/select_most_precise.hpp
View File

@@ -30,7 +30,6 @@ namespace boost { namespace geometry
namespace detail { namespace select_most_precise
{
// 0 - void
// 1 - integral
// 2 - floating point
@@ -55,6 +54,43 @@ struct type_priority
>
{};
/*
// 0 - void
// 1 - integral (int, long int)
// 2 - floating point (float, double)
// 3 - long long int
// 4 - long double
// 5 - non-fundamental
template <typename T>
struct type_priority
: std::conditional_t
<
std::is_void<T>::value,
std::integral_constant<int, 0>,
std::conditional_t
<
std::is_fundamental<T>::value,
std::conditional_t
<
std::is_same<T, long double>::value,
std::integral_constant<int, 4>,
std::conditional_t
<
std::is_same<T, long long int>::value,
std::integral_constant<int, 3>,
std::conditional_t
<
std::is_floating_point<T>::value,
std::integral_constant<int, 2>,
std::integral_constant<int, 1>
>
>
>,
std::integral_constant<int, 5>
>
>
{};
*/
template <typename T>
struct type_size
@@ -119,7 +155,7 @@ struct select_most_precise<T1, T2>
T2,
std::conditional_t // priority1 == priority2
<
(priority1 == 0 || priority1 == 3), // both void or non-fundamental
(priority1 == 0 || priority1 == 5), // both void or non-fundamental
T1,
std::conditional_t // both fundamental
<