intersection.hpp 30 KB

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