]>
Commit | Line | Data |
---|---|---|
b32b8144 FG |
1 | // Boost.Geometry |
2 | ||
11fdf7f2 TL |
3 | // Copyright (c) 2017 Adam Wulkiewicz, Lodz, Poland. |
4 | ||
92f5a8d4 | 5 | // Copyright (c) 2016-2019, 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> | |
16 | ||
17 | #include <boost/geometry/core/cs.hpp> | |
18 | #include <boost/geometry/core/access.hpp> | |
19 | #include <boost/geometry/core/radian_access.hpp> | |
b32b8144 FG |
20 | #include <boost/geometry/core/tags.hpp> |
21 | ||
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> | |
26 | ||
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> | |
31 | ||
32 | #include <boost/geometry/geometries/concepts/point_concept.hpp> | |
33 | #include <boost/geometry/geometries/concepts/segment_concept.hpp> | |
34 | ||
35 | #include <boost/geometry/policies/robustness/segment_ratio.hpp> | |
36 | ||
11fdf7f2 TL |
37 | #include <boost/geometry/srs/spheroid.hpp> |
38 | ||
b32b8144 | 39 | #include <boost/geometry/strategies/geographic/area.hpp> |
92f5a8d4 | 40 | #include <boost/geometry/strategies/geographic/disjoint_segment_box.hpp> |
b32b8144 | 41 | #include <boost/geometry/strategies/geographic/distance.hpp> |
92f5a8d4 | 42 | #include <boost/geometry/strategies/geographic/envelope.hpp> |
b32b8144 FG |
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> | |
92f5a8d4 TL |
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> | |
b32b8144 FG |
49 | #include <boost/geometry/strategies/intersection.hpp> |
50 | #include <boost/geometry/strategies/intersection_result.hpp> | |
51 | #include <boost/geometry/strategies/side_info.hpp> | |
52 | ||
53 | #include <boost/geometry/util/math.hpp> | |
54 | #include <boost/geometry/util/select_calculation_type.hpp> | |
55 | ||
56 | ||
57 | namespace boost { namespace geometry | |
58 | { | |
59 | ||
60 | namespace strategy { namespace intersection | |
61 | { | |
62 | ||
63 | // CONSIDER: Improvement of the robustness/accuracy/repeatability by | |
64 | // moving all segments to 0 longitude | |
65 | // picking latitudes closer to 0 | |
66 | // etc. | |
67 | ||
68 | template | |
69 | < | |
70 | typename FormulaPolicy = strategy::andoyer, | |
71 | unsigned int Order = strategy::default_order<FormulaPolicy>::value, | |
72 | typename Spheroid = srs::spheroid<double>, | |
73 | typename CalculationType = void | |
74 | > | |
75 | struct geographic_segments | |
76 | { | |
92f5a8d4 TL |
77 | typedef geographic_tag cs_tag; |
78 | ||
b32b8144 FG |
79 | typedef side::geographic |
80 | < | |
81 | FormulaPolicy, Spheroid, CalculationType | |
82 | > side_strategy_type; | |
83 | ||
84 | inline side_strategy_type get_side_strategy() const | |
85 | { | |
86 | return side_strategy_type(m_spheroid); | |
87 | } | |
88 | ||
89 | template <typename Geometry1, typename Geometry2> | |
90 | struct point_in_geometry_strategy | |
91 | { | |
92 | typedef strategy::within::geographic_winding | |
93 | < | |
94 | typename point_type<Geometry1>::type, | |
95 | typename point_type<Geometry2>::type, | |
96 | FormulaPolicy, | |
97 | Spheroid, | |
98 | CalculationType | |
99 | > type; | |
100 | }; | |
101 | ||
102 | template <typename Geometry1, typename Geometry2> | |
103 | inline typename point_in_geometry_strategy<Geometry1, Geometry2>::type | |
104 | get_point_in_geometry_strategy() const | |
105 | { | |
106 | typedef typename point_in_geometry_strategy | |
107 | < | |
108 | Geometry1, Geometry2 | |
109 | >::type strategy_type; | |
110 | return strategy_type(m_spheroid); | |
111 | } | |
112 | ||
113 | template <typename Geometry> | |
114 | struct area_strategy | |
115 | { | |
116 | typedef area::geographic | |
117 | < | |
b32b8144 FG |
118 | FormulaPolicy, |
119 | Order, | |
120 | Spheroid, | |
121 | CalculationType | |
122 | > type; | |
123 | }; | |
124 | ||
125 | template <typename Geometry> | |
126 | inline typename area_strategy<Geometry>::type get_area_strategy() const | |
127 | { | |
128 | typedef typename area_strategy<Geometry>::type strategy_type; | |
129 | return strategy_type(m_spheroid); | |
130 | } | |
131 | ||
132 | template <typename Geometry> | |
133 | struct distance_strategy | |
134 | { | |
135 | typedef distance::geographic | |
136 | < | |
137 | FormulaPolicy, | |
138 | Spheroid, | |
139 | CalculationType | |
140 | > type; | |
141 | }; | |
142 | ||
143 | template <typename Geometry> | |
144 | inline typename distance_strategy<Geometry>::type get_distance_strategy() const | |
145 | { | |
146 | typedef typename distance_strategy<Geometry>::type strategy_type; | |
147 | return strategy_type(m_spheroid); | |
148 | } | |
149 | ||
92f5a8d4 | 150 | typedef envelope::geographic<FormulaPolicy, Spheroid, CalculationType> |
b32b8144 FG |
151 | envelope_strategy_type; |
152 | ||
153 | inline envelope_strategy_type get_envelope_strategy() const | |
154 | { | |
155 | return envelope_strategy_type(m_spheroid); | |
156 | } | |
157 | ||
92f5a8d4 TL |
158 | typedef expand::geographic_segment<FormulaPolicy, Spheroid, CalculationType> |
159 | expand_strategy_type; | |
160 | ||
161 | inline expand_strategy_type get_expand_strategy() const | |
162 | { | |
163 | return expand_strategy_type(m_spheroid); | |
164 | } | |
165 | ||
166 | typedef within::spherical_point_point point_in_point_strategy_type; | |
167 | ||
168 | static inline point_in_point_strategy_type get_point_in_point_strategy() | |
169 | { | |
170 | return point_in_point_strategy_type(); | |
171 | } | |
172 | ||
173 | typedef within::spherical_point_point equals_point_point_strategy_type; | |
174 | ||
175 | static inline equals_point_point_strategy_type get_equals_point_point_strategy() | |
176 | { | |
177 | return equals_point_point_strategy_type(); | |
178 | } | |
179 | ||
180 | typedef disjoint::spherical_box_box disjoint_box_box_strategy_type; | |
181 | ||
182 | static inline disjoint_box_box_strategy_type get_disjoint_box_box_strategy() | |
183 | { | |
184 | return disjoint_box_box_strategy_type(); | |
185 | } | |
186 | ||
187 | typedef disjoint::segment_box_geographic | |
188 | < | |
189 | FormulaPolicy, Spheroid, CalculationType | |
190 | > disjoint_segment_box_strategy_type; | |
191 | ||
192 | inline disjoint_segment_box_strategy_type get_disjoint_segment_box_strategy() const | |
193 | { | |
194 | return disjoint_segment_box_strategy_type(m_spheroid); | |
195 | } | |
196 | ||
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; | |
202 | ||
b32b8144 FG |
203 | enum intersection_point_flag { ipi_inters = 0, ipi_at_a1, ipi_at_a2, ipi_at_b1, ipi_at_b2 }; |
204 | ||
205 | template <typename CoordinateType, typename SegmentRatio> | |
206 | struct segment_intersection_info | |
207 | { | |
b32b8144 | 208 | template <typename Point, typename Segment1, typename Segment2> |
11fdf7f2 | 209 | void calculate(Point& point, Segment1 const& a, Segment2 const& b) const |
b32b8144 FG |
210 | { |
211 | if (ip_flag == ipi_inters) | |
212 | { | |
213 | // TODO: assign the rest of coordinates | |
214 | set_from_radian<0>(point, lon); | |
215 | set_from_radian<1>(point, lat); | |
216 | } | |
217 | else if (ip_flag == ipi_at_a1) | |
218 | { | |
219 | detail::assign_point_from_index<0>(a, point); | |
220 | } | |
221 | else if (ip_flag == ipi_at_a2) | |
222 | { | |
223 | detail::assign_point_from_index<1>(a, point); | |
224 | } | |
225 | else if (ip_flag == ipi_at_b1) | |
226 | { | |
227 | detail::assign_point_from_index<0>(b, point); | |
228 | } | |
229 | else // ip_flag == ipi_at_b2 | |
230 | { | |
231 | detail::assign_point_from_index<1>(b, point); | |
232 | } | |
233 | } | |
234 | ||
235 | CoordinateType lon; | |
236 | CoordinateType lat; | |
237 | SegmentRatio robust_ra; | |
238 | SegmentRatio robust_rb; | |
239 | intersection_point_flag ip_flag; | |
240 | }; | |
241 | ||
242 | explicit geographic_segments(Spheroid const& spheroid = Spheroid()) | |
243 | : m_spheroid(spheroid) | |
244 | {} | |
245 | ||
246 | // Relate segments a and b | |
247 | template | |
248 | < | |
92f5a8d4 TL |
249 | typename UniqueSubRange1, |
250 | typename UniqueSubRange2, | |
251 | typename Policy | |
b32b8144 | 252 | > |
92f5a8d4 TL |
253 | inline typename Policy::return_type apply(UniqueSubRange1 const& range_p, |
254 | UniqueSubRange2 const& range_q, | |
255 | Policy const&) const | |
b32b8144 | 256 | { |
92f5a8d4 TL |
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; | |
261 | ||
262 | BOOST_CONCEPT_ASSERT( (concepts::ConstPoint<point1_type>) ); | |
263 | BOOST_CONCEPT_ASSERT( (concepts::ConstPoint<point2_type>) ); | |
264 | ||
265 | /* | |
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); | |
274 | */ | |
275 | ||
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)); | |
278 | ||
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); | |
b32b8144 FG |
287 | } |
288 | ||
289 | private: | |
290 | // Relate segments a and b | |
291 | template | |
292 | < | |
293 | typename Policy, | |
294 | typename Segment1, | |
295 | typename Segment2, | |
296 | typename Point1, | |
297 | typename Point2 | |
298 | > | |
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 | |
303 | { | |
304 | BOOST_CONCEPT_ASSERT( (concepts::ConstSegment<Segment1>) ); | |
305 | BOOST_CONCEPT_ASSERT( (concepts::ConstSegment<Segment2>) ); | |
306 | ||
307 | typedef typename select_calculation_type | |
308 | <Segment1, Segment2, CalculationType>::type calc_t; | |
309 | ||
310 | typedef srs::spheroid<calc_t> spheroid_type; | |
311 | ||
312 | static const calc_t c0 = 0; | |
313 | ||
314 | // normalized spheroid | |
315 | spheroid_type spheroid = formula::unit_spheroid<spheroid_type>(m_spheroid); | |
316 | ||
317 | // TODO: check only 2 first coordinates here? | |
b32b8144 FG |
318 | bool a_is_point = equals_point_point(a1, a2); |
319 | bool b_is_point = equals_point_point(b1, b2); | |
320 | ||
321 | if(a_is_point && b_is_point) | |
322 | { | |
323 | return equals_point_point(a1, b2) | |
324 | ? Policy::degenerate(a, true) | |
325 | : Policy::disjoint() | |
326 | ; | |
327 | } | |
328 | ||
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); | |
337 | ||
338 | side_info sides; | |
339 | ||
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; | |
346 | ||
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; | |
352 | ||
353 | inverse_result res_b1_b2, res_b1_a1, res_b1_a2; | |
354 | if (! b_is_point) | |
355 | { | |
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)) | |
358 | { | |
359 | b_is_point = true; | |
360 | degen_neq_coords = true; | |
361 | } | |
362 | else | |
363 | { | |
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)) | |
366 | { | |
367 | is_equal_a1_b1 = true; | |
368 | } | |
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)) | |
371 | { | |
372 | is_equal_a2_b1 = true; | |
373 | } | |
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)); | |
376 | if (sides.same<0>()) | |
377 | { | |
378 | // Both points are at the same side of other segment, we can leave | |
379 | return Policy::disjoint(); | |
380 | } | |
381 | } | |
382 | } | |
383 | ||
384 | bool is_equal_a1_b2 = equals_point_point(a1, b2); | |
385 | ||
386 | inverse_result res_a1_a2, res_a1_b1, res_a1_b2; | |
387 | if (! a_is_point) | |
388 | { | |
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)) | |
391 | { | |
392 | a_is_point = true; | |
393 | degen_neq_coords = true; | |
394 | } | |
395 | else | |
396 | { | |
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)) | |
399 | { | |
400 | is_equal_a1_b1 = true; | |
401 | } | |
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)) | |
404 | { | |
405 | is_equal_a1_b2 = true; | |
406 | } | |
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)); | |
409 | if (sides.same<1>()) | |
410 | { | |
411 | // Both points are at the same side of other segment, we can leave | |
412 | return Policy::disjoint(); | |
413 | } | |
414 | } | |
415 | } | |
416 | ||
417 | if(a_is_point && b_is_point) | |
418 | { | |
419 | return is_equal_a1_b2 | |
420 | ? Policy::degenerate(a, true) | |
421 | : Policy::disjoint() | |
422 | ; | |
423 | } | |
424 | ||
425 | // NOTE: at this point the segments may still be disjoint | |
426 | // NOTE: at this point one of the segments may be degenerated | |
427 | ||
428 | bool collinear = sides.collinear(); | |
429 | ||
430 | if (! collinear) | |
431 | { | |
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 | |
435 | ||
436 | if (sides.get<0, 0>() == 0 && sides.get<0, 1>() == 0) | |
437 | { | |
438 | collinear = true; | |
439 | sides.set<1>(0, 0); | |
440 | } | |
441 | else if (sides.get<1, 0>() == 0 && sides.get<1, 1>() == 0) | |
442 | { | |
443 | collinear = true; | |
444 | sides.set<0>(0, 0); | |
445 | } | |
446 | } | |
447 | ||
448 | if (collinear) | |
449 | { | |
450 | if (a_is_point) | |
451 | { | |
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); | |
453 | } | |
454 | else if (b_is_point) | |
455 | { | |
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); | |
457 | } | |
458 | else | |
459 | { | |
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) | |
464 | { | |
465 | calculate_collinear_data(a1, a2, b1, b2, res_a1_a2, res_a1_b1, res_a1_b2, dist_a1_a2, dist_a1_b1); | |
92f5a8d4 | 466 | calculate_collinear_data(a1, a2, b2, b1, res_a1_a2, res_a1_b2, res_a1_b1, dist_a1_a2, dist_a1_b2); |
b32b8144 FG |
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; | |
470 | } | |
471 | else | |
472 | { | |
473 | calculate_collinear_data(b1, b2, a1, a2, res_b1_b2, res_b1_a1, res_b1_a2, dist_b1_b2, dist_b1_a1); | |
92f5a8d4 | 474 | calculate_collinear_data(b1, b2, a2, a1, res_b1_b2, res_b1_a2, res_b1_a1, dist_b1_b2, dist_b1_a2); |
b32b8144 FG |
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; | |
478 | } | |
479 | ||
480 | // NOTE: this is probably not needed | |
b32b8144 FG |
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); | |
485 | ||
486 | if ((a1_on_b < 1 && a2_on_b < 1) || (a1_on_b > 3 && a2_on_b > 3)) | |
487 | { | |
488 | return Policy::disjoint(); | |
489 | } | |
490 | ||
491 | if (a1_on_b == 1) | |
492 | { | |
493 | dist_b1_a1 = 0; | |
494 | dist_a1_b1 = 0; | |
495 | } | |
496 | else if (a1_on_b == 3) | |
497 | { | |
498 | dist_b1_a1 = dist_b1_b2; | |
499 | dist_a1_b2 = 0; | |
500 | } | |
501 | ||
502 | if (a2_on_b == 1) | |
503 | { | |
504 | dist_b1_a2 = 0; | |
505 | dist_a1_b1 = dist_a1_a2; | |
506 | } | |
507 | else if (a2_on_b == 3) | |
508 | { | |
509 | dist_b1_a2 = dist_b1_b2; | |
510 | dist_a1_b2 = dist_a1_a2; | |
511 | } | |
512 | ||
513 | bool opposite = ! same_direction(res_a1_a2.azimuth, res_b1_b2.azimuth); | |
514 | ||
515 | // NOTE: If segment was reversed opposite, positions and segment ratios has to be altered | |
516 | if (is_a_reversed) | |
517 | { | |
518 | // opposite | |
519 | opposite = ! opposite; | |
520 | // positions | |
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; | |
528 | } | |
529 | if (is_b_reversed) | |
530 | { | |
531 | // opposite | |
532 | opposite = ! opposite; | |
533 | // positions | |
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); | |
541 | } | |
542 | ||
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); | |
547 | ||
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); | |
551 | } | |
552 | } | |
553 | else // crossing or touching | |
554 | { | |
555 | if (a_is_point || b_is_point) | |
556 | { | |
557 | return Policy::disjoint(); | |
558 | } | |
559 | ||
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, | |
568 | sides, spheroid, | |
569 | lon, lat, | |
570 | dist_a1_a2, dist_a1_i1, dist_b1_b2, dist_b1_i1, | |
571 | ip_flag)) | |
572 | { | |
573 | // NOTE: If segment was reversed sides and segment ratios has to be altered | |
574 | if (is_a_reversed) | |
575 | { | |
576 | // sides | |
577 | sides_reverse_segment<0>(sides); | |
578 | // distance for ratio | |
579 | dist_a1_i1 = dist_a1_a2 - dist_a1_i1; | |
580 | // ip flag | |
581 | ip_flag_reverse_segment(ip_flag, ipi_at_a1, ipi_at_a2); | |
582 | } | |
583 | if (is_b_reversed) | |
584 | { | |
585 | // sides | |
586 | sides_reverse_segment<1>(sides); | |
587 | // distance for ratio | |
588 | dist_b1_i1 = dist_b1_b2 - dist_b1_i1; | |
589 | // ip flag | |
590 | ip_flag_reverse_segment(ip_flag, ipi_at_b1, ipi_at_b2); | |
591 | } | |
592 | ||
593 | // intersects | |
594 | segment_intersection_info | |
595 | < | |
596 | calc_t, | |
597 | segment_ratio<calc_t> | |
598 | > sinfo; | |
599 | ||
600 | sinfo.lon = lon; | |
601 | sinfo.lat = lat; | |
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; | |
605 | ||
606 | return Policy::segments_crosses(sides, sinfo, a, b); | |
607 | } | |
608 | else | |
609 | { | |
610 | return Policy::disjoint(); | |
611 | } | |
612 | } | |
613 | } | |
614 | ||
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) | |
625 | { | |
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)) | |
628 | { | |
629 | return Policy::disjoint(); | |
630 | } | |
631 | ||
632 | // NOTE: If segment was reversed segment ratio has to be altered | |
633 | if (is_other_reversed) | |
634 | { | |
635 | // distance for ratio | |
636 | dist_1_o = dist_1_2 - dist_1_o; | |
637 | } | |
638 | ||
639 | return Policy::one_degenerate(segment, segment_ratio<CalcT>(dist_1_o, dist_1_2), degenerated_a); | |
640 | } | |
641 | ||
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 | |
92f5a8d4 | 646 | Point2 const& b1, Point2 const& /*b2*/, // in |
b32b8144 FG |
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 | |
92f5a8d4 | 651 | CalcT& dist_a1_b1, // out |
b32b8144 FG |
652 | bool degen_neq_coords = false) // in |
653 | { | |
654 | dist_a1_a2 = res_a1_a2.distance; | |
655 | ||
92f5a8d4 | 656 | dist_a1_b1 = res_a1_b1.distance; |
b32b8144 FG |
657 | if (! same_direction(res_a1_b1.azimuth, res_a1_a2.azimuth)) |
658 | { | |
92f5a8d4 | 659 | dist_a1_b1 = -dist_a1_b1; |
b32b8144 FG |
660 | } |
661 | ||
92f5a8d4 TL |
662 | // if b1 is close a1 |
663 | if (is_endpoint_equal(dist_a1_b1, a1, b1)) | |
b32b8144 | 664 | { |
92f5a8d4 | 665 | dist_a1_b1 = 0; |
b32b8144 FG |
666 | return true; |
667 | } | |
92f5a8d4 TL |
668 | // if b1 is close a2 |
669 | else if (is_endpoint_equal(dist_a1_a2 - dist_a1_b1, a2, b1)) | |
b32b8144 | 670 | { |
92f5a8d4 | 671 | dist_a1_b1 = dist_a1_a2; |
b32b8144 FG |
672 | return true; |
673 | } | |
674 | ||
92f5a8d4 | 675 | // check the other endpoint of degenerated segment near a pole |
b32b8144 FG |
676 | if (degen_neq_coords) |
677 | { | |
678 | static CalcT const c0 = 0; | |
679 | if (math::equals(res_a1_b2.distance, c0)) | |
680 | { | |
92f5a8d4 | 681 | dist_a1_b1 = 0; |
b32b8144 FG |
682 | return true; |
683 | } | |
684 | else if (math::equals(dist_a1_a2 - res_a1_b2.distance, c0)) | |
685 | { | |
92f5a8d4 | 686 | dist_a1_b1 = dist_a1_a2; |
b32b8144 FG |
687 | return true; |
688 | } | |
689 | } | |
690 | ||
691 | // or i1 is on b | |
92f5a8d4 | 692 | return segment_ratio<CalcT>(dist_a1_b1, dist_a1_a2).on_segment(); |
b32b8144 FG |
693 | } |
694 | ||
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 | |
714 | { | |
715 | dist_a1_a2 = res_a1_a2.distance; | |
716 | dist_b1_b2 = res_b1_b2.distance; | |
717 | ||
718 | // assign the IP if some endpoints overlap | |
b32b8144 FG |
719 | if (equals_point_point(a1, b1)) |
720 | { | |
721 | lon = a1_lon; | |
722 | lat = a1_lat; | |
723 | dist_a1_ip = 0; | |
724 | dist_b1_ip = 0; | |
725 | ip_flag = ipi_at_a1; | |
726 | return true; | |
727 | } | |
728 | else if (equals_point_point(a1, b2)) | |
729 | { | |
730 | lon = a1_lon; | |
731 | lat = a1_lat; | |
732 | dist_a1_ip = 0; | |
733 | dist_b1_ip = dist_b1_b2; | |
734 | ip_flag = ipi_at_a1; | |
735 | return true; | |
736 | } | |
737 | else if (equals_point_point(a2, b1)) | |
738 | { | |
739 | lon = a2_lon; | |
740 | lat = a2_lat; | |
741 | dist_a1_ip = dist_a1_a2; | |
742 | dist_b1_ip = 0; | |
743 | ip_flag = ipi_at_a2; | |
744 | return true; | |
745 | } | |
746 | else if (equals_point_point(a2, b2)) | |
747 | { | |
748 | lon = a2_lon; | |
749 | lat = a2_lat; | |
750 | dist_a1_ip = dist_a1_a2; | |
751 | dist_b1_ip = dist_b1_b2; | |
752 | ip_flag = ipi_at_a2; | |
753 | return true; | |
754 | } | |
755 | ||
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 | |
759 | { | |
760 | if (res_b1_a1.distance <= res_b1_b2.distance | |
761 | && same_direction(res_b1_a1.azimuth, res_b1_b2.azimuth)) | |
762 | { | |
763 | lon = a1_lon; | |
764 | lat = a1_lat; | |
765 | dist_a1_ip = 0; | |
766 | dist_b1_ip = res_b1_a1.distance; | |
767 | ip_flag = ipi_at_a1; | |
768 | return true; | |
769 | } | |
770 | else | |
771 | { | |
772 | return false; | |
773 | } | |
774 | } | |
775 | else if (sides.template get<0, 1>() == 0) // a2 wrt b | |
776 | { | |
777 | if (res_b1_a2.distance <= res_b1_b2.distance | |
778 | && same_direction(res_b1_a2.azimuth, res_b1_b2.azimuth)) | |
779 | { | |
780 | lon = a2_lon; | |
781 | lat = a2_lat; | |
782 | dist_a1_ip = res_a1_a2.distance; | |
783 | dist_b1_ip = res_b1_a2.distance; | |
784 | ip_flag = ipi_at_a2; | |
785 | return true; | |
786 | } | |
787 | else | |
788 | { | |
789 | return false; | |
790 | } | |
791 | } | |
792 | else if (sides.template get<1, 0>() == 0) // b1 wrt a | |
793 | { | |
794 | if (res_a1_b1.distance <= res_a1_a2.distance | |
795 | && same_direction(res_a1_b1.azimuth, res_a1_a2.azimuth)) | |
796 | { | |
797 | lon = b1_lon; | |
798 | lat = b1_lat; | |
799 | dist_a1_ip = res_a1_b1.distance; | |
800 | dist_b1_ip = 0; | |
801 | ip_flag = ipi_at_b1; | |
802 | return true; | |
803 | } | |
804 | else | |
805 | { | |
806 | return false; | |
807 | } | |
808 | } | |
809 | else if (sides.template get<1, 1>() == 0) // b2 wrt a | |
810 | { | |
811 | if (res_a1_b2.distance <= res_a1_a2.distance | |
812 | && same_direction(res_a1_b2.azimuth, res_a1_a2.azimuth)) | |
813 | { | |
814 | lon = b2_lon; | |
815 | lat = b2_lat; | |
816 | dist_a1_ip = res_a1_b2.distance; | |
817 | dist_b1_ip = res_b1_b2.distance; | |
818 | ip_flag = ipi_at_b2; | |
819 | return true; | |
820 | } | |
821 | else | |
822 | { | |
823 | return false; | |
824 | } | |
825 | } | |
826 | ||
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 | |
830 | ||
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, | |
834 | lon, lat, spheroid); | |
835 | ||
836 | if (! ok) | |
837 | { | |
838 | return false; | |
839 | } | |
840 | ||
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; | |
843 | ||
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)) | |
847 | { | |
848 | dist_a1_ip = -dist_a1_ip; | |
849 | } | |
850 | ||
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); | |
855 | ||
856 | if (! (is_on_a || is_on_a1 || is_on_a2)) | |
857 | { | |
858 | return false; | |
859 | } | |
860 | ||
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)) | |
864 | { | |
865 | dist_b1_ip = -dist_b1_ip; | |
866 | } | |
867 | ||
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); | |
872 | ||
873 | if (! (is_on_b || is_on_b1 || is_on_b2)) | |
874 | { | |
875 | return false; | |
876 | } | |
877 | ||
92f5a8d4 TL |
878 | typedef typename FormulaPolicy::template inverse<CalcT, true, false, false, false, false> inverse_dist; |
879 | ||
b32b8144 FG |
880 | ip_flag = ipi_inters; |
881 | ||
882 | if (is_on_b1) | |
883 | { | |
884 | lon = b1_lon; | |
885 | lat = b1_lat; | |
92f5a8d4 | 886 | dist_a1_ip = inverse_dist::apply(a1_lon, a1_lat, lon, lat, spheroid).distance; // for consistency |
b32b8144 FG |
887 | dist_b1_ip = 0; |
888 | ip_flag = ipi_at_b1; | |
889 | } | |
890 | else if (is_on_b2) | |
891 | { | |
892 | lon = b2_lon; | |
893 | lat = b2_lat; | |
92f5a8d4 | 894 | dist_a1_ip = inverse_dist::apply(a1_lon, a1_lat, lon, lat, spheroid).distance; // for consistency |
b32b8144 FG |
895 | dist_b1_ip = res_b1_b2.distance; |
896 | ip_flag = ipi_at_b2; | |
897 | } | |
898 | ||
899 | if (is_on_a1) | |
900 | { | |
901 | lon = a1_lon; | |
902 | lat = a1_lat; | |
903 | dist_a1_ip = 0; | |
92f5a8d4 | 904 | dist_b1_ip = inverse_dist::apply(b1_lon, b1_lat, lon, lat, spheroid).distance; // for consistency |
b32b8144 FG |
905 | ip_flag = ipi_at_a1; |
906 | } | |
907 | else if (is_on_a2) | |
908 | { | |
909 | lon = a2_lon; | |
910 | lat = a2_lat; | |
911 | dist_a1_ip = res_a1_a2.distance; | |
92f5a8d4 | 912 | dist_b1_ip = inverse_dist::apply(b1_lon, b1_lat, lon, lat, spheroid).distance; // for consistency |
b32b8144 FG |
913 | ip_flag = ipi_at_a2; |
914 | } | |
915 | ||
916 | return true; | |
917 | } | |
918 | ||
919 | template <typename CalcT, typename P1, typename P2> | |
920 | static inline bool is_endpoint_equal(CalcT const& dist, | |
92f5a8d4 | 921 | P1 const& ai, P2 const& b1) |
b32b8144 FG |
922 | { |
923 | static CalcT const c0 = 0; | |
92f5a8d4 | 924 | return is_near(dist) && (math::equals(dist, c0) || equals_point_point(ai, b1)); |
b32b8144 FG |
925 | } |
926 | ||
927 | template <typename CalcT> | |
928 | static inline bool is_near(CalcT const& dist) | |
929 | { | |
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; | |
933 | } | |
934 | ||
935 | template <typename ProjCoord1, typename ProjCoord2> | |
936 | static inline int position_value(ProjCoord1 const& ca1, | |
937 | ProjCoord2 const& cb1, | |
938 | ProjCoord2 const& cb2) | |
939 | { | |
940 | // S1x 0 1 2 3 4 | |
941 | // S2 |----------> | |
942 | return math::equals(ca1, cb1) ? 1 | |
943 | : math::equals(ca1, cb2) ? 3 | |
944 | : cb1 < cb2 ? | |
945 | ( ca1 < cb1 ? 0 | |
946 | : ca1 > cb2 ? 4 | |
947 | : 2 ) | |
948 | : ( ca1 > cb1 ? 0 | |
949 | : ca1 < cb2 ? 4 | |
950 | : 2 ); | |
951 | } | |
952 | ||
953 | template <typename CalcT> | |
954 | static inline bool same_direction(CalcT const& azimuth1, CalcT const& azimuth2) | |
955 | { | |
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>(); | |
959 | } | |
960 | ||
961 | template <int Which> | |
962 | static inline void sides_reverse_segment(side_info & sides) | |
963 | { | |
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); | |
972 | } | |
973 | ||
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) | |
977 | { | |
978 | ip_flag = ip_flag == ipi_at_p1 ? ipi_at_p2 : | |
979 | ip_flag == ipi_at_p2 ? ipi_at_p1 : | |
980 | ip_flag; | |
981 | } | |
982 | ||
92f5a8d4 TL |
983 | template <typename Point1, typename Point2> |
984 | static inline bool equals_point_point(Point1 const& point1, Point2 const& point2) | |
985 | { | |
986 | return detail::equals::equals_point_point(point1, point2, | |
987 | point_in_point_strategy_type()); | |
988 | } | |
989 | ||
b32b8144 FG |
990 | private: |
991 | Spheroid m_spheroid; | |
992 | }; | |
993 | ||
994 | ||
995 | }} // namespace strategy::intersection | |
996 | ||
997 | }} // namespace boost::geometry | |
998 | ||
999 | ||
1000 | #endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_INTERSECTION_HPP |