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