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