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