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, 2016, 2017, 2018, 2019.
   7 // Modifications copyright (c) 2014-2019, 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/cartesian/area.hpp>
  37 #include <boost/geometry/strategies/cartesian/disjoint_box_box.hpp>
  38 #include <boost/geometry/strategies/cartesian/disjoint_segment_box.hpp>
  39 #include <boost/geometry/strategies/cartesian/distance_pythagoras.hpp>
  40 #include <boost/geometry/strategies/cartesian/envelope.hpp>
  41 #include <boost/geometry/strategies/cartesian/expand_box.hpp>
  42 #include <boost/geometry/strategies/cartesian/expand_segment.hpp>
  43 #include <boost/geometry/strategies/cartesian/point_in_point.hpp>
  44 #include <boost/geometry/strategies/cartesian/point_in_poly_winding.hpp>
  45 #include <boost/geometry/strategies/cartesian/side_by_triangle.hpp>
  46 #include <boost/geometry/strategies/covered_by.hpp>
  47 #include <boost/geometry/strategies/intersection.hpp>
  48 #include <boost/geometry/strategies/intersection_result.hpp>
  49 #include <boost/geometry/strategies/side.hpp>
  50 #include <boost/geometry/strategies/side_info.hpp>
  51 #include <boost/geometry/strategies/within.hpp>
  52
  53 #include <boost/geometry/policies/robustness/rescale_policy_tags.hpp>
  54 #include <boost/geometry/policies/robustness/robust_point_type.hpp>
  55
  56 #if defined(BOOST_GEOMETRY_DEBUG_ROBUSTNESS)
  57 #  include <boost/geometry/io/wkt/write.hpp>
  58 #endif
  59
  60
  61 namespace boost { namespace geometry
  62 {
  63
  64
  65 namespace strategy { namespace intersection
  66 {
  67
  68
  69 /*!
  70     \see http://mathworld.wolfram.com/Line-LineIntersection.html
  71  */
  72 template
  73 <
  74     typename CalculationType = void
  75 >
  76 struct cartesian_segments
  77 {
  78     typedef cartesian_tag cs_tag;
  79
  80     typedef side::side_by_triangle<CalculationType> side_strategy_type;
  81
  82     static inline side_strategy_type get_side_strategy()
  83     {
  84         return side_strategy_type();
  85     }
  86
  87     template <typename Geometry1, typename Geometry2>
  88     struct point_in_geometry_strategy
  89     {
  90         typedef strategy::within::cartesian_winding
  91             <
  92                 typename point_type<Geometry1>::type,
  93                 typename point_type<Geometry2>::type,
  94                 CalculationType
  95             > type;
  96     };
  97
  98     template <typename Geometry1, typename Geometry2>
  99     static inline typename point_in_geometry_strategy<Geometry1, Geometry2>::type
 100         get_point_in_geometry_strategy()
 101     {
 102         typedef typename point_in_geometry_strategy
 103             <
 104                 Geometry1, Geometry2
 105             >::type strategy_type;
 106         return strategy_type();
 107     }
 108
 109     template <typename Geometry>
 110     struct area_strategy
 111     {
 112         typedef area::cartesian
 113             <
 114                 CalculationType
 115             > type;
 116     };
 117
 118     template <typename Geometry>
 119     static inline typename area_strategy<Geometry>::type get_area_strategy()
 120     {
 121         typedef typename area_strategy<Geometry>::type strategy_type;
 122         return strategy_type();
 123     }
 124
 125     template <typename Geometry>
 126     struct distance_strategy
 127     {
 128         typedef distance::pythagoras
 129             <
 130                 CalculationType
 131             > type;
 132     };
 133
 134     template <typename Geometry>
 135     static inline typename distance_strategy<Geometry>::type get_distance_strategy()
 136     {
 137         typedef typename distance_strategy<Geometry>::type strategy_type;
 138         return strategy_type();
 139     }
 140
 141     typedef envelope::cartesian<CalculationType> envelope_strategy_type;
 142
 143     static inline envelope_strategy_type get_envelope_strategy()
 144     {
 145         return envelope_strategy_type();
 146     }
 147
 148     typedef expand::cartesian_segment expand_strategy_type;
 149
 150     static inline expand_strategy_type get_expand_strategy()
 151     {
 152         return expand_strategy_type();
 153     }
 154
 155     typedef within::cartesian_point_point point_in_point_strategy_type;
 156
 157     static inline point_in_point_strategy_type get_point_in_point_strategy()
 158     {
 159         return point_in_point_strategy_type();
 160     }
 161
 162     typedef within::cartesian_point_point equals_point_point_strategy_type;
 163
 164     static inline equals_point_point_strategy_type get_equals_point_point_strategy()
 165     {
 166         return equals_point_point_strategy_type();
 167     }
 168
 169     typedef disjoint::cartesian_box_box disjoint_box_box_strategy_type;
 170
 171     static inline disjoint_box_box_strategy_type get_disjoint_box_box_strategy()
 172     {
 173         return disjoint_box_box_strategy_type();
 174     }
 175
 176     typedef disjoint::segment_box disjoint_segment_box_strategy_type;
 177
 178     static inline disjoint_segment_box_strategy_type get_disjoint_segment_box_strategy()
 179     {
 180         return disjoint_segment_box_strategy_type();
 181     }
 182
 183     typedef covered_by::cartesian_point_box disjoint_point_box_strategy_type;
 184     typedef covered_by::cartesian_point_box covered_by_point_box_strategy_type;
 185     typedef within::cartesian_point_box within_point_box_strategy_type;
 186     typedef envelope::cartesian_box envelope_box_strategy_type;
 187     typedef expand::cartesian_box expand_box_strategy_type;
 188
 189     template <typename CoordinateType, typename SegmentRatio>
 190     struct segment_intersection_info
 191     {
 192     private :
 193         typedef typename select_most_precise
 194             <
 195                 CoordinateType, double
 196             >::type promoted_type;
 197
 198         promoted_type comparable_length_a() const
 199         {
 200             return dx_a * dx_a + dy_a * dy_a;
 201         }
 202
 203         promoted_type comparable_length_b() const
 204         {
 205             return dx_b * dx_b + dy_b * dy_b;
 206         }
 207
 208         template <typename Point, typename Segment1, typename Segment2>
 209         void assign_a(Point& point, Segment1 const& a, Segment2 const& ) const
 210         {
 211             assign(point, a, dx_a, dy_a, robust_ra);
 212         }
 213         template <typename Point, typename Segment1, typename Segment2>
 214         void assign_b(Point& point, Segment1 const& , Segment2 const& b) const
 215         {
 216             assign(point, b, dx_b, dy_b, robust_rb);
 217         }
 218
 219         template <typename Point, typename Segment>
 220         void assign(Point& point, Segment const& segment, CoordinateType const& dx, CoordinateType const& dy, SegmentRatio const& ratio) const
 221         {
 222             // Calculate the intersection point based on segment_ratio
 223             // Up to now, division was postponed. Here we divide using numerator/
 224             // denominator. In case of integer this results in an integer
 225             // division.
 226             BOOST_GEOMETRY_ASSERT(ratio.denominator() != 0);
 227
 228             typedef typename promote_integral<CoordinateType>::type calc_type;
 229
 230             calc_type const numerator
 231                 = boost::numeric_cast<calc_type>(ratio.numerator());
 232             calc_type const denominator
 233                 = boost::numeric_cast<calc_type>(ratio.denominator());
 234             calc_type const dx_calc = boost::numeric_cast<calc_type>(dx);
 235             calc_type const dy_calc = boost::numeric_cast<calc_type>(dy);
 236
 237             set<0>(point, get<0, 0>(segment)
 238                    + boost::numeric_cast<CoordinateType>(numerator * dx_calc
 239                                                          / denominator));
 240             set<1>(point, get<0, 1>(segment)
 241                    + boost::numeric_cast<CoordinateType>(numerator * dy_calc
 242                                                          / denominator));
 243         }
 244
 245     public :
 246         template <typename Point, typename Segment1, typename Segment2>
 247         void calculate(Point& point, Segment1 const& a, Segment2 const& b) const
 248         {
 249             bool use_a = true;
 250
 251             // Prefer one segment if one is on or near an endpoint
 252             bool const a_near_end = robust_ra.near_end();
 253             bool const b_near_end = robust_rb.near_end();
 254             if (a_near_end && ! b_near_end)
 255             {
 256                 use_a = true;
 257             }
 258             else if (b_near_end && ! a_near_end)
 259             {
 260                 use_a = false;
 261             }
 262             else
 263             {
 264                 // Prefer shorter segment
 265                 promoted_type const len_a = comparable_length_a();
 266                 promoted_type const len_b = comparable_length_b();
 267                 if (len_b < len_a)
 268                 {
 269                     use_a = false;
 270                 }
 271                 // else use_a is true but was already assigned like that
 272             }
 273
 274             if (use_a)
 275             {
 276                 assign_a(point, a, b);
 277             }
 278             else
 279             {
 280                 assign_b(point, a, b);
 281             }
 282         }
 283
 284         CoordinateType dx_a, dy_a;
 285         CoordinateType dx_b, dy_b;
 286         SegmentRatio robust_ra;
 287         SegmentRatio robust_rb;
 288     };
 289
 290     template <typename D, typename W, typename ResultType>
 291     static inline void cramers_rule(D const& dx_a, D const& dy_a,
 292         D const& dx_b, D const& dy_b, W const& wx, W const& wy,
 293         // out:
 294         ResultType& d, ResultType& da)
 295     {
 296         // Cramers rule
 297         d = geometry::detail::determinant<ResultType>(dx_a, dy_a, dx_b, dy_b);
 298         da = geometry::detail::determinant<ResultType>(dx_b, dy_b, wx, wy);
 299         // Ratio is da/d , collinear if d == 0, intersecting if 0 <= r <= 1
 300         // IntersectionPoint = (x1 + r * dx_a, y1 + r * dy_a)
 301     }
 302
 303     // Version for non-rescaled policies
 304     template
 305     <
 306         typename UniqueSubRange1,
 307         typename UniqueSubRange2,
 308         typename Policy
 309     >
 310     static inline typename Policy::return_type
 311         apply(UniqueSubRange1 const& range_p,
 312               UniqueSubRange2 const& range_q,
 313               Policy const& policy)
 314     {
 315         // Pass the same ranges both as normal ranges and as robust ranges
 316         return apply(range_p, range_q, policy, range_p, range_q);
 317     }
 318
 319     // Main entry-routine, calculating intersections of segments p / q
 320     template
 321     <
 322         typename UniqueSubRange1,
 323         typename UniqueSubRange2,
 324         typename Policy,
 325         typename RobustUniqueSubRange1,
 326         typename RobustUniqueSubRange2
 327     >
 328     static inline typename Policy::return_type
 329         apply(UniqueSubRange1 const& range_p,
 330               UniqueSubRange2 const& range_q,
 331               Policy const&,
 332               RobustUniqueSubRange1 const& robust_range_p,
 333               RobustUniqueSubRange2 const& robust_range_q)
 334     {
 335         typedef typename UniqueSubRange1::point_type point1_type;
 336         typedef typename UniqueSubRange2::point_type point2_type;
 337
 338         BOOST_CONCEPT_ASSERT( (concepts::ConstPoint<point1_type>) );
 339         BOOST_CONCEPT_ASSERT( (concepts::ConstPoint<point2_type>) );
 340
 341         // Get robust points (to be omitted later)
 342         typedef typename RobustUniqueSubRange1::point_type robust_point1_type;
 343         typedef typename RobustUniqueSubRange2::point_type robust_point2_type;
 344
 345         point1_type const& p1 = range_p.at(0);
 346         point1_type const& p2 = range_p.at(1);
 347         point2_type const& q1 = range_q.at(0);
 348         point2_type const& q2 = range_q.at(1);
 349
 350         robust_point1_type const& robust_a1 = robust_range_p.at(0);
 351         robust_point1_type const& robust_a2 = robust_range_p.at(1);
 352         robust_point2_type const& robust_b1 = robust_range_q.at(0);
 353         robust_point2_type const& robust_b2 = robust_range_q.at(1);
 354
 355         using geometry::detail::equals::equals_point_point;
 356         bool const a_is_point = equals_point_point(robust_a1, robust_a2, point_in_point_strategy_type());
 357         bool const b_is_point = equals_point_point(robust_b1, robust_b2, point_in_point_strategy_type());
 358
 359         if(a_is_point && b_is_point)
 360         {
 361             // Take either a or b
 362             model::referring_segment<point1_type const> const d(p1, p2);
 363
 364             return equals_point_point(robust_a1, robust_b2, point_in_point_strategy_type())
 365                 ? Policy::degenerate(d, true)
 366                 : Policy::disjoint()
 367                 ;
 368         }
 369
 370         side_info sides;
 371         sides.set<0>(side_strategy_type::apply(robust_b1, robust_b2, robust_a1),
 372                      side_strategy_type::apply(robust_b1, robust_b2, robust_a2));
 373
 374         if (sides.same<0>())
 375         {
 376             // Both points are at same side of other segment, we can leave
 377             return Policy::disjoint();
 378         }
 379
 380         sides.set<1>(side_strategy_type::apply(robust_a1, robust_a2, robust_b1),
 381                      side_strategy_type::apply(robust_a1, robust_a2, robust_b2));
 382
 383         if (sides.same<1>())
 384         {
 385             // Both points are at same side of other segment, we can leave
 386             return Policy::disjoint();
 387         }
 388
 389         bool collinear = sides.collinear();
 390
 391         typedef typename select_most_precise
 392             <
 393                 typename geometry::coordinate_type<robust_point1_type>::type,
 394                 typename geometry::coordinate_type<robust_point2_type>::type
 395             >::type robust_coordinate_type;
 396
 397         typedef segment_ratio<robust_coordinate_type> ratio_type;
 398
 399         segment_intersection_info
 400             <
 401                 typename select_calculation_type<point1_type, point2_type, CalculationType>::type,
 402                 ratio_type
 403             > sinfo;
 404
 405         sinfo.dx_a = get<0>(p2) - get<0>(p1); // distance in x-dir
 406         sinfo.dx_b = get<0>(q2) - get<0>(q1);
 407         sinfo.dy_a = get<1>(p2) - get<1>(p1); // distance in y-dir
 408         sinfo.dy_b = get<1>(q2) - get<1>(q1);
 409
 410         robust_coordinate_type const robust_dx_a = get<0>(robust_a2) - get<0>(robust_a1);
 411         robust_coordinate_type const robust_dx_b = get<0>(robust_b2) - get<0>(robust_b1);
 412         robust_coordinate_type const robust_dy_a = get<1>(robust_a2) - get<1>(robust_a1);
 413         robust_coordinate_type const robust_dy_b = get<1>(robust_b2) - get<1>(robust_b1);
 414
 415         // r: ratio 0-1 where intersection divides A/B
 416         // (only calculated for non-collinear segments)
 417         if (! collinear)
 418         {
 419             robust_coordinate_type robust_da0, robust_da;
 420             robust_coordinate_type robust_db0, robust_db;
 421
 422             cramers_rule(robust_dx_a, robust_dy_a, robust_dx_b, robust_dy_b,
 423                 get<0>(robust_a1) - get<0>(robust_b1),
 424                 get<1>(robust_a1) - get<1>(robust_b1),
 425                 robust_da0, robust_da);
 426
 427             cramers_rule(robust_dx_b, robust_dy_b, robust_dx_a, robust_dy_a,
 428                 get<0>(robust_b1) - get<0>(robust_a1),
 429                 get<1>(robust_b1) - get<1>(robust_a1),
 430                 robust_db0, robust_db);
 431
 432             math::detail::equals_factor_policy<robust_coordinate_type>
 433                 policy(robust_dx_a, robust_dy_a, robust_dx_b, robust_dy_b);
 434             robust_coordinate_type const zero = 0;
 435             if (math::detail::equals_by_policy(robust_da0, zero, policy)
 436              || math::detail::equals_by_policy(robust_db0, zero, policy))
 437             {
 438                 // If this is the case, no rescaling is done for FP precision.
 439                 // We set it to collinear, but it indicates a robustness issue.
 440                 sides.set<0>(0,0);
 441                 sides.set<1>(0,0);
 442                 collinear = true;
 443             }
 444             else
 445             {
 446                 sinfo.robust_ra.assign(robust_da, robust_da0);
 447                 sinfo.robust_rb.assign(robust_db, robust_db0);
 448             }
 449         }
 450
 451         // Declare segments, currently necessary for the policies
 452         // (segment_crosses, segment_colinear, degenerate, one_degenerate, etc)
 453         model::referring_segment<point1_type const> const p(p1, p2);
 454         model::referring_segment<point2_type const> const q(q1, q2);
 455
 456         if (collinear)
 457         {
 458             std::pair<bool, bool> const collinear_use_first
 459                     = is_x_more_significant(geometry::math::abs(robust_dx_a),
 460                                             geometry::math::abs(robust_dy_a),
 461                                             geometry::math::abs(robust_dx_b),
 462                                             geometry::math::abs(robust_dy_b),
 463                                             a_is_point, b_is_point);
 464
 465             if (collinear_use_first.second)
 466             {
 467                 // Degenerate cases: segments of single point, lying on other segment, are not disjoint
 468                 // This situation is collinear too
 469
 470                 if (collinear_use_first.first)
 471                 {
 472                     return relate_collinear<0, Policy, ratio_type>(p, q,
 473                             robust_a1, robust_a2, robust_b1, robust_b2,
 474                             a_is_point, b_is_point);
 475                 }
 476                 else
 477                 {
 478                     // Y direction contains larger segments (maybe dx is zero)
 479                     return relate_collinear<1, Policy, ratio_type>(p, q,
 480                             robust_a1, robust_a2, robust_b1, robust_b2,
 481                             a_is_point, b_is_point);
 482                 }
 483             }
 484         }
 485
 486         return Policy::segments_crosses(sides, sinfo, p, q);
 487     }
 488
 489 private:
 490     // first is true if x is more significant
 491     // second is true if the more significant difference is not 0
 492     template <typename RobustCoordinateType>
 493     static inline std::pair<bool, bool>
 494         is_x_more_significant(RobustCoordinateType const& abs_robust_dx_a,
 495                               RobustCoordinateType const& abs_robust_dy_a,
 496                               RobustCoordinateType const& abs_robust_dx_b,
 497                               RobustCoordinateType const& abs_robust_dy_b,
 498                               bool const a_is_point,
 499                               bool const b_is_point)
 500     {
 501         //BOOST_GEOMETRY_ASSERT_MSG(!(a_is_point && b_is_point), "both segments shouldn't be degenerated");
 502
 503         // for degenerated segments the second is always true because this function
 504         // shouldn't be called if both segments were degenerated
 505
 506         if (a_is_point)
 507         {
 508             return std::make_pair(abs_robust_dx_b >= abs_robust_dy_b, true);
 509         }
 510         else if (b_is_point)
 511         {
 512             return std::make_pair(abs_robust_dx_a >= abs_robust_dy_a, true);
 513         }
 514         else
 515         {
 516             RobustCoordinateType const min_dx = (std::min)(abs_robust_dx_a, abs_robust_dx_b);
 517             RobustCoordinateType const min_dy = (std::min)(abs_robust_dy_a, abs_robust_dy_b);
 518             return min_dx == min_dy ?
 519                     std::make_pair(true, min_dx > RobustCoordinateType(0)) :
 520                     std::make_pair(min_dx > min_dy, true);
 521         }
 522     }
 523
 524     template
 525     <
 526         std::size_t Dimension,
 527         typename Policy,
 528         typename RatioType,
 529         typename Segment1,
 530         typename Segment2,
 531         typename RobustPoint1,
 532         typename RobustPoint2
 533     >
 534     static inline typename Policy::return_type
 535         relate_collinear(Segment1 const& a,
 536                          Segment2 const& b,
 537                          RobustPoint1 const& robust_a1, RobustPoint1 const& robust_a2,
 538                          RobustPoint2 const& robust_b1, RobustPoint2 const& robust_b2,
 539                          bool a_is_point, bool b_is_point)
 540     {
 541         if (a_is_point)
 542         {
 543             return relate_one_degenerate<Policy, RatioType>(a,
 544                 get<Dimension>(robust_a1),
 545                 get<Dimension>(robust_b1), get<Dimension>(robust_b2),
 546                 true);
 547         }
 548         if (b_is_point)
 549         {
 550             return relate_one_degenerate<Policy, RatioType>(b,
 551                 get<Dimension>(robust_b1),
 552                 get<Dimension>(robust_a1), get<Dimension>(robust_a2),
 553                 false);
 554         }
 555         return relate_collinear<Policy, RatioType>(a, b,
 556                                 get<Dimension>(robust_a1),
 557                                 get<Dimension>(robust_a2),
 558                                 get<Dimension>(robust_b1),
 559                                 get<Dimension>(robust_b2));
 560     }
 561
 562     /// Relate segments known collinear
 563     template
 564     <
 565         typename Policy,
 566         typename RatioType,
 567         typename Segment1,
 568         typename Segment2,
 569         typename RobustType1,
 570         typename RobustType2
 571     >
 572     static inline typename Policy::return_type
 573         relate_collinear(Segment1 const& a, Segment2 const& b,
 574                          RobustType1 oa_1, RobustType1 oa_2,
 575                          RobustType2 ob_1, RobustType2 ob_2)
 576     {
 577         // Calculate the ratios where a starts in b, b starts in a
 578         //         a1--------->a2         (2..7)
 579         //                b1----->b2      (5..8)
 580         // length_a: 7-2=5
 581         // length_b: 8-5=3
 582         // b1 is located w.r.t. a at ratio: (5-2)/5=3/5 (on a)
 583         // b2 is located w.r.t. a at ratio: (8-2)/5=6/5 (right of a)
 584         // a1 is located w.r.t. b at ratio: (2-5)/3=-3/3 (left of b)
 585         // a2 is located w.r.t. b at ratio: (7-5)/3=2/3 (on b)
 586         // A arrives (a2 on b), B departs (b1 on a)
 587
 588         // If both are reversed:
 589         //         a2<---------a1         (7..2)
 590         //                b2<-----b1      (8..5)
 591         // length_a: 2-7=-5
 592         // length_b: 5-8=-3
 593         // b1 is located w.r.t. a at ratio: (8-7)/-5=-1/5 (before a starts)
 594         // b2 is located w.r.t. a at ratio: (5-7)/-5=2/5 (on a)
 595         // a1 is located w.r.t. b at ratio: (7-8)/-3=1/3 (on b)
 596         // a2 is located w.r.t. b at ratio: (2-8)/-3=6/3 (after b ends)
 597
 598         // If both one is reversed:
 599         //         a1--------->a2         (2..7)
 600         //                b2<-----b1      (8..5)
 601         // length_a: 7-2=+5
 602         // length_b: 5-8=-3
 603         // b1 is located w.r.t. a at ratio: (8-2)/5=6/5 (after a ends)
 604         // b2 is located w.r.t. a at ratio: (5-2)/5=3/5 (on a)
 605         // a1 is located w.r.t. b at ratio: (2-8)/-3=6/3 (after b ends)
 606         // a2 is located w.r.t. b at ratio: (7-8)/-3=1/3 (on b)
 607         RobustType1 const length_a = oa_2 - oa_1; // no abs, see above
 608         RobustType2 const length_b = ob_2 - ob_1;
 609
 610         RatioType ra_from(oa_1 - ob_1, length_b);
 611         RatioType ra_to(oa_2 - ob_1, length_b);
 612         RatioType rb_from(ob_1 - oa_1, length_a);
 613         RatioType rb_to(ob_2 - oa_1, length_a);
 614
 615         // use absolute measure to detect endpoints intersection
 616         // NOTE: it'd be possible to calculate bx_wrt_a using ax_wrt_b values
 617         int const a1_wrt_b = position_value(oa_1, ob_1, ob_2);
 618         int const a2_wrt_b = position_value(oa_2, ob_1, ob_2);
 619         int const b1_wrt_a = position_value(ob_1, oa_1, oa_2);
 620         int const b2_wrt_a = position_value(ob_2, oa_1, oa_2);
 621
 622         // fix the ratios if necessary
 623         // CONSIDER: fixing ratios also in other cases, if they're inconsistent
 624         // e.g. if ratio == 1 or 0 (so IP at the endpoint)
 625         // but position value indicates that the IP is in the middle of the segment
 626         // because one of the segments is very long
 627         // In such case the ratios could be moved into the middle direction
 628         // by some small value (e.g. EPS+1ULP)
 629         if (a1_wrt_b == 1)
 630         {
 631             ra_from.assign(0, 1);
 632             rb_from.assign(0, 1);
 633         }
 634         else if (a1_wrt_b == 3)
 635         {
 636             ra_from.assign(1, 1);
 637             rb_to.assign(0, 1);
 638         }
 639
 640         if (a2_wrt_b == 1)
 641         {
 642             ra_to.assign(0, 1);
 643             rb_from.assign(1, 1);
 644         }
 645         else if (a2_wrt_b == 3)
 646         {
 647             ra_to.assign(1, 1);
 648             rb_to.assign(1, 1);
 649         }
 650
 651         if ((a1_wrt_b < 1 && a2_wrt_b < 1) || (a1_wrt_b > 3 && a2_wrt_b > 3))
 652         //if ((ra_from.left() && ra_to.left()) || (ra_from.right() && ra_to.right()))
 653         {
 654             return Policy::disjoint();
 655         }
 656
 657         bool const opposite = math::sign(length_a) != math::sign(length_b);
 658
 659         return Policy::segments_collinear(a, b, opposite,
 660                                           a1_wrt_b, a2_wrt_b, b1_wrt_a, b2_wrt_a,
 661                                           ra_from, ra_to, rb_from, rb_to);
 662     }
 663
 664     /// Relate segments where one is degenerate
 665     template
 666     <
 667         typename Policy,
 668         typename RatioType,
 669         typename DegenerateSegment,
 670         typename RobustType1,
 671         typename RobustType2
 672     >
 673     static inline typename Policy::return_type
 674         relate_one_degenerate(DegenerateSegment const& degenerate_segment,
 675                               RobustType1 d, RobustType2 s1, RobustType2 s2,
 676                               bool a_degenerate)
 677     {
 678         // Calculate the ratios where ds starts in s
 679         //         a1--------->a2         (2..6)
 680         //              b1/b2      (4..4)
 681         // Ratio: (4-2)/(6-2)
 682         RatioType const ratio(d - s1, s2 - s1);
 683
 684         if (!ratio.on_segment())
 685         {
 686             return Policy::disjoint();
 687         }
 688
 689         return Policy::one_degenerate(degenerate_segment, ratio, a_degenerate);
 690     }
 691
 692     template <typename ProjCoord1, typename ProjCoord2>
 693     static inline int position_value(ProjCoord1 const& ca1,
 694                                      ProjCoord2 const& cb1,
 695                                      ProjCoord2 const& cb2)
 696     {
 697         // S1x  0   1    2     3   4
 698         // S2       |---------->
 699         return math::equals(ca1, cb1) ? 1
 700              : math::equals(ca1, cb2) ? 3
 701              : cb1 < cb2 ?
 702                 ( ca1 < cb1 ? 0
 703                 : ca1 > cb2 ? 4
 704                 : 2 )
 705               : ( ca1 > cb1 ? 0
 706                 : ca1 < cb2 ? 4
 707                 : 2 );
 708     }
 709 };
 710
 711
 712 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
 713 namespace services
 714 {
 715
 716 template <typename CalculationType>
 717 struct default_strategy<cartesian_tag, CalculationType>
 718 {
 719     typedef cartesian_segments<CalculationType> type;
 720 };
 721
 722 } // namespace services
 723 #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
 724
 725
 726 }} // namespace strategy::intersection
 727
 728 namespace strategy
 729 {
 730
 731 namespace within { namespace services
 732 {
 733
 734 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
 735 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, linear_tag, cartesian_tag, cartesian_tag>
 736 {
 737     typedef strategy::intersection::cartesian_segments<> type;
 738 };
 739
 740 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
 741 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, polygonal_tag, cartesian_tag, cartesian_tag>
 742 {
 743     typedef strategy::intersection::cartesian_segments<> type;
 744 };
 745
 746 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
 747 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, linear_tag, cartesian_tag, cartesian_tag>
 748 {
 749     typedef strategy::intersection::cartesian_segments<> type;
 750 };
 751
 752 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
 753 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, polygonal_tag, cartesian_tag, cartesian_tag>
 754 {
 755     typedef strategy::intersection::cartesian_segments<> type;
 756 };
 757
 758 }} // within::services
 759
 760 namespace covered_by { namespace services
 761 {
 762
 763 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
 764 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, linear_tag, cartesian_tag, cartesian_tag>
 765 {
 766     typedef strategy::intersection::cartesian_segments<> type;
 767 };
 768
 769 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
 770 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, polygonal_tag, cartesian_tag, cartesian_tag>
 771 {
 772     typedef strategy::intersection::cartesian_segments<> type;
 773 };
 774
 775 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
 776 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, linear_tag, cartesian_tag, cartesian_tag>
 777 {
 778     typedef strategy::intersection::cartesian_segments<> type;
 779 };
 780
 781 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
 782 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, polygonal_tag, cartesian_tag, cartesian_tag>
 783 {
 784     typedef strategy::intersection::cartesian_segments<> type;
 785 };
 786
 787 }} // within::services
 788
 789 } // strategy
 790
 791 }} // namespace boost::geometry
 792
 793
 794 #endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_INTERSECTION_HPP