boost/geometry/strategies/agnostic/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_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
  49     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>
  87     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>
 124     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>
 149     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>
 224     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
 265     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
 272     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
 298     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
 354         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
 365         inline counter()
 366             : m_count(0)
 367             , m_touches(false)
 368         {}
 369
 370     };
 371
 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     winding()
 392     {}
 393
 394     explicit winding(SideStrategy const& side_strategy)
 395         : m_side_strategy(side_strategy)
 396     {}
 397
 398     // Typedefs and static methods to fulfill the concept
 399     typedef Point point_type;
 400     typedef PointOfSegment segment_point_type;
 401     typedef counter state_type;
 402
 403     inline bool apply(Point const& point,
 404                       PointOfSegment const& s1, PointOfSegment const& s2,
 405                       counter& state) const
 406     {
 407         typedef typename cs_tag<Point>::type cs_t;
 408
 409         bool eq1 = false;
 410         bool eq2 = false;
 411         boost::ignore_unused(eq2);
 412
 413         int count = check_segment(point, s1, s2, state, eq1, eq2);
 414         if (count != 0)
 415         {
 416             int side = 0;
 417             if (count == 1 || count == -1)
 418             {
 419                 side = winding_side_equal<cs_t>::apply(point, eq1 ? s1 : s2, count);
 420             }
 421             else // count == 2 || count == -2
 422             {
 423                 // 1 left, -1 right
 424                 side = m_side_strategy.apply(s1, s2, point);
 425             }
 426
 427             if (side == 0)
 428             {
 429                 // Point is lying on segment
 430                 state.m_touches = true;
 431                 state.m_count = 0;
 432                 return false;
 433             }
 434
 435             // Side is NEG for right, POS for left.
 436             // The count is -2 for down, 2 for up (or -1/1)
 437             // Side positive thus means UP and LEFTSIDE or DOWN and RIGHTSIDE
 438             // See accompagnying figure (TODO)
 439             if (side * count > 0)
 440             {
 441                 state.m_count += count;
 442             }
 443         }
 444         return ! state.m_touches;
 445     }
 446
 447     static inline int result(counter const& state)
 448     {
 449         return state.code();
 450     }
 451
 452 private:
 453     SideStrategy m_side_strategy;
 454 };
 455
 456
 457 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
 458
 459 namespace services
 460 {
 461
 462 template <typename Point, typename Geometry, typename AnyTag>
 463 struct default_strategy<Point, Geometry, point_tag, AnyTag, pointlike_tag, polygonal_tag, cartesian_tag, cartesian_tag>
 464 {
 465     typedef winding<Point, typename geometry::point_type<Geometry>::type> type;
 466 };
 467
 468 template <typename Point, typename Geometry, typename AnyTag>
 469 struct default_strategy<Point, Geometry, point_tag, AnyTag, pointlike_tag, polygonal_tag, spherical_tag, spherical_tag>
 470 {
 471     typedef winding<Point, typename geometry::point_type<Geometry>::type> type;
 472 };
 473
 474 template <typename Point, typename Geometry, typename AnyTag>
 475 struct default_strategy<Point, Geometry, point_tag, AnyTag, pointlike_tag, linear_tag, cartesian_tag, cartesian_tag>
 476 {
 477     typedef winding<Point, typename geometry::point_type<Geometry>::type> type;
 478 };
 479
 480 template <typename Point, typename Geometry, typename AnyTag>
 481 struct default_strategy<Point, Geometry, point_tag, AnyTag, pointlike_tag, linear_tag, spherical_tag, spherical_tag>
 482 {
 483     typedef winding<Point, typename geometry::point_type<Geometry>::type> type;
 484 };
 485
 486 } // namespace services
 487
 488 #endif
 489
 490
 491 }} // namespace strategy::within
 492
 493
 494 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
 495 namespace strategy { namespace covered_by { namespace services
 496 {
 497
 498 template <typename Point, typename Geometry, typename AnyTag>
 499 struct default_strategy<Point, Geometry, point_tag, AnyTag, pointlike_tag, polygonal_tag, cartesian_tag, cartesian_tag>
 500 {
 501     typedef strategy::within::winding<Point, typename geometry::point_type<Geometry>::type> type;
 502 };
 503
 504 template <typename Point, typename Geometry, typename AnyTag>
 505 struct default_strategy<Point, Geometry, point_tag, AnyTag, pointlike_tag, polygonal_tag, spherical_tag, spherical_tag>
 506 {
 507     typedef strategy::within::winding<Point, typename geometry::point_type<Geometry>::type> type;
 508 };
 509
 510 template <typename Point, typename Geometry, typename AnyTag>
 511 struct default_strategy<Point, Geometry, point_tag, AnyTag, pointlike_tag, linear_tag, cartesian_tag, cartesian_tag>
 512 {
 513     typedef strategy::within::winding<Point, typename geometry::point_type<Geometry>::type> type;
 514 };
 515
 516 template <typename Point, typename Geometry, typename AnyTag>
 517 struct default_strategy<Point, Geometry, point_tag, AnyTag, pointlike_tag, linear_tag, spherical_tag, spherical_tag>
 518 {
 519     typedef strategy::within::winding<Point, typename geometry::point_type<Geometry>::type> type;
 520 };
 521
 522 }}} // namespace strategy::covered_by::services
 523 #endif
 524
 525
 526 }} // namespace boost::geometry
 527
 528
 529 #endif // BOOST_GEOMETRY_STRATEGY_AGNOSTIC_POINT_IN_POLY_WINDING_HPP