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