1 // Boost.Geometry 2 3 // Copyright (c) 2017 Adam Wulkiewicz, Lodz, Poland. 4 5 // Copyright (c) 2016-2019, 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 17 #include <boost/geometry/core/cs.hpp> 18 #include <boost/geometry/core/access.hpp> 19 #include <boost/geometry/core/radian_access.hpp> 20 #include <boost/geometry/core/tags.hpp> 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 37 #include <boost/geometry/srs/spheroid.hpp> 38 39 #include <boost/geometry/strategies/geographic/area.hpp> 40 #include <boost/geometry/strategies/geographic/disjoint_segment_box.hpp> 41 #include <boost/geometry/strategies/geographic/distance.hpp> 42 #include <boost/geometry/strategies/geographic/envelope.hpp> 43 #include <boost/geometry/strategies/geographic/parameters.hpp> 44 #include <boost/geometry/strategies/geographic/point_in_poly_winding.hpp> 45 #include <boost/geometry/strategies/geographic/side.hpp> 46 #include <boost/geometry/strategies/spherical/expand_box.hpp> 47 #include <boost/geometry/strategies/spherical/disjoint_box_box.hpp> 48 #include <boost/geometry/strategies/spherical/point_in_point.hpp> 49 #include <boost/geometry/strategies/intersection.hpp> 50 #include <boost/geometry/strategies/intersection_result.hpp> 51 #include <boost/geometry/strategies/side_info.hpp> 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 { 77 typedef geographic_tag cs_tag; 78 79 typedef side::geographic 80 < 81 FormulaPolicy, Spheroid, CalculationType 82 > side_strategy_type; 83 get_side_strategyboost::geometry::strategy::intersection::geographic_segments84 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 get_point_in_geometry_strategyboost::geometry::strategy::intersection::geographic_segments104 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 < 118 FormulaPolicy, 119 Order, 120 Spheroid, 121 CalculationType 122 > type; 123 }; 124 125 template <typename Geometry> get_area_strategyboost::geometry::strategy::intersection::geographic_segments126 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> get_distance_strategyboost::geometry::strategy::intersection::geographic_segments144 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 150 typedef envelope::geographic<FormulaPolicy, Spheroid, CalculationType> 151 envelope_strategy_type; 152 get_envelope_strategyboost::geometry::strategy::intersection::geographic_segments153 inline envelope_strategy_type get_envelope_strategy() const 154 { 155 return envelope_strategy_type(m_spheroid); 156 } 157 158 typedef expand::geographic_segment<FormulaPolicy, Spheroid, CalculationType> 159 expand_strategy_type; 160 get_expand_strategyboost::geometry::strategy::intersection::geographic_segments161 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 get_point_in_point_strategyboost::geometry::strategy::intersection::geographic_segments168 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 get_equals_point_point_strategyboost::geometry::strategy::intersection::geographic_segments175 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 get_disjoint_box_box_strategyboost::geometry::strategy::intersection::geographic_segments182 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 get_disjoint_segment_box_strategyboost::geometry::strategy::intersection::geographic_segments192 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 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 { 208 template <typename Point, typename Segment1, typename Segment2> calculateboost::geometry::strategy::intersection::geographic_segments::segment_intersection_info209 void calculate(Point& point, Segment1 const& a, Segment2 const& b) const 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 geographic_segmentsboost::geometry::strategy::intersection::geographic_segments242 explicit geographic_segments(Spheroid const& spheroid = Spheroid()) 243 : m_spheroid(spheroid) 244 {} 245 246 // Relate segments a and b 247 template 248 < 249 typename UniqueSubRange1, 250 typename UniqueSubRange2, 251 typename Policy 252 > applyboost::geometry::strategy::intersection::geographic_segments253 inline typename Policy::return_type apply(UniqueSubRange1 const& range_p, 254 UniqueSubRange2 const& range_q, 255 Policy const&) const 256 { 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); 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 > applyboost::geometry::strategy::intersection::geographic_segments299 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? 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); 466 calculate_collinear_data(a1, a2, b2, b1, res_a1_a2, res_a1_b2, res_a1_b1, dist_a1_a2, dist_a1_b2); 467 dist_b1_b2 = dist_a1_b2 - dist_a1_b1; 468 dist_b1_a1 = -dist_a1_b1; 469 dist_b1_a2 = dist_a1_a2 - dist_a1_b1; 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); 474 calculate_collinear_data(b1, b2, a2, a1, res_b1_b2, res_b1_a2, res_b1_a1, dist_b1_b2, dist_b1_a2); 475 dist_a1_a2 = dist_b1_a2 - dist_b1_a1; 476 dist_a1_b1 = -dist_b1_a1; 477 dist_a1_b2 = dist_b1_b2 - dist_b1_a1; 478 } 479 480 // NOTE: this is probably not needed 481 int a1_on_b = position_value(c0, dist_a1_b1, dist_a1_b2); 482 int a2_on_b = position_value(dist_a1_a2, dist_a1_b1, dist_a1_b2); 483 int b1_on_a = position_value(c0, dist_b1_a1, dist_b1_a2); 484 int b2_on_a = position_value(dist_b1_b2, dist_b1_a1, dist_b1_a2); 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 collinear_one_degeneratedboost::geometry::strategy::intersection::geographic_segments617 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> calculate_collinear_databoost::geometry::strategy::intersection::geographic_segments645 static inline bool calculate_collinear_data(Point1 const& a1, Point1 const& a2, // in 646 Point2 const& b1, Point2 const& /*b2*/, // in 647 ResultInverse const& res_a1_a2, // in 648 ResultInverse const& res_a1_b1, // in 649 ResultInverse const& res_a1_b2, // in 650 CalcT& dist_a1_a2, // out 651 CalcT& dist_a1_b1, // out 652 bool degen_neq_coords = false) // in 653 { 654 dist_a1_a2 = res_a1_a2.distance; 655 656 dist_a1_b1 = res_a1_b1.distance; 657 if (! same_direction(res_a1_b1.azimuth, res_a1_a2.azimuth)) 658 { 659 dist_a1_b1 = -dist_a1_b1; 660 } 661 662 // if b1 is close a1 663 if (is_endpoint_equal(dist_a1_b1, a1, b1)) 664 { 665 dist_a1_b1 = 0; 666 return true; 667 } 668 // if b1 is close a2 669 else if (is_endpoint_equal(dist_a1_a2 - dist_a1_b1, a2, b1)) 670 { 671 dist_a1_b1 = dist_a1_a2; 672 return true; 673 } 674 675 // check the other endpoint of degenerated segment near a pole 676 if (degen_neq_coords) 677 { 678 static CalcT const c0 = 0; 679 if (math::equals(res_a1_b2.distance, c0)) 680 { 681 dist_a1_b1 = 0; 682 return true; 683 } 684 else if (math::equals(dist_a1_a2 - res_a1_b2.distance, c0)) 685 { 686 dist_a1_b1 = dist_a1_a2; 687 return true; 688 } 689 } 690 691 // or i1 is on b 692 return segment_ratio<CalcT>(dist_a1_b1, dist_a1_a2).on_segment(); 693 } 694 695 template <typename Point1, typename Point2, typename CalcT, typename ResultInverse, typename Spheroid_> calculate_ip_databoost::geometry::strategy::intersection::geographic_segments696 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 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 878 typedef typename FormulaPolicy::template inverse<CalcT, true, false, false, false, false> inverse_dist; 879 880 ip_flag = ipi_inters; 881 882 if (is_on_b1) 883 { 884 lon = b1_lon; 885 lat = b1_lat; 886 dist_a1_ip = inverse_dist::apply(a1_lon, a1_lat, lon, lat, spheroid).distance; // for consistency 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; 894 dist_a1_ip = inverse_dist::apply(a1_lon, a1_lat, lon, lat, spheroid).distance; // for consistency 895 dist_b1_ip = res_b1_b2.distance; 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; 904 dist_b1_ip = inverse_dist::apply(b1_lon, b1_lat, lon, lat, spheroid).distance; // for consistency 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; 912 dist_b1_ip = inverse_dist::apply(b1_lon, b1_lat, lon, lat, spheroid).distance; // for consistency 913 ip_flag = ipi_at_a2; 914 } 915 916 return true; 917 } 918 919 template <typename CalcT, typename P1, typename P2> is_endpoint_equalboost::geometry::strategy::intersection::geographic_segments920 static inline bool is_endpoint_equal(CalcT const& dist, 921 P1 const& ai, P2 const& b1) 922 { 923 static CalcT const c0 = 0; 924 return is_near(dist) && (math::equals(dist, c0) || equals_point_point(ai, b1)); 925 } 926 927 template <typename CalcT> is_nearboost::geometry::strategy::intersection::geographic_segments928 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> position_valueboost::geometry::strategy::intersection::geographic_segments936 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> same_directionboost::geometry::strategy::intersection::geographic_segments954 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> sides_reverse_segmentboost::geometry::strategy::intersection::geographic_segments962 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 ip_flag_reverse_segmentboost::geometry::strategy::intersection::geographic_segments974 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 983 template <typename Point1, typename Point2> equals_point_pointboost::geometry::strategy::intersection::geographic_segments984 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 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 1001