ceph/src/boost/boost/geometry/strategies/spherical/point_in_poly_winding.hpp

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