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