intersection.hpp 28 KB

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