ceph/src/boost/libs/geometry/include/boost/geometry/strategies/cartesian/cart_intersect.hpp

   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.
   7 // Modifications copyright (c) 2014-2016, 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 // Temporary / will be Strategy as template parameter
  37 #include <boost/geometry/strategies/side.hpp>
  38 #include <boost/geometry/strategies/cartesian/side_by_triangle.hpp>
  39
  40 #include <boost/geometry/strategies/side_info.hpp>
  41 #include <boost/geometry/strategies/intersection.hpp>
  42 #include <boost/geometry/strategies/intersection_result.hpp>
  43
  44 #include <boost/geometry/policies/robustness/robust_point_type.hpp>
  45 #include <boost/geometry/policies/robustness/segment_ratio_type.hpp>
  46
  47
  48 #if defined(BOOST_GEOMETRY_DEBUG_ROBUSTNESS)
  49 #  include <boost/geometry/io/wkt/write.hpp>
  50 #endif
  51
  52
  53 namespace boost { namespace geometry
  54 {
  55
  56
  57 namespace strategy { namespace intersection
  58 {
  59
  60
  61 /*!
  62     \see http://mathworld.wolfram.com/Line-LineIntersection.html
  63  */
  64 template <typename Policy, typename CalculationType = void>
  65 struct relate_cartesian_segments
  66 {
  67     typedef typename Policy::return_type return_type;
  68
  69     template <typename CoordinateType, typename SegmentRatio>
  70     struct segment_intersection_info
  71     {
  72         typedef typename select_most_precise
  73             <
  74                 CoordinateType, double
  75             >::type promoted_type;
  76
  77         promoted_type comparable_length_a() const
  78         {
  79             return dx_a * dx_a + dy_a * dy_a;
  80         }
  81
  82         promoted_type comparable_length_b() const
  83         {
  84             return dx_b * dx_b + dy_b * dy_b;
  85         }
  86
  87         template <typename Point, typename Segment1, typename Segment2>
  88         void assign_a(Point& point, Segment1 const& a, Segment2 const& ) const
  89         {
  90             assign(point, a, dx_a, dy_a, robust_ra);
  91         }
  92         template <typename Point, typename Segment1, typename Segment2>
  93         void assign_b(Point& point, Segment1 const& , Segment2 const& b) const
  94         {
  95             assign(point, b, dx_b, dy_b, robust_rb);
  96         }
  97
  98         template <typename Point, typename Segment>
  99         void assign(Point& point, Segment const& segment, CoordinateType const& dx, CoordinateType const& dy, SegmentRatio const& ratio) const
 100         {
 101             // Calculate the intersection point based on segment_ratio
 102             // Up to now, division was postponed. Here we divide using numerator/
 103             // denominator. In case of integer this results in an integer
 104             // division.
 105             BOOST_GEOMETRY_ASSERT(ratio.denominator() != 0);
 106
 107             typedef typename promote_integral<CoordinateType>::type promoted_type;
 108
 109             promoted_type const numerator
 110                 = boost::numeric_cast<promoted_type>(ratio.numerator());
 111             promoted_type const denominator
 112                 = boost::numeric_cast<promoted_type>(ratio.denominator());
 113             promoted_type const dx_promoted = boost::numeric_cast<promoted_type>(dx);
 114             promoted_type const dy_promoted = boost::numeric_cast<promoted_type>(dy);
 115
 116             set<0>(point, get<0, 0>(segment) + boost::numeric_cast
 117                 <
 118                     CoordinateType
 119                 >(numerator * dx_promoted / denominator));
 120             set<1>(point, get<0, 1>(segment) + boost::numeric_cast
 121                 <
 122                     CoordinateType
 123                 >(numerator * dy_promoted / denominator));
 124         }
 125
 126         CoordinateType dx_a, dy_a;
 127         CoordinateType dx_b, dy_b;
 128         SegmentRatio robust_ra;
 129         SegmentRatio robust_rb;
 130     };
 131
 132     template <typename D, typename W, typename ResultType>
 133     static inline void cramers_rule(D const& dx_a, D const& dy_a,
 134         D const& dx_b, D const& dy_b, W const& wx, W const& wy,
 135         // out:
 136         ResultType& d, ResultType& da)
 137     {
 138         // Cramers rule
 139         d = geometry::detail::determinant<ResultType>(dx_a, dy_a, dx_b, dy_b);
 140         da = geometry::detail::determinant<ResultType>(dx_b, dy_b, wx, wy);
 141         // Ratio is da/d , collinear if d == 0, intersecting if 0 <= r <= 1
 142         // IntersectionPoint = (x1 + r * dx_a, y1 + r * dy_a)
 143     }
 144
 145
 146     // Relate segments a and b
 147     template <typename Segment1, typename Segment2, typename RobustPolicy>
 148     static inline return_type apply(Segment1 const& a, Segment2 const& b,
 149                 RobustPolicy const& robust_policy)
 150     {
 151         // type them all as in Segment1 - TODO reconsider this, most precise?
 152         typedef typename geometry::point_type<Segment1>::type point_type;
 153
 154         typedef typename geometry::robust_point_type
 155             <
 156                 point_type, RobustPolicy
 157             >::type robust_point_type;
 158
 159         point_type a0, a1, b0, b1;
 160         robust_point_type a0_rob, a1_rob, b0_rob, b1_rob;
 161
 162         detail::assign_point_from_index<0>(a, a0);
 163         detail::assign_point_from_index<1>(a, a1);
 164         detail::assign_point_from_index<0>(b, b0);
 165         detail::assign_point_from_index<1>(b, b1);
 166
 167         geometry::recalculate(a0_rob, a0, robust_policy);
 168         geometry::recalculate(a1_rob, a1, robust_policy);
 169         geometry::recalculate(b0_rob, b0, robust_policy);
 170         geometry::recalculate(b1_rob, b1, robust_policy);
 171
 172         return apply(a, b, robust_policy, a0_rob, a1_rob, b0_rob, b1_rob);
 173     }
 174
 175     // The main entry-routine, calculating intersections of segments a / b
 176     // NOTE: Robust* types may be the same as Segments' point types
 177     template <typename Segment1, typename Segment2,
 178               typename RobustPolicy,
 179               typename RobustPoint1, typename RobustPoint2>
 180     static inline return_type apply(Segment1 const& a, Segment2 const& b,
 181             RobustPolicy const& /*robust_policy*/,
 182             RobustPoint1 const& robust_a1, RobustPoint1 const& robust_a2,
 183             RobustPoint2 const& robust_b1, RobustPoint2 const& robust_b2)
 184     {
 185         BOOST_CONCEPT_ASSERT( (concepts::ConstSegment<Segment1>) );
 186         BOOST_CONCEPT_ASSERT( (concepts::ConstSegment<Segment2>) );
 187
 188         using geometry::detail::equals::equals_point_point;
 189         bool const a_is_point = equals_point_point(robust_a1, robust_a2);
 190         bool const b_is_point = equals_point_point(robust_b1, robust_b2);
 191
 192         if(a_is_point && b_is_point)
 193         {
 194             return equals_point_point(robust_a1, robust_b2)
 195                 ? Policy::degenerate(a, true)
 196                 : Policy::disjoint()
 197                 ;
 198         }
 199
 200         typedef typename select_calculation_type
 201             <Segment1, Segment2, CalculationType>::type coordinate_type;
 202
 203         typedef side::side_by_triangle<coordinate_type> side;
 204
 205         side_info sides;
 206         sides.set<0>(side::apply(robust_b1, robust_b2, robust_a1),
 207                      side::apply(robust_b1, robust_b2, robust_a2));
 208
 209         if (sides.same<0>())
 210         {
 211             // Both points are at same side of other segment, we can leave
 212             return Policy::disjoint();
 213         }
 214
 215         sides.set<1>(side::apply(robust_a1, robust_a2, robust_b1),
 216                      side::apply(robust_a1, robust_a2, robust_b2));
 217
 218         if (sides.same<1>())
 219         {
 220             // Both points are at same side of other segment, we can leave
 221             return Policy::disjoint();
 222         }
 223
 224         bool collinear = sides.collinear();
 225
 226         typedef typename select_most_precise
 227             <
 228                 typename geometry::coordinate_type<RobustPoint1>::type,
 229                 typename geometry::coordinate_type<RobustPoint2>::type
 230             >::type robust_coordinate_type;
 231
 232         typedef typename segment_ratio_type
 233         <
 234             typename geometry::point_type<Segment1>::type, // TODO: most precise point?
 235             RobustPolicy
 236         >::type ratio_type;
 237
 238         segment_intersection_info
 239         <
 240             coordinate_type,
 241             ratio_type
 242         > sinfo;
 243
 244         sinfo.dx_a = get<1, 0>(a) - get<0, 0>(a); // distance in x-dir
 245         sinfo.dx_b = get<1, 0>(b) - get<0, 0>(b);
 246         sinfo.dy_a = get<1, 1>(a) - get<0, 1>(a); // distance in y-dir
 247         sinfo.dy_b = get<1, 1>(b) - get<0, 1>(b);
 248
 249         robust_coordinate_type const robust_dx_a = get<0>(robust_a2) - get<0>(robust_a1);
 250         robust_coordinate_type const robust_dx_b = get<0>(robust_b2) - get<0>(robust_b1);
 251         robust_coordinate_type const robust_dy_a = get<1>(robust_a2) - get<1>(robust_a1);
 252         robust_coordinate_type const robust_dy_b = get<1>(robust_b2) - get<1>(robust_b1);
 253
 254         // r: ratio 0-1 where intersection divides A/B
 255         // (only calculated for non-collinear segments)
 256         if (! collinear)
 257         {
 258             robust_coordinate_type robust_da0, robust_da;
 259             robust_coordinate_type robust_db0, robust_db;
 260
 261             cramers_rule(robust_dx_a, robust_dy_a, robust_dx_b, robust_dy_b,
 262                 get<0>(robust_a1) - get<0>(robust_b1),
 263                 get<1>(robust_a1) - get<1>(robust_b1),
 264                 robust_da0, robust_da);
 265
 266             cramers_rule(robust_dx_b, robust_dy_b, robust_dx_a, robust_dy_a,
 267                 get<0>(robust_b1) - get<0>(robust_a1),
 268                 get<1>(robust_b1) - get<1>(robust_a1),
 269                 robust_db0, robust_db);
 270
 271             math::detail::equals_factor_policy<robust_coordinate_type>
 272                 policy(robust_dx_a, robust_dy_a, robust_dx_b, robust_dy_b);
 273             robust_coordinate_type const zero = 0;
 274             if (math::detail::equals_by_policy(robust_da0, zero, policy)
 275              || math::detail::equals_by_policy(robust_db0, zero, policy))
 276             {
 277                 // If this is the case, no rescaling is done for FP precision.
 278                 // We set it to collinear, but it indicates a robustness issue.
 279                 sides.set<0>(0,0);
 280                 sides.set<1>(0,0);
 281                 collinear = true;
 282             }
 283             else
 284             {
 285                 sinfo.robust_ra.assign(robust_da, robust_da0);
 286                 sinfo.robust_rb.assign(robust_db, robust_db0);
 287             }
 288         }
 289
 290         if (collinear)
 291         {
 292             std::pair<bool, bool> const collinear_use_first
 293                     = is_x_more_significant(geometry::math::abs(robust_dx_a),
 294                                             geometry::math::abs(robust_dy_a),
 295                                             geometry::math::abs(robust_dx_b),
 296                                             geometry::math::abs(robust_dy_b),
 297                                             a_is_point, b_is_point);
 298
 299             if (collinear_use_first.second)
 300             {
 301                 // Degenerate cases: segments of single point, lying on other segment, are not disjoint
 302                 // This situation is collinear too
 303
 304                 if (collinear_use_first.first)
 305                 {
 306                     return relate_collinear<0, ratio_type>(a, b,
 307                             robust_a1, robust_a2, robust_b1, robust_b2,
 308                             a_is_point, b_is_point);
 309                 }
 310                 else
 311                 {
 312                     // Y direction contains larger segments (maybe dx is zero)
 313                     return relate_collinear<1, ratio_type>(a, b,
 314                             robust_a1, robust_a2, robust_b1, robust_b2,
 315                             a_is_point, b_is_point);
 316                 }
 317             }
 318         }
 319
 320         return Policy::segments_crosses(sides, sinfo, a, b);
 321     }
 322
 323 private:
 324     // first is true if x is more significant
 325     // second is true if the more significant difference is not 0
 326     template <typename RobustCoordinateType>
 327     static inline std::pair<bool, bool>
 328         is_x_more_significant(RobustCoordinateType const& abs_robust_dx_a,
 329                               RobustCoordinateType const& abs_robust_dy_a,
 330                               RobustCoordinateType const& abs_robust_dx_b,
 331                               RobustCoordinateType const& abs_robust_dy_b,
 332                               bool const a_is_point,
 333                               bool const b_is_point)
 334     {
 335         //BOOST_GEOMETRY_ASSERT_MSG(!(a_is_point && b_is_point), "both segments shouldn't be degenerated");
 336
 337         // for degenerated segments the second is always true because this function
 338         // shouldn't be called if both segments were degenerated
 339
 340         if (a_is_point)
 341         {
 342             return std::make_pair(abs_robust_dx_b >= abs_robust_dy_b, true);
 343         }
 344         else if (b_is_point)
 345         {
 346             return std::make_pair(abs_robust_dx_a >= abs_robust_dy_a, true);
 347         }
 348         else
 349         {
 350             RobustCoordinateType const min_dx = (std::min)(abs_robust_dx_a, abs_robust_dx_b);
 351             RobustCoordinateType const min_dy = (std::min)(abs_robust_dy_a, abs_robust_dy_b);
 352             return min_dx == min_dy ?
 353                     std::make_pair(true, min_dx > RobustCoordinateType(0)) :
 354                     std::make_pair(min_dx > min_dy, true);
 355         }
 356     }
 357
 358     template
 359     <
 360         std::size_t Dimension,
 361         typename RatioType,
 362         typename Segment1,
 363         typename Segment2,
 364         typename RobustPoint1,
 365         typename RobustPoint2
 366     >
 367     static inline return_type relate_collinear(Segment1 const& a,
 368             Segment2 const& b,
 369             RobustPoint1 const& robust_a1, RobustPoint1 const& robust_a2,
 370             RobustPoint2 const& robust_b1, RobustPoint2 const& robust_b2,
 371             bool a_is_point, bool b_is_point)
 372     {
 373         if (a_is_point)
 374         {
 375             return relate_one_degenerate<RatioType>(a,
 376                 get<Dimension>(robust_a1),
 377                 get<Dimension>(robust_b1), get<Dimension>(robust_b2),
 378                 true);
 379         }
 380         if (b_is_point)
 381         {
 382             return relate_one_degenerate<RatioType>(b,
 383                 get<Dimension>(robust_b1),
 384                 get<Dimension>(robust_a1), get<Dimension>(robust_a2),
 385                 false);
 386         }
 387         return relate_collinear<RatioType>(a, b,
 388                                 get<Dimension>(robust_a1),
 389                                 get<Dimension>(robust_a2),
 390                                 get<Dimension>(robust_b1),
 391                                 get<Dimension>(robust_b2));
 392     }
 393
 394     /// Relate segments known collinear
 395     template
 396     <
 397         typename RatioType,
 398         typename Segment1,
 399         typename Segment2,
 400         typename RobustType1,
 401         typename RobustType2
 402     >
 403     static inline return_type relate_collinear(Segment1 const& a
 404             , Segment2 const& b
 405             , RobustType1 oa_1, RobustType1 oa_2
 406             , RobustType2 ob_1, RobustType2 ob_2
 407             )
 408     {
 409         // Calculate the ratios where a starts in b, b starts in a
 410         //         a1--------->a2         (2..7)
 411         //                b1----->b2      (5..8)
 412         // length_a: 7-2=5
 413         // length_b: 8-5=3
 414         // b1 is located w.r.t. a at ratio: (5-2)/5=3/5 (on a)
 415         // b2 is located w.r.t. a at ratio: (8-2)/5=6/5 (right of a)
 416         // a1 is located w.r.t. b at ratio: (2-5)/3=-3/3 (left of b)
 417         // a2 is located w.r.t. b at ratio: (7-5)/3=2/3 (on b)
 418         // A arrives (a2 on b), B departs (b1 on a)
 419
 420         // If both are reversed:
 421         //         a2<---------a1         (7..2)
 422         //                b2<-----b1      (8..5)
 423         // length_a: 2-7=-5
 424         // length_b: 5-8=-3
 425         // b1 is located w.r.t. a at ratio: (8-7)/-5=-1/5 (before a starts)
 426         // b2 is located w.r.t. a at ratio: (5-7)/-5=2/5 (on a)
 427         // a1 is located w.r.t. b at ratio: (7-8)/-3=1/3 (on b)
 428         // a2 is located w.r.t. b at ratio: (2-8)/-3=6/3 (after b ends)
 429
 430         // If both one is reversed:
 431         //         a1--------->a2         (2..7)
 432         //                b2<-----b1      (8..5)
 433         // length_a: 7-2=+5
 434         // length_b: 5-8=-3
 435         // b1 is located w.r.t. a at ratio: (8-2)/5=6/5 (after a ends)
 436         // b2 is located w.r.t. a at ratio: (5-2)/5=3/5 (on a)
 437         // a1 is located w.r.t. b at ratio: (2-8)/-3=6/3 (after b ends)
 438         // a2 is located w.r.t. b at ratio: (7-8)/-3=1/3 (on b)
 439         RobustType1 const length_a = oa_2 - oa_1; // no abs, see above
 440         RobustType2 const length_b = ob_2 - ob_1;
 441
 442         RatioType ra_from(oa_1 - ob_1, length_b);
 443         RatioType ra_to(oa_2 - ob_1, length_b);
 444         RatioType rb_from(ob_1 - oa_1, length_a);
 445         RatioType rb_to(ob_2 - oa_1, length_a);
 446
 447         // use absolute measure to detect endpoints intersection
 448         // NOTE: it'd be possible to calculate bx_wrt_a using ax_wrt_b values
 449         int const a1_wrt_b = position_value(oa_1, ob_1, ob_2);
 450         int const a2_wrt_b = position_value(oa_2, ob_1, ob_2);
 451         int const b1_wrt_a = position_value(ob_1, oa_1, oa_2);
 452         int const b2_wrt_a = position_value(ob_2, oa_1, oa_2);
 453
 454         // fix the ratios if necessary
 455         // CONSIDER: fixing ratios also in other cases, if they're inconsistent
 456         // e.g. if ratio == 1 or 0 (so IP at the endpoint)
 457         // but position value indicates that the IP is in the middle of the segment
 458         // because one of the segments is very long
 459         // In such case the ratios could be moved into the middle direction
 460         // by some small value (e.g. EPS+1ULP)
 461         if (a1_wrt_b == 1)
 462         {
 463             ra_from.assign(0, 1);
 464             rb_from.assign(0, 1);
 465         }
 466         else if (a1_wrt_b == 3)
 467         {
 468             ra_from.assign(1, 1);
 469             rb_to.assign(0, 1);
 470         }
 471
 472         if (a2_wrt_b == 1)
 473         {
 474             ra_to.assign(0, 1);
 475             rb_from.assign(1, 1);
 476         }
 477         else if (a2_wrt_b == 3)
 478         {
 479             ra_to.assign(1, 1);
 480             rb_to.assign(1, 1);
 481         }
 482
 483         if ((a1_wrt_b < 1 && a2_wrt_b < 1) || (a1_wrt_b > 3 && a2_wrt_b > 3))
 484         //if ((ra_from.left() && ra_to.left()) || (ra_from.right() && ra_to.right()))
 485         {
 486             return Policy::disjoint();
 487         }
 488
 489         bool const opposite = math::sign(length_a) != math::sign(length_b);
 490
 491         return Policy::segments_collinear(a, b, opposite,
 492                                           a1_wrt_b, a2_wrt_b, b1_wrt_a, b2_wrt_a,
 493                                           ra_from, ra_to, rb_from, rb_to);
 494     }
 495
 496     /// Relate segments where one is degenerate
 497     template
 498     <
 499         typename RatioType,
 500         typename DegenerateSegment,
 501         typename RobustType1,
 502         typename RobustType2
 503     >
 504     static inline return_type relate_one_degenerate(
 505             DegenerateSegment const& degenerate_segment
 506             , RobustType1 d
 507             , RobustType2 s1, RobustType2 s2
 508             , bool a_degenerate
 509             )
 510     {
 511         // Calculate the ratios where ds starts in s
 512         //         a1--------->a2         (2..6)
 513         //              b1/b2      (4..4)
 514         // Ratio: (4-2)/(6-2)
 515         RatioType const ratio(d - s1, s2 - s1);
 516
 517         if (!ratio.on_segment())
 518         {
 519             return Policy::disjoint();
 520         }
 521
 522         return Policy::one_degenerate(degenerate_segment, ratio, a_degenerate);
 523     }
 524
 525     template <typename ProjCoord1, typename ProjCoord2>
 526     static inline int position_value(ProjCoord1 const& ca1,
 527                                      ProjCoord2 const& cb1,
 528                                      ProjCoord2 const& cb2)
 529     {
 530         // S1x  0   1    2     3   4
 531         // S2       |---------->
 532         return math::equals(ca1, cb1) ? 1
 533              : math::equals(ca1, cb2) ? 3
 534              : cb1 < cb2 ?
 535                 ( ca1 < cb1 ? 0
 536                 : ca1 > cb2 ? 4
 537                 : 2 )
 538               : ( ca1 > cb1 ? 0
 539                 : ca1 < cb2 ? 4
 540                 : 2 );
 541     }
 542 };
 543
 544
 545 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
 546 namespace services
 547 {
 548
 549 template <typename Policy, typename CalculationType>
 550 struct default_strategy<cartesian_tag, Policy, CalculationType>
 551 {
 552     typedef relate_cartesian_segments<Policy, CalculationType> type;
 553 };
 554
 555 } // namespace services
 556 #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
 557
 558
 559 }} // namespace strategy::intersection
 560
 561 }} // namespace boost::geometry
 562
 563
 564 #endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_INTERSECTION_HPP