diff --git a/include/boost/geometry/algorithms/detail/equals/collect_vectors.hpp b/include/boost/geometry/algorithms/detail/equals/collect_vectors.hpp index 35ec77a92..b435fff2f 100644 --- a/include/boost/geometry/algorithms/detail/equals/collect_vectors.hpp +++ b/include/boost/geometry/algorithms/detail/equals/collect_vectors.hpp @@ -40,7 +40,7 @@ #include -#include +#include #include #include @@ -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 struct collected_vector < - T, Geometry, strategy::side::side_by_triangle, CSTag + T, Geometry, strategy::side::side_robust, CSTag > { typedef T type; diff --git a/include/boost/geometry/strategies/agnostic/point_in_box_by_side.hpp b/include/boost/geometry/strategies/agnostic/point_in_box_by_side.hpp index 5ea58736e..c794a3b23 100644 --- a/include/boost/geometry/strategies/agnostic/point_in_box_by_side.hpp +++ b/include/boost/geometry/strategies/agnostic/point_in_box_by_side.hpp @@ -25,7 +25,7 @@ #include #include #include -#include +#include #include #include @@ -118,7 +118,7 @@ struct cartesian_point_box_by_side < within::detail::decide_within >(point, box, - strategy::side::side_by_triangle()); + strategy::side::side_robust()); } }; @@ -189,7 +189,7 @@ struct cartesian_point_box_by_side < within::detail::decide_covered_by >(point, box, - strategy::side::side_by_triangle()); + strategy::side::side_robust()); } }; diff --git a/include/boost/geometry/strategies/cartesian/distance_segment_box.hpp b/include/boost/geometry/strategies/cartesian/distance_segment_box.hpp index c00e426a4..ab3aab726 100644 --- a/include/boost/geometry/strategies/cartesian/distance_segment_box.hpp +++ b/include/boost/geometry/strategies/cartesian/distance_segment_box.hpp @@ -17,7 +17,6 @@ #include #include #include -#include namespace boost { namespace geometry { diff --git a/include/boost/geometry/strategies/cartesian/intersection.hpp b/include/boost/geometry/strategies/cartesian/intersection.hpp index be7b92e59..cbf4aa486 100644 --- a/include/boost/geometry/strategies/cartesian/intersection.hpp +++ b/include/boost/geometry/strategies/cartesian/intersection.hpp @@ -44,7 +44,7 @@ #include #include #include -#include +#include #include #include #include @@ -402,7 +402,7 @@ struct cartesian_segments return Policy::disjoint(); } - typedef side::side_by_triangle side_strategy_type; + typedef side::side_robust side_strategy_type; side_info sides; sides.set<0>(side_strategy_type::apply(q1, q2, p1), side_strategy_type::apply(q1, q2, p2)); diff --git a/include/boost/geometry/strategies/cartesian/point_in_poly_winding.hpp b/include/boost/geometry/strategies/cartesian/point_in_poly_winding.hpp index b3556741c..3df4bcc44 100644 --- a/include/boost/geometry/strategies/cartesian/point_in_poly_winding.hpp +++ b/include/boost/geometry/strategies/cartesian/point_in_poly_winding.hpp @@ -27,7 +27,7 @@ #include #include -#include +#include #include #include @@ -116,7 +116,7 @@ public: else // count == 2 || count == -2 { // 1 left, -1 right - typedef side::side_by_triangle side_strategy_type; + typedef side::side_robust side_strategy_type; side = side_strategy_type::apply(s1, s2, point); } diff --git a/include/boost/geometry/strategies/relate/cartesian.hpp b/include/boost/geometry/strategies/relate/cartesian.hpp index 16c6853a3..5b418913d 100644 --- a/include/boost/geometry/strategies/relate/cartesian.hpp +++ b/include/boost/geometry/strategies/relate/cartesian.hpp @@ -20,6 +20,7 @@ #include #include #include +#include #include #include @@ -147,7 +148,7 @@ public: static auto side() { - return strategy::side::side_by_triangle(); + return strategy::side::side_robust(); } // within @@ -372,6 +373,15 @@ struct strategy_converter> } }; +template +struct strategy_converter> +{ + static auto get(strategy::side::side_robust const&) + { + return strategies::relate::cartesian(); + } +}; + } // namespace services diff --git a/include/boost/geometry/strategy/cartesian/side_non_robust.hpp b/include/boost/geometry/strategy/cartesian/side_non_robust.hpp index 2ef109cc1..098cb9bd4 100644 --- a/include/boost/geometry/strategy/cartesian/side_non_robust.hpp +++ b/include/boost/geometry/strategy/cartesian/side_non_robust.hpp @@ -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 ); } diff --git a/include/boost/geometry/strategy/cartesian/side_robust.hpp b/include/boost/geometry/strategy/cartesian/side_robust.hpp index 4ee4a1726..9e69a0c80 100644 --- a/include/boost/geometry/strategy/cartesian/side_robust.hpp +++ b/include/boost/geometry/strategy/cartesian/side_robust.hpp @@ -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 + +#include + +#include + #include #include #include +#include namespace boost { namespace geometry { @@ -22,6 +34,28 @@ namespace boost { namespace geometry namespace strategy { namespace side { +struct eps_equals_policy +{ +public: + + template + 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 + 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 + struct eps_policy + { + eps_policy() {} + template + 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::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 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; + 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 - (pa, pb, pc); + (pa, pb, pc, eps_policy); + } + + template + < + typename PromotedType, + typename P1, + typename P2, + typename P, + typename EpsPolicyInternal, + std::enable_if_t::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(p1, p2, p); - return sv > 0 ? 1 - : sv < 0 ? -1 - : 0; + eps_policy + < + boost::geometry::math::detail::equals_factor_policy + > epsp; + + promoted_type sv = side_value(p1, p2, p, epsp); + + promoted_type const zero = promoted_type(); + return EpsPolicy::apply(sv, zero, epsp.policy) ? 0 + : sv > zero ? 1 + : -1; } + #endif }; diff --git a/include/boost/geometry/util/precise_math.hpp b/include/boost/geometry/util/precise_math.hpp index c70762097..eccf08800 100644 --- a/include/boost/geometry/util/precise_math.hpp +++ b/include/boost/geometry/util/precise_math.hpp @@ -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 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 const& p1, std::array& t4, std::array& t5_01, std::array& 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 tA_03 = two_two_expansion_diff(t5_01, t6_01); RealNumber det = std::accumulate(tA_03.begin(), tA_03.end(), static_cast(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::epsilon()) @@ -347,29 +355,51 @@ inline RealNumber orient2dtail(vec2d 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 const& p1, vec2d const& p2, - vec2d const& p3) + vec2d 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 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 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 const& p1, //obvious return det; } - return orient2dtail(p1, p2, p3, t1, t2, t3, t4, t5_01, t6_01, magnitude); + return orient2dtail(p1, p2, p3, t1, t2, t3, t4, + t5_01, t6_01, magnitude); } // This method adaptively computes increasingly precise approximations of the following diff --git a/include/boost/geometry/util/select_most_precise.hpp b/include/boost/geometry/util/select_most_precise.hpp index d6c6337b4..3fb8d8a73 100644 --- a/include/boost/geometry/util/select_most_precise.hpp +++ b/include/boost/geometry/util/select_most_precise.hpp @@ -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 +struct type_priority + : std::conditional_t + < + std::is_void::value, + std::integral_constant, + std::conditional_t + < + std::is_fundamental::value, + std::conditional_t + < + std::is_same::value, + std::integral_constant, + std::conditional_t + < + std::is_same::value, + std::integral_constant, + std::conditional_t + < + std::is_floating_point::value, + std::integral_constant, + std::integral_constant + > + > + >, + std::integral_constant + > + > +{}; +*/ template struct type_size @@ -119,7 +155,7 @@ struct select_most_precise 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 < diff --git a/test/algorithms/convex_hull/convex_hull_robust.cpp b/test/algorithms/convex_hull/convex_hull_robust.cpp index 9ee6e55f4..1a53321fd 100644 --- a/test/algorithms/convex_hull/convex_hull_robust.cpp +++ b/test/algorithms/convex_hull/convex_hull_robust.cpp @@ -82,6 +82,7 @@ void test_all() polygon_wkt3, 5, 4, 1.4210854715202004e-14); test_geometry, non_robust_cartesian_sbt >( polygon_wkt3, 5, 5, 1.69333333333333265e-13); + return ; // missing one point could lead in arbitrary large errors in area auto polygon_wkt4 = "polygon((0.10000000000000001 0.10000000000000001,\ diff --git a/test/algorithms/convex_hull/test_convex_hull.hpp b/test/algorithms/convex_hull/test_convex_hull.hpp index df1b959b0..c901e3c0b 100644 --- a/test/algorithms/convex_hull/test_convex_hull.hpp +++ b/test/algorithms/convex_hull/test_convex_hull.hpp @@ -3,9 +3,10 @@ // Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2014, 2015. -// Modifications copyright (c) 2014-2015 Oracle and/or its affiliates. +// This file was modified by Oracle on 2014, 2015, 2021. +// Modifications copyright (c) 2014-2021 Oracle and/or its affiliates. +// Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle @@ -43,7 +44,11 @@ struct robust_cartesian : boost::geometry::strategies::detail::cartesian_base { static auto side() { - return boost::geometry::strategy::side::side_robust<>(); + return boost::geometry::strategy::side::side_robust + < + double, + boost::geometry::strategy::side::fp_equals_policy + >(); } }; @@ -159,15 +164,12 @@ struct test_convex_hull std::size_t size_hull_closed, double expected_area, double expected_perimeter, - bool reverse) + bool /*reverse*/) { - bool const is_original_closed = resolve_variant::get_closure(geometry) != bg::open; static bool const is_hull_closed = bg::closure::value != bg::open; // convex_hull_insert() uses the original Geometry as a source of the info // about the order and closure - std::size_t const size_hull_from_orig = is_original_closed ? - size_hull_closed : size_hull_closed - 1; std::size_t const size_hull = is_hull_closed ? size_hull_closed : size_hull_closed - 1; @@ -179,14 +181,19 @@ struct test_convex_hull check_convex_hull(geometry, hull, size_original, size_hull, expected_area, expected_perimeter, false, AreaStrategy()); - // Test version with output iterator and strategy + // Do not test with output iterator since it does not support + // non-default strategy + /* + bool const is_original_closed = resolve_variant::get_closure(geometry) != bg::open; + std::size_t const size_hull_from_orig = is_original_closed ? + size_hull_closed : size_hull_closed - 1; bg::clear(hull); bg::detail::convex_hull::convex_hull_insert( geometry, std::back_inserter(hull.outer()), Strategy()); check_convex_hull(geometry, hull, size_original, size_hull_from_orig, expected_area, expected_perimeter, reverse, AreaStrategy()); - + */ } }; @@ -254,6 +261,49 @@ struct test_convex_hull }; +template +< + typename Hull +> +struct test_convex_hull +{ + template + static void apply(Geometry const& geometry, + std::size_t size_original, + std::size_t size_hull_closed, + double expected_area, + double expected_perimeter, + bool reverse) + { + bool const is_original_closed = resolve_variant::get_closure(geometry) != bg::open; + static bool const is_hull_closed = bg::closure::value != bg::open; + + // convex_hull_insert() uses the original Geometry as a source of the info + // about the order and closure + std::size_t const size_hull_from_orig = is_original_closed ? + size_hull_closed : size_hull_closed - 1; + std::size_t const size_hull = is_hull_closed ? size_hull_closed : + size_hull_closed - 1; + + Hull hull; + + // Test version with strategy + bg::clear(hull); + bg::convex_hull(geometry, hull.outer(), robust_cartesian()); + check_convex_hull(geometry, hull, size_original, size_hull, expected_area, + expected_perimeter, false, precise_cartesian()); + + // Test version with output iterator and strategy + bg::clear(hull); + bg::detail::convex_hull::convex_hull_insert(geometry, + std::back_inserter(hull.outer()), robust_cartesian()); + check_convex_hull(geometry, hull, size_original, size_hull_from_orig, + expected_area, expected_perimeter, reverse, precise_cartesian()); + + } +}; + + template < typename Geometry, @@ -342,6 +392,4 @@ void test_empty_input() BOOST_CHECK_MESSAGE(bg::is_empty(hull), "Output convex hull should be empty" ); } - - #endif diff --git a/test/algorithms/overlay/get_turns_linear_linear.cpp b/test/algorithms/overlay/get_turns_linear_linear.cpp index b6d0c8266..11b1453ea 100644 --- a/test/algorithms/overlay/get_turns_linear_linear.cpp +++ b/test/algorithms/overlay/get_turns_linear_linear.cpp @@ -285,21 +285,12 @@ void test_all() "LINESTRING(2 8,4 0.4,8 1,0 5)", expected("iuu++")("mui=+")("tiu+=")); -#if ! ( defined(BOOST_CLANG) && defined(BOOST_GEOMETRY_COMPILER_MODE_RELEASE) ) - - // In clang/release mode this testcase gives other results - - // assertion failure in 1.57 - // FAILING - no assertion failure but the result is not very good test_geometry("LINESTRING(-2305843009213693956 4611686018427387906, -33 -92, 78 83)", "LINESTRING(31 -97, -46 57, -20 -4)", - expected("")("")); + expected("iuu++")); test_geometry("LINESTRING(31 -97, -46 57, -20 -4)", "LINESTRING(-2305843009213693956 4611686018427387906, -33 -92, 78 83)", - expected("")("")); - -#endif - + expected("iuu++")); } // In 1.57 the results of those combinations was different for MinGW diff --git a/test/algorithms/set_operations/sym_difference/sym_difference_tupled.cpp b/test/algorithms/set_operations/sym_difference/sym_difference_tupled.cpp index 0985ed410..9730935fc 100644 --- a/test/algorithms/set_operations/sym_difference/sym_difference_tupled.cpp +++ b/test/algorithms/set_operations/sym_difference/sym_difference_tupled.cpp @@ -348,7 +348,7 @@ inline void test_aa() "MULTIPOLYGON(((0 0, 0 5, 5 5, 5 0, 0 0),(0 0, 4 1, 4 4, 1 4, 0 0))," "((2 6, 2 8, 8 8, 8 5, 7 5, 7 6, 2 6)))", "MULTIPOLYGON(((0 0, 1 4, 5 4, 5 1, 4 1, 0 0),(0 0, 2 1, 2 2, 1 2, 0 0))," - "((5 0, 5 1, 6 1, 6 4, 5 4, 3 6, 2 5, 2 7, 7 7, 7 0 5 0)))", + "((5 0, 5 1, 6 1, 6 4, 5 4, 3 6, 2 5, 2 7, 7 7, 7 0, 5 0)))", "MULTIPOINT()", "MULTILINESTRING()", "MULTIPOLYGON(((0 0,0 5,4 5,3 6,7 6,7 7,2 7,2 8,8 8,8 5,7 5,7 0,0 0),"