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