1 // Boost.Geometry (aka GGL, Generic Geometry Library) 2 3 // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. 4 // Copyright (c) 2013 Adam Wulkiewicz, Lodz, Poland. 5 6 // This file was modified by Oracle on 2013, 2014, 2016, 2017. 7 // Modifications copyright (c) 2013-2017 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_AGNOSTIC_POINT_IN_POLY_WINDING_HPP 18 #define BOOST_GEOMETRY_STRATEGY_AGNOSTIC_POINT_IN_POLY_WINDING_HPP 19 20 21 #include <boost/core/ignore_unused.hpp> 22 23 #include <boost/geometry/util/math.hpp> 24 #include <boost/geometry/util/select_calculation_type.hpp> 25 26 #include <boost/geometry/strategies/side.hpp> 27 #include <boost/geometry/strategies/covered_by.hpp> 28 #include <boost/geometry/strategies/within.hpp> 29 30 31 namespace boost { namespace geometry 32 { 33 34 namespace strategy { namespace within 35 { 36 37 // 1 deg or pi/180 rad 38 template <typename Point, 39 typename CalculationType = typename coordinate_type<Point>::type> 40 struct winding_small_angle 41 { 42 typedef typename coordinate_system<Point>::type cs_t; 43 typedef math::detail::constants_on_spheroid 44 < 45 CalculationType, 46 typename cs_t::units 47 > constants; 48 applyboost::geometry::strategy::within::winding_small_angle49 static inline CalculationType apply() 50 { 51 return constants::half_period() / CalculationType(180); 52 } 53 }; 54 55 56 // Fix for https://svn.boost.org/trac/boost/ticket/9628 57 // For floating point coordinates, the <D> coordinate of a point is compared 58 // with the segment's points using some EPS. If the coordinates are "equal" 59 // the sides are calculated. Therefore we can treat a segment as a long areal 60 // geometry having some width. There is a small ~triangular area somewhere 61 // between the segment's effective area and a segment's line used in sides 62 // calculation where the segment is on the one side of the line but on the 63 // other side of a segment (due to the width). 64 // Below picture assuming D = 1, if D = 0 horiz<->vert, E<->N, RIGHT<->UP. 65 // For the s1 of a segment going NE the real side is RIGHT but the point may 66 // be detected as LEFT, like this: 67 // RIGHT 68 // ___-----> 69 // ^ O Pt __ __ 70 // EPS __ __ 71 // v__ __ BUT DETECTED AS LEFT OF THIS LINE 72 // _____7 73 // _____/ 74 // _____/ 75 // In the code below actually D = 0, so segments are nearly-vertical 76 // Called when the point is on the same level as one of the segment's points 77 // but the point is not aligned with a vertical segment 78 template <typename CSTag> 79 struct winding_side_equal 80 { 81 typedef typename strategy::side::services::default_strategy 82 < 83 CSTag 84 >::type strategy_side_type; 85 86 template <typename Point, typename PointOfSegment> applyboost::geometry::strategy::within::winding_side_equal87 static inline int apply(Point const& point, 88 PointOfSegment const& se, 89 int count) 90 { 91 typedef typename coordinate_type<PointOfSegment>::type scoord_t; 92 typedef typename coordinate_system<PointOfSegment>::type::units units_t; 93 94 if (math::equals(get<1>(point), get<1>(se))) 95 return 0; 96 97 // Create a horizontal segment intersecting the original segment's endpoint 98 // equal to the point, with the derived direction (E/W). 99 PointOfSegment ss1, ss2; 100 set<1>(ss1, get<1>(se)); 101 set<0>(ss1, get<0>(se)); 102 set<1>(ss2, get<1>(se)); 103 scoord_t ss20 = get<0>(se); 104 if (count > 0) 105 { 106 ss20 += winding_small_angle<PointOfSegment>::apply(); 107 } 108 else 109 { 110 ss20 -= winding_small_angle<PointOfSegment>::apply(); 111 } 112 math::normalize_longitude<units_t>(ss20); 113 set<0>(ss2, ss20); 114 115 // Check the side using this vertical segment 116 return strategy_side_type::apply(ss1, ss2, point); 117 } 118 }; 119 // The optimization for cartesian 120 template <> 121 struct winding_side_equal<cartesian_tag> 122 { 123 template <typename Point, typename PointOfSegment> applyboost::geometry::strategy::within::winding_side_equal124 static inline int apply(Point const& point, 125 PointOfSegment const& se, 126 int count) 127 { 128 // NOTE: for D=0 the signs would be reversed 129 return math::equals(get<1>(point), get<1>(se)) ? 130 0 : 131 get<1>(point) < get<1>(se) ? 132 // assuming count is equal to 1 or -1 133 -count : // ( count > 0 ? -1 : 1) : 134 count; // ( count > 0 ? 1 : -1) ; 135 } 136 }; 137 138 139 template <typename Point, 140 typename CalculationType, 141 typename CSTag = typename cs_tag<Point>::type> 142 struct winding_check_touch 143 { 144 typedef CalculationType calc_t; 145 typedef typename coordinate_system<Point>::type::units units_t; 146 typedef math::detail::constants_on_spheroid<CalculationType, units_t> constants; 147 148 template <typename PointOfSegment, typename State> applyboost::geometry::strategy::within::winding_check_touch149 static inline int apply(Point const& point, 150 PointOfSegment const& seg1, 151 PointOfSegment const& seg2, 152 State& state, 153 bool& eq1, 154 bool& eq2) 155 { 156 calc_t const pi = constants::half_period(); 157 calc_t const pi2 = pi / calc_t(2); 158 159 calc_t const px = get<0>(point); 160 calc_t const s1x = get<0>(seg1); 161 calc_t const s2x = get<0>(seg2); 162 calc_t const py = get<1>(point); 163 calc_t const s1y = get<1>(seg1); 164 calc_t const s2y = get<1>(seg2); 165 166 // NOTE: lat in {-90, 90} and arbitrary lon 167 // it doesn't matter what lon it is if it's a pole 168 // so e.g. if one of the segment endpoints is a pole 169 // then only the other lon matters 170 171 bool eq1_strict = math::equals(s1x, px); 172 bool eq2_strict = math::equals(s2x, px); 173 174 eq1 = eq1_strict // lon strictly equal to s1 175 || math::equals(s1y, pi2) || math::equals(s1y, -pi2); // s1 is pole 176 eq2 = eq2_strict // lon strictly equal to s2 177 || math::equals(s2y, pi2) || math::equals(s2y, -pi2); // s2 is pole 178 179 // segment overlapping pole 180 calc_t s1x_anti = s1x + constants::half_period(); 181 math::normalize_longitude<units_t, calc_t>(s1x_anti); 182 bool antipodal = math::equals(s2x, s1x_anti); 183 if (antipodal) 184 { 185 eq1 = eq2 = eq1 || eq2; 186 187 // segment overlapping pole and point is pole 188 if (math::equals(py, pi2) || math::equals(py, -pi2)) 189 { 190 eq1 = eq2 = true; 191 } 192 } 193 194 // Both equal p -> segment vertical 195 // The only thing which has to be done is check if point is ON segment 196 if (eq1 && eq2) 197 { 198 // segment endpoints on the same sides of the globe 199 if (! antipodal 200 // p's lat between segment endpoints' lats 201 ? (s1y <= py && s2y >= py) || (s2y <= py && s1y >= py) 202 // going through north or south pole? 203 : (pi - s1y - s2y <= pi 204 ? (eq1_strict && s1y <= py) || (eq2_strict && s2y <= py) // north 205 || math::equals(py, pi2) // point on north pole 206 : (eq1_strict && s1y >= py) || (eq2_strict && s2y >= py)) // south 207 || math::equals(py, -pi2) // point on south pole 208 ) 209 { 210 state.m_touches = true; 211 } 212 return true; 213 } 214 return false; 215 } 216 }; 217 // The optimization for cartesian 218 template <typename Point, typename CalculationType> 219 struct winding_check_touch<Point, CalculationType, cartesian_tag> 220 { 221 typedef CalculationType calc_t; 222 223 template <typename PointOfSegment, typename State> applyboost::geometry::strategy::within::winding_check_touch224 static inline bool apply(Point const& point, 225 PointOfSegment const& seg1, 226 PointOfSegment const& seg2, 227 State& state, 228 bool& eq1, 229 bool& eq2) 230 { 231 calc_t const px = get<0>(point); 232 calc_t const s1x = get<0>(seg1); 233 calc_t const s2x = get<0>(seg2); 234 235 eq1 = math::equals(s1x, px); 236 eq2 = math::equals(s2x, px); 237 238 // Both equal p -> segment vertical 239 // The only thing which has to be done is check if point is ON segment 240 if (eq1 && eq2) 241 { 242 calc_t const py = get<1>(point); 243 calc_t const s1y = get<1>(seg1); 244 calc_t const s2y = get<1>(seg2); 245 if ((s1y <= py && s2y >= py) || (s2y <= py && s1y >= py)) 246 { 247 state.m_touches = true; 248 } 249 return true; 250 } 251 return false; 252 } 253 }; 254 255 256 // Called if point is not aligned with a vertical segment 257 template <typename Point, 258 typename CalculationType, 259 typename CSTag = typename cs_tag<Point>::type> 260 struct winding_calculate_count 261 { 262 typedef CalculationType calc_t; 263 typedef typename coordinate_system<Point>::type::units units_t; 264 greaterboost::geometry::strategy::within::winding_calculate_count265 static inline bool greater(calc_t const& l, calc_t const& r) 266 { 267 calc_t diff = l - r; 268 math::normalize_longitude<units_t, calc_t>(diff); 269 return diff > calc_t(0); 270 } 271 applyboost::geometry::strategy::within::winding_calculate_count272 static inline int apply(calc_t const& p, 273 calc_t const& s1, calc_t const& s2, 274 bool eq1, bool eq2) 275 { 276 // Probably could be optimized by avoiding normalization for some comparisons 277 // e.g. s1 > p could be calculated from p > s1 278 279 // If both segment endpoints were poles below checks wouldn't be enough 280 // but this means that either both are the same or that they are N/S poles 281 // and therefore the segment is not valid. 282 // If needed (eq1 && eq2 ? 0) could be returned 283 284 return 285 eq1 ? (greater(s2, p) ? 1 : -1) // Point on level s1, E/W depending on s2 286 : eq2 ? (greater(s1, p) ? -1 : 1) // idem 287 : greater(p, s1) && greater(s2, p) ? 2 // Point between s1 -> s2 --> E 288 : greater(p, s2) && greater(s1, p) ? -2 // Point between s2 -> s1 --> W 289 : 0; 290 } 291 }; 292 // The optimization for cartesian 293 template <typename Point, typename CalculationType> 294 struct winding_calculate_count<Point, CalculationType, cartesian_tag> 295 { 296 typedef CalculationType calc_t; 297 applyboost::geometry::strategy::within::winding_calculate_count298 static inline int apply(calc_t const& p, 299 calc_t const& s1, calc_t const& s2, 300 bool eq1, bool eq2) 301 { 302 return 303 eq1 ? (s2 > p ? 1 : -1) // Point on level s1, E/W depending on s2 304 : eq2 ? (s1 > p ? -1 : 1) // idem 305 : s1 < p && s2 > p ? 2 // Point between s1 -> s2 --> E 306 : s2 < p && s1 > p ? -2 // Point between s2 -> s1 --> W 307 : 0; 308 } 309 }; 310 311 312 /*! 313 \brief Within detection using winding rule 314 \ingroup strategies 315 \tparam Point \tparam_point 316 \tparam PointOfSegment \tparam_segment_point 317 \tparam SideStrategy Side strategy 318 \tparam CalculationType \tparam_calculation 319 \author Barend Gehrels 320 \note The implementation is inspired by terralib http://www.terralib.org (LGPL) 321 \note but totally revised afterwards, especially for cases on segments 322 \note Only dependant on "side", -> agnostic, suitable for spherical/latlong 323 324 \qbk{ 325 [heading See also] 326 [link geometry.reference.algorithms.within.within_3_with_strategy within (with strategy)] 327 } 328 */ 329 template 330 < 331 typename Point, 332 typename PointOfSegment = Point, 333 typename SideStrategy = typename strategy::side::services::default_strategy 334 < 335 typename cs_tag<Point>::type 336 >::type, 337 typename CalculationType = void 338 > 339 class winding 340 { 341 typedef typename select_calculation_type 342 < 343 Point, 344 PointOfSegment, 345 CalculationType 346 >::type calculation_type; 347 348 /*! subclass to keep state */ 349 class counter 350 { 351 int m_count; 352 bool m_touches; 353 code() const354 inline int code() const 355 { 356 return m_touches ? 0 : m_count == 0 ? -1 : 1; 357 } 358 359 public : 360 friend class winding; 361 362 template <typename P, typename CT, typename CST> 363 friend struct winding_check_touch; 364 counter()365 inline counter() 366 : m_count(0) 367 , m_touches(false) 368 {} 369 370 }; 371 check_segment(Point const & point,PointOfSegment const & seg1,PointOfSegment const & seg2,counter & state,bool & eq1,bool & eq2)372 static inline int check_segment(Point const& point, 373 PointOfSegment const& seg1, PointOfSegment const& seg2, 374 counter& state, bool& eq1, bool& eq2) 375 { 376 if (winding_check_touch<Point, calculation_type> 377 ::apply(point, seg1, seg2, state, eq1, eq2)) 378 { 379 return 0; 380 } 381 382 calculation_type const p = get<0>(point); 383 calculation_type const s1 = get<0>(seg1); 384 calculation_type const s2 = get<0>(seg2); 385 return winding_calculate_count<Point, calculation_type> 386 ::apply(p, s1, s2, eq1, eq2); 387 } 388 389 390 public: 391 typedef typename SideStrategy::envelope_strategy_type envelope_strategy_type; 392 get_envelope_strategy() const393 inline envelope_strategy_type get_envelope_strategy() const 394 { 395 return m_side_strategy.get_envelope_strategy(); 396 } 397 398 typedef typename SideStrategy::disjoint_strategy_type disjoint_strategy_type; 399 get_disjoint_strategy() const400 inline disjoint_strategy_type get_disjoint_strategy() const 401 { 402 return m_side_strategy.get_disjoint_strategy(); 403 } 404 winding()405 winding() 406 {} 407 winding(SideStrategy const & side_strategy)408 explicit winding(SideStrategy const& side_strategy) 409 : m_side_strategy(side_strategy) 410 {} 411 412 // Typedefs and static methods to fulfill the concept 413 typedef Point point_type; 414 typedef PointOfSegment segment_point_type; 415 typedef counter state_type; 416 apply(Point const & point,PointOfSegment const & s1,PointOfSegment const & s2,counter & state) const417 inline bool apply(Point const& point, 418 PointOfSegment const& s1, PointOfSegment const& s2, 419 counter& state) const 420 { 421 typedef typename cs_tag<Point>::type cs_t; 422 423 bool eq1 = false; 424 bool eq2 = false; 425 boost::ignore_unused(eq2); 426 427 int count = check_segment(point, s1, s2, state, eq1, eq2); 428 if (count != 0) 429 { 430 int side = 0; 431 if (count == 1 || count == -1) 432 { 433 side = winding_side_equal<cs_t>::apply(point, eq1 ? s1 : s2, count); 434 } 435 else // count == 2 || count == -2 436 { 437 // 1 left, -1 right 438 side = m_side_strategy.apply(s1, s2, point); 439 } 440 441 if (side == 0) 442 { 443 // Point is lying on segment 444 state.m_touches = true; 445 state.m_count = 0; 446 return false; 447 } 448 449 // Side is NEG for right, POS for left. 450 // The count is -2 for down, 2 for up (or -1/1) 451 // Side positive thus means UP and LEFTSIDE or DOWN and RIGHTSIDE 452 // See accompagnying figure (TODO) 453 if (side * count > 0) 454 { 455 state.m_count += count; 456 } 457 } 458 return ! state.m_touches; 459 } 460 result(counter const & state)461 static inline int result(counter const& state) 462 { 463 return state.code(); 464 } 465 466 private: 467 SideStrategy m_side_strategy; 468 }; 469 470 471 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS 472 473 namespace services 474 { 475 476 template <typename PointLike, typename Geometry, typename AnyTag1, typename AnyTag2> 477 struct default_strategy<PointLike, Geometry, AnyTag1, AnyTag2, pointlike_tag, polygonal_tag, cartesian_tag, cartesian_tag> 478 { 479 typedef winding 480 < 481 typename geometry::point_type<PointLike>::type, 482 typename geometry::point_type<Geometry>::type 483 > type; 484 }; 485 486 template <typename PointLike, typename Geometry, typename AnyTag1, typename AnyTag2> 487 struct default_strategy<PointLike, Geometry, AnyTag1, AnyTag2, pointlike_tag, polygonal_tag, spherical_tag, spherical_tag> 488 { 489 typedef winding 490 < 491 typename geometry::point_type<PointLike>::type, 492 typename geometry::point_type<Geometry>::type 493 > type; 494 }; 495 496 template <typename PointLike, typename Geometry, typename AnyTag1, typename AnyTag2> 497 struct default_strategy<PointLike, Geometry, AnyTag1, AnyTag2, pointlike_tag, linear_tag, cartesian_tag, cartesian_tag> 498 { 499 typedef winding 500 < 501 typename geometry::point_type<PointLike>::type, 502 typename geometry::point_type<Geometry>::type 503 > type; 504 }; 505 506 template <typename PointLike, typename Geometry, typename AnyTag1, typename AnyTag2> 507 struct default_strategy<PointLike, Geometry, AnyTag1, AnyTag2, pointlike_tag, linear_tag, spherical_tag, spherical_tag> 508 { 509 typedef winding 510 < 511 typename geometry::point_type<PointLike>::type, 512 typename geometry::point_type<Geometry>::type 513 > type; 514 }; 515 516 } // namespace services 517 518 #endif 519 520 521 }} // namespace strategy::within 522 523 524 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS 525 namespace strategy { namespace covered_by { namespace services 526 { 527 528 template <typename PointLike, typename Geometry, typename AnyTag1, typename AnyTag2> 529 struct default_strategy<PointLike, Geometry, AnyTag1, AnyTag2, pointlike_tag, polygonal_tag, cartesian_tag, cartesian_tag> 530 { 531 typedef within::winding 532 < 533 typename geometry::point_type<PointLike>::type, 534 typename geometry::point_type<Geometry>::type 535 > type; 536 }; 537 538 template <typename PointLike, typename Geometry, typename AnyTag1, typename AnyTag2> 539 struct default_strategy<PointLike, Geometry, AnyTag1, AnyTag2, pointlike_tag, polygonal_tag, spherical_tag, spherical_tag> 540 { 541 typedef within::winding 542 < 543 typename geometry::point_type<PointLike>::type, 544 typename geometry::point_type<Geometry>::type 545 > type; 546 }; 547 548 template <typename PointLike, typename Geometry, typename AnyTag1, typename AnyTag2> 549 struct default_strategy<PointLike, Geometry, AnyTag1, AnyTag2, pointlike_tag, linear_tag, cartesian_tag, cartesian_tag> 550 { 551 typedef within::winding 552 < 553 typename geometry::point_type<PointLike>::type, 554 typename geometry::point_type<Geometry>::type 555 > type; 556 }; 557 558 template <typename PointLike, typename Geometry, typename AnyTag1, typename AnyTag2> 559 struct default_strategy<PointLike, Geometry, AnyTag1, AnyTag2, pointlike_tag, linear_tag, spherical_tag, spherical_tag> 560 { 561 typedef within::winding 562 < 563 typename geometry::point_type<PointLike>::type, 564 typename geometry::point_type<Geometry>::type 565 > type; 566 }; 567 568 }}} // namespace strategy::covered_by::services 569 #endif 570 571 572 }} // namespace boost::geometry 573 574 575 #endif // BOOST_GEOMETRY_STRATEGY_AGNOSTIC_POINT_IN_POLY_WINDING_HPP 576