ceph/src/boost/boost/geometry/strategies/cartesian/intersection.hpp

   1 // Boost.Geometry (aka GGL, Generic Geometry Library)
   2
   3 // Copyright (c) 2007-2014 Barend Gehrels, Amsterdam, the Netherlands.
   4 // Copyright (c) 2013-2017 Adam Wulkiewicz, Lodz, Poland.
   5
   6 // This file was modified by Oracle on 2014-2021.
   7 // Modifications copyright (c) 2014-2021, 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 #include <boost/geometry/geometries/segment.hpp>
  26
  27 #include <boost/geometry/arithmetic/determinant.hpp>
  28 #include <boost/geometry/algorithms/detail/assign_values.hpp>
  29 #include <boost/geometry/algorithms/detail/assign_indexed_point.hpp>
  30 #include <boost/geometry/algorithms/detail/equals/point_point.hpp>
  31 #include <boost/geometry/algorithms/detail/recalculate.hpp>
  32
  33 #include <boost/geometry/util/math.hpp>
  34 #include <boost/geometry/util/promote_integral.hpp>
  35 #include <boost/geometry/util/select_calculation_type.hpp>
  36
  37 #include <boost/geometry/strategy/cartesian/area.hpp>
  38 #include <boost/geometry/strategy/cartesian/envelope.hpp>
  39 #include <boost/geometry/strategy/cartesian/expand_box.hpp>
  40 #include <boost/geometry/strategy/cartesian/expand_segment.hpp>
  41
  42 #include <boost/geometry/strategies/cartesian/disjoint_box_box.hpp>
  43 #include <boost/geometry/strategies/cartesian/disjoint_segment_box.hpp>
  44 #include <boost/geometry/strategies/cartesian/distance_pythagoras.hpp>
  45 #include <boost/geometry/strategies/cartesian/point_in_point.hpp>
  46 #include <boost/geometry/strategies/cartesian/point_in_poly_winding.hpp>
  47 #include <boost/geometry/strategies/covered_by.hpp>
  48 #include <boost/geometry/strategies/intersection.hpp>
  49 #include <boost/geometry/strategies/intersection_result.hpp>
  50 #include <boost/geometry/strategies/side.hpp>
  51 #include <boost/geometry/strategies/side_info.hpp>
  52 #include <boost/geometry/strategies/within.hpp>
  53
  54 #include <boost/geometry/policies/robustness/rescale_policy_tags.hpp>
  55 #include <boost/geometry/policies/robustness/robust_point_type.hpp>
  56
  57
  58 #if defined(BOOST_GEOMETRY_DEBUG_ROBUSTNESS)
  59 #  include <boost/geometry/io/wkt/write.hpp>
  60 #endif
  61
  62
  63 namespace boost { namespace geometry
  64 {
  65
  66
  67 namespace strategy { namespace intersection
  68 {
  69
  70 namespace detail_usage
  71 {
  72
  73 // When calculating the intersection, the information of "a" or "b" can be used.
  74 // Theoretically this gives equal results, but due to floating point precision
  75 // there might be tiny differences. These are edge cases.
  76 // This structure is to determine if "a" or "b" should be used.
  77 // Prefer the segment closer to the endpoint.
  78 // If both are about equally close, then prefer the longer segment
  79 // To avoid hard thresholds, behavior is made fluent.
  80 // Calculate comparable length indications,
  81 // the longer the segment (relatively), the lower the value
  82 // such that the shorter lengths are evaluated higher and will
  83 // be preferred.
  84 template <bool IsArithmetic>
  85 struct use_a
  86 {
  87   template <typename Ct, typename Ev>
  88   static bool apply(Ct const& cla, Ct const& clb, Ev const& eva, Ev const& evb)
  89   {
  90       auto const clm = (std::max)(cla, clb);
  91       if (clm <= 0)
  92       {
  93           return true;
  94       }
  95
  96       // Relative comparible length
  97       auto const rcla = Ct(1.0) - cla / clm;
  98       auto const rclb = Ct(1.0) - clb / clm;
  99
 100       // Multipliers for edgevalue (ev) and relative comparible length (rcl)
 101       // They determine the balance between edge value (should be larger)
 102       // and segment length. In 99.9xx% of the cases there is no difference
 103       // at all (if either a or b is used). Therefore the values of the
 104       // constants are not sensitive for the majority of the situations.
 105       // One known case is #mysql_23023665_6 (difference) which needs mev >= 2
 106       Ev const mev = 5;
 107       Ev const mrcl = 1;
 108
 109       return mev * eva + mrcl * rcla > mev * evb + mrcl * rclb;
 110   }
 111 };
 112
 113 // Specialization for non arithmetic types. They will always use "a"
 114 template <>
 115 struct use_a<false>
 116 {
 117     template <typename Ct, typename Ev>
 118     static bool apply(Ct const& , Ct const& , Ev const& , Ev const& )
 119     {
 120         return true;
 121     }
 122 };
 123
 124 }
 125
 126 /*!
 127     \see http://mathworld.wolfram.com/Line-LineIntersection.html
 128  */
 129 template
 130 <
 131     typename CalculationType = void
 132 >
 133 struct cartesian_segments
 134 {
 135     typedef cartesian_tag cs_tag;
 136
 137     template <typename CoordinateType, typename SegmentRatio>
 138     struct segment_intersection_info
 139     {
 140     private :
 141         typedef typename select_most_precise
 142             <
 143                 CoordinateType, double
 144             >::type promoted_type;
 145
 146         promoted_type comparable_length_a() const
 147         {
 148             return dx_a * dx_a + dy_a * dy_a;
 149         }
 150
 151         promoted_type comparable_length_b() const
 152         {
 153             return dx_b * dx_b + dy_b * dy_b;
 154         }
 155
 156         template <typename Point, typename Segment1, typename Segment2>
 157         void assign_a(Point& point, Segment1 const& a, Segment2 const& ) const
 158         {
 159             assign(point, a, dx_a, dy_a, robust_ra);
 160         }
 161         template <typename Point, typename Segment1, typename Segment2>
 162         void assign_b(Point& point, Segment1 const& , Segment2 const& b) const
 163         {
 164             assign(point, b, dx_b, dy_b, robust_rb);
 165         }
 166
 167         template <typename Point, typename Segment>
 168         void assign(Point& point, Segment const& segment,
 169                     CoordinateType const& dx, CoordinateType const& dy,
 170                     SegmentRatio const& ratio) const
 171         {
 172             // Calculate the intersection point based on segment_ratio
 173             // The division, postponed until here, is done now. In case of integer this
 174             // results in an integer which rounds to the nearest integer.
 175             BOOST_GEOMETRY_ASSERT(ratio.denominator() != typename SegmentRatio::int_type(0));
 176
 177             typedef typename promote_integral<CoordinateType>::type calc_type;
 178
 179             calc_type const numerator
 180                 = boost::numeric_cast<calc_type>(ratio.numerator());
 181             calc_type const denominator
 182                 = boost::numeric_cast<calc_type>(ratio.denominator());
 183             calc_type const dx_calc = boost::numeric_cast<calc_type>(dx);
 184             calc_type const dy_calc = boost::numeric_cast<calc_type>(dy);
 185
 186             set<0>(point, get<0, 0>(segment)
 187                    + boost::numeric_cast<CoordinateType>(
 188                          math::divide<calc_type>(numerator * dx_calc, denominator)));
 189             set<1>(point, get<0, 1>(segment)
 190                    + boost::numeric_cast<CoordinateType>(
 191                          math::divide<calc_type>(numerator * dy_calc, denominator)));
 192         }
 193
 194         template <int Index, int Dim, typename Point, typename Segment>
 195         static bool exceeds_side_in_dimension(Point& p, Segment const& s)
 196         {
 197             // Situation a (positive)
 198             //     0>-------------->1     segment
 199             // *                          point left of segment<I> in D x or y
 200             // Situation b (negative)
 201             //     1<--------------<0     segment
 202             // *                          point right of segment<I>
 203             // Situation c (degenerate), return false (check other dimension)
 204             auto const& c = get<Dim>(p);
 205             auto const& c0 = get<Index, Dim>(s);
 206             auto const& c1 = get<1 - Index, Dim>(s);
 207             return c0 < c1 ? math::smaller(c, c0)
 208                  : c0 > c1 ? math::larger(c, c0)
 209                  : false;
 210         }
 211
 212         template <int Index, typename Point, typename Segment>
 213         static bool exceeds_side_of_segment(Point& p, Segment const& s)
 214         {
 215             return exceeds_side_in_dimension<Index, 0>(p, s)
 216                 || exceeds_side_in_dimension<Index, 1>(p, s);
 217         }
 218
 219         template <typename Point, typename Segment>
 220         static void assign_if_exceeds(Point& point, Segment const& s)
 221         {
 222             if (exceeds_side_of_segment<0>(point, s))
 223             {
 224                 detail::assign_point_from_index<0>(s, point);
 225             }
 226             else if (exceeds_side_of_segment<1>(point, s))
 227             {
 228                 detail::assign_point_from_index<1>(s, point);
 229             }
 230         }
 231
 232     public :
 233         template <typename Point, typename Segment1, typename Segment2>
 234         void calculate(Point& point, Segment1 const& a, Segment2 const& b) const
 235         {
 236             bool const use_a
 237                 = detail_usage::use_a
 238                      <
 239                          std::is_arithmetic<CoordinateType>::value
 240                      >::apply(comparable_length_a(), comparable_length_b(),
 241                          robust_ra.edge_value(), robust_rb.edge_value());
 242
 243             if (use_a)
 244             {
 245                 assign_a(point, a, b);
 246             }
 247             else
 248             {
 249                 assign_b(point, a, b);
 250             }
 251
 252 #if defined(BOOST_GEOMETRY_USE_RESCALING)
 253             return;
 254 #endif
 255
 256             // Verify nearly collinear cases (the threshold is arbitrary
 257             // but influences performance). If the intersection is located
 258             // outside the segments, then it should be moved.
 259             if (robust_ra.possibly_collinear(1.0e-3)
 260                 && robust_rb.possibly_collinear(1.0e-3))
 261             {
 262                 // The segments are nearly collinear and because of the calculation
 263                 // method with very small denominator, the IP appears outside the
 264                 // segment(s). Correct it to the end point.
 265                 // Because they are nearly collinear, it doesn't really matter to
 266                 // to which endpoint (or it is corrected twice).
 267                 assign_if_exceeds(point, a);
 268                 assign_if_exceeds(point, b);
 269             }
 270         }
 271
 272         CoordinateType dx_a, dy_a;
 273         CoordinateType dx_b, dy_b;
 274         SegmentRatio robust_ra;
 275         SegmentRatio robust_rb;
 276     };
 277
 278     template <typename D, typename W, typename ResultType>
 279     static inline void cramers_rule(D const& dx_a, D const& dy_a,
 280         D const& dx_b, D const& dy_b, W const& wx, W const& wy,
 281         // out:
 282         ResultType& nominator, ResultType& denominator)
 283     {
 284         // Cramers rule
 285         nominator = geometry::detail::determinant<ResultType>(dx_b, dy_b, wx, wy);
 286         denominator = geometry::detail::determinant<ResultType>(dx_a, dy_a, dx_b, dy_b);
 287         // Ratio r = nominator/denominator
 288         // Collinear if denominator == 0, intersecting if 0 <= r <= 1
 289         // IntersectionPoint = (x1 + r * dx_a, y1 + r * dy_a)
 290     }
 291
 292     // Version for non-rescaled policies
 293     template
 294     <
 295         typename UniqueSubRange1,
 296         typename UniqueSubRange2,
 297         typename Policy
 298     >
 299     static inline typename Policy::return_type
 300         apply(UniqueSubRange1 const& range_p,
 301               UniqueSubRange2 const& range_q,
 302               Policy const& policy)
 303     {
 304         // Pass the same ranges both as normal ranges and as modelled ranges
 305         return apply(range_p, range_q, policy, range_p, range_q);
 306     }
 307
 308     // Version for non rescaled versions.
 309     // The "modelled" parameter might be rescaled (will be removed later)
 310     template
 311     <
 312         typename UniqueSubRange1,
 313         typename UniqueSubRange2,
 314         typename Policy,
 315         typename ModelledUniqueSubRange1,
 316         typename ModelledUniqueSubRange2
 317     >
 318     static inline typename Policy::return_type
 319         apply(UniqueSubRange1 const& range_p,
 320               UniqueSubRange2 const& range_q,
 321               Policy const& policy,
 322               ModelledUniqueSubRange1 const& modelled_range_p,
 323               ModelledUniqueSubRange2 const& modelled_range_q)
 324     {
 325         typedef typename UniqueSubRange1::point_type point1_type;
 326         typedef typename UniqueSubRange2::point_type point2_type;
 327
 328         BOOST_CONCEPT_ASSERT( (concepts::ConstPoint<point1_type>) );
 329         BOOST_CONCEPT_ASSERT( (concepts::ConstPoint<point2_type>) );
 330
 331         point1_type const& p1 = range_p.at(0);
 332         point1_type const& p2 = range_p.at(1);
 333         point2_type const& q1 = range_q.at(0);
 334         point2_type const& q2 = range_q.at(1);
 335
 336         // Declare segments, currently necessary for the policies
 337         // (segment_crosses, segment_colinear, degenerate, one_degenerate, etc)
 338         model::referring_segment<point1_type const> const p(p1, p2);
 339         model::referring_segment<point2_type const> const q(q1, q2);
 340
 341         typedef typename select_most_precise
 342             <
 343                 typename geometry::coordinate_type<typename ModelledUniqueSubRange1::point_type>::type,
 344                 typename geometry::coordinate_type<typename ModelledUniqueSubRange1::point_type>::type
 345             >::type modelled_coordinate_type;
 346
 347         typedef segment_ratio<modelled_coordinate_type> ratio_type;
 348         segment_intersection_info
 349             <
 350                 typename select_calculation_type<point1_type, point2_type, CalculationType>::type,
 351                 ratio_type
 352             > sinfo;
 353
 354         sinfo.dx_a = get<0>(p2) - get<0>(p1); // distance in x-dir
 355         sinfo.dx_b = get<0>(q2) - get<0>(q1);
 356         sinfo.dy_a = get<1>(p2) - get<1>(p1); // distance in y-dir
 357         sinfo.dy_b = get<1>(q2) - get<1>(q1);
 358
 359         return unified<ratio_type>(sinfo, p, q, policy, modelled_range_p, modelled_range_q);
 360     }
 361
 362     //! Returns true if two segments do not overlap.
 363     //! If not, then no further calculations need to be done.
 364     template
 365     <
 366         std::size_t Dimension,
 367         typename PointP,
 368         typename PointQ
 369     >
 370     static inline bool disjoint_by_range(PointP const& p1, PointP const& p2,
 371                                          PointQ const& q1, PointQ const& q2)
 372     {
 373         auto minp = get<Dimension>(p1);
 374         auto maxp = get<Dimension>(p2);
 375         auto minq = get<Dimension>(q1);
 376         auto maxq = get<Dimension>(q2);
 377         if (minp > maxp)
 378         {
 379             std::swap(minp, maxp);
 380         }
 381         if (minq > maxq)
 382         {
 383             std::swap(minq, maxq);
 384         }
 385
 386         // In this case, max(p) < min(q)
 387         //     P         Q
 388         // <-------> <------->
 389         // (and the space in between is not extremely small)
 390         return math::smaller(maxp, minq) || math::smaller(maxq, minp);
 391     }
 392
 393     // Implementation for either rescaled or non rescaled versions.
 394     template
 395     <
 396         typename RatioType,
 397         typename SegmentInfo,
 398         typename Segment1,
 399         typename Segment2,
 400         typename Policy,
 401         typename UniqueSubRange1,
 402         typename UniqueSubRange2
 403     >
 404     static inline typename Policy::return_type
 405         unified(SegmentInfo& sinfo,
 406                 Segment1 const& p, Segment2 const& q, Policy const&,
 407                 UniqueSubRange1 const& range_p,
 408                 UniqueSubRange2 const& range_q)
 409     {
 410         typedef typename UniqueSubRange1::point_type point1_type;
 411         typedef typename UniqueSubRange2::point_type point2_type;
 412         typedef typename select_most_precise
 413             <
 414                 typename geometry::coordinate_type<point1_type>::type,
 415                 typename geometry::coordinate_type<point2_type>::type
 416             >::type coordinate_type;
 417
 418         point1_type const& p1 = range_p.at(0);
 419         point1_type const& p2 = range_p.at(1);
 420         point2_type const& q1 = range_q.at(0);
 421         point2_type const& q2 = range_q.at(1);
 422
 423         bool const p_is_point = equals_point_point(p1, p2);
 424         bool const q_is_point = equals_point_point(q1, q2);
 425
 426         if (p_is_point && q_is_point)
 427         {
 428             return equals_point_point(p1, q2)
 429                 ? Policy::degenerate(p, true)
 430                 : Policy::disjoint()
 431                 ;
 432         }
 433
 434         if (disjoint_by_range<0>(p1, p2, q1, q2)
 435          || disjoint_by_range<1>(p1, p2, q1, q2))
 436         {
 437             return Policy::disjoint();
 438         }
 439
 440         using side_strategy_type
 441             = typename side::services::default_strategy
 442                 <cartesian_tag, CalculationType>::type;
 443         side_info sides;
 444         sides.set<0>(side_strategy_type::apply(q1, q2, p1),
 445                      side_strategy_type::apply(q1, q2, p2));
 446
 447         if (sides.same<0>())
 448         {
 449             // Both points are at same side of other segment, we can leave
 450             return Policy::disjoint();
 451         }
 452
 453         sides.set<1>(side_strategy_type::apply(p1, p2, q1),
 454                      side_strategy_type::apply(p1, p2, q2));
 455
 456         if (sides.same<1>())
 457         {
 458             // Both points are at same side of other segment, we can leave
 459             return Policy::disjoint();
 460         }
 461
 462         bool collinear = sides.collinear();
 463
 464         // Calculate the differences again
 465         // (for rescaled version, this is different from dx_p etc)
 466         coordinate_type const dx_p = get<0>(p2) - get<0>(p1);
 467         coordinate_type const dx_q = get<0>(q2) - get<0>(q1);
 468         coordinate_type const dy_p = get<1>(p2) - get<1>(p1);
 469         coordinate_type const dy_q = get<1>(q2) - get<1>(q1);
 470
 471         // r: ratio 0-1 where intersection divides A/B
 472         // (only calculated for non-collinear segments)
 473         if (! collinear)
 474         {
 475             coordinate_type denominator_a, nominator_a;
 476             coordinate_type denominator_b, nominator_b;
 477
 478             cramers_rule(dx_p, dy_p, dx_q, dy_q,
 479                 get<0>(p1) - get<0>(q1),
 480                 get<1>(p1) - get<1>(q1),
 481                 nominator_a, denominator_a);
 482
 483             cramers_rule(dx_q, dy_q, dx_p, dy_p,
 484                 get<0>(q1) - get<0>(p1),
 485                 get<1>(q1) - get<1>(p1),
 486                 nominator_b, denominator_b);
 487
 488             math::detail::equals_factor_policy<coordinate_type>
 489                 policy(dx_p, dy_p, dx_q, dy_q);
 490
 491             coordinate_type const zero = 0;
 492             if (math::detail::equals_by_policy(denominator_a, zero, policy)
 493              || math::detail::equals_by_policy(denominator_b, zero, policy))
 494             {
 495                 // If this is the case, no rescaling is done for FP precision.
 496                 // We set it to collinear, but it indicates a robustness issue.
 497                 sides.set<0>(0, 0);
 498                 sides.set<1>(0, 0);
 499                 collinear = true;
 500             }
 501             else
 502             {
 503                 sinfo.robust_ra.assign(nominator_a, denominator_a);
 504                 sinfo.robust_rb.assign(nominator_b, denominator_b);
 505             }
 506         }
 507
 508         if (collinear)
 509         {
 510             std::pair<bool, bool> const collinear_use_first
 511                     = is_x_more_significant(geometry::math::abs(dx_p),
 512                                             geometry::math::abs(dy_p),
 513                                             geometry::math::abs(dx_q),
 514                                             geometry::math::abs(dy_q),
 515                                             p_is_point, q_is_point);
 516
 517             if (collinear_use_first.second)
 518             {
 519                 // Degenerate cases: segments of single point, lying on other segment, are not disjoint
 520                 // This situation is collinear too
 521
 522                 if (collinear_use_first.first)
 523                 {
 524                     return relate_collinear<0, Policy, RatioType>(p, q,
 525                             p1, p2, q1, q2,
 526                             p_is_point, q_is_point);
 527                 }
 528                 else
 529                 {
 530                     // Y direction contains larger segments (maybe dx is zero)
 531                     return relate_collinear<1, Policy, RatioType>(p, q,
 532                             p1, p2, q1, q2,
 533                             p_is_point, q_is_point);
 534                 }
 535             }
 536         }
 537
 538         return Policy::segments_crosses(sides, sinfo, p, q);
 539     }
 540
 541 private:
 542     // first is true if x is more significant
 543     // second is true if the more significant difference is not 0
 544     template <typename CoordinateType>
 545     static inline std::pair<bool, bool>
 546         is_x_more_significant(CoordinateType const& abs_dx_a,
 547                               CoordinateType const& abs_dy_a,
 548                               CoordinateType const& abs_dx_b,
 549                               CoordinateType const& abs_dy_b,
 550                               bool const a_is_point,
 551                               bool const b_is_point)
 552     {
 553         //BOOST_GEOMETRY_ASSERT_MSG(!(a_is_point && b_is_point), "both segments shouldn't be degenerated");
 554
 555         // for degenerated segments the second is always true because this function
 556         // shouldn't be called if both segments were degenerated
 557
 558         if (a_is_point)
 559         {
 560             return std::make_pair(abs_dx_b >= abs_dy_b, true);
 561         }
 562         else if (b_is_point)
 563         {
 564             return std::make_pair(abs_dx_a >= abs_dy_a, true);
 565         }
 566         else
 567         {
 568             CoordinateType const min_dx = (std::min)(abs_dx_a, abs_dx_b);
 569             CoordinateType const min_dy = (std::min)(abs_dy_a, abs_dy_b);
 570             return min_dx == min_dy ?
 571                     std::make_pair(true, min_dx > CoordinateType(0)) :
 572                     std::make_pair(min_dx > min_dy, true);
 573         }
 574     }
 575
 576     template
 577     <
 578         std::size_t Dimension,
 579         typename Policy,
 580         typename RatioType,
 581         typename Segment1,
 582         typename Segment2,
 583         typename RobustPoint1,
 584         typename RobustPoint2
 585     >
 586     static inline typename Policy::return_type
 587         relate_collinear(Segment1 const& a,
 588                          Segment2 const& b,
 589                          RobustPoint1 const& robust_a1, RobustPoint1 const& robust_a2,
 590                          RobustPoint2 const& robust_b1, RobustPoint2 const& robust_b2,
 591                          bool a_is_point, bool b_is_point)
 592     {
 593         if (a_is_point)
 594         {
 595             return relate_one_degenerate<Policy, RatioType>(a,
 596                 get<Dimension>(robust_a1),
 597                 get<Dimension>(robust_b1), get<Dimension>(robust_b2),
 598                 true);
 599         }
 600         if (b_is_point)
 601         {
 602             return relate_one_degenerate<Policy, RatioType>(b,
 603                 get<Dimension>(robust_b1),
 604                 get<Dimension>(robust_a1), get<Dimension>(robust_a2),
 605                 false);
 606         }
 607         return relate_collinear<Policy, RatioType>(a, b,
 608                                 get<Dimension>(robust_a1),
 609                                 get<Dimension>(robust_a2),
 610                                 get<Dimension>(robust_b1),
 611                                 get<Dimension>(robust_b2));
 612     }
 613
 614     /// Relate segments known collinear
 615     template
 616     <
 617         typename Policy,
 618         typename RatioType,
 619         typename Segment1,
 620         typename Segment2,
 621         typename Type1,
 622         typename Type2
 623     >
 624     static inline typename Policy::return_type
 625         relate_collinear(Segment1 const& a, Segment2 const& b,
 626                          Type1 oa_1, Type1 oa_2,
 627                          Type2 ob_1, Type2 ob_2)
 628     {
 629         // Calculate the ratios where a starts in b, b starts in a
 630         //         a1--------->a2         (2..7)
 631         //                b1----->b2      (5..8)
 632         // length_a: 7-2=5
 633         // length_b: 8-5=3
 634         // b1 is located w.r.t. a at ratio: (5-2)/5=3/5 (on a)
 635         // b2 is located w.r.t. a at ratio: (8-2)/5=6/5 (right of a)
 636         // a1 is located w.r.t. b at ratio: (2-5)/3=-3/3 (left of b)
 637         // a2 is located w.r.t. b at ratio: (7-5)/3=2/3 (on b)
 638         // A arrives (a2 on b), B departs (b1 on a)
 639
 640         // If both are reversed:
 641         //         a2<---------a1         (7..2)
 642         //                b2<-----b1      (8..5)
 643         // length_a: 2-7=-5
 644         // length_b: 5-8=-3
 645         // b1 is located w.r.t. a at ratio: (8-7)/-5=-1/5 (before a starts)
 646         // b2 is located w.r.t. a at ratio: (5-7)/-5=2/5 (on a)
 647         // a1 is located w.r.t. b at ratio: (7-8)/-3=1/3 (on b)
 648         // a2 is located w.r.t. b at ratio: (2-8)/-3=6/3 (after b ends)
 649
 650         // If both one is reversed:
 651         //         a1--------->a2         (2..7)
 652         //                b2<-----b1      (8..5)
 653         // length_a: 7-2=+5
 654         // length_b: 5-8=-3
 655         // b1 is located w.r.t. a at ratio: (8-2)/5=6/5 (after a ends)
 656         // b2 is located w.r.t. a at ratio: (5-2)/5=3/5 (on a)
 657         // a1 is located w.r.t. b at ratio: (2-8)/-3=6/3 (after b ends)
 658         // a2 is located w.r.t. b at ratio: (7-8)/-3=1/3 (on b)
 659         Type1 const length_a = oa_2 - oa_1; // no abs, see above
 660         Type2 const length_b = ob_2 - ob_1;
 661
 662         RatioType ra_from(oa_1 - ob_1, length_b);
 663         RatioType ra_to(oa_2 - ob_1, length_b);
 664         RatioType rb_from(ob_1 - oa_1, length_a);
 665         RatioType rb_to(ob_2 - oa_1, length_a);
 666
 667         // use absolute measure to detect endpoints intersection
 668         // NOTE: it'd be possible to calculate bx_wrt_a using ax_wrt_b values
 669         int const a1_wrt_b = position_value(oa_1, ob_1, ob_2);
 670         int const a2_wrt_b = position_value(oa_2, ob_1, ob_2);
 671         int const b1_wrt_a = position_value(ob_1, oa_1, oa_2);
 672         int const b2_wrt_a = position_value(ob_2, oa_1, oa_2);
 673
 674         // fix the ratios if necessary
 675         // CONSIDER: fixing ratios also in other cases, if they're inconsistent
 676         // e.g. if ratio == 1 or 0 (so IP at the endpoint)
 677         // but position value indicates that the IP is in the middle of the segment
 678         // because one of the segments is very long
 679         // In such case the ratios could be moved into the middle direction
 680         // by some small value (e.g. EPS+1ULP)
 681         if (a1_wrt_b == 1)
 682         {
 683             ra_from.assign(0, 1);
 684             rb_from.assign(0, 1);
 685         }
 686         else if (a1_wrt_b == 3)
 687         {
 688             ra_from.assign(1, 1);
 689             rb_to.assign(0, 1);
 690         }
 691
 692         if (a2_wrt_b == 1)
 693         {
 694             ra_to.assign(0, 1);
 695             rb_from.assign(1, 1);
 696         }
 697         else if (a2_wrt_b == 3)
 698         {
 699             ra_to.assign(1, 1);
 700             rb_to.assign(1, 1);
 701         }
 702
 703         if ((a1_wrt_b < 1 && a2_wrt_b < 1) || (a1_wrt_b > 3 && a2_wrt_b > 3))
 704         //if ((ra_from.left() && ra_to.left()) || (ra_from.right() && ra_to.right()))
 705         {
 706             return Policy::disjoint();
 707         }
 708
 709         bool const opposite = math::sign(length_a) != math::sign(length_b);
 710
 711         return Policy::segments_collinear(a, b, opposite,
 712                                           a1_wrt_b, a2_wrt_b, b1_wrt_a, b2_wrt_a,
 713                                           ra_from, ra_to, rb_from, rb_to);
 714     }
 715
 716     /// Relate segments where one is degenerate
 717     template
 718     <
 719         typename Policy,
 720         typename RatioType,
 721         typename DegenerateSegment,
 722         typename Type1,
 723         typename Type2
 724     >
 725     static inline typename Policy::return_type
 726         relate_one_degenerate(DegenerateSegment const& degenerate_segment,
 727                               Type1 d, Type2 s1, Type2 s2,
 728                               bool a_degenerate)
 729     {
 730         // Calculate the ratios where ds starts in s
 731         //         a1--------->a2         (2..6)
 732         //              b1/b2      (4..4)
 733         // Ratio: (4-2)/(6-2)
 734         RatioType const ratio(d - s1, s2 - s1);
 735
 736         if (!ratio.on_segment())
 737         {
 738             return Policy::disjoint();
 739         }
 740
 741         return Policy::one_degenerate(degenerate_segment, ratio, a_degenerate);
 742     }
 743
 744     template <typename ProjCoord1, typename ProjCoord2>
 745     static inline int position_value(ProjCoord1 const& ca1,
 746                                      ProjCoord2 const& cb1,
 747                                      ProjCoord2 const& cb2)
 748     {
 749         // S1x  0   1    2     3   4
 750         // S2       |---------->
 751         return math::equals(ca1, cb1) ? 1
 752              : math::equals(ca1, cb2) ? 3
 753              : cb1 < cb2 ?
 754                 ( ca1 < cb1 ? 0
 755                 : ca1 > cb2 ? 4
 756                 : 2 )
 757               : ( ca1 > cb1 ? 0
 758                 : ca1 < cb2 ? 4
 759                 : 2 );
 760     }
 761
 762     template <typename Point1, typename Point2>
 763     static inline bool equals_point_point(Point1 const& point1, Point2 const& point2)
 764     {
 765         return strategy::within::cartesian_point_point::apply(point1, point2);
 766     }
 767 };
 768
 769
 770 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
 771 namespace services
 772 {
 773
 774 template <typename CalculationType>
 775 struct default_strategy<cartesian_tag, CalculationType>
 776 {
 777     typedef cartesian_segments<CalculationType> type;
 778 };
 779
 780 } // namespace services
 781 #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
 782
 783
 784 }} // namespace strategy::intersection
 785
 786 namespace strategy
 787 {
 788
 789 namespace within { namespace services
 790 {
 791
 792 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
 793 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, linear_tag, cartesian_tag, cartesian_tag>
 794 {
 795     typedef strategy::intersection::cartesian_segments<> type;
 796 };
 797
 798 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
 799 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, polygonal_tag, cartesian_tag, cartesian_tag>
 800 {
 801     typedef strategy::intersection::cartesian_segments<> type;
 802 };
 803
 804 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
 805 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, linear_tag, cartesian_tag, cartesian_tag>
 806 {
 807     typedef strategy::intersection::cartesian_segments<> type;
 808 };
 809
 810 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
 811 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, polygonal_tag, cartesian_tag, cartesian_tag>
 812 {
 813     typedef strategy::intersection::cartesian_segments<> type;
 814 };
 815
 816 }} // within::services
 817
 818 namespace covered_by { namespace services
 819 {
 820
 821 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
 822 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, linear_tag, cartesian_tag, cartesian_tag>
 823 {
 824     typedef strategy::intersection::cartesian_segments<> type;
 825 };
 826
 827 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
 828 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, polygonal_tag, cartesian_tag, cartesian_tag>
 829 {
 830     typedef strategy::intersection::cartesian_segments<> type;
 831 };
 832
 833 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
 834 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, linear_tag, cartesian_tag, cartesian_tag>
 835 {
 836     typedef strategy::intersection::cartesian_segments<> type;
 837 };
 838
 839 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
 840 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, polygonal_tag, cartesian_tag, cartesian_tag>
 841 {
 842     typedef strategy::intersection::cartesian_segments<> type;
 843 };
 844
 845 }} // within::services
 846
 847 } // strategy
 848
 849 }} // namespace boost::geometry
 850
 851
 852 #endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_INTERSECTION_HPP