diff options
Diffstat (limited to 'boost/geometry/strategies/cartesian/cart_intersect.hpp')
| -rw-r--r-- | boost/geometry/strategies/cartesian/cart_intersect.hpp | 862 |
1 files changed, 244 insertions, 618 deletions
diff --git a/boost/geometry/strategies/cartesian/cart_intersect.hpp b/boost/geometry/strategies/cartesian/cart_intersect.hpp index ea92cf37b0..66af2d2e9c 100644 --- a/boost/geometry/strategies/cartesian/cart_intersect.hpp +++ b/boost/geometry/strategies/cartesian/cart_intersect.hpp @@ -1,6 +1,7 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. +// Copyright (c) 2013 Adam Wulkiewicz, Lodz, Poland. // Use, modification and distribution is subject to the Boost Software License, // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at @@ -18,6 +19,9 @@ #include <boost/geometry/arithmetic/determinant.hpp> #include <boost/geometry/algorithms/detail/assign_values.hpp> +#include <boost/geometry/algorithms/detail/assign_indexed_point.hpp> +#include <boost/geometry/algorithms/detail/equals/point_point.hpp> +#include <boost/geometry/algorithms/detail/recalculate.hpp> #include <boost/geometry/util/math.hpp> #include <boost/geometry/util/select_calculation_type.hpp> @@ -27,59 +31,24 @@ #include <boost/geometry/strategies/cartesian/side_by_triangle.hpp> #include <boost/geometry/strategies/side_info.hpp> +#include <boost/geometry/strategies/intersection_result.hpp> - -namespace boost { namespace geometry -{ +#include <boost/geometry/policies/robustness/robust_point_type.hpp> +#include <boost/geometry/policies/robustness/segment_ratio_type.hpp> -namespace strategy { namespace intersection -{ +#if defined(BOOST_GEOMETRY_DEBUG_ROBUSTNESS) +# include <boost/geometry/io/wkt/write.hpp> +#endif -#ifndef DOXYGEN_NO_DETAIL -namespace detail +namespace boost { namespace geometry { -template <std::size_t Dimension, typename Segment, typename T> -static inline void segment_arrange(Segment const& s, T& s_1, T& s_2, bool& swapped) -{ - s_1 = get<0, Dimension>(s); - s_2 = get<1, Dimension>(s); - if (s_1 > s_2) - { - std::swap(s_1, s_2); - swapped = true; - } -} -template <std::size_t Index, typename Segment> -inline typename geometry::point_type<Segment>::type get_from_index( - Segment const& segment) +namespace strategy { namespace intersection { - typedef typename geometry::point_type<Segment>::type point_type; - point_type point; - geometry::detail::assign::assign_point_from_index - < - Segment, point_type, Index, 0, dimension<Segment>::type::value - >::apply(segment, point); - return point; -} -} -#endif - -/*** -template <typename T> -inline std::string rdebug(T const& value) -{ - if (math::equals(value, 0)) return "'0'"; - if (math::equals(value, 1)) return "'1'"; - if (value < 0) return "<0"; - if (value > 1) return ">1"; - return "<0..1>"; -} -***/ /*! \see http://mathworld.wolfram.com/Line-LineIntersection.html @@ -88,650 +57,307 @@ template <typename Policy, typename CalculationType = void> struct relate_cartesian_segments { typedef typename Policy::return_type return_type; - typedef typename Policy::segment_type1 segment_type1; - typedef typename Policy::segment_type2 segment_type2; - - //typedef typename point_type<segment_type1>::type point_type; - //BOOST_CONCEPT_ASSERT( (concept::Point<point_type>) ); - - BOOST_CONCEPT_ASSERT( (concept::ConstSegment<segment_type1>) ); - BOOST_CONCEPT_ASSERT( (concept::ConstSegment<segment_type2>) ); - - typedef typename select_calculation_type - <segment_type1, segment_type2, CalculationType>::type coordinate_type; - /// Relate segments a and b - static inline return_type apply(segment_type1 const& a, segment_type2 const& b) + template <typename D, typename W, typename ResultType> + static inline void cramers_rule(D const& dx_a, D const& dy_a, + D const& dx_b, D const& dy_b, W const& wx, W const& wy, + // out: + ResultType& d, ResultType& da) { - coordinate_type const dx_a = get<1, 0>(a) - get<0, 0>(a); // distance in x-dir - coordinate_type const dx_b = get<1, 0>(b) - get<0, 0>(b); - coordinate_type const dy_a = get<1, 1>(a) - get<0, 1>(a); // distance in y-dir - coordinate_type const dy_b = get<1, 1>(b) - get<0, 1>(b); - return apply(a, b, dx_a, dy_a, dx_b, dy_b); + // Cramers rule + d = geometry::detail::determinant<ResultType>(dx_a, dy_a, dx_b, dy_b); + da = geometry::detail::determinant<ResultType>(dx_b, dy_b, wx, wy); + // Ratio is da/d , collinear if d == 0, intersecting if 0 <= r <= 1 + // IntersectionPoint = (x1 + r * dx_a, y1 + r * dy_a) } - // Relate segments a and b using precalculated differences. - // This can save two or four subtractions in many cases - static inline return_type apply(segment_type1 const& a, segment_type2 const& b, - coordinate_type const& dx_a, coordinate_type const& dy_a, - coordinate_type const& dx_b, coordinate_type const& dy_b) + // Relate segments a and b + template <typename Segment1, typename Segment2, typename RobustPolicy> + static inline return_type apply(Segment1 const& a, Segment2 const& b, + RobustPolicy const& robust_policy) { - typedef side::side_by_triangle<coordinate_type> side; - side_info sides; - - bool collinear_use_first = math::abs(dx_a) + math::abs(dx_b) >= math::abs(dy_a) + math::abs(dy_b); - - sides.set<0> - ( - side::apply(detail::get_from_index<0>(b) - , detail::get_from_index<1>(b) - , detail::get_from_index<0>(a)), - side::apply(detail::get_from_index<0>(b) - , detail::get_from_index<1>(b) - , detail::get_from_index<1>(a)) - ); - sides.set<1> - ( - side::apply(detail::get_from_index<0>(a) - , detail::get_from_index<1>(a) - , detail::get_from_index<0>(b)), - side::apply(detail::get_from_index<0>(a) - , detail::get_from_index<1>(a) - , detail::get_from_index<1>(b)) - ); + // type them all as in Segment1 - TODO reconsider this, most precise? + typedef typename geometry::point_type<Segment1>::type point_type; - bool collinear = sides.collinear(); + typedef typename geometry::robust_point_type + < + point_type, RobustPolicy + >::type robust_point_type; - robustness_verify_collinear(a, b, sides, collinear); - robustness_verify_meeting(a, b, sides, collinear, collinear_use_first); + point_type a0, a1, b0, b1; + robust_point_type a0_rob, a1_rob, b0_rob, b1_rob; - if (sides.same<0>() || sides.same<1>()) - { - // Both points are at same side of other segment, we can leave - if (robustness_verify_same_side(a, b, sides)) - { - return Policy::disjoint(); - } - } + detail::assign_point_from_index<0>(a, a0); + detail::assign_point_from_index<1>(a, a1); + detail::assign_point_from_index<0>(b, b0); + detail::assign_point_from_index<1>(b, b1); - // Degenerate cases: segments of single point, lying on other segment, non disjoint - coordinate_type const zero = 0; - if (math::equals(dx_a, zero) && math::equals(dy_a, zero)) - { - return Policy::degenerate(a, true); - } - if (math::equals(dx_b, zero) && math::equals(dy_b, zero)) - { - return Policy::degenerate(b, false); - } + geometry::recalculate(a0_rob, a0, robust_policy); + geometry::recalculate(a1_rob, a1, robust_policy); + geometry::recalculate(b0_rob, b0, robust_policy); + geometry::recalculate(b1_rob, b1, robust_policy); - typedef typename select_most_precise - < - coordinate_type, double - >::type promoted_type; + return apply(a, b, robust_policy, a0_rob, a1_rob, b0_rob, b1_rob); + } - // r: ratio 0-1 where intersection divides A/B - // (only calculated for non-collinear segments) - promoted_type r; - if (! collinear) - { - // Calculate determinants - Cramers rule - coordinate_type const wx = get<0, 0>(a) - get<0, 0>(b); - coordinate_type const wy = get<0, 1>(a) - get<0, 1>(b); - coordinate_type const d = geometry::detail::determinant<coordinate_type>(dx_a, dy_a, dx_b, dy_b); - coordinate_type const da = geometry::detail::determinant<coordinate_type>(dx_b, dy_b, wx, wy); - - coordinate_type const zero = coordinate_type(); - if (math::equals(d, zero)) - { - // This is still a collinear case (because of FP imprecision this can occur here) - // sides.debug(); - sides.set<0>(0,0); - sides.set<1>(0,0); - collinear = true; - } - else - { - r = promoted_type(da) / promoted_type(d); + // The main entry-routine, calculating intersections of segments a / b + template <typename Segment1, typename Segment2, typename RobustPolicy, typename RobustPoint> + static inline return_type apply(Segment1 const& a, Segment2 const& b, + RobustPolicy const& robust_policy, + RobustPoint const& robust_a1, RobustPoint const& robust_a2, + RobustPoint const& robust_b1, RobustPoint const& robust_b2) + { + BOOST_CONCEPT_ASSERT( (concept::ConstSegment<Segment1>) ); + BOOST_CONCEPT_ASSERT( (concept::ConstSegment<Segment2>) ); - if (! robustness_verify_r(a, b, r)) - { - return Policy::disjoint(); - } + boost::ignore_unused_variable_warning(robust_policy); - robustness_handle_meeting(a, b, sides, dx_a, dy_a, wx, wy, d, r); + typedef typename select_calculation_type + <Segment1, Segment2, CalculationType>::type coordinate_type; - if (robustness_verify_disjoint_at_one_collinear(a, b, sides)) - { - return Policy::disjoint(); - } + using geometry::detail::equals::equals_point_point; + bool const a_is_point = equals_point_point(robust_a1, robust_a2); + bool const b_is_point = equals_point_point(robust_b1, robust_b2); - } - } + typedef side::side_by_triangle<coordinate_type> side; - if(collinear) + if(a_is_point && b_is_point) { - if (collinear_use_first) - { - return relate_collinear<0>(a, b); - } - else - { - // Y direction contains larger segments (maybe dx is zero) - return relate_collinear<1>(a, b); - } + return equals_point_point(robust_a1, robust_b2) + ? Policy::degenerate(a, true) + : Policy::disjoint() + ; } - return Policy::segments_intersect(sides, r, - dx_a, dy_a, dx_b, dy_b, - a, b); - } - -private : + side_info sides; + sides.set<0>(side::apply(robust_b1, robust_b2, robust_a1), + side::apply(robust_b1, robust_b2, robust_a2)); + sides.set<1>(side::apply(robust_a1, robust_a2, robust_b1), + side::apply(robust_a1, robust_a2, robust_b2)); + bool collinear = sides.collinear(); - // Ratio should lie between 0 and 1 - // Also these three conditions might be of FP imprecision, the segments were actually (nearly) collinear - template <typename T> - static inline bool robustness_verify_r( - segment_type1 const& a, segment_type2 const& b, - T& r) - { - T const zero = 0; - T const one = 1; - if (r < zero || r > one) + if (sides.same<0>() || sides.same<1>()) { - if (verify_disjoint<0>(a, b) || verify_disjoint<1>(a, b)) - { - // Can still be disjoint (even if not one is left or right from another) - // This is e.g. in case #snake4 of buffer test. - return false; - } - - //std::cout << "ROBUSTNESS: correction of r " << r << std::endl; - // sides.debug(); + // Both points are at same side of other segment, we can leave + return Policy::disjoint(); + } - // ROBUSTNESS: the r value can in epsilon-cases much larger than 1, while (with perfect arithmetic) - // it should be one. It can be 1.14 or even 1.98049 or 2 (while still intersecting) + typedef typename select_most_precise + < + coordinate_type, double + >::type promoted_type; - // If segments are crossing (we can see that with the sides) - // and one is inside the other, there must be an intersection point. - // We correct for that. - // This is (only) in case #ggl_list_20110820_christophe in unit tests + typedef typename geometry::coordinate_type + < + RobustPoint + >::type robust_coordinate_type; - // If segments are touching (two sides zero), of course they should intersect - // This is (only) in case #buffer_rt_i in the unit tests) + typedef typename segment_ratio_type + < + typename geometry::point_type<Segment1>::type, // TODO: most precise point? + RobustPolicy + >::type ratio_type; - // If one touches in the middle, they also should intersect (#buffer_rt_j) + segment_intersection_info + < + coordinate_type, + promoted_type, + ratio_type + > sinfo; - // Note that even for ttmath r is occasionally > 1, e.g. 1.0000000000000000000000036191231203575 + sinfo.dx_a = get<1, 0>(a) - get<0, 0>(a); // distance in x-dir + sinfo.dx_b = get<1, 0>(b) - get<0, 0>(b); + sinfo.dy_a = get<1, 1>(a) - get<0, 1>(a); // distance in y-dir + sinfo.dy_b = get<1, 1>(b) - get<0, 1>(b); - if (r > one) - { - r = one; - } - else if (r < zero) - { - r = zero; - } - } - return true; - } + robust_coordinate_type const robust_dx_a = get<0>(robust_a2) - get<0>(robust_a1); + robust_coordinate_type const robust_dx_b = get<0>(robust_b2) - get<0>(robust_b1); + robust_coordinate_type const robust_dy_a = get<1>(robust_a2) - get<1>(robust_a1); + robust_coordinate_type const robust_dy_b = get<1>(robust_b2) - get<1>(robust_b1); - static inline void robustness_verify_collinear( - segment_type1 const& a, segment_type2 const& b, - side_info& sides, - bool& collinear) - { - if ((sides.zero<0>() && ! sides.zero<1>()) || (sides.zero<1>() && ! sides.zero<0>())) + // r: ratio 0-1 where intersection divides A/B + // (only calculated for non-collinear segments) + if (! collinear) { - // If one of the segments is collinear, the other must be as well. - // So handle it as collinear. - // (In float/double epsilon margins it can easily occur that one or two of them are -1/1) - // sides.debug(); - sides.set<0>(0,0); - sides.set<1>(0,0); - collinear = true; - } - } + robust_coordinate_type robust_da0, robust_da; + robust_coordinate_type robust_db0, robust_db; - static inline void robustness_verify_meeting( - segment_type1 const& a, segment_type2 const& b, - side_info& sides, - bool& collinear, bool& collinear_use_first) - { - if (sides.meeting()) - { - // If two segments meet each other at their segment-points, two sides are zero, - // the other two are not (unless collinear but we don't mean those here). - // However, in near-epsilon ranges it can happen that two sides are zero - // but they do not meet at their segment-points. - // In that case they are nearly collinear and handled as such. - if (! point_equals - ( - select(sides.zero_index<0>(), a), - select(sides.zero_index<1>(), b) - ) - ) + cramers_rule(robust_dx_a, robust_dy_a, robust_dx_b, robust_dy_b, + get<0>(robust_a1) - get<0>(robust_b1), + get<1>(robust_a1) - get<1>(robust_b1), + robust_da0, robust_da); + + cramers_rule(robust_dx_b, robust_dy_b, robust_dx_a, robust_dy_a, + get<0>(robust_b1) - get<0>(robust_a1), + get<1>(robust_b1) - get<1>(robust_a1), + robust_db0, robust_db); + + if (robust_da0 == 0) { + // If this is the case, no rescaling is done for FP precision. + // We set it to collinear, but it indicates a robustness issue. sides.set<0>(0,0); sides.set<1>(0,0); collinear = true; - - if (collinear_use_first && analyse_equal<0>(a, b)) - { - collinear_use_first = false; - } - else if (! collinear_use_first && analyse_equal<1>(a, b)) - { - collinear_use_first = true; - } - } - } - } - - // Verifies and if necessary correct missed touch because of robustness - // This is the case at multi_polygon_buffer unittest #rt_m - static inline bool robustness_verify_same_side( - segment_type1 const& a, segment_type2 const& b, - side_info& sides) - { - int corrected = 0; - if (sides.one_touching<0>()) - { - if (point_equals( - select(sides.zero_index<0>(), a), - select(0, b) - )) - { - sides.correct_to_zero<1, 0>(); - corrected = 1; - } - if (point_equals - ( - select(sides.zero_index<0>(), a), - select(1, b) - )) + else { - sides.correct_to_zero<1, 1>(); - corrected = 2; + sinfo.robust_ra.assign(robust_da, robust_da0); + sinfo.robust_rb.assign(robust_db, robust_db0); } } - else if (sides.one_touching<1>()) + + if (collinear) { - if (point_equals( - select(sides.zero_index<1>(), b), - select(0, a) - )) + bool const collinear_use_first + = geometry::math::abs(robust_dx_a) + geometry::math::abs(robust_dx_b) + >= geometry::math::abs(robust_dy_a) + geometry::math::abs(robust_dy_b); + + // Degenerate cases: segments of single point, lying on other segment, are not disjoint + // This situation is collinear too + + if (collinear_use_first) { - sides.correct_to_zero<0, 0>(); - corrected = 3; + return relate_collinear<0, ratio_type>(a, b, + robust_a1, robust_a2, robust_b1, robust_b2, + a_is_point, b_is_point); } - if (point_equals - ( - select(sides.zero_index<1>(), b), - select(1, a) - )) + else { - sides.correct_to_zero<0, 1>(); - corrected = 4; + // Y direction contains larger segments (maybe dx is zero) + return relate_collinear<1, ratio_type>(a, b, + robust_a1, robust_a2, robust_b1, robust_b2, + a_is_point, b_is_point); } } - return corrected == 0; + return Policy::segments_crosses(sides, sinfo, a, b); } - static inline bool robustness_verify_disjoint_at_one_collinear( - segment_type1 const& a, segment_type2 const& b, - side_info const& sides) +private: + template + < + std::size_t Dimension, + typename RatioType, + typename Segment1, + typename Segment2, + typename RobustPoint + > + static inline return_type relate_collinear(Segment1 const& a, + Segment2 const& b, + RobustPoint const& robust_a1, RobustPoint const& robust_a2, + RobustPoint const& robust_b1, RobustPoint const& robust_b2, + bool a_is_point, bool b_is_point) { - if (sides.one_of_all_zero()) + if (a_is_point) { - if (verify_disjoint<0>(a, b) || verify_disjoint<1>(a, b)) - { - return true; - } + return relate_one_degenerate<RatioType>(a, + get<Dimension>(robust_a1), + get<Dimension>(robust_b1), get<Dimension>(robust_b2), + true); } - return false; - } - - - // If r is one, or zero, segments should meet and their endpoints. - // Robustness issue: check if this is really the case. - // It turns out to be no problem, see buffer test #rt_s1 (and there are many cases generated) - // It generates an "ends in the middle" situation which is correct. - template <typename T, typename R> - static inline void robustness_handle_meeting(segment_type1 const& a, segment_type2 const& b, - side_info& sides, - T const& dx_a, T const& dy_a, T const& wx, T const& wy, - T const& d, R const& r) - { - return; - - T const db = geometry::detail::determinant<T>(dx_a, dy_a, wx, wy); - - R const zero = 0; - R const one = 1; - if (math::equals(r, zero) || math::equals(r, one)) + if (b_is_point) { - R rb = db / d; - if (rb <= 0 || rb >= 1 || math::equals(rb, 0) || math::equals(rb, 1)) - { - if (sides.one_zero<0>() && ! sides.one_zero<1>()) // or vice versa - { -#if defined(BOOST_GEOMETRY_COUNT_INTERSECTION_EQUAL) - extern int g_count_intersection_equal; - g_count_intersection_equal++; -#endif - sides.debug(); - std::cout << "E r=" << r << " r.b=" << rb << " "; - } - } + return relate_one_degenerate<RatioType>(b, + get<Dimension>(robust_b1), + get<Dimension>(robust_a1), get<Dimension>(robust_a2), + false); } + return relate_collinear<RatioType>(a, b, + get<Dimension>(robust_a1), + get<Dimension>(robust_a2), + get<Dimension>(robust_b1), + get<Dimension>(robust_b2)); } - template <std::size_t Dimension> - static inline bool verify_disjoint(segment_type1 const& a, - segment_type2 const& b) - { - coordinate_type a_1, a_2, b_1, b_2; - bool a_swapped = false, b_swapped = false; - detail::segment_arrange<Dimension>(a, a_1, a_2, a_swapped); - detail::segment_arrange<Dimension>(b, b_1, b_2, b_swapped); - return math::smaller(a_2, b_1) || math::larger(a_1, b_2); - } - - template <typename Segment> - static inline typename point_type<Segment>::type select(int index, Segment const& segment) - { - return index == 0 - ? detail::get_from_index<0>(segment) - : detail::get_from_index<1>(segment) - ; - } - - // We cannot use geometry::equals here. Besides that this will be changed - // to compare segment-coordinate-values directly (not necessary to retrieve point first) - template <typename Point1, typename Point2> - static inline bool point_equals(Point1 const& point1, Point2 const& point2) - { - return math::equals(get<0>(point1), get<0>(point2)) - && math::equals(get<1>(point1), get<1>(point2)) - ; - } - - // We cannot use geometry::equals here. Besides that this will be changed - // to compare segment-coordinate-values directly (not necessary to retrieve point first) - template <typename Point1, typename Point2> - static inline bool point_equality(Point1 const& point1, Point2 const& point2, - bool& equals_0, bool& equals_1) - { - equals_0 = math::equals(get<0>(point1), get<0>(point2)); - equals_1 = math::equals(get<1>(point1), get<1>(point2)); - return equals_0 && equals_1; - } - - template <std::size_t Dimension> - static inline bool analyse_equal(segment_type1 const& a, segment_type2 const& b) - { - coordinate_type const a_1 = geometry::get<0, Dimension>(a); - coordinate_type const a_2 = geometry::get<1, Dimension>(a); - coordinate_type const b_1 = geometry::get<0, Dimension>(b); - coordinate_type const b_2 = geometry::get<1, Dimension>(b); - return math::equals(a_1, b_1) - || math::equals(a_2, b_1) - || math::equals(a_1, b_2) - || math::equals(a_2, b_2) - ; - } - - template <std::size_t Dimension> - static inline return_type relate_collinear(segment_type1 const& a, - segment_type2 const& b) + /// Relate segments known collinear + template + < + typename RatioType, + typename Segment1, + typename Segment2, + typename RobustType + > + static inline return_type relate_collinear(Segment1 const& a + , Segment2 const& b + , RobustType oa_1, RobustType oa_2 + , RobustType ob_1, RobustType ob_2 + ) { - coordinate_type a_1, a_2, b_1, b_2; - bool a_swapped = false, b_swapped = false; - detail::segment_arrange<Dimension>(a, a_1, a_2, a_swapped); - detail::segment_arrange<Dimension>(b, b_1, b_2, b_swapped); - if (math::smaller(a_2, b_1) || math::larger(a_1, b_2)) - //if (a_2 < b_1 || a_1 > b_2) + // Calculate the ratios where a starts in b, b starts in a + // a1--------->a2 (2..7) + // b1----->b2 (5..8) + // length_a: 7-2=5 + // length_b: 8-5=3 + // b1 is located w.r.t. a at ratio: (5-2)/5=3/5 (on a) + // b2 is located w.r.t. a at ratio: (8-2)/5=6/5 (right of a) + // a1 is located w.r.t. b at ratio: (2-5)/3=-3/3 (left of b) + // a2 is located w.r.t. b at ratio: (7-5)/3=2/3 (on b) + // A arrives (a2 on b), B departs (b1 on a) + + // If both are reversed: + // a2<---------a1 (7..2) + // b2<-----b1 (8..5) + // length_a: 2-7=-5 + // length_b: 5-8=-3 + // b1 is located w.r.t. a at ratio: (8-7)/-5=-1/5 (before a starts) + // b2 is located w.r.t. a at ratio: (5-7)/-5=2/5 (on a) + // a1 is located w.r.t. b at ratio: (7-8)/-3=1/3 (on b) + // a2 is located w.r.t. b at ratio: (2-8)/-3=6/3 (after b ends) + + // If both one is reversed: + // a1--------->a2 (2..7) + // b2<-----b1 (8..5) + // length_a: 7-2=+5 + // length_b: 5-8=-3 + // b1 is located w.r.t. a at ratio: (8-2)/5=6/5 (after a ends) + // b2 is located w.r.t. a at ratio: (5-2)/5=3/5 (on a) + // a1 is located w.r.t. b at ratio: (2-8)/-3=6/3 (after b ends) + // a2 is located w.r.t. b at ratio: (7-8)/-3=1/3 (on b) + RobustType const length_a = oa_2 - oa_1; // no abs, see above + RobustType const length_b = ob_2 - ob_1; + + RatioType const ra_from(oa_1 - ob_1, length_b); + RatioType const ra_to(oa_2 - ob_1, length_b); + RatioType const rb_from(ob_1 - oa_1, length_a); + RatioType const rb_to(ob_2 - oa_1, length_a); + + if ((ra_from.left() && ra_to.left()) || (ra_from.right() && ra_to.right())) { return Policy::disjoint(); } - return relate_collinear(a, b, a_1, a_2, b_1, b_2, a_swapped, b_swapped); + + return Policy::segments_collinear(a, b, ra_from, ra_to, rb_from, rb_to); } - /// Relate segments known collinear - static inline return_type relate_collinear(segment_type1 const& a - , segment_type2 const& b - , coordinate_type a_1, coordinate_type a_2 - , coordinate_type b_1, coordinate_type b_2 - , bool a_swapped, bool b_swapped) + /// Relate segments where one is degenerate + template + < + typename RatioType, + typename DegenerateSegment, + typename RobustType + > + static inline return_type relate_one_degenerate( + DegenerateSegment const& degenerate_segment + , RobustType d + , RobustType s1, RobustType s2 + , bool a_degenerate + ) { - // All ca. 150 lines are about collinear rays - // The intersections, if any, are always boundary points of the segments. No need to calculate anything. - // However we want to find out HOW they intersect, there are many cases. - // Most sources only provide the intersection (above) or that there is a collinearity (but not the points) - // or some spare sources give the intersection points (calculated) but not how they align. - // This source tries to give everything and still be efficient. - // It is therefore (and because of the extensive clarification comments) rather long... - - // \see http://mpa.itc.it/radim/g50history/CMP/4.2.1-CERL-beta-libes/file475.txt - // \see http://docs.codehaus.org/display/GEOTDOC/Point+Set+Theory+and+the+DE-9IM+Matrix - // \see http://mathworld.wolfram.com/Line-LineIntersection.html - - // Because of collinearity the case is now one-dimensional and can be checked using intervals - // This function is called either horizontally or vertically - // We get then two intervals: - // a_1-------------a_2 where a_1 < a_2 - // b_1-------------b_2 where b_1 < b_2 - // In all figures below a_1/a_2 denotes arranged intervals, a1-a2 or a2-a1 are still unarranged - - // Handle "equal", in polygon neighbourhood comparisons a common case - - bool const opposite = a_swapped ^ b_swapped; - bool const both_swapped = a_swapped && b_swapped; - - // Check if segments are equal or opposite equal... - bool const swapped_a1_eq_b1 = math::equals(a_1, b_1); - bool const swapped_a2_eq_b2 = math::equals(a_2, b_2); - - if (swapped_a1_eq_b1 && swapped_a2_eq_b2) - { - return Policy::segment_equal(a, opposite); - } - - bool const swapped_a2_eq_b1 = math::equals(a_2, b_1); - bool const swapped_a1_eq_b2 = math::equals(a_1, b_2); - - bool const a1_eq_b1 = both_swapped ? swapped_a2_eq_b2 : a_swapped ? swapped_a2_eq_b1 : b_swapped ? swapped_a1_eq_b2 : swapped_a1_eq_b1; - bool const a2_eq_b2 = both_swapped ? swapped_a1_eq_b1 : a_swapped ? swapped_a1_eq_b2 : b_swapped ? swapped_a2_eq_b1 : swapped_a2_eq_b2; - - bool const a1_eq_b2 = both_swapped ? swapped_a2_eq_b1 : a_swapped ? swapped_a2_eq_b2 : b_swapped ? swapped_a1_eq_b1 : swapped_a1_eq_b2; - bool const a2_eq_b1 = both_swapped ? swapped_a1_eq_b2 : a_swapped ? swapped_a1_eq_b1 : b_swapped ? swapped_a2_eq_b2 : swapped_a2_eq_b1; - - - - - // The rest below will return one or two intersections. - // The delegated class can decide which is the intersection point, or two, build the Intersection Matrix (IM) - // For IM it is important to know which relates to which. So this information is given, - // without performance penalties to intersection calculation - - bool const has_common_points = swapped_a1_eq_b1 || swapped_a1_eq_b2 || swapped_a2_eq_b1 || swapped_a2_eq_b2; - - - // "Touch" -> one intersection point -> one but not two common points - // --------> A (or B) - // <---------- B (or A) - // a_2==b_1 (b_2==a_1 or a_2==b1) - - // The check a_2/b_1 is necessary because it excludes cases like - // -------> - // ---> - // ... which are handled lateron - - // Corresponds to 4 cases, of which the equal points are determined above - // #1: a1---->a2 b1--->b2 (a arrives at b's border) - // #2: a2<----a1 b2<---b1 (b arrives at a's border) - // #3: a1---->a2 b2<---b1 (both arrive at each others border) - // #4: a2<----a1 b1--->b2 (no arrival at all) - // Where the arranged forms have two forms: - // a_1-----a_2/b_1-------b_2 or reverse (B left of A) - if ((swapped_a2_eq_b1 || swapped_a1_eq_b2) && ! swapped_a1_eq_b1 && ! swapped_a2_eq_b2) - { - if (a2_eq_b1) return Policy::collinear_touch(get<1, 0>(a), get<1, 1>(a), 0, -1); - if (a1_eq_b2) return Policy::collinear_touch(get<0, 0>(a), get<0, 1>(a), -1, 0); - if (a2_eq_b2) return Policy::collinear_touch(get<1, 0>(a), get<1, 1>(a), 0, 0); - if (a1_eq_b1) return Policy::collinear_touch(get<0, 0>(a), get<0, 1>(a), -1, -1); - } - - - // "Touch/within" -> there are common points and also an intersection of interiors: - // Corresponds to many cases: - // #1a: a1------->a2 #1b: a1-->a2 - // b1--->b2 b1------->b2 - // #2a: a2<-------a1 #2b: a2<--a1 - // b1--->b2 b1------->b2 - // #3a: a1------->a2 #3b: a1-->a2 - // b2<---b1 b2<-------b1 - // #4a: a2<-------a1 #4b: a2<--a1 - // b2<---b1 b2<-------b1 - - // Note: next cases are similar and handled by the code - // #4c: a1--->a2 - // b1-------->b2 - // #4d: a1-------->a2 - // b1-->b2 - - // For case 1-4: a_1 < (b_1 or b_2) < a_2, two intersections are equal to segment B - // For case 5-8: b_1 < (a_1 or a_2) < b_2, two intersections are equal to segment A - if (has_common_points) - { - // Either A is in B, or B is in A, or (in case of robustness/equals) - // both are true, see below - bool a_in_b = (b_1 < a_1 && a_1 < b_2) || (b_1 < a_2 && a_2 < b_2); - bool b_in_a = (a_1 < b_1 && b_1 < a_2) || (a_1 < b_2 && b_2 < a_2); - - if (a_in_b && b_in_a) - { - // testcase "ggl_list_20110306_javier" - // In robustness it can occur that a point of A is inside B AND a point of B is inside A, - // still while has_common_points is true (so one point equals the other). - // If that is the case we select on length. - coordinate_type const length_a = geometry::math::abs(a_1 - a_2); - coordinate_type const length_b = geometry::math::abs(b_1 - b_2); - if (length_a > length_b) - { - a_in_b = false; - } - else - { - b_in_a = false; - } - } - - int const arrival_a = a_in_b ? 1 : -1; - if (a2_eq_b2) return Policy::collinear_interior_boundary_intersect(a_in_b ? a : b, a_in_b, 0, 0, false); - if (a1_eq_b2) return Policy::collinear_interior_boundary_intersect(a_in_b ? a : b, a_in_b, arrival_a, 0, true); - if (a2_eq_b1) return Policy::collinear_interior_boundary_intersect(a_in_b ? a : b, a_in_b, 0, -arrival_a, true); - if (a1_eq_b1) return Policy::collinear_interior_boundary_intersect(a_in_b ? a : b, a_in_b, arrival_a, -arrival_a, false); - } - - - - // "Inside", a completely within b or b completely within a - // 2 cases: - // case 1: - // a_1---a_2 -> take A's points as intersection points - // b_1------------b_2 - // case 2: - // a_1------------a_2 - // b_1---b_2 -> take B's points - if (a_1 > b_1 && a_2 < b_2) - { - // A within B - return Policy::collinear_a_in_b(a, opposite); - } - if (b_1 > a_1 && b_2 < a_2) - { - // B within A - return Policy::collinear_b_in_a(b, opposite); - } - - - /* - - Now that all cases with equal,touch,inside,disjoint, - degenerate are handled the only thing left is an overlap - - Either a1 is between b1,b2 - or a2 is between b1,b2 (a2 arrives) - - Next table gives an overview. - The IP's are ordered following the line A1->A2 - - | | - | a_2 in between | a_1 in between - | | - -----+---------------------------------+-------------------------- - | a1--------->a2 | a1--------->a2 - | b1----->b2 | b1----->b2 - | (b1,a2), a arrives | (a1,b2), b arrives - | | - -----+---------------------------------+-------------------------- - a sw.| a2<---------a1* | a2<---------a1* - | b1----->b2 | b1----->b2 - | (a1,b1), no arrival | (b2,a2), a and b arrive - | | - -----+---------------------------------+-------------------------- - | a1--------->a2 | a1--------->a2 - b sw.| b2<-----b1 | b2<-----b1 - | (b2,a2), a and b arrive | (a1,b1), no arrival - | | - -----+---------------------------------+-------------------------- - a sw.| a2<---------a1* | a2<---------a1* - b sw.| b2<-----b1 | b2<-----b1 - | (a1,b2), b arrives | (b1,a2), a arrives - | | - -----+---------------------------------+-------------------------- - * Note that a_1 < a_2, and a1 <> a_1; if a is swapped, - the picture might seem wrong but it (supposed to be) is right. - */ - - if (b_1 < a_2 && a_2 < b_2) - { - // Left column, from bottom to top - return - both_swapped ? Policy::collinear_overlaps(get<0, 0>(a), get<0, 1>(a), get<1, 0>(b), get<1, 1>(b), -1, 1, opposite) - : b_swapped ? Policy::collinear_overlaps(get<1, 0>(b), get<1, 1>(b), get<1, 0>(a), get<1, 1>(a), 1, 1, opposite) - : a_swapped ? Policy::collinear_overlaps(get<0, 0>(a), get<0, 1>(a), get<0, 0>(b), get<0, 1>(b), -1, -1, opposite) - : Policy::collinear_overlaps(get<0, 0>(b), get<0, 1>(b), get<1, 0>(a), get<1, 1>(a), 1, -1, opposite) - ; - } - if (b_1 < a_1 && a_1 < b_2) - { - // Right column, from bottom to top - return - both_swapped ? Policy::collinear_overlaps(get<0, 0>(b), get<0, 1>(b), get<1, 0>(a), get<1, 1>(a), 1, -1, opposite) - : b_swapped ? Policy::collinear_overlaps(get<0, 0>(a), get<0, 1>(a), get<0, 0>(b), get<0, 1>(b), -1, -1, opposite) - : a_swapped ? Policy::collinear_overlaps(get<1, 0>(b), get<1, 1>(b), get<1, 0>(a), get<1, 1>(a), 1, 1, opposite) - : Policy::collinear_overlaps(get<0, 0>(a), get<0, 1>(a), get<1, 0>(b), get<1, 1>(b), -1, 1, opposite) - ; - } - // Nothing should goes through. If any we have made an error - // std::cout << "Robustness issue, non-logical behaviour" << std::endl; - return Policy::error("Robustness issue, non-logical behaviour"); + // Calculate the ratios where ds starts in s + // a1--------->a2 (2..6) + // b1/b2 (4..4) + // Ratio: (4-2)/(6-2) + RatioType const ratio(d - s1, s2 - s1); + return Policy::one_degenerate(degenerate_segment, ratio, a_degenerate); } }; }} // namespace strategy::intersection - - }} // namespace boost::geometry |