]>
Commit | Line | Data |
---|---|---|
b32b8144 FG |
1 | // Boost.Geometry |
2 | ||
11fdf7f2 TL |
3 | // Copyright (c) 2017 Adam Wulkiewicz, Lodz, Poland. |
4 | ||
20effc67 | 5 | // Copyright (c) 2016-2020, Oracle and/or its affiliates. |
b32b8144 FG |
6 | // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle |
7 | ||
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 | ||
12 | #ifndef BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_INTERSECTION_HPP | |
13 | #define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_INTERSECTION_HPP | |
14 | ||
15 | #include <algorithm> | |
20effc67 | 16 | #include <type_traits> |
b32b8144 FG |
17 | |
18 | #include <boost/geometry/core/cs.hpp> | |
19 | #include <boost/geometry/core/access.hpp> | |
20 | #include <boost/geometry/core/radian_access.hpp> | |
b32b8144 FG |
21 | #include <boost/geometry/core/tags.hpp> |
22 | ||
23 | #include <boost/geometry/algorithms/detail/assign_values.hpp> | |
24 | #include <boost/geometry/algorithms/detail/assign_indexed_point.hpp> | |
25 | #include <boost/geometry/algorithms/detail/equals/point_point.hpp> | |
26 | #include <boost/geometry/algorithms/detail/recalculate.hpp> | |
27 | ||
28 | #include <boost/geometry/formulas/andoyer_inverse.hpp> | |
29 | #include <boost/geometry/formulas/sjoberg_intersection.hpp> | |
30 | #include <boost/geometry/formulas/spherical.hpp> | |
31 | #include <boost/geometry/formulas/unit_spheroid.hpp> | |
32 | ||
33 | #include <boost/geometry/geometries/concepts/point_concept.hpp> | |
34 | #include <boost/geometry/geometries/concepts/segment_concept.hpp> | |
20effc67 | 35 | #include <boost/geometry/geometries/segment.hpp> |
b32b8144 FG |
36 | |
37 | #include <boost/geometry/policies/robustness/segment_ratio.hpp> | |
38 | ||
11fdf7f2 TL |
39 | #include <boost/geometry/srs/spheroid.hpp> |
40 | ||
20effc67 TL |
41 | #include <boost/geometry/strategy/geographic/area.hpp> |
42 | #include <boost/geometry/strategy/geographic/envelope.hpp> | |
43 | #include <boost/geometry/strategy/geographic/expand_segment.hpp> | |
44 | #include <boost/geometry/strategy/spherical/expand_box.hpp> | |
45 | ||
92f5a8d4 | 46 | #include <boost/geometry/strategies/geographic/disjoint_segment_box.hpp> |
b32b8144 | 47 | #include <boost/geometry/strategies/geographic/distance.hpp> |
b32b8144 FG |
48 | #include <boost/geometry/strategies/geographic/parameters.hpp> |
49 | #include <boost/geometry/strategies/geographic/point_in_poly_winding.hpp> | |
50 | #include <boost/geometry/strategies/geographic/side.hpp> | |
92f5a8d4 TL |
51 | #include <boost/geometry/strategies/spherical/disjoint_box_box.hpp> |
52 | #include <boost/geometry/strategies/spherical/point_in_point.hpp> | |
b32b8144 FG |
53 | #include <boost/geometry/strategies/intersection.hpp> |
54 | #include <boost/geometry/strategies/intersection_result.hpp> | |
55 | #include <boost/geometry/strategies/side_info.hpp> | |
56 | ||
57 | #include <boost/geometry/util/math.hpp> | |
58 | #include <boost/geometry/util/select_calculation_type.hpp> | |
59 | ||
60 | ||
61 | namespace boost { namespace geometry | |
62 | { | |
63 | ||
64 | namespace strategy { namespace intersection | |
65 | { | |
66 | ||
67 | // CONSIDER: Improvement of the robustness/accuracy/repeatability by | |
68 | // moving all segments to 0 longitude | |
69 | // picking latitudes closer to 0 | |
70 | // etc. | |
71 | ||
72 | template | |
73 | < | |
74 | typename FormulaPolicy = strategy::andoyer, | |
75 | unsigned int Order = strategy::default_order<FormulaPolicy>::value, | |
76 | typename Spheroid = srs::spheroid<double>, | |
77 | typename CalculationType = void | |
78 | > | |
79 | struct geographic_segments | |
80 | { | |
92f5a8d4 TL |
81 | typedef geographic_tag cs_tag; |
82 | ||
b32b8144 FG |
83 | typedef side::geographic |
84 | < | |
85 | FormulaPolicy, Spheroid, CalculationType | |
86 | > side_strategy_type; | |
87 | ||
88 | inline side_strategy_type get_side_strategy() const | |
89 | { | |
90 | return side_strategy_type(m_spheroid); | |
91 | } | |
92 | ||
93 | template <typename Geometry1, typename Geometry2> | |
94 | struct point_in_geometry_strategy | |
95 | { | |
96 | typedef strategy::within::geographic_winding | |
97 | < | |
98 | typename point_type<Geometry1>::type, | |
99 | typename point_type<Geometry2>::type, | |
100 | FormulaPolicy, | |
101 | Spheroid, | |
102 | CalculationType | |
103 | > type; | |
104 | }; | |
105 | ||
106 | template <typename Geometry1, typename Geometry2> | |
107 | inline typename point_in_geometry_strategy<Geometry1, Geometry2>::type | |
108 | get_point_in_geometry_strategy() const | |
109 | { | |
110 | typedef typename point_in_geometry_strategy | |
111 | < | |
112 | Geometry1, Geometry2 | |
113 | >::type strategy_type; | |
114 | return strategy_type(m_spheroid); | |
115 | } | |
116 | ||
117 | template <typename Geometry> | |
118 | struct area_strategy | |
119 | { | |
120 | typedef area::geographic | |
121 | < | |
b32b8144 FG |
122 | FormulaPolicy, |
123 | Order, | |
124 | Spheroid, | |
125 | CalculationType | |
126 | > type; | |
127 | }; | |
128 | ||
129 | template <typename Geometry> | |
130 | inline typename area_strategy<Geometry>::type get_area_strategy() const | |
131 | { | |
132 | typedef typename area_strategy<Geometry>::type strategy_type; | |
133 | return strategy_type(m_spheroid); | |
134 | } | |
135 | ||
136 | template <typename Geometry> | |
137 | struct distance_strategy | |
138 | { | |
139 | typedef distance::geographic | |
140 | < | |
141 | FormulaPolicy, | |
142 | Spheroid, | |
143 | CalculationType | |
144 | > type; | |
145 | }; | |
146 | ||
147 | template <typename Geometry> | |
148 | inline typename distance_strategy<Geometry>::type get_distance_strategy() const | |
149 | { | |
150 | typedef typename distance_strategy<Geometry>::type strategy_type; | |
151 | return strategy_type(m_spheroid); | |
152 | } | |
153 | ||
92f5a8d4 | 154 | typedef envelope::geographic<FormulaPolicy, Spheroid, CalculationType> |
b32b8144 FG |
155 | envelope_strategy_type; |
156 | ||
157 | inline envelope_strategy_type get_envelope_strategy() const | |
158 | { | |
159 | return envelope_strategy_type(m_spheroid); | |
160 | } | |
161 | ||
92f5a8d4 TL |
162 | typedef expand::geographic_segment<FormulaPolicy, Spheroid, CalculationType> |
163 | expand_strategy_type; | |
164 | ||
165 | inline expand_strategy_type get_expand_strategy() const | |
166 | { | |
167 | return expand_strategy_type(m_spheroid); | |
168 | } | |
169 | ||
170 | typedef within::spherical_point_point point_in_point_strategy_type; | |
171 | ||
172 | static inline point_in_point_strategy_type get_point_in_point_strategy() | |
173 | { | |
174 | return point_in_point_strategy_type(); | |
175 | } | |
176 | ||
177 | typedef within::spherical_point_point equals_point_point_strategy_type; | |
178 | ||
179 | static inline equals_point_point_strategy_type get_equals_point_point_strategy() | |
180 | { | |
181 | return equals_point_point_strategy_type(); | |
182 | } | |
183 | ||
184 | typedef disjoint::spherical_box_box disjoint_box_box_strategy_type; | |
185 | ||
186 | static inline disjoint_box_box_strategy_type get_disjoint_box_box_strategy() | |
187 | { | |
188 | return disjoint_box_box_strategy_type(); | |
189 | } | |
190 | ||
191 | typedef disjoint::segment_box_geographic | |
192 | < | |
193 | FormulaPolicy, Spheroid, CalculationType | |
194 | > disjoint_segment_box_strategy_type; | |
195 | ||
196 | inline disjoint_segment_box_strategy_type get_disjoint_segment_box_strategy() const | |
197 | { | |
198 | return disjoint_segment_box_strategy_type(m_spheroid); | |
199 | } | |
200 | ||
201 | typedef covered_by::spherical_point_box disjoint_point_box_strategy_type; | |
202 | typedef covered_by::spherical_point_box covered_by_point_box_strategy_type; | |
203 | typedef within::spherical_point_box within_point_box_strategy_type; | |
204 | typedef envelope::spherical_box envelope_box_strategy_type; | |
205 | typedef expand::spherical_box expand_box_strategy_type; | |
206 | ||
b32b8144 FG |
207 | enum intersection_point_flag { ipi_inters = 0, ipi_at_a1, ipi_at_a2, ipi_at_b1, ipi_at_b2 }; |
208 | ||
209 | template <typename CoordinateType, typename SegmentRatio> | |
210 | struct segment_intersection_info | |
211 | { | |
b32b8144 | 212 | template <typename Point, typename Segment1, typename Segment2> |
11fdf7f2 | 213 | void calculate(Point& point, Segment1 const& a, Segment2 const& b) const |
b32b8144 FG |
214 | { |
215 | if (ip_flag == ipi_inters) | |
216 | { | |
217 | // TODO: assign the rest of coordinates | |
218 | set_from_radian<0>(point, lon); | |
219 | set_from_radian<1>(point, lat); | |
220 | } | |
221 | else if (ip_flag == ipi_at_a1) | |
222 | { | |
223 | detail::assign_point_from_index<0>(a, point); | |
224 | } | |
225 | else if (ip_flag == ipi_at_a2) | |
226 | { | |
227 | detail::assign_point_from_index<1>(a, point); | |
228 | } | |
229 | else if (ip_flag == ipi_at_b1) | |
230 | { | |
231 | detail::assign_point_from_index<0>(b, point); | |
232 | } | |
233 | else // ip_flag == ipi_at_b2 | |
234 | { | |
235 | detail::assign_point_from_index<1>(b, point); | |
236 | } | |
237 | } | |
238 | ||
239 | CoordinateType lon; | |
240 | CoordinateType lat; | |
241 | SegmentRatio robust_ra; | |
242 | SegmentRatio robust_rb; | |
243 | intersection_point_flag ip_flag; | |
244 | }; | |
245 | ||
246 | explicit geographic_segments(Spheroid const& spheroid = Spheroid()) | |
247 | : m_spheroid(spheroid) | |
248 | {} | |
249 | ||
250 | // Relate segments a and b | |
251 | template | |
252 | < | |
92f5a8d4 TL |
253 | typename UniqueSubRange1, |
254 | typename UniqueSubRange2, | |
255 | typename Policy | |
b32b8144 | 256 | > |
92f5a8d4 TL |
257 | inline typename Policy::return_type apply(UniqueSubRange1 const& range_p, |
258 | UniqueSubRange2 const& range_q, | |
259 | Policy const&) const | |
b32b8144 | 260 | { |
92f5a8d4 TL |
261 | typedef typename UniqueSubRange1::point_type point1_type; |
262 | typedef typename UniqueSubRange2::point_type point2_type; | |
263 | typedef model::referring_segment<point1_type const> segment_type1; | |
264 | typedef model::referring_segment<point2_type const> segment_type2; | |
265 | ||
266 | BOOST_CONCEPT_ASSERT( (concepts::ConstPoint<point1_type>) ); | |
267 | BOOST_CONCEPT_ASSERT( (concepts::ConstPoint<point2_type>) ); | |
268 | ||
269 | /* | |
270 | typename coordinate_type<Point1>::type | |
271 | const a1_lon = get<0>(a1), | |
272 | const a2_lon = get<0>(a2); | |
273 | typename coordinate_type<Point2>::type | |
274 | const b1_lon = get<0>(b1), | |
275 | const b2_lon = get<0>(b2); | |
276 | bool is_a_reversed = a1_lon > a2_lon || a1_lon == a2_lon && get<1>(a1) > get<1>(a2); | |
277 | bool is_b_reversed = b1_lon > b2_lon || b1_lon == b2_lon && get<1>(b1) > get<1>(b2); | |
278 | */ | |
279 | ||
20effc67 TL |
280 | point1_type const& p0 = range_p.at(0); |
281 | point1_type const& p1 = range_p.at(1); | |
282 | point2_type const& q0 = range_q.at(0); | |
283 | point2_type const& q1 = range_q.at(1); | |
284 | ||
285 | bool const is_p_reversed = get<1>(p0) > get<1>(p1); | |
286 | bool const is_q_reversed = get<1>(q0) > get<1>(q1); | |
92f5a8d4 TL |
287 | |
288 | // Call apply with original segments and ordered points | |
20effc67 TL |
289 | return apply<Policy>(segment_type1(p0, p1), |
290 | segment_type2(q0, q1), | |
291 | (is_p_reversed ? p1 : p0), | |
292 | (is_p_reversed ? p0 : p1), | |
293 | (is_q_reversed ? q1 : q0), | |
294 | (is_q_reversed ? q0 : q1), | |
92f5a8d4 | 295 | is_p_reversed, is_q_reversed); |
b32b8144 FG |
296 | } |
297 | ||
298 | private: | |
299 | // Relate segments a and b | |
300 | template | |
301 | < | |
302 | typename Policy, | |
303 | typename Segment1, | |
304 | typename Segment2, | |
305 | typename Point1, | |
306 | typename Point2 | |
307 | > | |
308 | inline typename Policy::return_type apply(Segment1 const& a, Segment2 const& b, | |
309 | Point1 const& a1, Point1 const& a2, | |
310 | Point2 const& b1, Point2 const& b2, | |
311 | bool is_a_reversed, bool is_b_reversed) const | |
312 | { | |
313 | BOOST_CONCEPT_ASSERT( (concepts::ConstSegment<Segment1>) ); | |
314 | BOOST_CONCEPT_ASSERT( (concepts::ConstSegment<Segment2>) ); | |
315 | ||
316 | typedef typename select_calculation_type | |
317 | <Segment1, Segment2, CalculationType>::type calc_t; | |
318 | ||
319 | typedef srs::spheroid<calc_t> spheroid_type; | |
320 | ||
321 | static const calc_t c0 = 0; | |
322 | ||
323 | // normalized spheroid | |
324 | spheroid_type spheroid = formula::unit_spheroid<spheroid_type>(m_spheroid); | |
325 | ||
326 | // TODO: check only 2 first coordinates here? | |
b32b8144 FG |
327 | bool a_is_point = equals_point_point(a1, a2); |
328 | bool b_is_point = equals_point_point(b1, b2); | |
329 | ||
330 | if(a_is_point && b_is_point) | |
331 | { | |
332 | return equals_point_point(a1, b2) | |
333 | ? Policy::degenerate(a, true) | |
334 | : Policy::disjoint() | |
335 | ; | |
336 | } | |
337 | ||
338 | calc_t const a1_lon = get_as_radian<0>(a1); | |
339 | calc_t const a1_lat = get_as_radian<1>(a1); | |
340 | calc_t const a2_lon = get_as_radian<0>(a2); | |
341 | calc_t const a2_lat = get_as_radian<1>(a2); | |
342 | calc_t const b1_lon = get_as_radian<0>(b1); | |
343 | calc_t const b1_lat = get_as_radian<1>(b1); | |
344 | calc_t const b2_lon = get_as_radian<0>(b2); | |
345 | calc_t const b2_lat = get_as_radian<1>(b2); | |
346 | ||
347 | side_info sides; | |
348 | ||
349 | // NOTE: potential optimization, don't calculate distance at this point | |
350 | // this would require to reimplement inverse strategy to allow | |
351 | // calculation of distance if needed, probably also storing intermediate | |
352 | // results somehow inside an object. | |
353 | typedef typename FormulaPolicy::template inverse<calc_t, true, true, false, false, false> inverse_dist_azi; | |
354 | typedef typename inverse_dist_azi::result_type inverse_result; | |
355 | ||
356 | // TODO: no need to call inverse formula if we know that the points are equal | |
357 | // distance can be set to 0 in this case and azimuth may be not calculated | |
358 | bool is_equal_a1_b1 = equals_point_point(a1, b1); | |
359 | bool is_equal_a2_b1 = equals_point_point(a2, b1); | |
360 | bool degen_neq_coords = false; | |
361 | ||
362 | inverse_result res_b1_b2, res_b1_a1, res_b1_a2; | |
363 | if (! b_is_point) | |
364 | { | |
365 | res_b1_b2 = inverse_dist_azi::apply(b1_lon, b1_lat, b2_lon, b2_lat, spheroid); | |
366 | if (math::equals(res_b1_b2.distance, c0)) | |
367 | { | |
368 | b_is_point = true; | |
369 | degen_neq_coords = true; | |
370 | } | |
371 | else | |
372 | { | |
373 | res_b1_a1 = inverse_dist_azi::apply(b1_lon, b1_lat, a1_lon, a1_lat, spheroid); | |
374 | if (math::equals(res_b1_a1.distance, c0)) | |
375 | { | |
376 | is_equal_a1_b1 = true; | |
377 | } | |
378 | res_b1_a2 = inverse_dist_azi::apply(b1_lon, b1_lat, a2_lon, a2_lat, spheroid); | |
379 | if (math::equals(res_b1_a2.distance, c0)) | |
380 | { | |
381 | is_equal_a2_b1 = true; | |
382 | } | |
383 | sides.set<0>(is_equal_a1_b1 ? 0 : formula::azimuth_side_value(res_b1_a1.azimuth, res_b1_b2.azimuth), | |
384 | is_equal_a2_b1 ? 0 : formula::azimuth_side_value(res_b1_a2.azimuth, res_b1_b2.azimuth)); | |
385 | if (sides.same<0>()) | |
386 | { | |
387 | // Both points are at the same side of other segment, we can leave | |
388 | return Policy::disjoint(); | |
389 | } | |
390 | } | |
391 | } | |
392 | ||
393 | bool is_equal_a1_b2 = equals_point_point(a1, b2); | |
394 | ||
395 | inverse_result res_a1_a2, res_a1_b1, res_a1_b2; | |
396 | if (! a_is_point) | |
397 | { | |
398 | res_a1_a2 = inverse_dist_azi::apply(a1_lon, a1_lat, a2_lon, a2_lat, spheroid); | |
399 | if (math::equals(res_a1_a2.distance, c0)) | |
400 | { | |
401 | a_is_point = true; | |
402 | degen_neq_coords = true; | |
403 | } | |
404 | else | |
405 | { | |
406 | res_a1_b1 = inverse_dist_azi::apply(a1_lon, a1_lat, b1_lon, b1_lat, spheroid); | |
407 | if (math::equals(res_a1_b1.distance, c0)) | |
408 | { | |
409 | is_equal_a1_b1 = true; | |
410 | } | |
411 | res_a1_b2 = inverse_dist_azi::apply(a1_lon, a1_lat, b2_lon, b2_lat, spheroid); | |
412 | if (math::equals(res_a1_b2.distance, c0)) | |
413 | { | |
414 | is_equal_a1_b2 = true; | |
415 | } | |
416 | sides.set<1>(is_equal_a1_b1 ? 0 : formula::azimuth_side_value(res_a1_b1.azimuth, res_a1_a2.azimuth), | |
417 | is_equal_a1_b2 ? 0 : formula::azimuth_side_value(res_a1_b2.azimuth, res_a1_a2.azimuth)); | |
418 | if (sides.same<1>()) | |
419 | { | |
420 | // Both points are at the same side of other segment, we can leave | |
421 | return Policy::disjoint(); | |
422 | } | |
423 | } | |
424 | } | |
425 | ||
426 | if(a_is_point && b_is_point) | |
427 | { | |
428 | return is_equal_a1_b2 | |
429 | ? Policy::degenerate(a, true) | |
430 | : Policy::disjoint() | |
431 | ; | |
432 | } | |
433 | ||
434 | // NOTE: at this point the segments may still be disjoint | |
435 | // NOTE: at this point one of the segments may be degenerated | |
436 | ||
437 | bool collinear = sides.collinear(); | |
438 | ||
439 | if (! collinear) | |
440 | { | |
441 | // WARNING: the side strategy doesn't have the info about the other | |
442 | // segment so it may return results inconsistent with this intersection | |
443 | // strategy, as it checks both segments for consistency | |
444 | ||
445 | if (sides.get<0, 0>() == 0 && sides.get<0, 1>() == 0) | |
446 | { | |
447 | collinear = true; | |
448 | sides.set<1>(0, 0); | |
449 | } | |
450 | else if (sides.get<1, 0>() == 0 && sides.get<1, 1>() == 0) | |
451 | { | |
452 | collinear = true; | |
453 | sides.set<0>(0, 0); | |
454 | } | |
455 | } | |
456 | ||
457 | if (collinear) | |
458 | { | |
459 | if (a_is_point) | |
460 | { | |
461 | 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); | |
462 | } | |
463 | else if (b_is_point) | |
464 | { | |
465 | 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); | |
466 | } | |
467 | else | |
468 | { | |
469 | calc_t dist_a1_a2, dist_a1_b1, dist_a1_b2; | |
470 | calc_t dist_b1_b2, dist_b1_a1, dist_b1_a2; | |
471 | // use shorter segment | |
472 | if (res_a1_a2.distance <= res_b1_b2.distance) | |
473 | { | |
474 | calculate_collinear_data(a1, a2, b1, b2, res_a1_a2, res_a1_b1, res_a1_b2, dist_a1_a2, dist_a1_b1); | |
92f5a8d4 | 475 | calculate_collinear_data(a1, a2, b2, b1, res_a1_a2, res_a1_b2, res_a1_b1, dist_a1_a2, dist_a1_b2); |
b32b8144 FG |
476 | dist_b1_b2 = dist_a1_b2 - dist_a1_b1; |
477 | dist_b1_a1 = -dist_a1_b1; | |
478 | dist_b1_a2 = dist_a1_a2 - dist_a1_b1; | |
479 | } | |
480 | else | |
481 | { | |
482 | calculate_collinear_data(b1, b2, a1, a2, res_b1_b2, res_b1_a1, res_b1_a2, dist_b1_b2, dist_b1_a1); | |
92f5a8d4 | 483 | calculate_collinear_data(b1, b2, a2, a1, res_b1_b2, res_b1_a2, res_b1_a1, dist_b1_b2, dist_b1_a2); |
b32b8144 FG |
484 | dist_a1_a2 = dist_b1_a2 - dist_b1_a1; |
485 | dist_a1_b1 = -dist_b1_a1; | |
486 | dist_a1_b2 = dist_b1_b2 - dist_b1_a1; | |
487 | } | |
488 | ||
489 | // NOTE: this is probably not needed | |
b32b8144 FG |
490 | int a1_on_b = position_value(c0, dist_a1_b1, dist_a1_b2); |
491 | int a2_on_b = position_value(dist_a1_a2, dist_a1_b1, dist_a1_b2); | |
492 | int b1_on_a = position_value(c0, dist_b1_a1, dist_b1_a2); | |
493 | int b2_on_a = position_value(dist_b1_b2, dist_b1_a1, dist_b1_a2); | |
494 | ||
495 | if ((a1_on_b < 1 && a2_on_b < 1) || (a1_on_b > 3 && a2_on_b > 3)) | |
496 | { | |
497 | return Policy::disjoint(); | |
498 | } | |
499 | ||
500 | if (a1_on_b == 1) | |
501 | { | |
502 | dist_b1_a1 = 0; | |
503 | dist_a1_b1 = 0; | |
504 | } | |
505 | else if (a1_on_b == 3) | |
506 | { | |
507 | dist_b1_a1 = dist_b1_b2; | |
508 | dist_a1_b2 = 0; | |
509 | } | |
510 | ||
511 | if (a2_on_b == 1) | |
512 | { | |
513 | dist_b1_a2 = 0; | |
514 | dist_a1_b1 = dist_a1_a2; | |
515 | } | |
516 | else if (a2_on_b == 3) | |
517 | { | |
518 | dist_b1_a2 = dist_b1_b2; | |
519 | dist_a1_b2 = dist_a1_a2; | |
520 | } | |
521 | ||
522 | bool opposite = ! same_direction(res_a1_a2.azimuth, res_b1_b2.azimuth); | |
523 | ||
524 | // NOTE: If segment was reversed opposite, positions and segment ratios has to be altered | |
525 | if (is_a_reversed) | |
526 | { | |
527 | // opposite | |
528 | opposite = ! opposite; | |
529 | // positions | |
530 | std::swap(a1_on_b, a2_on_b); | |
531 | b1_on_a = 4 - b1_on_a; | |
532 | b2_on_a = 4 - b2_on_a; | |
533 | // distances for ratios | |
534 | std::swap(dist_b1_a1, dist_b1_a2); | |
535 | dist_a1_b1 = dist_a1_a2 - dist_a1_b1; | |
536 | dist_a1_b2 = dist_a1_a2 - dist_a1_b2; | |
537 | } | |
538 | if (is_b_reversed) | |
539 | { | |
540 | // opposite | |
541 | opposite = ! opposite; | |
542 | // positions | |
543 | a1_on_b = 4 - a1_on_b; | |
544 | a2_on_b = 4 - a2_on_b; | |
545 | std::swap(b1_on_a, b2_on_a); | |
546 | // distances for ratios | |
547 | dist_b1_a1 = dist_b1_b2 - dist_b1_a1; | |
548 | dist_b1_a2 = dist_b1_b2 - dist_b1_a2; | |
549 | std::swap(dist_a1_b1, dist_a1_b2); | |
550 | } | |
551 | ||
552 | segment_ratio<calc_t> ra_from(dist_b1_a1, dist_b1_b2); | |
553 | segment_ratio<calc_t> ra_to(dist_b1_a2, dist_b1_b2); | |
554 | segment_ratio<calc_t> rb_from(dist_a1_b1, dist_a1_a2); | |
555 | segment_ratio<calc_t> rb_to(dist_a1_b2, dist_a1_a2); | |
556 | ||
557 | return Policy::segments_collinear(a, b, opposite, | |
558 | a1_on_b, a2_on_b, b1_on_a, b2_on_a, | |
559 | ra_from, ra_to, rb_from, rb_to); | |
560 | } | |
561 | } | |
562 | else // crossing or touching | |
563 | { | |
564 | if (a_is_point || b_is_point) | |
565 | { | |
566 | return Policy::disjoint(); | |
567 | } | |
568 | ||
569 | calc_t lon = 0, lat = 0; | |
570 | intersection_point_flag ip_flag; | |
571 | calc_t dist_a1_a2, dist_a1_i1, dist_b1_b2, dist_b1_i1; | |
572 | if (calculate_ip_data(a1, a2, b1, b2, | |
573 | a1_lon, a1_lat, a2_lon, a2_lat, | |
574 | b1_lon, b1_lat, b2_lon, b2_lat, | |
575 | res_a1_a2, res_a1_b1, res_a1_b2, | |
576 | res_b1_b2, res_b1_a1, res_b1_a2, | |
577 | sides, spheroid, | |
578 | lon, lat, | |
579 | dist_a1_a2, dist_a1_i1, dist_b1_b2, dist_b1_i1, | |
580 | ip_flag)) | |
581 | { | |
582 | // NOTE: If segment was reversed sides and segment ratios has to be altered | |
583 | if (is_a_reversed) | |
584 | { | |
585 | // sides | |
586 | sides_reverse_segment<0>(sides); | |
587 | // distance for ratio | |
588 | dist_a1_i1 = dist_a1_a2 - dist_a1_i1; | |
589 | // ip flag | |
590 | ip_flag_reverse_segment(ip_flag, ipi_at_a1, ipi_at_a2); | |
591 | } | |
592 | if (is_b_reversed) | |
593 | { | |
594 | // sides | |
595 | sides_reverse_segment<1>(sides); | |
596 | // distance for ratio | |
597 | dist_b1_i1 = dist_b1_b2 - dist_b1_i1; | |
598 | // ip flag | |
599 | ip_flag_reverse_segment(ip_flag, ipi_at_b1, ipi_at_b2); | |
600 | } | |
601 | ||
602 | // intersects | |
603 | segment_intersection_info | |
604 | < | |
605 | calc_t, | |
606 | segment_ratio<calc_t> | |
607 | > sinfo; | |
608 | ||
609 | sinfo.lon = lon; | |
610 | sinfo.lat = lat; | |
611 | sinfo.robust_ra.assign(dist_a1_i1, dist_a1_a2); | |
612 | sinfo.robust_rb.assign(dist_b1_i1, dist_b1_b2); | |
613 | sinfo.ip_flag = ip_flag; | |
614 | ||
615 | return Policy::segments_crosses(sides, sinfo, a, b); | |
616 | } | |
617 | else | |
618 | { | |
619 | return Policy::disjoint(); | |
620 | } | |
621 | } | |
622 | } | |
623 | ||
624 | template <typename Policy, typename CalcT, typename Segment, typename Point1, typename Point2, typename ResultInverse> | |
625 | static inline typename Policy::return_type | |
626 | collinear_one_degenerated(Segment const& segment, bool degenerated_a, | |
627 | Point1 const& a1, Point1 const& a2, | |
628 | Point2 const& b1, Point2 const& b2, | |
629 | ResultInverse const& res_a1_a2, | |
630 | ResultInverse const& res_a1_b1, | |
631 | ResultInverse const& res_a1_b2, | |
632 | bool is_other_reversed, | |
633 | bool degen_neq_coords) | |
634 | { | |
635 | CalcT dist_1_2, dist_1_o; | |
636 | 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)) | |
637 | { | |
638 | return Policy::disjoint(); | |
639 | } | |
640 | ||
641 | // NOTE: If segment was reversed segment ratio has to be altered | |
642 | if (is_other_reversed) | |
643 | { | |
644 | // distance for ratio | |
645 | dist_1_o = dist_1_2 - dist_1_o; | |
646 | } | |
647 | ||
648 | return Policy::one_degenerate(segment, segment_ratio<CalcT>(dist_1_o, dist_1_2), degenerated_a); | |
649 | } | |
650 | ||
651 | // TODO: instead of checks below test bi against a1 and a2 here? | |
652 | // in order to make this independent from is_near() | |
653 | template <typename Point1, typename Point2, typename ResultInverse, typename CalcT> | |
654 | static inline bool calculate_collinear_data(Point1 const& a1, Point1 const& a2, // in | |
92f5a8d4 | 655 | Point2 const& b1, Point2 const& /*b2*/, // in |
b32b8144 FG |
656 | ResultInverse const& res_a1_a2, // in |
657 | ResultInverse const& res_a1_b1, // in | |
658 | ResultInverse const& res_a1_b2, // in | |
659 | CalcT& dist_a1_a2, // out | |
92f5a8d4 | 660 | CalcT& dist_a1_b1, // out |
b32b8144 FG |
661 | bool degen_neq_coords = false) // in |
662 | { | |
663 | dist_a1_a2 = res_a1_a2.distance; | |
664 | ||
92f5a8d4 | 665 | dist_a1_b1 = res_a1_b1.distance; |
b32b8144 FG |
666 | if (! same_direction(res_a1_b1.azimuth, res_a1_a2.azimuth)) |
667 | { | |
92f5a8d4 | 668 | dist_a1_b1 = -dist_a1_b1; |
b32b8144 FG |
669 | } |
670 | ||
92f5a8d4 TL |
671 | // if b1 is close a1 |
672 | if (is_endpoint_equal(dist_a1_b1, a1, b1)) | |
b32b8144 | 673 | { |
92f5a8d4 | 674 | dist_a1_b1 = 0; |
b32b8144 FG |
675 | return true; |
676 | } | |
92f5a8d4 TL |
677 | // if b1 is close a2 |
678 | else if (is_endpoint_equal(dist_a1_a2 - dist_a1_b1, a2, b1)) | |
b32b8144 | 679 | { |
92f5a8d4 | 680 | dist_a1_b1 = dist_a1_a2; |
b32b8144 FG |
681 | return true; |
682 | } | |
683 | ||
92f5a8d4 | 684 | // check the other endpoint of degenerated segment near a pole |
b32b8144 FG |
685 | if (degen_neq_coords) |
686 | { | |
687 | static CalcT const c0 = 0; | |
688 | if (math::equals(res_a1_b2.distance, c0)) | |
689 | { | |
92f5a8d4 | 690 | dist_a1_b1 = 0; |
b32b8144 FG |
691 | return true; |
692 | } | |
693 | else if (math::equals(dist_a1_a2 - res_a1_b2.distance, c0)) | |
694 | { | |
92f5a8d4 | 695 | dist_a1_b1 = dist_a1_a2; |
b32b8144 FG |
696 | return true; |
697 | } | |
698 | } | |
699 | ||
700 | // or i1 is on b | |
92f5a8d4 | 701 | return segment_ratio<CalcT>(dist_a1_b1, dist_a1_a2).on_segment(); |
b32b8144 FG |
702 | } |
703 | ||
704 | template <typename Point1, typename Point2, typename CalcT, typename ResultInverse, typename Spheroid_> | |
705 | static inline bool calculate_ip_data(Point1 const& a1, Point1 const& a2, // in | |
706 | Point2 const& b1, Point2 const& b2, // in | |
707 | CalcT const& a1_lon, CalcT const& a1_lat, // in | |
708 | CalcT const& a2_lon, CalcT const& a2_lat, // in | |
709 | CalcT const& b1_lon, CalcT const& b1_lat, // in | |
710 | CalcT const& b2_lon, CalcT const& b2_lat, // in | |
711 | ResultInverse const& res_a1_a2, // in | |
712 | ResultInverse const& res_a1_b1, // in | |
713 | ResultInverse const& res_a1_b2, // in | |
714 | ResultInverse const& res_b1_b2, // in | |
715 | ResultInverse const& res_b1_a1, // in | |
716 | ResultInverse const& res_b1_a2, // in | |
717 | side_info const& sides, // in | |
718 | Spheroid_ const& spheroid, // in | |
719 | CalcT & lon, CalcT & lat, // out | |
720 | CalcT& dist_a1_a2, CalcT& dist_a1_ip, // out | |
721 | CalcT& dist_b1_b2, CalcT& dist_b1_ip, // out | |
722 | intersection_point_flag& ip_flag) // out | |
723 | { | |
724 | dist_a1_a2 = res_a1_a2.distance; | |
725 | dist_b1_b2 = res_b1_b2.distance; | |
726 | ||
727 | // assign the IP if some endpoints overlap | |
b32b8144 FG |
728 | if (equals_point_point(a1, b1)) |
729 | { | |
730 | lon = a1_lon; | |
731 | lat = a1_lat; | |
732 | dist_a1_ip = 0; | |
733 | dist_b1_ip = 0; | |
734 | ip_flag = ipi_at_a1; | |
735 | return true; | |
736 | } | |
737 | else if (equals_point_point(a1, b2)) | |
738 | { | |
739 | lon = a1_lon; | |
740 | lat = a1_lat; | |
741 | dist_a1_ip = 0; | |
742 | dist_b1_ip = dist_b1_b2; | |
743 | ip_flag = ipi_at_a1; | |
744 | return true; | |
745 | } | |
746 | else if (equals_point_point(a2, b1)) | |
747 | { | |
748 | lon = a2_lon; | |
749 | lat = a2_lat; | |
750 | dist_a1_ip = dist_a1_a2; | |
751 | dist_b1_ip = 0; | |
752 | ip_flag = ipi_at_a2; | |
753 | return true; | |
754 | } | |
755 | else if (equals_point_point(a2, b2)) | |
756 | { | |
757 | lon = a2_lon; | |
758 | lat = a2_lat; | |
759 | dist_a1_ip = dist_a1_a2; | |
760 | dist_b1_ip = dist_b1_b2; | |
761 | ip_flag = ipi_at_a2; | |
762 | return true; | |
763 | } | |
764 | ||
765 | // at this point we know that the endpoints doesn't overlap | |
766 | // check cases when an endpoint lies on the other geodesic | |
767 | if (sides.template get<0, 0>() == 0) // a1 wrt b | |
768 | { | |
769 | if (res_b1_a1.distance <= res_b1_b2.distance | |
770 | && same_direction(res_b1_a1.azimuth, res_b1_b2.azimuth)) | |
771 | { | |
772 | lon = a1_lon; | |
773 | lat = a1_lat; | |
774 | dist_a1_ip = 0; | |
775 | dist_b1_ip = res_b1_a1.distance; | |
776 | ip_flag = ipi_at_a1; | |
777 | return true; | |
778 | } | |
779 | else | |
780 | { | |
781 | return false; | |
782 | } | |
783 | } | |
784 | else if (sides.template get<0, 1>() == 0) // a2 wrt b | |
785 | { | |
786 | if (res_b1_a2.distance <= res_b1_b2.distance | |
787 | && same_direction(res_b1_a2.azimuth, res_b1_b2.azimuth)) | |
788 | { | |
789 | lon = a2_lon; | |
790 | lat = a2_lat; | |
791 | dist_a1_ip = res_a1_a2.distance; | |
792 | dist_b1_ip = res_b1_a2.distance; | |
793 | ip_flag = ipi_at_a2; | |
794 | return true; | |
795 | } | |
796 | else | |
797 | { | |
798 | return false; | |
799 | } | |
800 | } | |
801 | else if (sides.template get<1, 0>() == 0) // b1 wrt a | |
802 | { | |
803 | if (res_a1_b1.distance <= res_a1_a2.distance | |
804 | && same_direction(res_a1_b1.azimuth, res_a1_a2.azimuth)) | |
805 | { | |
806 | lon = b1_lon; | |
807 | lat = b1_lat; | |
808 | dist_a1_ip = res_a1_b1.distance; | |
809 | dist_b1_ip = 0; | |
810 | ip_flag = ipi_at_b1; | |
811 | return true; | |
812 | } | |
813 | else | |
814 | { | |
815 | return false; | |
816 | } | |
817 | } | |
818 | else if (sides.template get<1, 1>() == 0) // b2 wrt a | |
819 | { | |
820 | if (res_a1_b2.distance <= res_a1_a2.distance | |
821 | && same_direction(res_a1_b2.azimuth, res_a1_a2.azimuth)) | |
822 | { | |
823 | lon = b2_lon; | |
824 | lat = b2_lat; | |
825 | dist_a1_ip = res_a1_b2.distance; | |
826 | dist_b1_ip = res_b1_b2.distance; | |
827 | ip_flag = ipi_at_b2; | |
828 | return true; | |
829 | } | |
830 | else | |
831 | { | |
832 | return false; | |
833 | } | |
834 | } | |
835 | ||
836 | // At this point neither the endpoints overlaps | |
837 | // nor any andpoint lies on the other geodesic | |
838 | // So the endpoints should lie on the opposite sides of both geodesics | |
839 | ||
840 | bool const ok = formula::sjoberg_intersection<CalcT, FormulaPolicy::template inverse, Order> | |
841 | ::apply(a1_lon, a1_lat, a2_lon, a2_lat, res_a1_a2.azimuth, | |
842 | b1_lon, b1_lat, b2_lon, b2_lat, res_b1_b2.azimuth, | |
843 | lon, lat, spheroid); | |
844 | ||
845 | if (! ok) | |
846 | { | |
847 | return false; | |
848 | } | |
849 | ||
850 | typedef typename FormulaPolicy::template inverse<CalcT, true, true, false, false, false> inverse_dist_azi; | |
851 | typedef typename inverse_dist_azi::result_type inverse_result; | |
852 | ||
853 | inverse_result const res_a1_ip = inverse_dist_azi::apply(a1_lon, a1_lat, lon, lat, spheroid); | |
854 | dist_a1_ip = res_a1_ip.distance; | |
855 | if (! same_direction(res_a1_ip.azimuth, res_a1_a2.azimuth)) | |
856 | { | |
857 | dist_a1_ip = -dist_a1_ip; | |
858 | } | |
859 | ||
860 | bool is_on_a = segment_ratio<CalcT>(dist_a1_ip, dist_a1_a2).on_segment(); | |
861 | // NOTE: not fully consistent with equals_point_point() since radians are always used. | |
862 | bool is_on_a1 = math::equals(lon, a1_lon) && math::equals(lat, a1_lat); | |
863 | bool is_on_a2 = math::equals(lon, a2_lon) && math::equals(lat, a2_lat); | |
864 | ||
865 | if (! (is_on_a || is_on_a1 || is_on_a2)) | |
866 | { | |
867 | return false; | |
868 | } | |
869 | ||
870 | inverse_result const res_b1_ip = inverse_dist_azi::apply(b1_lon, b1_lat, lon, lat, spheroid); | |
871 | dist_b1_ip = res_b1_ip.distance; | |
872 | if (! same_direction(res_b1_ip.azimuth, res_b1_b2.azimuth)) | |
873 | { | |
874 | dist_b1_ip = -dist_b1_ip; | |
875 | } | |
876 | ||
877 | bool is_on_b = segment_ratio<CalcT>(dist_b1_ip, dist_b1_b2).on_segment(); | |
878 | // NOTE: not fully consistent with equals_point_point() since radians are always used. | |
879 | bool is_on_b1 = math::equals(lon, b1_lon) && math::equals(lat, b1_lat); | |
880 | bool is_on_b2 = math::equals(lon, b2_lon) && math::equals(lat, b2_lat); | |
881 | ||
882 | if (! (is_on_b || is_on_b1 || is_on_b2)) | |
883 | { | |
884 | return false; | |
885 | } | |
886 | ||
92f5a8d4 TL |
887 | typedef typename FormulaPolicy::template inverse<CalcT, true, false, false, false, false> inverse_dist; |
888 | ||
b32b8144 FG |
889 | ip_flag = ipi_inters; |
890 | ||
891 | if (is_on_b1) | |
892 | { | |
893 | lon = b1_lon; | |
894 | lat = b1_lat; | |
92f5a8d4 | 895 | dist_a1_ip = inverse_dist::apply(a1_lon, a1_lat, lon, lat, spheroid).distance; // for consistency |
b32b8144 FG |
896 | dist_b1_ip = 0; |
897 | ip_flag = ipi_at_b1; | |
898 | } | |
899 | else if (is_on_b2) | |
900 | { | |
901 | lon = b2_lon; | |
902 | lat = b2_lat; | |
92f5a8d4 | 903 | dist_a1_ip = inverse_dist::apply(a1_lon, a1_lat, lon, lat, spheroid).distance; // for consistency |
b32b8144 FG |
904 | dist_b1_ip = res_b1_b2.distance; |
905 | ip_flag = ipi_at_b2; | |
906 | } | |
907 | ||
908 | if (is_on_a1) | |
909 | { | |
910 | lon = a1_lon; | |
911 | lat = a1_lat; | |
912 | dist_a1_ip = 0; | |
92f5a8d4 | 913 | dist_b1_ip = inverse_dist::apply(b1_lon, b1_lat, lon, lat, spheroid).distance; // for consistency |
b32b8144 FG |
914 | ip_flag = ipi_at_a1; |
915 | } | |
916 | else if (is_on_a2) | |
917 | { | |
918 | lon = a2_lon; | |
919 | lat = a2_lat; | |
920 | dist_a1_ip = res_a1_a2.distance; | |
92f5a8d4 | 921 | dist_b1_ip = inverse_dist::apply(b1_lon, b1_lat, lon, lat, spheroid).distance; // for consistency |
b32b8144 FG |
922 | ip_flag = ipi_at_a2; |
923 | } | |
924 | ||
925 | return true; | |
926 | } | |
927 | ||
928 | template <typename CalcT, typename P1, typename P2> | |
929 | static inline bool is_endpoint_equal(CalcT const& dist, | |
92f5a8d4 | 930 | P1 const& ai, P2 const& b1) |
b32b8144 FG |
931 | { |
932 | static CalcT const c0 = 0; | |
92f5a8d4 | 933 | return is_near(dist) && (math::equals(dist, c0) || equals_point_point(ai, b1)); |
b32b8144 FG |
934 | } |
935 | ||
936 | template <typename CalcT> | |
937 | static inline bool is_near(CalcT const& dist) | |
938 | { | |
939 | // NOTE: This strongly depends on the Inverse method | |
20effc67 | 940 | CalcT const small_number = CalcT(std::is_same<CalcT, float>::value ? 0.0001 : 0.00000001); |
b32b8144 FG |
941 | return math::abs(dist) <= small_number; |
942 | } | |
943 | ||
944 | template <typename ProjCoord1, typename ProjCoord2> | |
945 | static inline int position_value(ProjCoord1 const& ca1, | |
946 | ProjCoord2 const& cb1, | |
947 | ProjCoord2 const& cb2) | |
948 | { | |
949 | // S1x 0 1 2 3 4 | |
950 | // S2 |----------> | |
951 | return math::equals(ca1, cb1) ? 1 | |
952 | : math::equals(ca1, cb2) ? 3 | |
953 | : cb1 < cb2 ? | |
954 | ( ca1 < cb1 ? 0 | |
955 | : ca1 > cb2 ? 4 | |
956 | : 2 ) | |
957 | : ( ca1 > cb1 ? 0 | |
958 | : ca1 < cb2 ? 4 | |
959 | : 2 ); | |
960 | } | |
961 | ||
962 | template <typename CalcT> | |
963 | static inline bool same_direction(CalcT const& azimuth1, CalcT const& azimuth2) | |
964 | { | |
965 | // distance between two angles normalized to (-180, 180] | |
966 | CalcT const angle_diff = math::longitude_distance_signed<radian>(azimuth1, azimuth2); | |
967 | return math::abs(angle_diff) <= math::half_pi<CalcT>(); | |
968 | } | |
969 | ||
970 | template <int Which> | |
971 | static inline void sides_reverse_segment(side_info & sides) | |
972 | { | |
973 | // names assuming segment A is reversed (Which == 0) | |
974 | int a1_wrt_b = sides.template get<Which, 0>(); | |
975 | int a2_wrt_b = sides.template get<Which, 1>(); | |
976 | std::swap(a1_wrt_b, a2_wrt_b); | |
977 | sides.template set<Which>(a1_wrt_b, a2_wrt_b); | |
978 | int b1_wrt_a = sides.template get<1 - Which, 0>(); | |
979 | int b2_wrt_a = sides.template get<1 - Which, 1>(); | |
980 | sides.template set<1 - Which>(-b1_wrt_a, -b2_wrt_a); | |
981 | } | |
982 | ||
983 | static inline void ip_flag_reverse_segment(intersection_point_flag & ip_flag, | |
984 | intersection_point_flag const& ipi_at_p1, | |
985 | intersection_point_flag const& ipi_at_p2) | |
986 | { | |
987 | ip_flag = ip_flag == ipi_at_p1 ? ipi_at_p2 : | |
988 | ip_flag == ipi_at_p2 ? ipi_at_p1 : | |
989 | ip_flag; | |
990 | } | |
991 | ||
92f5a8d4 TL |
992 | template <typename Point1, typename Point2> |
993 | static inline bool equals_point_point(Point1 const& point1, Point2 const& point2) | |
994 | { | |
995 | return detail::equals::equals_point_point(point1, point2, | |
996 | point_in_point_strategy_type()); | |
997 | } | |
998 | ||
b32b8144 FG |
999 | private: |
1000 | Spheroid m_spheroid; | |
1001 | }; | |
1002 | ||
1003 | ||
1004 | }} // namespace strategy::intersection | |
1005 | ||
1006 | }} // namespace boost::geometry | |
1007 | ||
1008 | ||
1009 | #endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_INTERSECTION_HPP |