strategies/spherical/point_in_poly_winding.hpp

0001 // Boost.Geometry (aka GGL, Generic Geometry Library)
0002
0003 // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
0004 // Copyright (c) 2013-2017 Adam Wulkiewicz, Lodz, Poland.
0005
0006 // This file was modified by Oracle on 2013-2023.
0007 // Modifications copyright (c) 2013-2023 Oracle and/or its affiliates.
0008 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
0009 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0010
0011 // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
0012 // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
0013
0014 // Use, modification and distribution is subject to the Boost Software License,
0015 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0016 // http://www.boost.org/LICENSE_1_0.txt)
0017
0018 #ifndef BOOST_GEOMETRY_STRATEGY_SPHERICAL_POINT_IN_POLY_WINDING_HPP
0019 #define BOOST_GEOMETRY_STRATEGY_SPHERICAL_POINT_IN_POLY_WINDING_HPP
0020
0021
0022 #include <boost/geometry/core/access.hpp>
0023 #include <boost/geometry/core/coordinate_system.hpp>
0024 #include <boost/geometry/core/cs.hpp>
0025 #include <boost/geometry/core/tags.hpp>
0026
0027 #include <boost/geometry/util/math.hpp>
0028 #include <boost/geometry/util/select_calculation_type.hpp>
0029 #include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>
0030
0031 #include <boost/geometry/strategy/spherical/expand_point.hpp>
0032
0033 #include <boost/geometry/strategies/cartesian/point_in_box.hpp>
0034 #include <boost/geometry/strategies/covered_by.hpp>
0035 #include <boost/geometry/strategies/side.hpp>
0036 #include <boost/geometry/strategies/spherical/disjoint_box_box.hpp>
0037 #include <boost/geometry/strategies/spherical/ssf.hpp>
0038 #include <boost/geometry/strategies/within.hpp>
0039
0040
0041 namespace boost { namespace geometry
0042 {
0043
0044 namespace strategy { namespace within
0045 {
0046
0047
0048 #ifndef DOXYGEN_NO_DETAIL
0049 namespace detail
0050 {
0051
0052 template <typename SideStrategy, typename CalculationType>
0053 class spherical_winding_base
0054 {
0055     template <typename Point, typename PointOfSegment>
0056     struct calculation_type
0057         : select_calculation_type
0058             <
0059                 Point,
0060                 PointOfSegment,
0061                 CalculationType
0062             >
0063     {};
0064
0065     /*! subclass to keep state */
0066     class counter
0067     {
0068         int m_count;
0069         //int m_count_n;
0070         int m_count_s;
0071         int m_raw_count;
0072         int m_raw_count_anti;
0073         bool m_touches;
0074
0075         inline int code() const
0076         {
0077             if (m_touches)
0078             {
0079                 return 0;
0080             }
0081
0082             if (m_raw_count != 0 && m_raw_count_anti != 0)
0083             {
0084                 if (m_raw_count > 0) // right, wrap around south pole
0085                 {
0086                     return (m_count + m_count_s) == 0 ? -1 : 1;
0087                 }
0088                 else // left, wrap around north pole
0089                 {
0090                     //return (m_count + m_count_n) == 0 ? -1 : 1;
0091                     // m_count_n is 0
0092                     return m_count == 0 ? -1 : 1;
0093                 }
0094             }
0095
0096             return m_count == 0 ? -1 : 1;
0097         }
0098
0099     public :
0100         friend class spherical_winding_base;
0101
0102         inline counter()
0103             : m_count(0)
0104             //, m_count_n(0)
0105             , m_count_s(0)
0106             , m_raw_count(0)
0107             , m_raw_count_anti(0)
0108             , m_touches(false)
0109         {}
0110
0111     };
0112
0113     struct count_info
0114     {
0115         explicit count_info(int c = 0, bool ia = false)
0116             : count(c)
0117             , is_anti(ia)
0118         {}
0119
0120         int count;
0121         bool is_anti;
0122     };
0123
0124 public:
0125     typedef typename SideStrategy::cs_tag cs_tag;
0126
0127     spherical_winding_base() = default;
0128
0129     template <typename Model>
0130     explicit spherical_winding_base(Model const& model)
0131         : m_side_strategy(model)
0132     {}
0133
0134     // Typedefs and static methods to fulfill the concept
0135     typedef counter state_type;
0136
0137     template <typename Point, typename PointOfSegment>
0138     inline bool apply(Point const& point,
0139                       PointOfSegment const& s1, PointOfSegment const& s2,
0140                       counter& state) const
0141     {
0142         typedef typename calculation_type<Point, PointOfSegment>::type calc_t;
0143         typedef typename geometry::detail::cs_angular_units<Point>::type units_t;
0144         typedef math::detail::constants_on_spheroid<calc_t, units_t> constants;
0145
0146         bool eq1 = false;
0147         bool eq2 = false;
0148         bool s_antipodal = false;
0149
0150         count_info ci = check_segment(point, s1, s2, state, eq1, eq2, s_antipodal);
0151         if (ci.count != 0)
0152         {
0153             if (! ci.is_anti)
0154             {
0155                 int side = 0;
0156                 if (ci.count == 1 || ci.count == -1)
0157                 {
0158                     side = side_equal(point, eq1 ? s1 : s2, ci);
0159                 }
0160                 else // count == 2 || count == -2
0161                 {
0162                     if (! s_antipodal)
0163                     {
0164                         // 1 left, -1 right
0165                         side = m_side_strategy.apply(s1, s2, point);
0166                     }
0167                     else
0168                     {
0169                         calc_t const pi = constants::half_period();
0170                         calc_t const s1_lat = get<1>(s1);
0171                         calc_t const s2_lat = get<1>(s2);
0172
0173                         side = math::sign(ci.count)
0174                              * (pi - s1_lat - s2_lat <= pi // segment goes through north pole
0175                                 ? -1 // going right all points will be on right side
0176                                 : 1); // going right all points will be on left side
0177                     }
0178                 }
0179
0180                 if (side == 0)
0181                 {
0182                     // Point is lying on segment
0183                     state.m_touches = true;
0184                     state.m_count = 0;
0185                     return false;
0186                 }
0187
0188                 // Side is NEG for right, POS for left.
0189                 // The count is -2 for left, 2 for right (or -1/1)
0190                 // Side positive thus means RIGHT and LEFTSIDE or LEFT and RIGHTSIDE
0191                 // See accompagnying figure (TODO)
0192                 if (side * ci.count > 0)
0193                 {
0194                     state.m_count += ci.count;
0195                 }
0196
0197                 state.m_raw_count += ci.count;
0198             }
0199             else
0200             {
0201                 // Count negated because the segment is on the other side of the globe
0202                 // so it is reversed to match this side of the globe
0203
0204                 // Assuming geometry wraps around north pole, for segments on the other side of the globe
0205                 // the point will always be RIGHT+RIGHTSIDE or LEFT+LEFTSIDE, so side*-count always < 0
0206                 //state.m_count_n -= 0;
0207
0208                 // Assuming geometry wraps around south pole, for segments on the other side of the globe
0209                 // the point will always be RIGHT+LEFTSIDE or LEFT+RIGHTSIDE, so side*-count always > 0
0210                 state.m_count_s -= ci.count;
0211
0212                 state.m_raw_count_anti -= ci.count;
0213             }
0214         }
0215         return ! state.m_touches;
0216     }
0217
0218     static inline int result(counter const& state)
0219     {
0220         return state.code();
0221     }
0222
0223 protected:
0224     template <typename Point, typename PointOfSegment>
0225     static inline count_info check_segment(Point const& point,
0226                                     PointOfSegment const& seg1,
0227                                     PointOfSegment const& seg2,
0228                                     counter& state,
0229                                     bool& eq1, bool& eq2, bool& s_antipodal)
0230     {
0231         if (check_touch(point, seg1, seg2, state, eq1, eq2, s_antipodal))
0232         {
0233             return count_info(0, false);
0234         }
0235
0236         return calculate_count(point, seg1, seg2, eq1, eq2, s_antipodal);
0237     }
0238
0239     template <typename Point, typename PointOfSegment>
0240     static inline int check_touch(Point const& point,
0241                                   PointOfSegment const& seg1,
0242                                   PointOfSegment const& seg2,
0243                                   counter& state,
0244                                   bool& eq1,
0245                                   bool& eq2,
0246                                   bool& s_antipodal)
0247     {
0248         typedef typename calculation_type<Point, PointOfSegment>::type calc_t;
0249         typedef typename geometry::detail::cs_angular_units<Point>::type units_t;
0250         typedef math::detail::constants_on_spheroid<calc_t, units_t> constants;
0251
0252         calc_t const c0 = 0;
0253         calc_t const c2 = 2;
0254         calc_t const pi = constants::half_period();
0255         calc_t const half_pi = pi / c2;
0256
0257         calc_t const p_lon = get<0>(point);
0258         calc_t const s1_lon = get<0>(seg1);
0259         calc_t const s2_lon = get<0>(seg2);
0260         calc_t const p_lat = get<1>(point);
0261         calc_t const s1_lat = get<1>(seg1);
0262         calc_t const s2_lat = get<1>(seg2);
0263
0264         // NOTE: lat in {-90, 90} and arbitrary lon
0265         //  it doesn't matter what lon it is if it's a pole
0266         //  so e.g. if one of the segment endpoints is a pole
0267         //  then only the other lon matters
0268
0269         bool eq1_strict = longitudes_equal<units_t>(s1_lon, p_lon);
0270         bool eq2_strict = longitudes_equal<units_t>(s2_lon, p_lon);
0271         bool eq1_anti = false;
0272         bool eq2_anti = false;
0273
0274         calc_t const anti_p_lon = p_lon + (p_lon <= c0 ? pi : -pi);
0275
0276         eq1 = eq1_strict // lon strictly equal to s1
0277             || (eq1_anti = longitudes_equal<units_t>(s1_lon, anti_p_lon)); // anti-lon strictly equal to s1
0278         eq2 = eq2_strict // lon strictly equal to s2
0279             || (eq2_anti = longitudes_equal<units_t>(s2_lon, anti_p_lon)); // anti-lon strictly equal to s2
0280
0281         // segment overlapping pole
0282         calc_t const s_lon_diff = math::longitude_distance_signed<units_t>(s1_lon, s2_lon);
0283         s_antipodal = math::equals(s_lon_diff, pi);
0284         if (s_antipodal)
0285         {
0286             eq1 = eq2 = eq1 || eq2;
0287
0288             // segment overlapping pole and point is pole
0289             if (math::equals(math::abs(p_lat), half_pi))
0290             {
0291                 eq1 = eq2 = true;
0292             }
0293         }
0294
0295         // check whether point is on a segment with a pole endpoint
0296         if (math::longitude_distance_signed<units_t>(s2_lon, p_lon) == c0)
0297         {
0298             bool const s1_north = math::equals(get<1>(seg1), half_pi);
0299             bool const s1_south = math::equals(get<1>(seg1), -half_pi);
0300             if (s1_north || s1_south)
0301             {
0302                 state.m_touches = s1_south ? s2_lat > p_lat : s2_lat < p_lat;
0303                 return state.m_touches;
0304             }
0305         }
0306         if (math::longitude_distance_signed<units_t>(s1_lon, p_lon) == c0)
0307         {
0308             bool const s2_north = math::equals(get<1>(seg2), half_pi);
0309             bool const s2_south = math::equals(get<1>(seg2), -half_pi);
0310             if (s2_north || s2_south)
0311             {
0312                 state.m_touches = s2_south ? s1_lat > p_lat : s1_lat < p_lat;
0313                 return state.m_touches;
0314             }
0315         }
0316
0317         // Both equal p -> segment vertical
0318         // The only thing which has to be done is check if point is ON segment
0319         if (eq1 && eq2)
0320         {
0321             // segment endpoints on the same sides of the globe
0322             if (! s_antipodal)
0323             {
0324                 // p's lat between segment endpoints' lats
0325                 if ( (s1_lat <= p_lat && s2_lat >= p_lat) || (s2_lat <= p_lat && s1_lat >= p_lat) )
0326                 {
0327                     if (!eq1_anti || !eq2_anti)
0328                     {
0329                         state.m_touches = true;
0330                     }
0331                 }
0332             }
0333             else
0334             {
0335                 // going through north or south pole?
0336                 if (pi - s1_lat - s2_lat <= pi)
0337                 {
0338                     if ( (eq1_strict && s1_lat <= p_lat) || (eq2_strict && s2_lat <= p_lat) // north
0339                             || math::equals(p_lat, half_pi) ) // point on north pole
0340                     {
0341                         state.m_touches = true;
0342                     }
0343                     else if (! eq1_strict && ! eq2_strict && math::equals(p_lat, -half_pi) ) // point on south pole
0344                     {
0345                         return false;
0346                     }
0347                 }
0348                 else // south pole
0349                 {
0350                     if ( (eq1_strict && s1_lat >= p_lat) || (eq2_strict && s2_lat >= p_lat) // south
0351                             || math::equals(p_lat, -half_pi) ) // point on south pole
0352                     {
0353                         state.m_touches = true;
0354                     }
0355                     else if (! eq1_strict && ! eq2_strict && math::equals(p_lat, half_pi) ) // point on north pole
0356                     {
0357                         return false;
0358                     }
0359                 }
0360             }
0361
0362             return true;
0363         }
0364
0365         return false;
0366     }
0367
0368     // Called if point is not aligned with a vertical segment
0369     template <typename Point, typename PointOfSegment>
0370     static inline count_info calculate_count(Point const& point,
0371                                              PointOfSegment const& seg1,
0372                                              PointOfSegment const& seg2,
0373                                              bool eq1, bool eq2, bool s_antipodal)
0374     {
0375         typedef typename calculation_type<Point, PointOfSegment>::type calc_t;
0376         typedef typename geometry::detail::cs_angular_units<Point>::type units_t;
0377         typedef math::detail::constants_on_spheroid<calc_t, units_t> constants;
0378
0379         // If both segment endpoints were poles below checks wouldn't be enough
0380         // but this means that either both are the same or that they are N/S poles
0381         // and therefore the segment is not valid.
0382         // If needed (eq1 && eq2 ? 0) could be returned
0383
0384         calc_t const c0 = 0;
0385         calc_t const c2 = 2;
0386         calc_t const pi = constants::half_period();
0387         calc_t const half_pi = pi / c2;
0388
0389         bool const s1_is_pole = math::equals(std::abs(get<1>(seg1)), half_pi);
0390         bool const s2_is_pole = math::equals(std::abs(get<1>(seg2)), half_pi);
0391
0392         if (s1_is_pole && s2_is_pole)
0393         {
0394             return count_info(0, false);
0395         }
0396
0397         calc_t const p = get<0>(point);
0398         calc_t const s1 = get<0>(seg1);
0399         calc_t const s2 = get<0>(seg2);
0400
0401         calc_t const s1_p = math::longitude_distance_signed<units_t>(s1, p);
0402
0403         if (s_antipodal)
0404         {
0405             return count_info(s1_p < c0 ? -2 : 2, false); // choose W/E
0406         }
0407
0408         calc_t const s1_s2 = math::longitude_distance_signed<units_t>(s1, s2);
0409
0410         if (eq1 || eq2) // Point on level s1 or s2
0411         {
0412             return count_info(s1_s2 < c0 ? -1 : 1, // choose W/E
0413                               longitudes_equal<units_t>(p + pi, (eq1 ? s1 : s2)));
0414         }
0415
0416         // Point between s1 and s2
0417         if ( math::sign(s1_p) == math::sign(s1_s2)
0418           && math::abs(s1_p) < math::abs(s1_s2) )
0419         {
0420             return count_info(s1_s2 < c0 ? -2 : 2, false); // choose W/E
0421         }
0422
0423         calc_t const s1_p_anti = math::longitude_distance_signed<units_t>(s1, p + pi);
0424
0425         // Anti-Point between s1 and s2
0426         if ( math::sign(s1_p_anti) == math::sign(s1_s2)
0427           && math::abs(s1_p_anti) < math::abs(s1_s2) )
0428         {
0429             return count_info(s1_s2 < c0 ? -2 : 2, true); // choose W/E
0430         }
0431
0432         return count_info(0, false);
0433     }
0434
0435
0436     // Fix for https://svn.boost.org/trac/boost/ticket/9628
0437     // For floating point coordinates, the <D> coordinate of a point is compared
0438     // with the segment's points using some EPS. If the coordinates are "equal"
0439     // the sides are calculated. Therefore we can treat a segment as a long areal
0440     // geometry having some width. There is a small ~triangular area somewhere
0441     // between the segment's effective area and a segment's line used in sides
0442     // calculation where the segment is on the one side of the line but on the
0443     // other side of a segment (due to the width).
0444     // Below picture assuming D = 1, if D = 0 horiz<->vert, E<->N, RIGHT<->UP.
0445     // For the s1 of a segment going NE the real side is RIGHT but the point may
0446     // be detected as LEFT, like this:
0447     //                     RIGHT
0448     //                 ___----->
0449     //                  ^      O Pt  __ __
0450     //                 EPS     __ __
0451     //                  v__ __ BUT DETECTED AS LEFT OF THIS LINE
0452     //             _____7
0453     //       _____/
0454     // _____/
0455     // In the code below actually D = 0, so segments are nearly-vertical
0456     // Called when the point is on the same level as one of the segment's points
0457     // but the point is not aligned with a vertical segment
0458     template <typename Point, typename PointOfSegment>
0459     inline int side_equal(Point const& point,
0460                           PointOfSegment const& se,
0461                           count_info const& ci) const
0462     {
0463         typedef typename coordinate_type<PointOfSegment>::type scoord_t;
0464         typedef typename geometry::detail::cs_angular_units<Point>::type units_t;
0465
0466         if (math::equals(get<1>(point), get<1>(se)))
0467         {
0468             return 0;
0469         }
0470
0471         // Create a horizontal segment intersecting the original segment's endpoint
0472         // equal to the point, with the derived direction (E/W).
0473         PointOfSegment ss1, ss2;
0474         set<1>(ss1, get<1>(se));
0475         set<0>(ss1, get<0>(se));
0476         set<1>(ss2, get<1>(se));
0477         scoord_t ss20 = get<0>(se);
0478         if (ci.count > 0)
0479         {
0480             ss20 += small_angle<scoord_t, units_t>();
0481         }
0482         else
0483         {
0484             ss20 -= small_angle<scoord_t, units_t>();
0485         }
0486         math::normalize_longitude<units_t>(ss20);
0487         set<0>(ss2, ss20);
0488
0489         // Check the side using this vertical segment
0490         return m_side_strategy.apply(ss1, ss2, point);
0491     }
0492
0493     // 1 deg or pi/180 rad
0494     template <typename CalcT, typename Units>
0495     static inline CalcT small_angle()
0496     {
0497         typedef math::detail::constants_on_spheroid<CalcT, Units> constants;
0498
0499         return constants::half_period() / CalcT(180);
0500     }
0501
0502     template <typename Units, typename CalcT>
0503     static inline bool longitudes_equal(CalcT const& lon1, CalcT const& lon2)
0504     {
0505         return math::equals(
0506                 math::longitude_distance_signed<Units>(lon1, lon2),
0507                 CalcT(0));
0508     }
0509
0510     SideStrategy m_side_strategy;
0511 };
0512
0513
0514 } // namespace detail
0515 #endif // DOXYGEN_NO_DETAIL
0516
0517
0518 /*!
0519 \brief Within detection using winding rule in spherical coordinate system.
0520 \ingroup strategies
0521 \tparam Point \tparam_point
0522 \tparam PointOfSegment \tparam_segment_point
0523 \tparam CalculationType \tparam_calculation
0524
0525 \qbk{
0526 [heading See also]
0527 [link geometry.reference.algorithms.within.within_3_with_strategy within (with strategy)]
0528 }
0529  */
0530 template
0531 <
0532     typename Point = void, // for backward compatibility
0533     typename PointOfSegment = Point, // for backward compatibility
0534     typename CalculationType = void
0535 >
0536 class spherical_winding
0537     : public within::detail::spherical_winding_base
0538         <
0539             side::spherical_side_formula<CalculationType>,
0540             CalculationType
0541         >
0542 {};
0543
0544
0545 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
0546
0547 namespace services
0548 {
0549
0550 template <typename PointLike, typename Geometry, typename AnyTag1, typename AnyTag2>
0551 struct default_strategy<PointLike, Geometry, AnyTag1, AnyTag2, pointlike_tag, polygonal_tag, spherical_tag, spherical_tag>
0552 {
0553     typedef within::spherical_winding<> type;
0554 };
0555
0556 template <typename PointLike, typename Geometry, typename AnyTag1, typename AnyTag2>
0557 struct default_strategy<PointLike, Geometry, AnyTag1, AnyTag2, pointlike_tag, linear_tag, spherical_tag, spherical_tag>
0558 {
0559     typedef within::spherical_winding<> type;
0560 };
0561
0562 } // namespace services
0563
0564 #endif
0565
0566
0567 }} // namespace strategy::within
0568
0569
0570 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
0571 namespace strategy { namespace covered_by { namespace services
0572 {
0573
0574 template <typename PointLike, typename Geometry, typename AnyTag1, typename AnyTag2>
0575 struct default_strategy<PointLike, Geometry, AnyTag1, AnyTag2, pointlike_tag, polygonal_tag, spherical_tag, spherical_tag>
0576 {
0577     typedef within::spherical_winding<> type;
0578 };
0579
0580 template <typename PointLike, typename Geometry, typename AnyTag1, typename AnyTag2>
0581 struct default_strategy<PointLike, Geometry, AnyTag1, AnyTag2, pointlike_tag, linear_tag, spherical_tag, spherical_tag>
0582 {
0583     typedef within::spherical_winding<> type;
0584 };
0585
0586 }}} // namespace strategy::covered_by::services
0587 #endif
0588
0589
0590 }} // namespace boost::geometry
0591
0592
0593 #endif // BOOST_GEOMETRY_STRATEGY_SPHERICAL_POINT_IN_POLY_WINDING_HPP