| 1 | // Boost.Geometry (aka GGL, Generic Geometry Library) |
| 2 | |
| 3 | // Copyright (c) 2007-2014 Barend Gehrels, Amsterdam, the Netherlands. |
| 4 | // Copyright (c) 2013-2014 Adam Wulkiewicz, Lodz, Poland. |
| 5 | |
| 6 | // This file was modified by Oracle on 2014, 2016, 2017. |
| 7 | // Modifications copyright (c) 2014-2017, Oracle and/or its affiliates. |
| 8 | |
| 9 | // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle |
| 10 | // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle |
| 11 | |
| 12 | // Use, modification and distribution is subject to the Boost Software License, |
| 13 | // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at |
| 14 | // http://www.boost.org/LICENSE_1_0.txt) |
| 15 | |
| 16 | #ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_INTERSECTION_HPP |
| 17 | #define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_INTERSECTION_HPP |
| 18 | |
| 19 | #include <algorithm> |
| 20 | |
| 21 | #include <boost/geometry/core/exception.hpp> |
| 22 | |
| 23 | #include <boost/geometry/geometries/concepts/point_concept.hpp> |
| 24 | #include <boost/geometry/geometries/concepts/segment_concept.hpp> |
| 25 | |
| 26 | #include <boost/geometry/arithmetic/determinant.hpp> |
| 27 | #include <boost/geometry/algorithms/detail/assign_values.hpp> |
| 28 | #include <boost/geometry/algorithms/detail/assign_indexed_point.hpp> |
| 29 | #include <boost/geometry/algorithms/detail/equals/point_point.hpp> |
| 30 | #include <boost/geometry/algorithms/detail/recalculate.hpp> |
| 31 | |
| 32 | #include <boost/geometry/util/math.hpp> |
| 33 | #include <boost/geometry/util/promote_integral.hpp> |
| 34 | #include <boost/geometry/util/select_calculation_type.hpp> |
| 35 | |
| 36 | #include <boost/geometry/strategies/agnostic/point_in_poly_winding.hpp> |
| 37 | #include <boost/geometry/strategies/cartesian/area_surveyor.hpp> |
| 38 | #include <boost/geometry/strategies/cartesian/distance_pythagoras.hpp> |
| 39 | #include <boost/geometry/strategies/cartesian/envelope_segment.hpp> |
| 40 | #include <boost/geometry/strategies/cartesian/side_by_triangle.hpp> |
| 41 | #include <boost/geometry/strategies/covered_by.hpp> |
| 42 | #include <boost/geometry/strategies/intersection.hpp> |
| 43 | #include <boost/geometry/strategies/intersection_result.hpp> |
| 44 | #include <boost/geometry/strategies/side.hpp> |
| 45 | #include <boost/geometry/strategies/side_info.hpp> |
| 46 | #include <boost/geometry/strategies/within.hpp> |
| 47 | |
| 48 | #include <boost/geometry/policies/robustness/robust_point_type.hpp> |
| 49 | #include <boost/geometry/policies/robustness/segment_ratio_type.hpp> |
| 50 | |
| 51 | |
| 52 | #if defined(BOOST_GEOMETRY_DEBUG_ROBUSTNESS) |
| 53 | # include <boost/geometry/io/wkt/write.hpp> |
| 54 | #endif |
| 55 | |
| 56 | |
| 57 | namespace boost { namespace geometry |
| 58 | { |
| 59 | |
| 60 | |
| 61 | namespace strategy { namespace intersection |
| 62 | { |
| 63 | |
| 64 | |
| 65 | /*! |
| 66 | \see http://mathworld.wolfram.com/Line-LineIntersection.html |
| 67 | */ |
| 68 | template |
| 69 | < |
| 70 | typename CalculationType = void |
| 71 | > |
| 72 | struct cartesian_segments |
| 73 | { |
| 74 | typedef side::side_by_triangle<CalculationType> side_strategy_type; |
| 75 | |
| 76 | static inline side_strategy_type get_side_strategy() |
| 77 | { |
| 78 | return side_strategy_type(); |
| 79 | } |
| 80 | |
| 81 | template <typename Geometry1, typename Geometry2> |
| 82 | struct point_in_geometry_strategy |
| 83 | { |
| 84 | typedef strategy::within::winding |
| 85 | < |
| 86 | typename point_type<Geometry1>::type, |
| 87 | typename point_type<Geometry2>::type, |
| 88 | side_strategy_type, |
| 89 | CalculationType |
| 90 | > type; |
| 91 | }; |
| 92 | |
| 93 | template <typename Geometry1, typename Geometry2> |
| 94 | static inline typename point_in_geometry_strategy<Geometry1, Geometry2>::type |
| 95 | get_point_in_geometry_strategy() |
| 96 | { |
| 97 | typedef typename point_in_geometry_strategy |
| 98 | < |
| 99 | Geometry1, Geometry2 |
| 100 | >::type strategy_type; |
| 101 | return strategy_type(); |
| 102 | } |
| 103 | |
| 104 | template <typename Geometry> |
| 105 | struct area_strategy |
| 106 | { |
| 107 | typedef area::surveyor |
| 108 | < |
| 109 | typename point_type<Geometry>::type, |
| 110 | CalculationType |
| 111 | > type; |
| 112 | }; |
| 113 | |
| 114 | template <typename Geometry> |
| 115 | static inline typename area_strategy<Geometry>::type get_area_strategy() |
| 116 | { |
| 117 | typedef typename area_strategy<Geometry>::type strategy_type; |
| 118 | return strategy_type(); |
| 119 | } |
| 120 | |
| 121 | template <typename Geometry> |
| 122 | struct distance_strategy |
| 123 | { |
| 124 | typedef distance::pythagoras |
| 125 | < |
| 126 | CalculationType |
| 127 | > type; |
| 128 | }; |
| 129 | |
| 130 | template <typename Geometry> |
| 131 | static inline typename distance_strategy<Geometry>::type get_distance_strategy() |
| 132 | { |
| 133 | typedef typename distance_strategy<Geometry>::type strategy_type; |
| 134 | return strategy_type(); |
| 135 | } |
| 136 | |
| 137 | typedef envelope::cartesian_segment<CalculationType> |
| 138 | envelope_strategy_type; |
| 139 | |
| 140 | static inline envelope_strategy_type get_envelope_strategy() |
| 141 | { |
| 142 | return envelope_strategy_type(); |
| 143 | } |
| 144 | |
| 145 | template <typename CoordinateType, typename SegmentRatio> |
| 146 | struct segment_intersection_info |
| 147 | { |
| 148 | typedef typename select_most_precise |
| 149 | < |
| 150 | CoordinateType, double |
| 151 | >::type promoted_type; |
| 152 | |
| 153 | promoted_type comparable_length_a() const |
| 154 | { |
| 155 | return dx_a * dx_a + dy_a * dy_a; |
| 156 | } |
| 157 | |
| 158 | promoted_type comparable_length_b() const |
| 159 | { |
| 160 | return dx_b * dx_b + dy_b * dy_b; |
| 161 | } |
| 162 | |
| 163 | template <typename Point, typename Segment1, typename Segment2> |
| 164 | void assign_a(Point& point, Segment1 const& a, Segment2 const& ) const |
| 165 | { |
| 166 | assign(point, a, dx_a, dy_a, robust_ra); |
| 167 | } |
| 168 | template <typename Point, typename Segment1, typename Segment2> |
| 169 | void assign_b(Point& point, Segment1 const& , Segment2 const& b) const |
| 170 | { |
| 171 | assign(point, b, dx_b, dy_b, robust_rb); |
| 172 | } |
| 173 | |
| 174 | template <typename Point, typename Segment> |
| 175 | void assign(Point& point, Segment const& segment, CoordinateType const& dx, CoordinateType const& dy, SegmentRatio const& ratio) const |
| 176 | { |
| 177 | // Calculate the intersection point based on segment_ratio |
| 178 | // Up to now, division was postponed. Here we divide using numerator/ |
| 179 | // denominator. In case of integer this results in an integer |
| 180 | // division. |
| 181 | BOOST_GEOMETRY_ASSERT(ratio.denominator() != 0); |
| 182 | |
| 183 | typedef typename promote_integral<CoordinateType>::type promoted_type; |
| 184 | |
| 185 | promoted_type const numerator |
| 186 | = boost::numeric_cast<promoted_type>(ratio.numerator()); |
| 187 | promoted_type const denominator |
| 188 | = boost::numeric_cast<promoted_type>(ratio.denominator()); |
| 189 | promoted_type const dx_promoted = boost::numeric_cast<promoted_type>(dx); |
| 190 | promoted_type const dy_promoted = boost::numeric_cast<promoted_type>(dy); |
| 191 | |
| 192 | set<0>(point, get<0, 0>(segment) + boost::numeric_cast |
| 193 | < |
| 194 | CoordinateType |
| 195 | >(numerator * dx_promoted / denominator)); |
| 196 | set<1>(point, get<0, 1>(segment) + boost::numeric_cast |
| 197 | < |
| 198 | CoordinateType |
| 199 | >(numerator * dy_promoted / denominator)); |
| 200 | } |
| 201 | |
| 202 | CoordinateType dx_a, dy_a; |
| 203 | CoordinateType dx_b, dy_b; |
| 204 | SegmentRatio robust_ra; |
| 205 | SegmentRatio robust_rb; |
| 206 | }; |
| 207 | |
| 208 | template <typename D, typename W, typename ResultType> |
| 209 | static inline void cramers_rule(D const& dx_a, D const& dy_a, |
| 210 | D const& dx_b, D const& dy_b, W const& wx, W const& wy, |
| 211 | // out: |
| 212 | ResultType& d, ResultType& da) |
| 213 | { |
| 214 | // Cramers rule |
| 215 | d = geometry::detail::determinant<ResultType>(dx_a, dy_a, dx_b, dy_b); |
| 216 | da = geometry::detail::determinant<ResultType>(dx_b, dy_b, wx, wy); |
| 217 | // Ratio is da/d , collinear if d == 0, intersecting if 0 <= r <= 1 |
| 218 | // IntersectionPoint = (x1 + r * dx_a, y1 + r * dy_a) |
| 219 | } |
| 220 | |
| 221 | |
| 222 | // Relate segments a and b |
| 223 | template |
| 224 | < |
| 225 | typename Segment1, |
| 226 | typename Segment2, |
| 227 | typename Policy, |
| 228 | typename RobustPolicy |
| 229 | > |
| 230 | static inline typename Policy::return_type |
| 231 | apply(Segment1 const& a, Segment2 const& b, |
| 232 | Policy const& policy, RobustPolicy const& robust_policy) |
| 233 | { |
| 234 | // type them all as in Segment1 - TODO reconsider this, most precise? |
| 235 | typedef typename geometry::point_type<Segment1>::type point_type; |
| 236 | |
| 237 | typedef typename geometry::robust_point_type |
| 238 | < |
| 239 | point_type, RobustPolicy |
| 240 | >::type robust_point_type; |
| 241 | |
| 242 | point_type a0, a1, b0, b1; |
| 243 | robust_point_type a0_rob, a1_rob, b0_rob, b1_rob; |
| 244 | |
| 245 | detail::assign_point_from_index<0>(a, a0); |
| 246 | detail::assign_point_from_index<1>(a, a1); |
| 247 | detail::assign_point_from_index<0>(b, b0); |
| 248 | detail::assign_point_from_index<1>(b, b1); |
| 249 | |
| 250 | geometry::recalculate(a0_rob, a0, robust_policy); |
| 251 | geometry::recalculate(a1_rob, a1, robust_policy); |
| 252 | geometry::recalculate(b0_rob, b0, robust_policy); |
| 253 | geometry::recalculate(b1_rob, b1, robust_policy); |
| 254 | |
| 255 | return apply(a, b, policy, robust_policy, a0_rob, a1_rob, b0_rob, b1_rob); |
| 256 | } |
| 257 | |
| 258 | // The main entry-routine, calculating intersections of segments a / b |
| 259 | // NOTE: Robust* types may be the same as Segments' point types |
| 260 | template |
| 261 | < |
| 262 | typename Segment1, |
| 263 | typename Segment2, |
| 264 | typename Policy, |
| 265 | typename RobustPolicy, |
| 266 | typename RobustPoint1, |
| 267 | typename RobustPoint2 |
| 268 | > |
| 269 | static inline typename Policy::return_type |
| 270 | apply(Segment1 const& a, Segment2 const& b, |
| 271 | Policy const&, RobustPolicy const& /*robust_policy*/, |
| 272 | RobustPoint1 const& robust_a1, RobustPoint1 const& robust_a2, |
| 273 | RobustPoint2 const& robust_b1, RobustPoint2 const& robust_b2) |
| 274 | { |
| 275 | BOOST_CONCEPT_ASSERT( (concepts::ConstSegment<Segment1>) ); |
| 276 | BOOST_CONCEPT_ASSERT( (concepts::ConstSegment<Segment2>) ); |
| 277 | |
| 278 | using geometry::detail::equals::equals_point_point; |
| 279 | bool const a_is_point = equals_point_point(robust_a1, robust_a2); |
| 280 | bool const b_is_point = equals_point_point(robust_b1, robust_b2); |
| 281 | |
| 282 | if(a_is_point && b_is_point) |
| 283 | { |
| 284 | return equals_point_point(robust_a1, robust_b2) |
| 285 | ? Policy::degenerate(a, true) |
| 286 | : Policy::disjoint() |
| 287 | ; |
| 288 | } |
| 289 | |
| 290 | side_info sides; |
| 291 | sides.set<0>(side_strategy_type::apply(robust_b1, robust_b2, robust_a1), |
| 292 | side_strategy_type::apply(robust_b1, robust_b2, robust_a2)); |
| 293 | |
| 294 | if (sides.same<0>()) |
| 295 | { |
| 296 | // Both points are at same side of other segment, we can leave |
| 297 | return Policy::disjoint(); |
| 298 | } |
| 299 | |
| 300 | sides.set<1>(side_strategy_type::apply(robust_a1, robust_a2, robust_b1), |
| 301 | side_strategy_type::apply(robust_a1, robust_a2, robust_b2)); |
| 302 | |
| 303 | if (sides.same<1>()) |
| 304 | { |
| 305 | // Both points are at same side of other segment, we can leave |
| 306 | return Policy::disjoint(); |
| 307 | } |
| 308 | |
| 309 | bool collinear = sides.collinear(); |
| 310 | |
| 311 | typedef typename select_most_precise |
| 312 | < |
| 313 | typename geometry::coordinate_type<RobustPoint1>::type, |
| 314 | typename geometry::coordinate_type<RobustPoint2>::type |
| 315 | >::type robust_coordinate_type; |
| 316 | |
| 317 | typedef typename segment_ratio_type |
| 318 | < |
| 319 | typename geometry::point_type<Segment1>::type, // TODO: most precise point? |
| 320 | RobustPolicy |
| 321 | >::type ratio_type; |
| 322 | |
| 323 | segment_intersection_info |
| 324 | < |
| 325 | typename select_calculation_type<Segment1, Segment2, CalculationType>::type, |
| 326 | ratio_type |
| 327 | > sinfo; |
| 328 | |
| 329 | sinfo.dx_a = get<1, 0>(a) - get<0, 0>(a); // distance in x-dir |
| 330 | sinfo.dx_b = get<1, 0>(b) - get<0, 0>(b); |
| 331 | sinfo.dy_a = get<1, 1>(a) - get<0, 1>(a); // distance in y-dir |
| 332 | sinfo.dy_b = get<1, 1>(b) - get<0, 1>(b); |
| 333 | |
| 334 | robust_coordinate_type const robust_dx_a = get<0>(robust_a2) - get<0>(robust_a1); |
| 335 | robust_coordinate_type const robust_dx_b = get<0>(robust_b2) - get<0>(robust_b1); |
| 336 | robust_coordinate_type const robust_dy_a = get<1>(robust_a2) - get<1>(robust_a1); |
| 337 | robust_coordinate_type const robust_dy_b = get<1>(robust_b2) - get<1>(robust_b1); |
| 338 | |
| 339 | // r: ratio 0-1 where intersection divides A/B |
| 340 | // (only calculated for non-collinear segments) |
| 341 | if (! collinear) |
| 342 | { |
| 343 | robust_coordinate_type robust_da0, robust_da; |
| 344 | robust_coordinate_type robust_db0, robust_db; |
| 345 | |
| 346 | cramers_rule(robust_dx_a, robust_dy_a, robust_dx_b, robust_dy_b, |
| 347 | get<0>(robust_a1) - get<0>(robust_b1), |
| 348 | get<1>(robust_a1) - get<1>(robust_b1), |
| 349 | robust_da0, robust_da); |
| 350 | |
| 351 | cramers_rule(robust_dx_b, robust_dy_b, robust_dx_a, robust_dy_a, |
| 352 | get<0>(robust_b1) - get<0>(robust_a1), |
| 353 | get<1>(robust_b1) - get<1>(robust_a1), |
| 354 | robust_db0, robust_db); |
| 355 | |
| 356 | math::detail::equals_factor_policy<robust_coordinate_type> |
| 357 | policy(robust_dx_a, robust_dy_a, robust_dx_b, robust_dy_b); |
| 358 | robust_coordinate_type const zero = 0; |
| 359 | if (math::detail::equals_by_policy(robust_da0, zero, policy) |
| 360 | || math::detail::equals_by_policy(robust_db0, zero, policy)) |
| 361 | { |
| 362 | // If this is the case, no rescaling is done for FP precision. |
| 363 | // We set it to collinear, but it indicates a robustness issue. |
| 364 | sides.set<0>(first: 0,second: 0); |
| 365 | sides.set<1>(first: 0,second: 0); |
| 366 | collinear = true; |
| 367 | } |
| 368 | else |
| 369 | { |
| 370 | sinfo.robust_ra.assign(robust_da, robust_da0); |
| 371 | sinfo.robust_rb.assign(robust_db, robust_db0); |
| 372 | } |
| 373 | } |
| 374 | |
| 375 | if (collinear) |
| 376 | { |
| 377 | std::pair<bool, bool> const collinear_use_first |
| 378 | = is_x_more_significant(geometry::math::abs(robust_dx_a), |
| 379 | geometry::math::abs(robust_dy_a), |
| 380 | geometry::math::abs(robust_dx_b), |
| 381 | geometry::math::abs(robust_dy_b), |
| 382 | a_is_point, b_is_point); |
| 383 | |
| 384 | if (collinear_use_first.second) |
| 385 | { |
| 386 | // Degenerate cases: segments of single point, lying on other segment, are not disjoint |
| 387 | // This situation is collinear too |
| 388 | |
| 389 | if (collinear_use_first.first) |
| 390 | { |
| 391 | return relate_collinear<0, Policy, ratio_type>(a, b, |
| 392 | robust_a1, robust_a2, robust_b1, robust_b2, |
| 393 | a_is_point, b_is_point); |
| 394 | } |
| 395 | else |
| 396 | { |
| 397 | // Y direction contains larger segments (maybe dx is zero) |
| 398 | return relate_collinear<1, Policy, ratio_type>(a, b, |
| 399 | robust_a1, robust_a2, robust_b1, robust_b2, |
| 400 | a_is_point, b_is_point); |
| 401 | } |
| 402 | } |
| 403 | } |
| 404 | |
| 405 | return Policy::segments_crosses(sides, sinfo, a, b); |
| 406 | } |
| 407 | |
| 408 | private: |
| 409 | // first is true if x is more significant |
| 410 | // second is true if the more significant difference is not 0 |
| 411 | template <typename RobustCoordinateType> |
| 412 | static inline std::pair<bool, bool> |
| 413 | is_x_more_significant(RobustCoordinateType const& abs_robust_dx_a, |
| 414 | RobustCoordinateType const& abs_robust_dy_a, |
| 415 | RobustCoordinateType const& abs_robust_dx_b, |
| 416 | RobustCoordinateType const& abs_robust_dy_b, |
| 417 | bool const a_is_point, |
| 418 | bool const b_is_point) |
| 419 | { |
| 420 | //BOOST_GEOMETRY_ASSERT_MSG(!(a_is_point && b_is_point), "both segments shouldn't be degenerated"); |
| 421 | |
| 422 | // for degenerated segments the second is always true because this function |
| 423 | // shouldn't be called if both segments were degenerated |
| 424 | |
| 425 | if (a_is_point) |
| 426 | { |
| 427 | return std::make_pair(abs_robust_dx_b >= abs_robust_dy_b, true); |
| 428 | } |
| 429 | else if (b_is_point) |
| 430 | { |
| 431 | return std::make_pair(abs_robust_dx_a >= abs_robust_dy_a, true); |
| 432 | } |
| 433 | else |
| 434 | { |
| 435 | RobustCoordinateType const min_dx = (std::min)(abs_robust_dx_a, abs_robust_dx_b); |
| 436 | RobustCoordinateType const min_dy = (std::min)(abs_robust_dy_a, abs_robust_dy_b); |
| 437 | return min_dx == min_dy ? |
| 438 | std::make_pair(true, min_dx > RobustCoordinateType(0)) : |
| 439 | std::make_pair(min_dx > min_dy, true); |
| 440 | } |
| 441 | } |
| 442 | |
| 443 | template |
| 444 | < |
| 445 | std::size_t Dimension, |
| 446 | typename Policy, |
| 447 | typename RatioType, |
| 448 | typename Segment1, |
| 449 | typename Segment2, |
| 450 | typename RobustPoint1, |
| 451 | typename RobustPoint2 |
| 452 | > |
| 453 | static inline typename Policy::return_type |
| 454 | relate_collinear(Segment1 const& a, |
| 455 | Segment2 const& b, |
| 456 | RobustPoint1 const& robust_a1, RobustPoint1 const& robust_a2, |
| 457 | RobustPoint2 const& robust_b1, RobustPoint2 const& robust_b2, |
| 458 | bool a_is_point, bool b_is_point) |
| 459 | { |
| 460 | if (a_is_point) |
| 461 | { |
| 462 | return relate_one_degenerate<Policy, RatioType>(a, |
| 463 | get<Dimension>(robust_a1), |
| 464 | get<Dimension>(robust_b1), get<Dimension>(robust_b2), |
| 465 | true); |
| 466 | } |
| 467 | if (b_is_point) |
| 468 | { |
| 469 | return relate_one_degenerate<Policy, RatioType>(b, |
| 470 | get<Dimension>(robust_b1), |
| 471 | get<Dimension>(robust_a1), get<Dimension>(robust_a2), |
| 472 | false); |
| 473 | } |
| 474 | return relate_collinear<Policy, RatioType>(a, b, |
| 475 | get<Dimension>(robust_a1), |
| 476 | get<Dimension>(robust_a2), |
| 477 | get<Dimension>(robust_b1), |
| 478 | get<Dimension>(robust_b2)); |
| 479 | } |
| 480 | |
| 481 | /// Relate segments known collinear |
| 482 | template |
| 483 | < |
| 484 | typename Policy, |
| 485 | typename RatioType, |
| 486 | typename Segment1, |
| 487 | typename Segment2, |
| 488 | typename RobustType1, |
| 489 | typename RobustType2 |
| 490 | > |
| 491 | static inline typename Policy::return_type |
| 492 | relate_collinear(Segment1 const& a, Segment2 const& b, |
| 493 | RobustType1 oa_1, RobustType1 oa_2, |
| 494 | RobustType2 ob_1, RobustType2 ob_2) |
| 495 | { |
| 496 | // Calculate the ratios where a starts in b, b starts in a |
| 497 | // a1--------->a2 (2..7) |
| 498 | // b1----->b2 (5..8) |
| 499 | // length_a: 7-2=5 |
| 500 | // length_b: 8-5=3 |
| 501 | // b1 is located w.r.t. a at ratio: (5-2)/5=3/5 (on a) |
| 502 | // b2 is located w.r.t. a at ratio: (8-2)/5=6/5 (right of a) |
| 503 | // a1 is located w.r.t. b at ratio: (2-5)/3=-3/3 (left of b) |
| 504 | // a2 is located w.r.t. b at ratio: (7-5)/3=2/3 (on b) |
| 505 | // A arrives (a2 on b), B departs (b1 on a) |
| 506 | |
| 507 | // If both are reversed: |
| 508 | // a2<---------a1 (7..2) |
| 509 | // b2<-----b1 (8..5) |
| 510 | // length_a: 2-7=-5 |
| 511 | // length_b: 5-8=-3 |
| 512 | // b1 is located w.r.t. a at ratio: (8-7)/-5=-1/5 (before a starts) |
| 513 | // b2 is located w.r.t. a at ratio: (5-7)/-5=2/5 (on a) |
| 514 | // a1 is located w.r.t. b at ratio: (7-8)/-3=1/3 (on b) |
| 515 | // a2 is located w.r.t. b at ratio: (2-8)/-3=6/3 (after b ends) |
| 516 | |
| 517 | // If both one is reversed: |
| 518 | // a1--------->a2 (2..7) |
| 519 | // b2<-----b1 (8..5) |
| 520 | // length_a: 7-2=+5 |
| 521 | // length_b: 5-8=-3 |
| 522 | // b1 is located w.r.t. a at ratio: (8-2)/5=6/5 (after a ends) |
| 523 | // b2 is located w.r.t. a at ratio: (5-2)/5=3/5 (on a) |
| 524 | // a1 is located w.r.t. b at ratio: (2-8)/-3=6/3 (after b ends) |
| 525 | // a2 is located w.r.t. b at ratio: (7-8)/-3=1/3 (on b) |
| 526 | RobustType1 const length_a = oa_2 - oa_1; // no abs, see above |
| 527 | RobustType2 const length_b = ob_2 - ob_1; |
| 528 | |
| 529 | RatioType ra_from(oa_1 - ob_1, length_b); |
| 530 | RatioType ra_to(oa_2 - ob_1, length_b); |
| 531 | RatioType rb_from(ob_1 - oa_1, length_a); |
| 532 | RatioType rb_to(ob_2 - oa_1, length_a); |
| 533 | |
| 534 | // use absolute measure to detect endpoints intersection |
| 535 | // NOTE: it'd be possible to calculate bx_wrt_a using ax_wrt_b values |
| 536 | int const a1_wrt_b = position_value(oa_1, ob_1, ob_2); |
| 537 | int const a2_wrt_b = position_value(oa_2, ob_1, ob_2); |
| 538 | int const b1_wrt_a = position_value(ob_1, oa_1, oa_2); |
| 539 | int const b2_wrt_a = position_value(ob_2, oa_1, oa_2); |
| 540 | |
| 541 | // fix the ratios if necessary |
| 542 | // CONSIDER: fixing ratios also in other cases, if they're inconsistent |
| 543 | // e.g. if ratio == 1 or 0 (so IP at the endpoint) |
| 544 | // but position value indicates that the IP is in the middle of the segment |
| 545 | // because one of the segments is very long |
| 546 | // In such case the ratios could be moved into the middle direction |
| 547 | // by some small value (e.g. EPS+1ULP) |
| 548 | if (a1_wrt_b == 1) |
| 549 | { |
| 550 | ra_from.assign(0, 1); |
| 551 | rb_from.assign(0, 1); |
| 552 | } |
| 553 | else if (a1_wrt_b == 3) |
| 554 | { |
| 555 | ra_from.assign(1, 1); |
| 556 | rb_to.assign(0, 1); |
| 557 | } |
| 558 | |
| 559 | if (a2_wrt_b == 1) |
| 560 | { |
| 561 | ra_to.assign(0, 1); |
| 562 | rb_from.assign(1, 1); |
| 563 | } |
| 564 | else if (a2_wrt_b == 3) |
| 565 | { |
| 566 | ra_to.assign(1, 1); |
| 567 | rb_to.assign(1, 1); |
| 568 | } |
| 569 | |
| 570 | if ((a1_wrt_b < 1 && a2_wrt_b < 1) || (a1_wrt_b > 3 && a2_wrt_b > 3)) |
| 571 | //if ((ra_from.left() && ra_to.left()) || (ra_from.right() && ra_to.right())) |
| 572 | { |
| 573 | return Policy::disjoint(); |
| 574 | } |
| 575 | |
| 576 | bool const opposite = math::sign(length_a) != math::sign(length_b); |
| 577 | |
| 578 | return Policy::segments_collinear(a, b, opposite, |
| 579 | a1_wrt_b, a2_wrt_b, b1_wrt_a, b2_wrt_a, |
| 580 | ra_from, ra_to, rb_from, rb_to); |
| 581 | } |
| 582 | |
| 583 | /// Relate segments where one is degenerate |
| 584 | template |
| 585 | < |
| 586 | typename Policy, |
| 587 | typename RatioType, |
| 588 | typename DegenerateSegment, |
| 589 | typename RobustType1, |
| 590 | typename RobustType2 |
| 591 | > |
| 592 | static inline typename Policy::return_type |
| 593 | relate_one_degenerate(DegenerateSegment const& degenerate_segment, |
| 594 | RobustType1 d, RobustType2 s1, RobustType2 s2, |
| 595 | bool a_degenerate) |
| 596 | { |
| 597 | // Calculate the ratios where ds starts in s |
| 598 | // a1--------->a2 (2..6) |
| 599 | // b1/b2 (4..4) |
| 600 | // Ratio: (4-2)/(6-2) |
| 601 | RatioType const ratio(d - s1, s2 - s1); |
| 602 | |
| 603 | if (!ratio.on_segment()) |
| 604 | { |
| 605 | return Policy::disjoint(); |
| 606 | } |
| 607 | |
| 608 | return Policy::one_degenerate(degenerate_segment, ratio, a_degenerate); |
| 609 | } |
| 610 | |
| 611 | template <typename ProjCoord1, typename ProjCoord2> |
| 612 | static inline int position_value(ProjCoord1 const& ca1, |
| 613 | ProjCoord2 const& cb1, |
| 614 | ProjCoord2 const& cb2) |
| 615 | { |
| 616 | // S1x 0 1 2 3 4 |
| 617 | // S2 |----------> |
| 618 | return math::equals(ca1, cb1) ? 1 |
| 619 | : math::equals(ca1, cb2) ? 3 |
| 620 | : cb1 < cb2 ? |
| 621 | ( ca1 < cb1 ? 0 |
| 622 | : ca1 > cb2 ? 4 |
| 623 | : 2 ) |
| 624 | : ( ca1 > cb1 ? 0 |
| 625 | : ca1 < cb2 ? 4 |
| 626 | : 2 ); |
| 627 | } |
| 628 | }; |
| 629 | |
| 630 | |
| 631 | #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS |
| 632 | namespace services |
| 633 | { |
| 634 | |
| 635 | template <typename CalculationType> |
| 636 | struct default_strategy<cartesian_tag, CalculationType> |
| 637 | { |
| 638 | typedef cartesian_segments<CalculationType> type; |
| 639 | }; |
| 640 | |
| 641 | } // namespace services |
| 642 | #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS |
| 643 | |
| 644 | |
| 645 | }} // namespace strategy::intersection |
| 646 | |
| 647 | namespace strategy |
| 648 | { |
| 649 | |
| 650 | namespace within { namespace services |
| 651 | { |
| 652 | |
| 653 | template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> |
| 654 | struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, linear_tag, cartesian_tag, cartesian_tag> |
| 655 | { |
| 656 | typedef strategy::intersection::cartesian_segments<> type; |
| 657 | }; |
| 658 | |
| 659 | template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> |
| 660 | struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, polygonal_tag, cartesian_tag, cartesian_tag> |
| 661 | { |
| 662 | typedef strategy::intersection::cartesian_segments<> type; |
| 663 | }; |
| 664 | |
| 665 | template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> |
| 666 | struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, linear_tag, cartesian_tag, cartesian_tag> |
| 667 | { |
| 668 | typedef strategy::intersection::cartesian_segments<> type; |
| 669 | }; |
| 670 | |
| 671 | template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> |
| 672 | struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, polygonal_tag, cartesian_tag, cartesian_tag> |
| 673 | { |
| 674 | typedef strategy::intersection::cartesian_segments<> type; |
| 675 | }; |
| 676 | |
| 677 | }} // within::services |
| 678 | |
| 679 | namespace covered_by { namespace services |
| 680 | { |
| 681 | |
| 682 | template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> |
| 683 | struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, linear_tag, cartesian_tag, cartesian_tag> |
| 684 | { |
| 685 | typedef strategy::intersection::cartesian_segments<> type; |
| 686 | }; |
| 687 | |
| 688 | template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> |
| 689 | struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, polygonal_tag, cartesian_tag, cartesian_tag> |
| 690 | { |
| 691 | typedef strategy::intersection::cartesian_segments<> type; |
| 692 | }; |
| 693 | |
| 694 | template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> |
| 695 | struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, linear_tag, cartesian_tag, cartesian_tag> |
| 696 | { |
| 697 | typedef strategy::intersection::cartesian_segments<> type; |
| 698 | }; |
| 699 | |
| 700 | template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> |
| 701 | struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, polygonal_tag, cartesian_tag, cartesian_tag> |
| 702 | { |
| 703 | typedef strategy::intersection::cartesian_segments<> type; |
| 704 | }; |
| 705 | |
| 706 | }} // within::services |
| 707 | |
| 708 | } // strategy |
| 709 | |
| 710 | }} // namespace boost::geometry |
| 711 | |
| 712 | |
| 713 | #endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_INTERSECTION_HPP |
| 714 | |