diff options
Diffstat (limited to 'boost/geometry/strategies/cartesian/intersection.hpp')
| -rw-r--r-- | boost/geometry/strategies/cartesian/intersection.hpp | 119 |
1 files changed, 77 insertions, 42 deletions
diff --git a/boost/geometry/strategies/cartesian/intersection.hpp b/boost/geometry/strategies/cartesian/intersection.hpp index be7b92e59e..94d33401c5 100644 --- a/boost/geometry/strategies/cartesian/intersection.hpp +++ b/boost/geometry/strategies/cartesian/intersection.hpp @@ -1,7 +1,7 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) // Copyright (c) 2007-2014 Barend Gehrels, Amsterdam, the Netherlands. -// Copyright (c) 2013-2017 Adam Wulkiewicz, Lodz, Poland. +// Copyright (c) 2013-2023 Adam Wulkiewicz, Lodz, Poland. // This file was modified by Oracle on 2014-2021. // Modifications copyright (c) 2014-2021, Oracle and/or its affiliates. @@ -44,7 +44,6 @@ #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/strategies/covered_by.hpp> #include <boost/geometry/strategies/intersection.hpp> #include <boost/geometry/strategies/intersection_result.hpp> @@ -68,6 +67,61 @@ namespace boost { namespace geometry namespace strategy { namespace intersection { +namespace detail_usage +{ + +// When calculating the intersection, the information of "a" or "b" can be used. +// Theoretically this gives equal results, but due to floating point precision +// there might be tiny differences. These are edge cases. +// This structure is to determine if "a" or "b" should be used. +// Prefer the segment closer to the endpoint. +// If both are about equally close, then prefer the longer segment +// To avoid hard thresholds, behavior is made fluent. +// Calculate comparable length indications, +// the longer the segment (relatively), the lower the value +// such that the shorter lengths are evaluated higher and will +// be preferred. +template <bool IsArithmetic> +struct use_a +{ + template <typename Ct, typename Ev> + static bool apply(Ct const& cla, Ct const& clb, Ev const& eva, Ev const& evb) + { + auto const clm = (std::max)(cla, clb); + if (clm <= 0) + { + return true; + } + + // Relative comparible length + auto const rcla = Ct(1.0) - cla / clm; + auto const rclb = Ct(1.0) - clb / clm; + + // Multipliers for edgevalue (ev) and relative comparible length (rcl) + // They determine the balance between edge value (should be larger) + // and segment length. In 99.9xx% of the cases there is no difference + // at all (if either a or b is used). Therefore the values of the + // constants are not sensitive for the majority of the situations. + // One known case is #mysql_23023665_6 (difference) which needs mev >= 2 + Ev const mev = 5; + Ev const mrcl = 1; + + return mev * eva + mrcl * rcla > mev * evb + mrcl * rclb; + } +}; + +// Specialization for non arithmetic types. They will always use "a" +template <> +struct use_a<false> +{ + template <typename Ct, typename Ev> + static bool apply(Ct const& , Ct const& , Ev const& , Ev const& ) + { + return true; + } +}; + +} /*! \see http://mathworld.wolfram.com/Line-LineIntersection.html @@ -116,10 +170,9 @@ struct cartesian_segments SegmentRatio const& ratio) const { // Calculate the intersection point based on segment_ratio - // Up to now, division was postponed. Here we divide using numerator/ - // denominator. In case of integer this results in an integer - // division. - BOOST_GEOMETRY_ASSERT(ratio.denominator() != 0); + // The division, postponed until here, is done now. In case of integer this + // results in an integer which rounds to the nearest integer. + BOOST_GEOMETRY_ASSERT(ratio.denominator() != typename SegmentRatio::int_type(0)); typedef typename promote_integral<CoordinateType>::type calc_type; @@ -131,11 +184,11 @@ struct cartesian_segments calc_type const dy_calc = boost::numeric_cast<calc_type>(dy); set<0>(point, get<0, 0>(segment) - + boost::numeric_cast<CoordinateType>(numerator * dx_calc - / denominator)); + + boost::numeric_cast<CoordinateType>( + math::divide<calc_type>(numerator * dx_calc, denominator))); set<1>(point, get<0, 1>(segment) - + boost::numeric_cast<CoordinateType>(numerator * dy_calc - / denominator)); + + boost::numeric_cast<CoordinateType>( + math::divide<calc_type>(numerator * dy_calc, denominator))); } template <int Index, int Dim, typename Point, typename Segment> @@ -180,30 +233,12 @@ struct cartesian_segments template <typename Point, typename Segment1, typename Segment2> void calculate(Point& point, Segment1 const& a, Segment2 const& b) const { - bool use_a = true; - - // Prefer one segment if one is on or near an endpoint - bool const a_near_end = robust_ra.near_end(); - bool const b_near_end = robust_rb.near_end(); - if (a_near_end && ! b_near_end) - { - use_a = true; - } - else if (b_near_end && ! a_near_end) - { - use_a = false; - } - else - { - // Prefer shorter segment - promoted_type const len_a = comparable_length_a(); - promoted_type const len_b = comparable_length_b(); - if (len_b < len_a) - { - use_a = false; - } - // else use_a is true but was already assigned like that - } + bool const use_a + = detail_usage::use_a + < + std::is_arithmetic<CoordinateType>::value + >::apply(comparable_length_a(), comparable_length_b(), + robust_ra.edge_value(), robust_rb.edge_value()); if (use_a) { @@ -214,10 +249,7 @@ struct cartesian_segments assign_b(point, a, b); } -#if defined(BOOST_GEOMETRY_USE_RESCALING) - return; -#endif - +#ifndef BOOST_GEOMETRY_USE_RESCALING // Verify nearly collinear cases (the threshold is arbitrary // but influences performance). If the intersection is located // outside the segments, then it should be moved. @@ -232,6 +264,7 @@ struct cartesian_segments assign_if_exceeds(point, a); assign_if_exceeds(point, b); } +#endif } CoordinateType dx_a, dy_a; @@ -402,7 +435,9 @@ struct cartesian_segments return Policy::disjoint(); } - typedef side::side_by_triangle<CalculationType> side_strategy_type; + using side_strategy_type + = typename side::services::default_strategy + <cartesian_tag, CalculationType>::type; side_info sides; sides.set<0>(side_strategy_type::apply(q1, q2, p1), side_strategy_type::apply(q1, q2, p2)); @@ -415,7 +450,7 @@ struct cartesian_segments sides.set<1>(side_strategy_type::apply(p1, p2, q1), side_strategy_type::apply(p1, p2, q2)); - + if (sides.same<1>()) { // Both points are at same side of other segment, we can leave @@ -633,7 +668,7 @@ private: int const a2_wrt_b = position_value(oa_2, ob_1, ob_2); int const b1_wrt_a = position_value(ob_1, oa_1, oa_2); int const b2_wrt_a = position_value(ob_2, oa_1, oa_2); - + // fix the ratios if necessary // CONSIDER: fixing ratios also in other cases, if they're inconsistent // e.g. if ratio == 1 or 0 (so IP at the endpoint) @@ -650,7 +685,7 @@ private: { ra_from.assign(1, 1); rb_to.assign(0, 1); - } + } if (a2_wrt_b == 1) { |