1 // Boost.Geometry (aka GGL, Generic Geometry Library) 2 3 // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. 4 // Copyright (c) 2013-2017 Adam Wulkiewicz, Lodz, Poland. 5 6 // This file was modified by Oracle on 2013, 2014, 2016, 2017, 2018, 2019. 7 // Modifications copyright (c) 2013-2019 Oracle and/or its affiliates. 8 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle 9 10 // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library 11 // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands. 12 13 // Use, modification and distribution is subject to the Boost Software License, 14 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at 15 // http://www.boost.org/LICENSE_1_0.txt) 16 17 #ifndef BOOST_GEOMETRY_STRATEGY_SPHERICAL_POINT_IN_POLY_WINDING_HPP 18 #define BOOST_GEOMETRY_STRATEGY_SPHERICAL_POINT_IN_POLY_WINDING_HPP 19 20 21 #include <boost/geometry/core/access.hpp> 22 #include <boost/geometry/core/coordinate_system.hpp> 23 #include <boost/geometry/core/cs.hpp> 24 #include <boost/geometry/core/tags.hpp> 25 26 #include <boost/geometry/util/math.hpp> 27 #include <boost/geometry/util/select_calculation_type.hpp> 28 #include <boost/geometry/util/normalize_spheroidal_coordinates.hpp> 29 30 #include <boost/geometry/strategies/cartesian/point_in_box.hpp> 31 #include <boost/geometry/strategies/covered_by.hpp> 32 #include <boost/geometry/strategies/side.hpp> 33 #include <boost/geometry/strategies/spherical/disjoint_box_box.hpp> 34 #include <boost/geometry/strategies/spherical/ssf.hpp> 35 #include <boost/geometry/strategies/within.hpp> 36 37 38 namespace boost { namespace geometry 39 { 40 41 namespace strategy { namespace within 42 { 43 44 45 #ifndef DOXYGEN_NO_DETAIL 46 namespace detail 47 { 48 49 template <typename SideStrategy, typename CalculationType> 50 class spherical_winding_base 51 { 52 template <typename Point, typename PointOfSegment> 53 struct calculation_type 54 : select_calculation_type 55 < 56 Point, 57 PointOfSegment, 58 CalculationType 59 > 60 {}; 61 62 /*! subclass to keep state */ 63 class counter 64 { 65 int m_count; 66 //int m_count_n; 67 int m_count_s; 68 int m_raw_count; 69 int m_raw_count_anti; 70 bool m_touches; 71 code() const72 inline int code() const 73 { 74 if (m_touches) 75 { 76 return 0; 77 } 78 79 if (m_raw_count != 0 && m_raw_count_anti != 0) 80 { 81 if (m_raw_count > 0) // right, wrap around south pole 82 { 83 return (m_count + m_count_s) == 0 ? -1 : 1; 84 } 85 else // left, wrap around north pole 86 { 87 //return (m_count + m_count_n) == 0 ? -1 : 1; 88 // m_count_n is 0 89 return m_count == 0 ? -1 : 1; 90 } 91 } 92 93 return m_count == 0 ? -1 : 1; 94 } 95 96 public : 97 friend class spherical_winding_base; 98 counter()99 inline counter() 100 : m_count(0) 101 //, m_count_n(0) 102 , m_count_s(0) 103 , m_raw_count(0) 104 , m_raw_count_anti(0) 105 , m_touches(false) 106 {} 107 108 }; 109 110 struct count_info 111 { count_infoboost::geometry::strategy::within::detail::spherical_winding_base::count_info112 explicit count_info(int c = 0, bool ia = false) 113 : count(c) 114 , is_anti(ia) 115 {} 116 117 int count; 118 bool is_anti; 119 }; 120 121 public: 122 typedef typename SideStrategy::cs_tag cs_tag; 123 124 typedef SideStrategy side_strategy_type; 125 get_side_strategy() const126 inline side_strategy_type get_side_strategy() const 127 { 128 return m_side_strategy; 129 } 130 131 typedef expand::spherical_point expand_point_strategy_type; 132 133 typedef typename SideStrategy::envelope_strategy_type envelope_strategy_type; 134 get_envelope_strategy() const135 inline envelope_strategy_type get_envelope_strategy() const 136 { 137 return m_side_strategy.get_envelope_strategy(); 138 } 139 140 typedef typename SideStrategy::disjoint_strategy_type disjoint_strategy_type; 141 get_disjoint_strategy() const142 inline disjoint_strategy_type get_disjoint_strategy() const 143 { 144 return m_side_strategy.get_disjoint_strategy(); 145 } 146 147 typedef typename SideStrategy::equals_point_point_strategy_type equals_point_point_strategy_type; get_equals_point_point_strategy() const148 inline equals_point_point_strategy_type get_equals_point_point_strategy() const 149 { 150 return m_side_strategy.get_equals_point_point_strategy(); 151 } 152 153 typedef disjoint::spherical_box_box disjoint_box_box_strategy_type; get_disjoint_box_box_strategy()154 static inline disjoint_box_box_strategy_type get_disjoint_box_box_strategy() 155 { 156 return disjoint_box_box_strategy_type(); 157 } 158 159 typedef covered_by::spherical_point_box disjoint_point_box_strategy_type; 160 spherical_winding_base()161 spherical_winding_base() 162 {} 163 164 template <typename Model> spherical_winding_base(Model const & model)165 explicit spherical_winding_base(Model const& model) 166 : m_side_strategy(model) 167 {} 168 169 // Typedefs and static methods to fulfill the concept 170 typedef counter state_type; 171 172 template <typename Point, typename PointOfSegment> apply(Point const & point,PointOfSegment const & s1,PointOfSegment const & s2,counter & state) const173 inline bool apply(Point const& point, 174 PointOfSegment const& s1, PointOfSegment const& s2, 175 counter& state) const 176 { 177 typedef typename calculation_type<Point, PointOfSegment>::type calc_t; 178 typedef typename geometry::detail::cs_angular_units<Point>::type units_t; 179 typedef math::detail::constants_on_spheroid<calc_t, units_t> constants; 180 181 bool eq1 = false; 182 bool eq2 = false; 183 bool s_antipodal = false; 184 185 count_info ci = check_segment(point, s1, s2, state, eq1, eq2, s_antipodal); 186 if (ci.count != 0) 187 { 188 if (! ci.is_anti) 189 { 190 int side = 0; 191 if (ci.count == 1 || ci.count == -1) 192 { 193 side = side_equal(point, eq1 ? s1 : s2, ci); 194 } 195 else // count == 2 || count == -2 196 { 197 if (! s_antipodal) 198 { 199 // 1 left, -1 right 200 side = m_side_strategy.apply(s1, s2, point); 201 } 202 else 203 { 204 calc_t const pi = constants::half_period(); 205 calc_t const s1_lat = get<1>(s1); 206 calc_t const s2_lat = get<1>(s2); 207 208 side = math::sign(ci.count) 209 * (pi - s1_lat - s2_lat <= pi // segment goes through north pole 210 ? -1 // going right all points will be on right side 211 : 1); // going right all points will be on left side 212 } 213 } 214 215 if (side == 0) 216 { 217 // Point is lying on segment 218 state.m_touches = true; 219 state.m_count = 0; 220 return false; 221 } 222 223 // Side is NEG for right, POS for left. 224 // The count is -2 for left, 2 for right (or -1/1) 225 // Side positive thus means RIGHT and LEFTSIDE or LEFT and RIGHTSIDE 226 // See accompagnying figure (TODO) 227 if (side * ci.count > 0) 228 { 229 state.m_count += ci.count; 230 } 231 232 state.m_raw_count += ci.count; 233 } 234 else 235 { 236 // Count negated because the segment is on the other side of the globe 237 // so it is reversed to match this side of the globe 238 239 // Assuming geometry wraps around north pole, for segments on the other side of the globe 240 // the point will always be RIGHT+RIGHTSIDE or LEFT+LEFTSIDE, so side*-count always < 0 241 //state.m_count_n -= 0; 242 243 // Assuming geometry wraps around south pole, for segments on the other side of the globe 244 // the point will always be RIGHT+LEFTSIDE or LEFT+RIGHTSIDE, so side*-count always > 0 245 state.m_count_s -= ci.count; 246 247 state.m_raw_count_anti -= ci.count; 248 } 249 } 250 return ! state.m_touches; 251 } 252 result(counter const & state)253 static inline int result(counter const& state) 254 { 255 return state.code(); 256 } 257 258 protected: 259 template <typename Point, typename PointOfSegment> check_segment(Point const & point,PointOfSegment const & seg1,PointOfSegment const & seg2,counter & state,bool & eq1,bool & eq2,bool & s_antipodal)260 static inline count_info check_segment(Point const& point, 261 PointOfSegment const& seg1, 262 PointOfSegment const& seg2, 263 counter& state, 264 bool& eq1, bool& eq2, bool& s_antipodal) 265 { 266 if (check_touch(point, seg1, seg2, state, eq1, eq2, s_antipodal)) 267 { 268 return count_info(0, false); 269 } 270 271 return calculate_count(point, seg1, seg2, eq1, eq2, s_antipodal); 272 } 273 274 template <typename Point, typename PointOfSegment> check_touch(Point const & point,PointOfSegment const & seg1,PointOfSegment const & seg2,counter & state,bool & eq1,bool & eq2,bool & s_antipodal)275 static inline int check_touch(Point const& point, 276 PointOfSegment const& seg1, 277 PointOfSegment const& seg2, 278 counter& state, 279 bool& eq1, 280 bool& eq2, 281 bool& s_antipodal) 282 { 283 typedef typename calculation_type<Point, PointOfSegment>::type calc_t; 284 typedef typename geometry::detail::cs_angular_units<Point>::type units_t; 285 typedef math::detail::constants_on_spheroid<calc_t, units_t> constants; 286 287 calc_t const c0 = 0; 288 calc_t const c2 = 2; 289 calc_t const pi = constants::half_period(); 290 calc_t const half_pi = pi / c2; 291 292 calc_t const p_lon = get<0>(point); 293 calc_t const s1_lon = get<0>(seg1); 294 calc_t const s2_lon = get<0>(seg2); 295 calc_t const p_lat = get<1>(point); 296 calc_t const s1_lat = get<1>(seg1); 297 calc_t const s2_lat = get<1>(seg2); 298 299 // NOTE: lat in {-90, 90} and arbitrary lon 300 // it doesn't matter what lon it is if it's a pole 301 // so e.g. if one of the segment endpoints is a pole 302 // then only the other lon matters 303 304 bool eq1_strict = longitudes_equal<units_t>(s1_lon, p_lon); 305 bool eq2_strict = longitudes_equal<units_t>(s2_lon, p_lon); 306 bool eq1_anti = false; 307 bool eq2_anti = false; 308 309 calc_t const anti_p_lon = p_lon + (p_lon <= c0 ? pi : -pi); 310 311 eq1 = eq1_strict // lon strictly equal to s1 312 || (eq1_anti = longitudes_equal<units_t>(s1_lon, anti_p_lon)) // anti-lon strictly equal to s1 313 || math::equals(math::abs(s1_lat), half_pi); // s1 is pole 314 eq2 = eq2_strict // lon strictly equal to s2 315 || (eq2_anti = longitudes_equal<units_t>(s2_lon, anti_p_lon)) // anti-lon strictly equal to s2 316 || math::equals(math::abs(s2_lat), half_pi); // s2 is pole 317 318 // segment overlapping pole 319 calc_t const s_lon_diff = math::longitude_distance_signed<units_t>(s1_lon, s2_lon); 320 s_antipodal = math::equals(s_lon_diff, pi); 321 if (s_antipodal) 322 { 323 eq1 = eq2 = eq1 || eq2; 324 325 // segment overlapping pole and point is pole 326 if (math::equals(math::abs(p_lat), half_pi)) 327 { 328 eq1 = eq2 = true; 329 } 330 } 331 332 // Both equal p -> segment vertical 333 // The only thing which has to be done is check if point is ON segment 334 if (eq1 && eq2) 335 { 336 // segment endpoints on the same sides of the globe 337 if (! s_antipodal) 338 { 339 // p's lat between segment endpoints' lats 340 if ( (s1_lat <= p_lat && s2_lat >= p_lat) || (s2_lat <= p_lat && s1_lat >= p_lat) ) 341 { 342 if (!eq1_anti || !eq2_anti) 343 { 344 state.m_touches = true; 345 } 346 } 347 } 348 else 349 { 350 // going through north or south pole? 351 if (pi - s1_lat - s2_lat <= pi) 352 { 353 if ( (eq1_strict && s1_lat <= p_lat) || (eq2_strict && s2_lat <= p_lat) // north 354 || math::equals(p_lat, half_pi) ) // point on north pole 355 { 356 state.m_touches = true; 357 } 358 else if (! eq1_strict && ! eq2_strict && math::equals(p_lat, -half_pi) ) // point on south pole 359 { 360 return false; 361 } 362 } 363 else // south pole 364 { 365 if ( (eq1_strict && s1_lat >= p_lat) || (eq2_strict && s2_lat >= p_lat) // south 366 || math::equals(p_lat, -half_pi) ) // point on south pole 367 { 368 state.m_touches = true; 369 } 370 else if (! eq1_strict && ! eq2_strict && math::equals(p_lat, half_pi) ) // point on north pole 371 { 372 return false; 373 } 374 } 375 } 376 377 return true; 378 } 379 380 return false; 381 } 382 383 // Called if point is not aligned with a vertical segment 384 template <typename Point, typename PointOfSegment> calculate_count(Point const & point,PointOfSegment const & seg1,PointOfSegment const & seg2,bool eq1,bool eq2,bool s_antipodal)385 static inline count_info calculate_count(Point const& point, 386 PointOfSegment const& seg1, 387 PointOfSegment const& seg2, 388 bool eq1, bool eq2, bool s_antipodal) 389 { 390 typedef typename calculation_type<Point, PointOfSegment>::type calc_t; 391 typedef typename geometry::detail::cs_angular_units<Point>::type units_t; 392 typedef math::detail::constants_on_spheroid<calc_t, units_t> constants; 393 394 // If both segment endpoints were poles below checks wouldn't be enough 395 // but this means that either both are the same or that they are N/S poles 396 // and therefore the segment is not valid. 397 // If needed (eq1 && eq2 ? 0) could be returned 398 399 calc_t const c0 = 0; 400 calc_t const pi = constants::half_period(); 401 402 calc_t const p = get<0>(point); 403 calc_t const s1 = get<0>(seg1); 404 calc_t const s2 = get<0>(seg2); 405 406 calc_t const s1_p = math::longitude_distance_signed<units_t>(s1, p); 407 408 if (s_antipodal) 409 { 410 return count_info(s1_p < c0 ? -2 : 2, false); // choose W/E 411 } 412 413 calc_t const s1_s2 = math::longitude_distance_signed<units_t>(s1, s2); 414 415 if (eq1 || eq2) // Point on level s1 or s2 416 { 417 return count_info(s1_s2 < c0 ? -1 : 1, // choose W/E 418 longitudes_equal<units_t>(p + pi, (eq1 ? s1 : s2))); 419 } 420 421 // Point between s1 and s2 422 if ( math::sign(s1_p) == math::sign(s1_s2) 423 && math::abs(s1_p) < math::abs(s1_s2) ) 424 { 425 return count_info(s1_s2 < c0 ? -2 : 2, false); // choose W/E 426 } 427 428 calc_t const s1_p_anti = math::longitude_distance_signed<units_t>(s1, p + pi); 429 430 // Anti-Point between s1 and s2 431 if ( math::sign(s1_p_anti) == math::sign(s1_s2) 432 && math::abs(s1_p_anti) < math::abs(s1_s2) ) 433 { 434 return count_info(s1_s2 < c0 ? -2 : 2, true); // choose W/E 435 } 436 437 return count_info(0, false); 438 } 439 440 441 // Fix for https://svn.boost.org/trac/boost/ticket/9628 442 // For floating point coordinates, the <D> coordinate of a point is compared 443 // with the segment's points using some EPS. If the coordinates are "equal" 444 // the sides are calculated. Therefore we can treat a segment as a long areal 445 // geometry having some width. There is a small ~triangular area somewhere 446 // between the segment's effective area and a segment's line used in sides 447 // calculation where the segment is on the one side of the line but on the 448 // other side of a segment (due to the width). 449 // Below picture assuming D = 1, if D = 0 horiz<->vert, E<->N, RIGHT<->UP. 450 // For the s1 of a segment going NE the real side is RIGHT but the point may 451 // be detected as LEFT, like this: 452 // RIGHT 453 // ___-----> 454 // ^ O Pt __ __ 455 // EPS __ __ 456 // v__ __ BUT DETECTED AS LEFT OF THIS LINE 457 // _____7 458 // _____/ 459 // _____/ 460 // In the code below actually D = 0, so segments are nearly-vertical 461 // Called when the point is on the same level as one of the segment's points 462 // but the point is not aligned with a vertical segment 463 template <typename Point, typename PointOfSegment> side_equal(Point const & point,PointOfSegment const & se,count_info const & ci) const464 inline int side_equal(Point const& point, 465 PointOfSegment const& se, 466 count_info const& ci) const 467 { 468 typedef typename coordinate_type<PointOfSegment>::type scoord_t; 469 typedef typename geometry::detail::cs_angular_units<Point>::type units_t; 470 471 if (math::equals(get<1>(point), get<1>(se))) 472 { 473 return 0; 474 } 475 476 // Create a horizontal segment intersecting the original segment's endpoint 477 // equal to the point, with the derived direction (E/W). 478 PointOfSegment ss1, ss2; 479 set<1>(ss1, get<1>(se)); 480 set<0>(ss1, get<0>(se)); 481 set<1>(ss2, get<1>(se)); 482 scoord_t ss20 = get<0>(se); 483 if (ci.count > 0) 484 { 485 ss20 += small_angle<scoord_t, units_t>(); 486 } 487 else 488 { 489 ss20 -= small_angle<scoord_t, units_t>(); 490 } 491 math::normalize_longitude<units_t>(ss20); 492 set<0>(ss2, ss20); 493 494 // Check the side using this vertical segment 495 return m_side_strategy.apply(ss1, ss2, point); 496 } 497 498 // 1 deg or pi/180 rad 499 template <typename CalcT, typename Units> small_angle()500 static inline CalcT small_angle() 501 { 502 typedef math::detail::constants_on_spheroid<CalcT, Units> constants; 503 504 return constants::half_period() / CalcT(180); 505 } 506 507 template <typename Units, typename CalcT> longitudes_equal(CalcT const & lon1,CalcT const & lon2)508 static inline bool longitudes_equal(CalcT const& lon1, CalcT const& lon2) 509 { 510 return math::equals( 511 math::longitude_distance_signed<Units>(lon1, lon2), 512 CalcT(0)); 513 } 514 515 SideStrategy m_side_strategy; 516 }; 517 518 519 } // namespace detail 520 #endif // DOXYGEN_NO_DETAIL 521 522 523 /*! 524 \brief Within detection using winding rule in spherical coordinate system. 525 \ingroup strategies 526 \tparam Point \tparam_point 527 \tparam PointOfSegment \tparam_segment_point 528 \tparam CalculationType \tparam_calculation 529 530 \qbk{ 531 [heading See also] 532 [link geometry.reference.algorithms.within.within_3_with_strategy within (with strategy)] 533 } 534 */ 535 template 536 < 537 typename Point = void, // for backward compatibility 538 typename PointOfSegment = Point, // for backward compatibility 539 typename CalculationType = void 540 > 541 class spherical_winding 542 : public within::detail::spherical_winding_base 543 < 544 side::spherical_side_formula<CalculationType>, 545 CalculationType 546 > 547 {}; 548 549 550 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS 551 552 namespace services 553 { 554 555 template <typename PointLike, typename Geometry, typename AnyTag1, typename AnyTag2> 556 struct default_strategy<PointLike, Geometry, AnyTag1, AnyTag2, pointlike_tag, polygonal_tag, spherical_tag, spherical_tag> 557 { 558 typedef within::detail::spherical_winding_base 559 < 560 typename strategy::side::services::default_strategy 561 < 562 typename cs_tag<PointLike>::type 563 >::type, 564 void 565 > type; 566 }; 567 568 template <typename PointLike, typename Geometry, typename AnyTag1, typename AnyTag2> 569 struct default_strategy<PointLike, Geometry, AnyTag1, AnyTag2, pointlike_tag, linear_tag, spherical_tag, spherical_tag> 570 { 571 typedef within::detail::spherical_winding_base 572 < 573 typename strategy::side::services::default_strategy 574 < 575 typename cs_tag<PointLike>::type 576 >::type, 577 void 578 > type; 579 }; 580 581 } // namespace services 582 583 #endif 584 585 586 }} // namespace strategy::within 587 588 589 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS 590 namespace strategy { namespace covered_by { namespace services 591 { 592 593 template <typename PointLike, typename Geometry, typename AnyTag1, typename AnyTag2> 594 struct default_strategy<PointLike, Geometry, AnyTag1, AnyTag2, pointlike_tag, polygonal_tag, spherical_tag, spherical_tag> 595 { 596 typedef within::detail::spherical_winding_base 597 < 598 typename strategy::side::services::default_strategy 599 < 600 typename cs_tag<PointLike>::type 601 >::type, 602 void 603 > type; 604 }; 605 606 template <typename PointLike, typename Geometry, typename AnyTag1, typename AnyTag2> 607 struct default_strategy<PointLike, Geometry, AnyTag1, AnyTag2, pointlike_tag, linear_tag, spherical_tag, spherical_tag> 608 { 609 typedef within::detail::spherical_winding_base 610 < 611 typename strategy::side::services::default_strategy 612 < 613 typename cs_tag<PointLike>::type 614 >::type, 615 void 616 > type; 617 }; 618 619 }}} // namespace strategy::covered_by::services 620 #endif 621 622 623 }} // namespace boost::geometry 624 625 626 #endif // BOOST_GEOMETRY_STRATEGY_SPHERICAL_POINT_IN_POLY_WINDING_HPP 627