3 // Copyright (c) 2017 Adam Wulkiewicz, Lodz, Poland.
5 // Copyright (c) 2016-2019, Oracle and/or its affiliates.
6 // 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)
12 #ifndef BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_INTERSECTION_HPP
13 #define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_INTERSECTION_HPP
17 #include <boost/geometry/core/cs.hpp>
18 #include <boost/geometry/core/access.hpp>
19 #include <boost/geometry/core/radian_access.hpp>
20 #include <boost/geometry/core/tags.hpp>
22 #include <boost/geometry/algorithms/detail/assign_values.hpp>
23 #include <boost/geometry/algorithms/detail/assign_indexed_point.hpp>
24 #include <boost/geometry/algorithms/detail/equals/point_point.hpp>
25 #include <boost/geometry/algorithms/detail/recalculate.hpp>
27 #include <boost/geometry/formulas/andoyer_inverse.hpp>
28 #include <boost/geometry/formulas/sjoberg_intersection.hpp>
29 #include <boost/geometry/formulas/spherical.hpp>
30 #include <boost/geometry/formulas/unit_spheroid.hpp>
32 #include <boost/geometry/geometries/concepts/point_concept.hpp>
33 #include <boost/geometry/geometries/concepts/segment_concept.hpp>
35 #include <boost/geometry/policies/robustness/segment_ratio.hpp>
37 #include <boost/geometry/srs/spheroid.hpp>
39 #include <boost/geometry/strategies/geographic/area.hpp>
40 #include <boost/geometry/strategies/geographic/disjoint_segment_box.hpp>
41 #include <boost/geometry/strategies/geographic/distance.hpp>
42 #include <boost/geometry/strategies/geographic/envelope.hpp>
43 #include <boost/geometry/strategies/geographic/parameters.hpp>
44 #include <boost/geometry/strategies/geographic/point_in_poly_winding.hpp>
45 #include <boost/geometry/strategies/geographic/side.hpp>
46 #include <boost/geometry/strategies/spherical/expand_box.hpp>
47 #include <boost/geometry/strategies/spherical/disjoint_box_box.hpp>
48 #include <boost/geometry/strategies/spherical/point_in_point.hpp>
49 #include <boost/geometry/strategies/intersection.hpp>
50 #include <boost/geometry/strategies/intersection_result.hpp>
51 #include <boost/geometry/strategies/side_info.hpp>
53 #include <boost/geometry/util/math.hpp>
54 #include <boost/geometry/util/select_calculation_type.hpp>
57 namespace boost { namespace geometry
60 namespace strategy { namespace intersection
63 // CONSIDER: Improvement of the robustness/accuracy/repeatability by
64 // moving all segments to 0 longitude
65 // picking latitudes closer to 0
70 typename FormulaPolicy = strategy::andoyer,
71 unsigned int Order = strategy::default_order<FormulaPolicy>::value,
72 typename Spheroid = srs::spheroid<double>,
73 typename CalculationType = void
75 struct geographic_segments
77 typedef geographic_tag cs_tag;
79 typedef side::geographic
81 FormulaPolicy, Spheroid, CalculationType
84 inline side_strategy_type get_side_strategy() const
86 return side_strategy_type(m_spheroid);
89 template <typename Geometry1, typename Geometry2>
90 struct point_in_geometry_strategy
92 typedef strategy::within::geographic_winding
94 typename point_type<Geometry1>::type,
95 typename point_type<Geometry2>::type,
102 template <typename Geometry1, typename Geometry2>
103 inline typename point_in_geometry_strategy<Geometry1, Geometry2>::type
104 get_point_in_geometry_strategy() const
106 typedef typename point_in_geometry_strategy
109 >::type strategy_type;
110 return strategy_type(m_spheroid);
113 template <typename Geometry>
116 typedef area::geographic
125 template <typename Geometry>
126 inline typename area_strategy<Geometry>::type get_area_strategy() const
128 typedef typename area_strategy<Geometry>::type strategy_type;
129 return strategy_type(m_spheroid);
132 template <typename Geometry>
133 struct distance_strategy
135 typedef distance::geographic
143 template <typename Geometry>
144 inline typename distance_strategy<Geometry>::type get_distance_strategy() const
146 typedef typename distance_strategy<Geometry>::type strategy_type;
147 return strategy_type(m_spheroid);
150 typedef envelope::geographic<FormulaPolicy, Spheroid, CalculationType>
151 envelope_strategy_type;
153 inline envelope_strategy_type get_envelope_strategy() const
155 return envelope_strategy_type(m_spheroid);
158 typedef expand::geographic_segment<FormulaPolicy, Spheroid, CalculationType>
159 expand_strategy_type;
161 inline expand_strategy_type get_expand_strategy() const
163 return expand_strategy_type(m_spheroid);
166 typedef within::spherical_point_point point_in_point_strategy_type;
168 static inline point_in_point_strategy_type get_point_in_point_strategy()
170 return point_in_point_strategy_type();
173 typedef within::spherical_point_point equals_point_point_strategy_type;
175 static inline equals_point_point_strategy_type get_equals_point_point_strategy()
177 return equals_point_point_strategy_type();
180 typedef disjoint::spherical_box_box disjoint_box_box_strategy_type;
182 static inline disjoint_box_box_strategy_type get_disjoint_box_box_strategy()
184 return disjoint_box_box_strategy_type();
187 typedef disjoint::segment_box_geographic
189 FormulaPolicy, Spheroid, CalculationType
190 > disjoint_segment_box_strategy_type;
192 inline disjoint_segment_box_strategy_type get_disjoint_segment_box_strategy() const
194 return disjoint_segment_box_strategy_type(m_spheroid);
197 typedef covered_by::spherical_point_box disjoint_point_box_strategy_type;
198 typedef covered_by::spherical_point_box covered_by_point_box_strategy_type;
199 typedef within::spherical_point_box within_point_box_strategy_type;
200 typedef envelope::spherical_box envelope_box_strategy_type;
201 typedef expand::spherical_box expand_box_strategy_type;
203 enum intersection_point_flag { ipi_inters = 0, ipi_at_a1, ipi_at_a2, ipi_at_b1, ipi_at_b2 };
205 template <typename CoordinateType, typename SegmentRatio>
206 struct segment_intersection_info
208 template <typename Point, typename Segment1, typename Segment2>
209 void calculate(Point& point, Segment1 const& a, Segment2 const& b) const
211 if (ip_flag == ipi_inters)
213 // TODO: assign the rest of coordinates
214 set_from_radian<0>(point, lon);
215 set_from_radian<1>(point, lat);
217 else if (ip_flag == ipi_at_a1)
219 detail::assign_point_from_index<0>(a, point);
221 else if (ip_flag == ipi_at_a2)
223 detail::assign_point_from_index<1>(a, point);
225 else if (ip_flag == ipi_at_b1)
227 detail::assign_point_from_index<0>(b, point);
229 else // ip_flag == ipi_at_b2
231 detail::assign_point_from_index<1>(b, point);
237 SegmentRatio robust_ra;
238 SegmentRatio robust_rb;
239 intersection_point_flag ip_flag;
242 explicit geographic_segments(Spheroid const& spheroid = Spheroid())
243 : m_spheroid(spheroid)
246 // Relate segments a and b
249 typename UniqueSubRange1,
250 typename UniqueSubRange2,
253 inline typename Policy::return_type apply(UniqueSubRange1 const& range_p,
254 UniqueSubRange2 const& range_q,
257 typedef typename UniqueSubRange1::point_type point1_type;
258 typedef typename UniqueSubRange2::point_type point2_type;
259 typedef model::referring_segment<point1_type const> segment_type1;
260 typedef model::referring_segment<point2_type const> segment_type2;
262 BOOST_CONCEPT_ASSERT( (concepts::ConstPoint<point1_type>) );
263 BOOST_CONCEPT_ASSERT( (concepts::ConstPoint<point2_type>) );
266 typename coordinate_type<Point1>::type
267 const a1_lon = get<0>(a1),
268 const a2_lon = get<0>(a2);
269 typename coordinate_type<Point2>::type
270 const b1_lon = get<0>(b1),
271 const b2_lon = get<0>(b2);
272 bool is_a_reversed = a1_lon > a2_lon || a1_lon == a2_lon && get<1>(a1) > get<1>(a2);
273 bool is_b_reversed = b1_lon > b2_lon || b1_lon == b2_lon && get<1>(b1) > get<1>(b2);
276 bool const is_p_reversed = get<1>(range_p.at(0)) > get<1>(range_p.at(1));
277 bool const is_q_reversed = get<1>(range_q.at(0)) > get<1>(range_q.at(1));
279 // Call apply with original segments and ordered points
280 return apply<Policy>(segment_type1(range_p.at(0), range_p.at(1)),
281 segment_type2(range_q.at(0), range_q.at(1)),
282 range_p.at(is_p_reversed ? 1 : 0),
283 range_p.at(is_p_reversed ? 0 : 1),
284 range_q.at(is_q_reversed ? 1 : 0),
285 range_q.at(is_q_reversed ? 0 : 1),
286 is_p_reversed, is_q_reversed);
290 // Relate segments a and b
299 inline typename Policy::return_type apply(Segment1 const& a, Segment2 const& b,
300 Point1 const& a1, Point1 const& a2,
301 Point2 const& b1, Point2 const& b2,
302 bool is_a_reversed, bool is_b_reversed) const
304 BOOST_CONCEPT_ASSERT( (concepts::ConstSegment<Segment1>) );
305 BOOST_CONCEPT_ASSERT( (concepts::ConstSegment<Segment2>) );
307 typedef typename select_calculation_type
308 <Segment1, Segment2, CalculationType>::type calc_t;
310 typedef srs::spheroid<calc_t> spheroid_type;
312 static const calc_t c0 = 0;
314 // normalized spheroid
315 spheroid_type spheroid = formula::unit_spheroid<spheroid_type>(m_spheroid);
317 // TODO: check only 2 first coordinates here?
318 bool a_is_point = equals_point_point(a1, a2);
319 bool b_is_point = equals_point_point(b1, b2);
321 if(a_is_point && b_is_point)
323 return equals_point_point(a1, b2)
324 ? Policy::degenerate(a, true)
329 calc_t const a1_lon = get_as_radian<0>(a1);
330 calc_t const a1_lat = get_as_radian<1>(a1);
331 calc_t const a2_lon = get_as_radian<0>(a2);
332 calc_t const a2_lat = get_as_radian<1>(a2);
333 calc_t const b1_lon = get_as_radian<0>(b1);
334 calc_t const b1_lat = get_as_radian<1>(b1);
335 calc_t const b2_lon = get_as_radian<0>(b2);
336 calc_t const b2_lat = get_as_radian<1>(b2);
340 // NOTE: potential optimization, don't calculate distance at this point
341 // this would require to reimplement inverse strategy to allow
342 // calculation of distance if needed, probably also storing intermediate
343 // results somehow inside an object.
344 typedef typename FormulaPolicy::template inverse<calc_t, true, true, false, false, false> inverse_dist_azi;
345 typedef typename inverse_dist_azi::result_type inverse_result;
347 // TODO: no need to call inverse formula if we know that the points are equal
348 // distance can be set to 0 in this case and azimuth may be not calculated
349 bool is_equal_a1_b1 = equals_point_point(a1, b1);
350 bool is_equal_a2_b1 = equals_point_point(a2, b1);
351 bool degen_neq_coords = false;
353 inverse_result res_b1_b2, res_b1_a1, res_b1_a2;
356 res_b1_b2 = inverse_dist_azi::apply(b1_lon, b1_lat, b2_lon, b2_lat, spheroid);
357 if (math::equals(res_b1_b2.distance, c0))
360 degen_neq_coords = true;
364 res_b1_a1 = inverse_dist_azi::apply(b1_lon, b1_lat, a1_lon, a1_lat, spheroid);
365 if (math::equals(res_b1_a1.distance, c0))
367 is_equal_a1_b1 = true;
369 res_b1_a2 = inverse_dist_azi::apply(b1_lon, b1_lat, a2_lon, a2_lat, spheroid);
370 if (math::equals(res_b1_a2.distance, c0))
372 is_equal_a2_b1 = true;
374 sides.set<0>(is_equal_a1_b1 ? 0 : formula::azimuth_side_value(res_b1_a1.azimuth, res_b1_b2.azimuth),
375 is_equal_a2_b1 ? 0 : formula::azimuth_side_value(res_b1_a2.azimuth, res_b1_b2.azimuth));
378 // Both points are at the same side of other segment, we can leave
379 return Policy::disjoint();
384 bool is_equal_a1_b2 = equals_point_point(a1, b2);
386 inverse_result res_a1_a2, res_a1_b1, res_a1_b2;
389 res_a1_a2 = inverse_dist_azi::apply(a1_lon, a1_lat, a2_lon, a2_lat, spheroid);
390 if (math::equals(res_a1_a2.distance, c0))
393 degen_neq_coords = true;
397 res_a1_b1 = inverse_dist_azi::apply(a1_lon, a1_lat, b1_lon, b1_lat, spheroid);
398 if (math::equals(res_a1_b1.distance, c0))
400 is_equal_a1_b1 = true;
402 res_a1_b2 = inverse_dist_azi::apply(a1_lon, a1_lat, b2_lon, b2_lat, spheroid);
403 if (math::equals(res_a1_b2.distance, c0))
405 is_equal_a1_b2 = true;
407 sides.set<1>(is_equal_a1_b1 ? 0 : formula::azimuth_side_value(res_a1_b1.azimuth, res_a1_a2.azimuth),
408 is_equal_a1_b2 ? 0 : formula::azimuth_side_value(res_a1_b2.azimuth, res_a1_a2.azimuth));
411 // Both points are at the same side of other segment, we can leave
412 return Policy::disjoint();
417 if(a_is_point && b_is_point)
419 return is_equal_a1_b2
420 ? Policy::degenerate(a, true)
425 // NOTE: at this point the segments may still be disjoint
426 // NOTE: at this point one of the segments may be degenerated
428 bool collinear = sides.collinear();
432 // WARNING: the side strategy doesn't have the info about the other
433 // segment so it may return results inconsistent with this intersection
434 // strategy, as it checks both segments for consistency
436 if (sides.get<0, 0>() == 0 && sides.get<0, 1>() == 0)
441 else if (sides.get<1, 0>() == 0 && sides.get<1, 1>() == 0)
452 return collinear_one_degenerated<Policy, calc_t>(a, true, b1, b2, a1, a2, res_b1_b2, res_b1_a1, res_b1_a2, is_b_reversed, degen_neq_coords);
456 return collinear_one_degenerated<Policy, calc_t>(b, false, a1, a2, b1, b2, res_a1_a2, res_a1_b1, res_a1_b2, is_a_reversed, degen_neq_coords);
460 calc_t dist_a1_a2, dist_a1_b1, dist_a1_b2;
461 calc_t dist_b1_b2, dist_b1_a1, dist_b1_a2;
462 // use shorter segment
463 if (res_a1_a2.distance <= res_b1_b2.distance)
465 calculate_collinear_data(a1, a2, b1, b2, res_a1_a2, res_a1_b1, res_a1_b2, dist_a1_a2, dist_a1_b1);
466 calculate_collinear_data(a1, a2, b2, b1, res_a1_a2, res_a1_b2, res_a1_b1, dist_a1_a2, dist_a1_b2);
467 dist_b1_b2 = dist_a1_b2 - dist_a1_b1;
468 dist_b1_a1 = -dist_a1_b1;
469 dist_b1_a2 = dist_a1_a2 - dist_a1_b1;
473 calculate_collinear_data(b1, b2, a1, a2, res_b1_b2, res_b1_a1, res_b1_a2, dist_b1_b2, dist_b1_a1);
474 calculate_collinear_data(b1, b2, a2, a1, res_b1_b2, res_b1_a2, res_b1_a1, dist_b1_b2, dist_b1_a2);
475 dist_a1_a2 = dist_b1_a2 - dist_b1_a1;
476 dist_a1_b1 = -dist_b1_a1;
477 dist_a1_b2 = dist_b1_b2 - dist_b1_a1;
480 // NOTE: this is probably not needed
481 int a1_on_b = position_value(c0, dist_a1_b1, dist_a1_b2);
482 int a2_on_b = position_value(dist_a1_a2, dist_a1_b1, dist_a1_b2);
483 int b1_on_a = position_value(c0, dist_b1_a1, dist_b1_a2);
484 int b2_on_a = position_value(dist_b1_b2, dist_b1_a1, dist_b1_a2);
486 if ((a1_on_b < 1 && a2_on_b < 1) || (a1_on_b > 3 && a2_on_b > 3))
488 return Policy::disjoint();
496 else if (a1_on_b == 3)
498 dist_b1_a1 = dist_b1_b2;
505 dist_a1_b1 = dist_a1_a2;
507 else if (a2_on_b == 3)
509 dist_b1_a2 = dist_b1_b2;
510 dist_a1_b2 = dist_a1_a2;
513 bool opposite = ! same_direction(res_a1_a2.azimuth, res_b1_b2.azimuth);
515 // NOTE: If segment was reversed opposite, positions and segment ratios has to be altered
519 opposite = ! opposite;
521 std::swap(a1_on_b, a2_on_b);
522 b1_on_a = 4 - b1_on_a;
523 b2_on_a = 4 - b2_on_a;
524 // distances for ratios
525 std::swap(dist_b1_a1, dist_b1_a2);
526 dist_a1_b1 = dist_a1_a2 - dist_a1_b1;
527 dist_a1_b2 = dist_a1_a2 - dist_a1_b2;
532 opposite = ! opposite;
534 a1_on_b = 4 - a1_on_b;
535 a2_on_b = 4 - a2_on_b;
536 std::swap(b1_on_a, b2_on_a);
537 // distances for ratios
538 dist_b1_a1 = dist_b1_b2 - dist_b1_a1;
539 dist_b1_a2 = dist_b1_b2 - dist_b1_a2;
540 std::swap(dist_a1_b1, dist_a1_b2);
543 segment_ratio<calc_t> ra_from(dist_b1_a1, dist_b1_b2);
544 segment_ratio<calc_t> ra_to(dist_b1_a2, dist_b1_b2);
545 segment_ratio<calc_t> rb_from(dist_a1_b1, dist_a1_a2);
546 segment_ratio<calc_t> rb_to(dist_a1_b2, dist_a1_a2);
548 return Policy::segments_collinear(a, b, opposite,
549 a1_on_b, a2_on_b, b1_on_a, b2_on_a,
550 ra_from, ra_to, rb_from, rb_to);
553 else // crossing or touching
555 if (a_is_point || b_is_point)
557 return Policy::disjoint();
560 calc_t lon = 0, lat = 0;
561 intersection_point_flag ip_flag;
562 calc_t dist_a1_a2, dist_a1_i1, dist_b1_b2, dist_b1_i1;
563 if (calculate_ip_data(a1, a2, b1, b2,
564 a1_lon, a1_lat, a2_lon, a2_lat,
565 b1_lon, b1_lat, b2_lon, b2_lat,
566 res_a1_a2, res_a1_b1, res_a1_b2,
567 res_b1_b2, res_b1_a1, res_b1_a2,
570 dist_a1_a2, dist_a1_i1, dist_b1_b2, dist_b1_i1,
573 // NOTE: If segment was reversed sides and segment ratios has to be altered
577 sides_reverse_segment<0>(sides);
578 // distance for ratio
579 dist_a1_i1 = dist_a1_a2 - dist_a1_i1;
581 ip_flag_reverse_segment(ip_flag, ipi_at_a1, ipi_at_a2);
586 sides_reverse_segment<1>(sides);
587 // distance for ratio
588 dist_b1_i1 = dist_b1_b2 - dist_b1_i1;
590 ip_flag_reverse_segment(ip_flag, ipi_at_b1, ipi_at_b2);
594 segment_intersection_info
597 segment_ratio<calc_t>
602 sinfo.robust_ra.assign(dist_a1_i1, dist_a1_a2);
603 sinfo.robust_rb.assign(dist_b1_i1, dist_b1_b2);
604 sinfo.ip_flag = ip_flag;
606 return Policy::segments_crosses(sides, sinfo, a, b);
610 return Policy::disjoint();
615 template <typename Policy, typename CalcT, typename Segment, typename Point1, typename Point2, typename ResultInverse>
616 static inline typename Policy::return_type
617 collinear_one_degenerated(Segment const& segment, bool degenerated_a,
618 Point1 const& a1, Point1 const& a2,
619 Point2 const& b1, Point2 const& b2,
620 ResultInverse const& res_a1_a2,
621 ResultInverse const& res_a1_b1,
622 ResultInverse const& res_a1_b2,
623 bool is_other_reversed,
624 bool degen_neq_coords)
626 CalcT dist_1_2, dist_1_o;
627 if (! calculate_collinear_data(a1, a2, b1, b2, res_a1_a2, res_a1_b1, res_a1_b2, dist_1_2, dist_1_o, degen_neq_coords))
629 return Policy::disjoint();
632 // NOTE: If segment was reversed segment ratio has to be altered
633 if (is_other_reversed)
635 // distance for ratio
636 dist_1_o = dist_1_2 - dist_1_o;
639 return Policy::one_degenerate(segment, segment_ratio<CalcT>(dist_1_o, dist_1_2), degenerated_a);
642 // TODO: instead of checks below test bi against a1 and a2 here?
643 // in order to make this independent from is_near()
644 template <typename Point1, typename Point2, typename ResultInverse, typename CalcT>
645 static inline bool calculate_collinear_data(Point1 const& a1, Point1 const& a2, // in
646 Point2 const& b1, Point2 const& /*b2*/, // in
647 ResultInverse const& res_a1_a2, // in
648 ResultInverse const& res_a1_b1, // in
649 ResultInverse const& res_a1_b2, // in
650 CalcT& dist_a1_a2, // out
651 CalcT& dist_a1_b1, // out
652 bool degen_neq_coords = false) // in
654 dist_a1_a2 = res_a1_a2.distance;
656 dist_a1_b1 = res_a1_b1.distance;
657 if (! same_direction(res_a1_b1.azimuth, res_a1_a2.azimuth))
659 dist_a1_b1 = -dist_a1_b1;
663 if (is_endpoint_equal(dist_a1_b1, a1, b1))
669 else if (is_endpoint_equal(dist_a1_a2 - dist_a1_b1, a2, b1))
671 dist_a1_b1 = dist_a1_a2;
675 // check the other endpoint of degenerated segment near a pole
676 if (degen_neq_coords)
678 static CalcT const c0 = 0;
679 if (math::equals(res_a1_b2.distance, c0))
684 else if (math::equals(dist_a1_a2 - res_a1_b2.distance, c0))
686 dist_a1_b1 = dist_a1_a2;
692 return segment_ratio<CalcT>(dist_a1_b1, dist_a1_a2).on_segment();
695 template <typename Point1, typename Point2, typename CalcT, typename ResultInverse, typename Spheroid_>
696 static inline bool calculate_ip_data(Point1 const& a1, Point1 const& a2, // in
697 Point2 const& b1, Point2 const& b2, // in
698 CalcT const& a1_lon, CalcT const& a1_lat, // in
699 CalcT const& a2_lon, CalcT const& a2_lat, // in
700 CalcT const& b1_lon, CalcT const& b1_lat, // in
701 CalcT const& b2_lon, CalcT const& b2_lat, // in
702 ResultInverse const& res_a1_a2, // in
703 ResultInverse const& res_a1_b1, // in
704 ResultInverse const& res_a1_b2, // in
705 ResultInverse const& res_b1_b2, // in
706 ResultInverse const& res_b1_a1, // in
707 ResultInverse const& res_b1_a2, // in
708 side_info const& sides, // in
709 Spheroid_ const& spheroid, // in
710 CalcT & lon, CalcT & lat, // out
711 CalcT& dist_a1_a2, CalcT& dist_a1_ip, // out
712 CalcT& dist_b1_b2, CalcT& dist_b1_ip, // out
713 intersection_point_flag& ip_flag) // out
715 dist_a1_a2 = res_a1_a2.distance;
716 dist_b1_b2 = res_b1_b2.distance;
718 // assign the IP if some endpoints overlap
719 if (equals_point_point(a1, b1))
728 else if (equals_point_point(a1, b2))
733 dist_b1_ip = dist_b1_b2;
737 else if (equals_point_point(a2, b1))
741 dist_a1_ip = dist_a1_a2;
746 else if (equals_point_point(a2, b2))
750 dist_a1_ip = dist_a1_a2;
751 dist_b1_ip = dist_b1_b2;
756 // at this point we know that the endpoints doesn't overlap
757 // check cases when an endpoint lies on the other geodesic
758 if (sides.template get<0, 0>() == 0) // a1 wrt b
760 if (res_b1_a1.distance <= res_b1_b2.distance
761 && same_direction(res_b1_a1.azimuth, res_b1_b2.azimuth))
766 dist_b1_ip = res_b1_a1.distance;
775 else if (sides.template get<0, 1>() == 0) // a2 wrt b
777 if (res_b1_a2.distance <= res_b1_b2.distance
778 && same_direction(res_b1_a2.azimuth, res_b1_b2.azimuth))
782 dist_a1_ip = res_a1_a2.distance;
783 dist_b1_ip = res_b1_a2.distance;
792 else if (sides.template get<1, 0>() == 0) // b1 wrt a
794 if (res_a1_b1.distance <= res_a1_a2.distance
795 && same_direction(res_a1_b1.azimuth, res_a1_a2.azimuth))
799 dist_a1_ip = res_a1_b1.distance;
809 else if (sides.template get<1, 1>() == 0) // b2 wrt a
811 if (res_a1_b2.distance <= res_a1_a2.distance
812 && same_direction(res_a1_b2.azimuth, res_a1_a2.azimuth))
816 dist_a1_ip = res_a1_b2.distance;
817 dist_b1_ip = res_b1_b2.distance;
827 // At this point neither the endpoints overlaps
828 // nor any andpoint lies on the other geodesic
829 // So the endpoints should lie on the opposite sides of both geodesics
831 bool const ok = formula::sjoberg_intersection<CalcT, FormulaPolicy::template inverse, Order>
832 ::apply(a1_lon, a1_lat, a2_lon, a2_lat, res_a1_a2.azimuth,
833 b1_lon, b1_lat, b2_lon, b2_lat, res_b1_b2.azimuth,
841 typedef typename FormulaPolicy::template inverse<CalcT, true, true, false, false, false> inverse_dist_azi;
842 typedef typename inverse_dist_azi::result_type inverse_result;
844 inverse_result const res_a1_ip = inverse_dist_azi::apply(a1_lon, a1_lat, lon, lat, spheroid);
845 dist_a1_ip = res_a1_ip.distance;
846 if (! same_direction(res_a1_ip.azimuth, res_a1_a2.azimuth))
848 dist_a1_ip = -dist_a1_ip;
851 bool is_on_a = segment_ratio<CalcT>(dist_a1_ip, dist_a1_a2).on_segment();
852 // NOTE: not fully consistent with equals_point_point() since radians are always used.
853 bool is_on_a1 = math::equals(lon, a1_lon) && math::equals(lat, a1_lat);
854 bool is_on_a2 = math::equals(lon, a2_lon) && math::equals(lat, a2_lat);
856 if (! (is_on_a || is_on_a1 || is_on_a2))
861 inverse_result const res_b1_ip = inverse_dist_azi::apply(b1_lon, b1_lat, lon, lat, spheroid);
862 dist_b1_ip = res_b1_ip.distance;
863 if (! same_direction(res_b1_ip.azimuth, res_b1_b2.azimuth))
865 dist_b1_ip = -dist_b1_ip;
868 bool is_on_b = segment_ratio<CalcT>(dist_b1_ip, dist_b1_b2).on_segment();
869 // NOTE: not fully consistent with equals_point_point() since radians are always used.
870 bool is_on_b1 = math::equals(lon, b1_lon) && math::equals(lat, b1_lat);
871 bool is_on_b2 = math::equals(lon, b2_lon) && math::equals(lat, b2_lat);
873 if (! (is_on_b || is_on_b1 || is_on_b2))
878 typedef typename FormulaPolicy::template inverse<CalcT, true, false, false, false, false> inverse_dist;
880 ip_flag = ipi_inters;
886 dist_a1_ip = inverse_dist::apply(a1_lon, a1_lat, lon, lat, spheroid).distance; // for consistency
894 dist_a1_ip = inverse_dist::apply(a1_lon, a1_lat, lon, lat, spheroid).distance; // for consistency
895 dist_b1_ip = res_b1_b2.distance;
904 dist_b1_ip = inverse_dist::apply(b1_lon, b1_lat, lon, lat, spheroid).distance; // for consistency
911 dist_a1_ip = res_a1_a2.distance;
912 dist_b1_ip = inverse_dist::apply(b1_lon, b1_lat, lon, lat, spheroid).distance; // for consistency
919 template <typename CalcT, typename P1, typename P2>
920 static inline bool is_endpoint_equal(CalcT const& dist,
921 P1 const& ai, P2 const& b1)
923 static CalcT const c0 = 0;
924 return is_near(dist) && (math::equals(dist, c0) || equals_point_point(ai, b1));
927 template <typename CalcT>
928 static inline bool is_near(CalcT const& dist)
930 // NOTE: This strongly depends on the Inverse method
931 CalcT const small_number = CalcT(boost::is_same<CalcT, float>::value ? 0.0001 : 0.00000001);
932 return math::abs(dist) <= small_number;
935 template <typename ProjCoord1, typename ProjCoord2>
936 static inline int position_value(ProjCoord1 const& ca1,
937 ProjCoord2 const& cb1,
938 ProjCoord2 const& cb2)
942 return math::equals(ca1, cb1) ? 1
943 : math::equals(ca1, cb2) ? 3
953 template <typename CalcT>
954 static inline bool same_direction(CalcT const& azimuth1, CalcT const& azimuth2)
956 // distance between two angles normalized to (-180, 180]
957 CalcT const angle_diff = math::longitude_distance_signed<radian>(azimuth1, azimuth2);
958 return math::abs(angle_diff) <= math::half_pi<CalcT>();
962 static inline void sides_reverse_segment(side_info & sides)
964 // names assuming segment A is reversed (Which == 0)
965 int a1_wrt_b = sides.template get<Which, 0>();
966 int a2_wrt_b = sides.template get<Which, 1>();
967 std::swap(a1_wrt_b, a2_wrt_b);
968 sides.template set<Which>(a1_wrt_b, a2_wrt_b);
969 int b1_wrt_a = sides.template get<1 - Which, 0>();
970 int b2_wrt_a = sides.template get<1 - Which, 1>();
971 sides.template set<1 - Which>(-b1_wrt_a, -b2_wrt_a);
974 static inline void ip_flag_reverse_segment(intersection_point_flag & ip_flag,
975 intersection_point_flag const& ipi_at_p1,
976 intersection_point_flag const& ipi_at_p2)
978 ip_flag = ip_flag == ipi_at_p1 ? ipi_at_p2 :
979 ip_flag == ipi_at_p2 ? ipi_at_p1 :
983 template <typename Point1, typename Point2>
984 static inline bool equals_point_point(Point1 const& point1, Point2 const& point2)
986 return detail::equals::equals_point_point(point1, point2,
987 point_in_point_strategy_type());
995 }} // namespace strategy::intersection
997 }} // namespace boost::geometry
1000 #endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_INTERSECTION_HPP