strategies/cartesian/intersection.hpp

0001 // Boost.Geometry (aka GGL, Generic Geometry Library)
0002
0003 // Copyright (c) 2007-2014 Barend Gehrels, Amsterdam, the Netherlands.
0004 // Copyright (c) 2013-2023 Adam Wulkiewicz, Lodz, Poland.
0005
0006 // This file was modified by Oracle on 2014-2024.
0007 // Modifications copyright (c) 2014-2024, Oracle and/or its affiliates.
0008 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
0009 // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
0010 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0011
0012 // Use, modification and distribution is subject to the Boost Software License,
0013 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0014 // http://www.boost.org/LICENSE_1_0.txt)
0015
0016 #ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_INTERSECTION_HPP
0017 #define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_INTERSECTION_HPP
0018
0019 #include <algorithm>
0020
0021 #include <boost/geometry/core/exception.hpp>
0022
0023 #include <boost/geometry/geometries/concepts/point_concept.hpp>
0024 #include <boost/geometry/geometries/concepts/segment_concept.hpp>
0025 #include <boost/geometry/geometries/segment.hpp>
0026
0027 #include <boost/geometry/arithmetic/determinant.hpp>
0028 #include <boost/geometry/algorithms/detail/assign_values.hpp>
0029 #include <boost/geometry/algorithms/detail/assign_indexed_point.hpp>
0030 #include <boost/geometry/algorithms/detail/equals/point_point.hpp>
0031
0032 #include <boost/geometry/util/math.hpp>
0033 #include <boost/geometry/util/numeric_cast.hpp>
0034 #include <boost/geometry/util/promote_integral.hpp>
0035 #include <boost/geometry/util/select_calculation_type.hpp>
0036
0037 #include <boost/geometry/strategy/cartesian/area.hpp>
0038 #include <boost/geometry/strategy/cartesian/envelope.hpp>
0039 #include <boost/geometry/strategy/cartesian/expand_box.hpp>
0040 #include <boost/geometry/strategy/cartesian/expand_segment.hpp>
0041
0042 #include <boost/geometry/strategies/cartesian/disjoint_box_box.hpp>
0043 #include <boost/geometry/strategies/cartesian/disjoint_segment_box.hpp>
0044 #include <boost/geometry/strategies/cartesian/distance_pythagoras.hpp>
0045 #include <boost/geometry/strategies/cartesian/point_in_point.hpp>
0046 #include <boost/geometry/strategies/cartesian/point_in_poly_winding.hpp>
0047 #include <boost/geometry/strategies/covered_by.hpp>
0048 #include <boost/geometry/strategies/intersection.hpp>
0049 #include <boost/geometry/strategies/intersection_result.hpp>
0050 #include <boost/geometry/strategies/side.hpp>
0051 #include <boost/geometry/strategies/side_info.hpp>
0052 #include <boost/geometry/strategies/within.hpp>
0053
0054 namespace boost { namespace geometry
0055 {
0056
0057
0058 namespace strategy { namespace intersection
0059 {
0060
0061 namespace detail_usage
0062 {
0063
0064 // When calculating the intersection, the information of "a" or "b" can be used.
0065 // Theoretically this gives equal results, but due to floating point precision
0066 // there might be tiny differences. These are edge cases.
0067 // This structure is to determine if "a" or "b" should be used.
0068 // Prefer the segment closer to the endpoint.
0069 // If both are about equally close, then prefer the longer segment
0070 // To avoid hard thresholds, behavior is made fluent.
0071 // Calculate comparable length indications,
0072 // the longer the segment (relatively), the lower the value
0073 // such that the shorter lengths are evaluated higher and will
0074 // be preferred.
0075 template <bool IsArithmetic>
0076 struct use_a
0077 {
0078   template <typename Ct, typename Ev>
0079   static bool apply(Ct const& cla, Ct const& clb, Ev const& eva, Ev const& evb)
0080   {
0081       auto const clm = (std::max)(cla, clb);
0082       if (clm <= 0)
0083       {
0084           return true;
0085       }
0086
0087       // Relative comparible length
0088       auto const rcla = Ct(1.0) - cla / clm;
0089       auto const rclb = Ct(1.0) - clb / clm;
0090
0091       // Multipliers for edgevalue (ev) and relative comparible length (rcl)
0092       // They determine the balance between edge value (should be larger)
0093       // and segment length. In 99.9xx% of the cases there is no difference
0094       // at all (if either a or b is used). Therefore the values of the
0095       // constants are not sensitive for the majority of the situations.
0096       // One known case is #mysql_23023665_6 (difference) which needs mev >= 2
0097       Ev const mev = 5;
0098       Ev const mrcl = 1;
0099
0100       return mev * eva + mrcl * rcla > mev * evb + mrcl * rclb;
0101   }
0102 };
0103
0104 // Specialization for non arithmetic types. They will always use "a"
0105 template <>
0106 struct use_a<false>
0107 {
0108     template <typename Ct, typename Ev>
0109     static bool apply(Ct const& , Ct const& , Ev const& , Ev const& )
0110     {
0111         return true;
0112     }
0113 };
0114
0115 }
0116
0117 /*!
0118     \see http://mathworld.wolfram.com/Line-LineIntersection.html
0119  */
0120 template
0121 <
0122     typename CalculationType = void
0123 >
0124 struct cartesian_segments
0125 {
0126     using cs_tag = cartesian_tag;
0127
0128     template <typename CoordinateType, typename SegmentRatio>
0129     struct segment_intersection_info
0130     {
0131     private :
0132         typedef typename select_most_precise
0133             <
0134                 CoordinateType, double
0135             >::type promoted_type;
0136
0137         promoted_type comparable_length_a() const
0138         {
0139             return dx_a * dx_a + dy_a * dy_a;
0140         }
0141
0142         promoted_type comparable_length_b() const
0143         {
0144             return dx_b * dx_b + dy_b * dy_b;
0145         }
0146
0147         template <typename Point, typename Segment1, typename Segment2>
0148         void assign_a(Point& point, Segment1 const& a, Segment2 const& ) const
0149         {
0150             assign(point, a, dx_a, dy_a, ra);
0151         }
0152         template <typename Point, typename Segment1, typename Segment2>
0153         void assign_b(Point& point, Segment1 const& , Segment2 const& b) const
0154         {
0155             assign(point, b, dx_b, dy_b, rb);
0156         }
0157
0158         template <typename Point, typename Segment>
0159         void assign(Point& point, Segment const& segment,
0160                     CoordinateType const& dx, CoordinateType const& dy,
0161                     SegmentRatio const& ratio) const
0162         {
0163             // Calculate the intersection point based on segment_ratio
0164             // The division, postponed until here, is done now. In case of integer this
0165             // results in an integer which rounds to the nearest integer.
0166             BOOST_GEOMETRY_ASSERT(ratio.denominator() != typename SegmentRatio::int_type(0));
0167
0168             typedef typename promote_integral<CoordinateType>::type calc_type;
0169
0170             calc_type const numerator
0171                 = util::numeric_cast<calc_type>(ratio.numerator());
0172             calc_type const denominator
0173                 = util::numeric_cast<calc_type>(ratio.denominator());
0174             calc_type const dx_calc = util::numeric_cast<calc_type>(dx);
0175             calc_type const dy_calc = util::numeric_cast<calc_type>(dy);
0176
0177             set<0>(point, get<0, 0>(segment)
0178                    + util::numeric_cast<CoordinateType>(
0179                          math::divide<calc_type>(numerator * dx_calc, denominator)));
0180             set<1>(point, get<0, 1>(segment)
0181                    + util::numeric_cast<CoordinateType>(
0182                          math::divide<calc_type>(numerator * dy_calc, denominator)));
0183         }
0184
0185         template <int Index, int Dim, typename Point, typename Segment>
0186         static bool exceeds_side_in_dimension(Point& p, Segment const& s)
0187         {
0188             // Situation a (positive)
0189             //     0>-------------->1     segment
0190             // *                          point left of segment<I> in D x or y
0191             // Situation b (negative)
0192             //     1<--------------<0     segment
0193             // *                          point right of segment<I>
0194             // Situation c (degenerate), return false (check other dimension)
0195             auto const& c = get<Dim>(p);
0196             auto const& c0 = get<Index, Dim>(s);
0197             auto const& c1 = get<1 - Index, Dim>(s);
0198             return c0 < c1 ? math::smaller(c, c0)
0199                  : c0 > c1 ? math::larger(c, c0)
0200                  : false;
0201         }
0202
0203         template <int Index, typename Point, typename Segment>
0204         static bool exceeds_side_of_segment(Point& p, Segment const& s)
0205         {
0206             return exceeds_side_in_dimension<Index, 0>(p, s)
0207                 || exceeds_side_in_dimension<Index, 1>(p, s);
0208         }
0209
0210         template <typename Point, typename Segment>
0211         static void assign_if_exceeds(Point& point, Segment const& s)
0212         {
0213             if (exceeds_side_of_segment<0>(point, s))
0214             {
0215                 detail::assign_point_from_index<0>(s, point);
0216             }
0217             else if (exceeds_side_of_segment<1>(point, s))
0218             {
0219                 detail::assign_point_from_index<1>(s, point);
0220             }
0221         }
0222
0223     public :
0224         template <typename Point, typename Segment1, typename Segment2>
0225         void calculate(Point& point, Segment1 const& a, Segment2 const& b) const
0226         {
0227             bool const use_a
0228                 = detail_usage::use_a
0229                      <
0230                          std::is_arithmetic<CoordinateType>::value
0231                      >::apply(comparable_length_a(), comparable_length_b(),
0232                          ra.edge_value(), rb.edge_value());
0233
0234             if (use_a)
0235             {
0236                 assign_a(point, a, b);
0237             }
0238             else
0239             {
0240                 assign_b(point, a, b);
0241             }
0242             // Verify nearly collinear cases (the threshold is arbitrary
0243             // but influences performance). If the intersection is located
0244             // outside the segments, then it should be moved.
0245             if (ra.possibly_collinear(1.0e-3) && rb.possibly_collinear(1.0e-3))
0246             {
0247                 // The segments are nearly collinear and because of the calculation
0248                 // method with very small denominator, the IP appears outside the
0249                 // segment(s). Correct it to the end point.
0250                 // Because they are nearly collinear, it doesn't really matter to
0251                 // to which endpoint (or it is corrected twice).
0252                 assign_if_exceeds(point, a);
0253                 assign_if_exceeds(point, b);
0254             }
0255         }
0256
0257         CoordinateType dx_a, dy_a;
0258         CoordinateType dx_b, dy_b;
0259         SegmentRatio ra;
0260         SegmentRatio rb;
0261     };
0262
0263     template <typename D, typename W, typename ResultType>
0264     static inline void cramers_rule(D const& dx_a, D const& dy_a,
0265         D const& dx_b, D const& dy_b, W const& wx, W const& wy,
0266         // out:
0267         ResultType& nominator, ResultType& denominator)
0268     {
0269         // Cramers rule
0270         nominator = geometry::detail::determinant<ResultType>(dx_b, dy_b, wx, wy);
0271         denominator = geometry::detail::determinant<ResultType>(dx_a, dy_a, dx_b, dy_b);
0272         // Ratio r = nominator/denominator
0273         // Collinear if denominator == 0, intersecting if 0 <= r <= 1
0274         // IntersectionPoint = (x1 + r * dx_a, y1 + r * dy_a)
0275     }
0276
0277     template
0278     <
0279         typename UniqueSubRange1,
0280         typename UniqueSubRange2,
0281         typename Policy
0282     >
0283     static inline typename Policy::return_type apply(UniqueSubRange1 const& range_p,
0284                                                      UniqueSubRange2 const& range_q,
0285                                                      Policy const& policy)
0286     {
0287         typedef typename UniqueSubRange1::point_type point1_type;
0288         typedef typename UniqueSubRange2::point_type point2_type;
0289
0290         BOOST_CONCEPT_ASSERT( (concepts::ConstPoint<point1_type>) );
0291         BOOST_CONCEPT_ASSERT( (concepts::ConstPoint<point2_type>) );
0292
0293         point1_type const& p1 = range_p.at(0);
0294         point1_type const& p2 = range_p.at(1);
0295         point2_type const& q1 = range_q.at(0);
0296         point2_type const& q2 = range_q.at(1);
0297
0298         using coordinate_type = typename select_most_precise
0299             <
0300                 geometry::coordinate_type_t<point1_type>,
0301                 geometry::coordinate_type_t<point2_type>
0302             >::type;
0303
0304         using ratio_type = segment_ratio<coordinate_type>;
0305         segment_intersection_info
0306             <
0307                 typename select_calculation_type<point1_type, point2_type, CalculationType>::type,
0308                 ratio_type
0309             > sinfo;
0310
0311         sinfo.dx_a = get<0>(p2) - get<0>(p1); // distance in x-dir
0312         sinfo.dx_b = get<0>(q2) - get<0>(q1);
0313         sinfo.dy_a = get<1>(p2) - get<1>(p1); // distance in y-dir
0314         sinfo.dy_b = get<1>(q2) - get<1>(q1);
0315
0316         return unified<ratio_type>(sinfo, range_p, range_q, policy);
0317     }
0318
0319     //! Returns true if two segments do not overlap.
0320     //! If not, then no further calculations need to be done.
0321     template
0322     <
0323         std::size_t Dimension,
0324         typename PointP,
0325         typename PointQ
0326     >
0327     static inline bool disjoint_by_range(PointP const& p1, PointP const& p2,
0328                                          PointQ const& q1, PointQ const& q2)
0329     {
0330         auto minp = get<Dimension>(p1);
0331         auto maxp = get<Dimension>(p2);
0332         auto minq = get<Dimension>(q1);
0333         auto maxq = get<Dimension>(q2);
0334         if (minp > maxp)
0335         {
0336             std::swap(minp, maxp);
0337         }
0338         if (minq > maxq)
0339         {
0340             std::swap(minq, maxq);
0341         }
0342
0343         // In this case, max(p) < min(q)
0344         //     P         Q
0345         // <-------> <------->
0346         // (and the space in between is not extremely small)
0347         return math::smaller(maxp, minq) || math::smaller(maxq, minp);
0348     }
0349
0350     template
0351     <
0352         typename RatioType,
0353         typename SegmentInfo,
0354         typename UniqueSubRange1,
0355         typename UniqueSubRange2,
0356         typename Policy
0357     >
0358     static inline typename Policy::return_type
0359         unified(SegmentInfo& sinfo,
0360                 UniqueSubRange1 const& range_p,
0361                 UniqueSubRange2 const& range_q,
0362                 Policy const&)
0363     {
0364         using point1_type = typename UniqueSubRange1::point_type;
0365         using point2_type = typename UniqueSubRange2::point_type;
0366
0367         point1_type const& p1 = range_p.at(0);
0368         point1_type const& p2 = range_p.at(1);
0369         point2_type const& q1 = range_q.at(0);
0370         point2_type const& q2 = range_q.at(1);
0371
0372         bool const p_is_point = equals_point_point(p1, p2);
0373         bool const q_is_point = equals_point_point(q1, q2);
0374
0375         // Declare segments, currently necessary for the policies
0376         // (segment_crosses, segment_colinear, degenerate, one_degenerate, etc)
0377         model::referring_segment<point1_type const> const p(p1, p2);
0378         model::referring_segment<point2_type const> const q(q1, q2);
0379
0380         if (p_is_point && q_is_point)
0381         {
0382             return equals_point_point(p1, q2)
0383                 ? Policy::degenerate(p, true)
0384                 : Policy::disjoint()
0385                 ;
0386         }
0387
0388         if (disjoint_by_range<0>(p1, p2, q1, q2)
0389          || disjoint_by_range<1>(p1, p2, q1, q2))
0390         {
0391             return Policy::disjoint();
0392         }
0393
0394         using side_strategy_type
0395             = typename side::services::default_strategy
0396                 <cartesian_tag, CalculationType>::type;
0397         side_info sides;
0398         sides.set<0>(side_strategy_type::apply(q1, q2, p1),
0399                      side_strategy_type::apply(q1, q2, p2));
0400
0401         if (sides.same<0>())
0402         {
0403             // Both points are at same side of other segment, we can leave
0404             return Policy::disjoint();
0405         }
0406
0407         sides.set<1>(side_strategy_type::apply(p1, p2, q1),
0408                      side_strategy_type::apply(p1, p2, q2));
0409
0410         if (sides.same<1>())
0411         {
0412             // Both points are at same side of other segment, we can leave
0413             return Policy::disjoint();
0414         }
0415
0416         bool collinear = sides.collinear();
0417
0418         // r: ratio 0-1 where intersection divides A/B
0419         // (only calculated for non-collinear segments)
0420         if (! collinear)
0421         {
0422             using coordinate_type = typename select_most_precise
0423                 <
0424                     geometry::coordinate_type_t<point1_type>,
0425                     geometry::coordinate_type_t<point2_type>
0426                 >::type;
0427
0428             coordinate_type denominator_a;
0429             coordinate_type nominator_a;
0430             coordinate_type denominator_b;
0431             coordinate_type nominator_b;
0432
0433             cramers_rule(sinfo.dx_a, sinfo.dy_a, sinfo.dx_b, sinfo.dy_b,
0434                 get<0>(p1) - get<0>(q1),
0435                 get<1>(p1) - get<1>(q1),
0436                 nominator_a, denominator_a);
0437
0438             cramers_rule(sinfo.dx_b, sinfo.dy_b, sinfo.dx_a, sinfo.dy_a,
0439                 get<0>(q1) - get<0>(p1),
0440                 get<1>(q1) - get<1>(p1),
0441                 nominator_b, denominator_b);
0442
0443             math::detail::equals_factor_policy<coordinate_type>
0444                 policy(sinfo.dx_a, sinfo.dy_a, sinfo.dx_b, sinfo.dy_b);
0445
0446             coordinate_type const zero = 0;
0447             if (math::detail::equals_by_policy(denominator_a, zero, policy)
0448              || math::detail::equals_by_policy(denominator_b, zero, policy))
0449             {
0450                 // We set it to collinear, but it indicates a robustness issue.
0451                 sides.set<0>(0, 0);
0452                 sides.set<1>(0, 0);
0453                 collinear = true;
0454             }
0455             else
0456             {
0457                 sinfo.ra.assign(nominator_a, denominator_a);
0458                 sinfo.rb.assign(nominator_b, denominator_b);
0459             }
0460         }
0461
0462         if (collinear)
0463         {
0464             std::pair<bool, bool> const collinear_use_first
0465                     = is_x_more_significant(geometry::math::abs(sinfo.dx_a),
0466                                             geometry::math::abs(sinfo.dy_a),
0467                                             geometry::math::abs(sinfo.dx_b),
0468                                             geometry::math::abs(sinfo.dy_b),
0469                                             p_is_point, q_is_point);
0470
0471             if (collinear_use_first.second)
0472             {
0473                 // Degenerate cases: segments of single point, lying on other segment, are not disjoint
0474                 // This situation is collinear too
0475
0476                 if (collinear_use_first.first)
0477                 {
0478                     return relate_collinear<0, Policy, RatioType>(p, q,
0479                             p1, p2, q1, q2,
0480                             p_is_point, q_is_point);
0481                 }
0482                 else
0483                 {
0484                     // Y direction contains larger segments (maybe dx is zero)
0485                     return relate_collinear<1, Policy, RatioType>(p, q,
0486                             p1, p2, q1, q2,
0487                             p_is_point, q_is_point);
0488                 }
0489             }
0490         }
0491
0492         if (equals_point_point(p1, q1) || equals_point_point(p1, q2))
0493         {
0494             return Policy::segments_share_common_point(sides, sinfo, p1);
0495         }
0496
0497         if (equals_point_point(p2, q1) || equals_point_point(p2, q2))
0498         {
0499             return Policy::segments_share_common_point(sides, sinfo, p2);
0500         }
0501
0502         return Policy::segments_crosses(sides, sinfo, p, q);
0503     }
0504
0505 private:
0506     // first is true if x is more significant
0507     // second is true if the more significant difference is not 0
0508     template <typename CoordinateType>
0509     static inline std::pair<bool, bool>
0510         is_x_more_significant(CoordinateType const& abs_dx_a,
0511                               CoordinateType const& abs_dy_a,
0512                               CoordinateType const& abs_dx_b,
0513                               CoordinateType const& abs_dy_b,
0514                               bool const a_is_point,
0515                               bool const b_is_point)
0516     {
0517         //BOOST_GEOMETRY_ASSERT_MSG(!(a_is_point && b_is_point), "both segments shouldn't be degenerated");
0518
0519         // for degenerated segments the second is always true because this function
0520         // shouldn't be called if both segments were degenerated
0521
0522         if (a_is_point)
0523         {
0524             return std::make_pair(abs_dx_b >= abs_dy_b, true);
0525         }
0526         else if (b_is_point)
0527         {
0528             return std::make_pair(abs_dx_a >= abs_dy_a, true);
0529         }
0530         else
0531         {
0532             CoordinateType const min_dx = (std::min)(abs_dx_a, abs_dx_b);
0533             CoordinateType const min_dy = (std::min)(abs_dy_a, abs_dy_b);
0534             return min_dx == min_dy ?
0535                     std::make_pair(true, min_dx > CoordinateType(0)) :
0536                     std::make_pair(min_dx > min_dy, true);
0537         }
0538     }
0539
0540     template
0541     <
0542         std::size_t Dimension,
0543         typename Policy,
0544         typename RatioType,
0545         typename Segment1,
0546         typename Segment2,
0547         typename Point1,
0548         typename Point2
0549     >
0550     static inline typename Policy::return_type
0551         relate_collinear(Segment1 const& a,
0552                          Segment2 const& b,
0553                          Point1 const& a1, Point1 const& a2,
0554                          Point2 const& b1, Point2 const& b2,
0555                          bool a_is_point, bool b_is_point)
0556     {
0557         if (a_is_point)
0558         {
0559             return relate_one_degenerate<Policy, RatioType>(a,
0560                 get<Dimension>(a1),
0561                 get<Dimension>(b1), get<Dimension>(b2),
0562                 true);
0563         }
0564         if (b_is_point)
0565         {
0566             return relate_one_degenerate<Policy, RatioType>(b,
0567                 get<Dimension>(b1),
0568                 get<Dimension>(a1), get<Dimension>(a2),
0569                 false);
0570         }
0571         return relate_collinear<Policy, RatioType>(a, b,
0572                                 get<Dimension>(a1),
0573                                 get<Dimension>(a2),
0574                                 get<Dimension>(b1),
0575                                 get<Dimension>(b2));
0576     }
0577
0578     /// Relate segments known collinear
0579     template
0580     <
0581         typename Policy,
0582         typename RatioType,
0583         typename Segment1,
0584         typename Segment2,
0585         typename Type1,
0586         typename Type2
0587     >
0588     static inline typename Policy::return_type
0589         relate_collinear(Segment1 const& a, Segment2 const& b,
0590                          Type1 oa_1, Type1 oa_2,
0591                          Type2 ob_1, Type2 ob_2)
0592     {
0593         // Calculate the ratios where a starts in b, b starts in a
0594         //         a1--------->a2         (2..7)
0595         //                b1----->b2      (5..8)
0596         // length_a: 7-2=5
0597         // length_b: 8-5=3
0598         // b1 is located w.r.t. a at ratio: (5-2)/5=3/5 (on a)
0599         // b2 is located w.r.t. a at ratio: (8-2)/5=6/5 (right of a)
0600         // a1 is located w.r.t. b at ratio: (2-5)/3=-3/3 (left of b)
0601         // a2 is located w.r.t. b at ratio: (7-5)/3=2/3 (on b)
0602         // A arrives (a2 on b), B departs (b1 on a)
0603
0604         // If both are reversed:
0605         //         a2<---------a1         (7..2)
0606         //                b2<-----b1      (8..5)
0607         // length_a: 2-7=-5
0608         // length_b: 5-8=-3
0609         // b1 is located w.r.t. a at ratio: (8-7)/-5=-1/5 (before a starts)
0610         // b2 is located w.r.t. a at ratio: (5-7)/-5=2/5 (on a)
0611         // a1 is located w.r.t. b at ratio: (7-8)/-3=1/3 (on b)
0612         // a2 is located w.r.t. b at ratio: (2-8)/-3=6/3 (after b ends)
0613
0614         // If both one is reversed:
0615         //         a1--------->a2         (2..7)
0616         //                b2<-----b1      (8..5)
0617         // length_a: 7-2=+5
0618         // length_b: 5-8=-3
0619         // b1 is located w.r.t. a at ratio: (8-2)/5=6/5 (after a ends)
0620         // b2 is located w.r.t. a at ratio: (5-2)/5=3/5 (on a)
0621         // a1 is located w.r.t. b at ratio: (2-8)/-3=6/3 (after b ends)
0622         // a2 is located w.r.t. b at ratio: (7-8)/-3=1/3 (on b)
0623         Type1 const length_a = oa_2 - oa_1; // no abs, see above
0624         Type2 const length_b = ob_2 - ob_1;
0625
0626         RatioType ra_from(oa_1 - ob_1, length_b);
0627         RatioType ra_to(oa_2 - ob_1, length_b);
0628         RatioType rb_from(ob_1 - oa_1, length_a);
0629         RatioType rb_to(ob_2 - oa_1, length_a);
0630
0631         // use absolute measure to detect endpoints intersection
0632         // NOTE: it'd be possible to calculate bx_wrt_a using ax_wrt_b values
0633         int const a1_wrt_b = position_value(oa_1, ob_1, ob_2);
0634         int const a2_wrt_b = position_value(oa_2, ob_1, ob_2);
0635         int const b1_wrt_a = position_value(ob_1, oa_1, oa_2);
0636         int const b2_wrt_a = position_value(ob_2, oa_1, oa_2);
0637
0638         // fix the ratios if necessary
0639         // CONSIDER: fixing ratios also in other cases, if they're inconsistent
0640         // e.g. if ratio == 1 or 0 (so IP at the endpoint)
0641         // but position value indicates that the IP is in the middle of the segment
0642         // because one of the segments is very long
0643         // In such case the ratios could be moved into the middle direction
0644         // by some small value (e.g. EPS+1ULP)
0645         if (a1_wrt_b == 1)
0646         {
0647             ra_from.assign(0, 1);
0648             rb_from.assign(0, 1);
0649         }
0650         else if (a1_wrt_b == 3)
0651         {
0652             ra_from.assign(1, 1);
0653             rb_to.assign(0, 1);
0654         }
0655
0656         if (a2_wrt_b == 1)
0657         {
0658             ra_to.assign(0, 1);
0659             rb_from.assign(1, 1);
0660         }
0661         else if (a2_wrt_b == 3)
0662         {
0663             ra_to.assign(1, 1);
0664             rb_to.assign(1, 1);
0665         }
0666
0667         if ((a1_wrt_b < 1 && a2_wrt_b < 1) || (a1_wrt_b > 3 && a2_wrt_b > 3))
0668         //if ((ra_from.left() && ra_to.left()) || (ra_from.right() && ra_to.right()))
0669         {
0670             return Policy::disjoint();
0671         }
0672
0673         bool const opposite = math::sign(length_a) != math::sign(length_b);
0674
0675         return Policy::segments_collinear(a, b, opposite,
0676                                           a1_wrt_b, a2_wrt_b, b1_wrt_a, b2_wrt_a,
0677                                           ra_from, ra_to, rb_from, rb_to);
0678     }
0679
0680     /// Relate segments where one is degenerate
0681     template
0682     <
0683         typename Policy,
0684         typename RatioType,
0685         typename DegenerateSegment,
0686         typename Type1,
0687         typename Type2
0688     >
0689     static inline typename Policy::return_type
0690         relate_one_degenerate(DegenerateSegment const& degenerate_segment,
0691                               Type1 d, Type2 s1, Type2 s2,
0692                               bool a_degenerate)
0693     {
0694         // Calculate the ratios where ds starts in s
0695         //         a1--------->a2         (2..6)
0696         //              b1/b2      (4..4)
0697         // Ratio: (4-2)/(6-2)
0698         RatioType const ratio(d - s1, s2 - s1);
0699
0700         if (!ratio.on_segment())
0701         {
0702             return Policy::disjoint();
0703         }
0704
0705         return Policy::one_degenerate(degenerate_segment, ratio, a_degenerate);
0706     }
0707
0708     template <typename ProjCoord1, typename ProjCoord2>
0709     static inline int position_value(ProjCoord1 const& ca1,
0710                                      ProjCoord2 const& cb1,
0711                                      ProjCoord2 const& cb2)
0712     {
0713         // S1x  0   1    2     3   4
0714         // S2       |---------->
0715         return math::equals(ca1, cb1) ? 1
0716              : math::equals(ca1, cb2) ? 3
0717              : cb1 < cb2 ?
0718                 ( ca1 < cb1 ? 0
0719                 : ca1 > cb2 ? 4
0720                 : 2 )
0721               : ( ca1 > cb1 ? 0
0722                 : ca1 < cb2 ? 4
0723                 : 2 );
0724     }
0725
0726     template <typename Point1, typename Point2>
0727     static inline bool equals_point_point(Point1 const& point1, Point2 const& point2)
0728     {
0729         return strategy::within::cartesian_point_point::apply(point1, point2);
0730     }
0731 };
0732
0733
0734 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
0735 namespace services
0736 {
0737
0738 template <typename CalculationType>
0739 struct default_strategy<cartesian_tag, CalculationType>
0740 {
0741     typedef cartesian_segments<CalculationType> type;
0742 };
0743
0744 } // namespace services
0745 #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
0746
0747
0748 }} // namespace strategy::intersection
0749
0750 namespace strategy
0751 {
0752
0753 namespace within { namespace services
0754 {
0755
0756 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
0757 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, linear_tag, cartesian_tag, cartesian_tag>
0758 {
0759     typedef strategy::intersection::cartesian_segments<> type;
0760 };
0761
0762 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
0763 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, polygonal_tag, cartesian_tag, cartesian_tag>
0764 {
0765     typedef strategy::intersection::cartesian_segments<> type;
0766 };
0767
0768 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
0769 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, linear_tag, cartesian_tag, cartesian_tag>
0770 {
0771     typedef strategy::intersection::cartesian_segments<> type;
0772 };
0773
0774 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
0775 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, polygonal_tag, cartesian_tag, cartesian_tag>
0776 {
0777     typedef strategy::intersection::cartesian_segments<> type;
0778 };
0779
0780 }} // within::services
0781
0782 namespace covered_by { namespace services
0783 {
0784
0785 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
0786 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, linear_tag, cartesian_tag, cartesian_tag>
0787 {
0788     typedef strategy::intersection::cartesian_segments<> type;
0789 };
0790
0791 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
0792 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, polygonal_tag, cartesian_tag, cartesian_tag>
0793 {
0794     typedef strategy::intersection::cartesian_segments<> type;
0795 };
0796
0797 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
0798 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, linear_tag, cartesian_tag, cartesian_tag>
0799 {
0800     typedef strategy::intersection::cartesian_segments<> type;
0801 };
0802
0803 template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2>
0804 struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, polygonal_tag, cartesian_tag, cartesian_tag>
0805 {
0806     typedef strategy::intersection::cartesian_segments<> type;
0807 };
0808
0809 }} // within::services
0810
0811 } // strategy
0812
0813 }} // namespace boost::geometry
0814
0815
0816 #endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_INTERSECTION_HPP