ceph/src/boost/libs/geometry/include/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.
   7 // Modifications copyright (c) 2013-2016 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 CalculationType \tparam_calculation
 318 \author Barend Gehrels
 319 \note The implementation is inspired by terralib http://www.terralib.org (LGPL)
 320 \note but totally revised afterwards, especially for cases on segments
 321 \note Only dependant on "side", -> agnostic, suitable for spherical/latlong
 322
 323 \qbk{
 324 [heading See also]
 325 [link geometry.reference.algorithms.within.within_3_with_strategy within (with strategy)]
 326 }
 327  */
 328 template
 329 <
 330     typename Point,
 331     typename PointOfSegment = Point,
 332     typename CalculationType = void
 333 >
 334 class winding
 335 {
 336     typedef typename select_calculation_type
 337         <
 338             Point,
 339             PointOfSegment,
 340             CalculationType
 341         >::type calculation_type;
 342
 343
 344     typedef typename strategy::side::services::default_strategy
 345         <
 346             typename cs_tag<Point>::type
 347         >::type strategy_side_type;
 348
 349
 350     /*! subclass to keep state */
 351     class counter
 352     {
 353         int m_count;
 354         bool m_touches;
 355
 356         inline int code() const
 357         {
 358             return m_touches ? 0 : m_count == 0 ? -1 : 1;
 359         }
 360
 361     public :
 362         friend class winding;
 363
 364         template <typename P, typename CT, typename CST>
 365         friend struct winding_check_touch;
 366
 367         inline counter()
 368             : m_count(0)
 369             , m_touches(false)
 370         {}
 371
 372     };
 373
 374
 375     static inline int check_segment(Point const& point,
 376                 PointOfSegment const& seg1, PointOfSegment const& seg2,
 377                 counter& state, bool& eq1, bool& eq2)
 378     {
 379         if (winding_check_touch<Point, calculation_type>
 380                 ::apply(point, seg1, seg2, state, eq1, eq2))
 381         {
 382             return 0;
 383         }
 384
 385         calculation_type const p = get<0>(point);
 386         calculation_type const s1 = get<0>(seg1);
 387         calculation_type const s2 = get<0>(seg2);
 388         return winding_calculate_count<Point, calculation_type>
 389                     ::apply(p, s1, s2, eq1, eq2);
 390     }
 391
 392
 393 public :
 394
 395     // Typedefs and static methods to fulfill the concept
 396     typedef Point point_type;
 397     typedef PointOfSegment segment_point_type;
 398     typedef counter state_type;
 399
 400     static inline bool apply(Point const& point,
 401                 PointOfSegment const& s1, PointOfSegment const& s2,
 402                 counter& state)
 403     {
 404         typedef typename cs_tag<Point>::type cs_t;
 405
 406         bool eq1 = false;
 407         bool eq2 = false;
 408         boost::ignore_unused(eq2);
 409
 410         int count = check_segment(point, s1, s2, state, eq1, eq2);
 411         if (count != 0)
 412         {
 413             int side = 0;
 414             if (count == 1 || count == -1)
 415             {
 416                 side = winding_side_equal<cs_t>::apply(point, eq1 ? s1 : s2, count);
 417             }
 418             else // count == 2 || count == -2
 419             {
 420                 // 1 left, -1 right
 421                 side = strategy_side_type::apply(s1, s2, point);
 422             }
 423
 424             if (side == 0)
 425             {
 426                 // Point is lying on segment
 427                 state.m_touches = true;
 428                 state.m_count = 0;
 429                 return false;
 430             }
 431
 432             // Side is NEG for right, POS for left.
 433             // The count is -2 for down, 2 for up (or -1/1)
 434             // Side positive thus means UP and LEFTSIDE or DOWN and RIGHTSIDE
 435             // See accompagnying figure (TODO)
 436             if (side * count > 0)
 437             {
 438                 state.m_count += count;
 439             }
 440         }
 441         return ! state.m_touches;
 442     }
 443
 444     static inline int result(counter const& state)
 445     {
 446         return state.code();
 447     }
 448 };
 449
 450
 451 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
 452
 453 namespace services
 454 {
 455
 456 // Register using "areal_tag" for ring, polygon, multi-polygon
 457 template <typename AnyTag, typename Point, typename Geometry>
 458 struct default_strategy<point_tag, AnyTag, point_tag, areal_tag, cartesian_tag, cartesian_tag, Point, Geometry>
 459 {
 460     typedef winding<Point, typename geometry::point_type<Geometry>::type> type;
 461 };
 462
 463 template <typename AnyTag, typename Point, typename Geometry>
 464 struct default_strategy<point_tag, AnyTag, point_tag, areal_tag, spherical_tag, spherical_tag, Point, Geometry>
 465 {
 466     typedef winding<Point, typename geometry::point_type<Geometry>::type> type;
 467 };
 468
 469 // TODO: use linear_tag and pointlike_tag the same way how areal_tag is used
 470
 471 template <typename Point, typename Geometry, typename AnyTag>
 472 struct default_strategy<point_tag, AnyTag, point_tag, AnyTag, cartesian_tag, cartesian_tag, Point, Geometry>
 473 {
 474     typedef winding<Point, typename geometry::point_type<Geometry>::type> type;
 475 };
 476
 477 template <typename Point, typename Geometry, typename AnyTag>
 478 struct default_strategy<point_tag, AnyTag, point_tag, AnyTag, spherical_tag, spherical_tag, Point, Geometry>
 479 {
 480     typedef winding<Point, typename geometry::point_type<Geometry>::type> type;
 481 };
 482
 483 } // namespace services
 484
 485 #endif
 486
 487
 488 }} // namespace strategy::within
 489
 490
 491
 492 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
 493 namespace strategy { namespace covered_by { namespace services
 494 {
 495
 496 // Register using "areal_tag" for ring, polygon, multi-polygon
 497 template <typename AnyTag, typename Point, typename Geometry>
 498 struct default_strategy<point_tag, AnyTag, point_tag, areal_tag, cartesian_tag, cartesian_tag, Point, Geometry>
 499 {
 500     typedef strategy::within::winding<Point, typename geometry::point_type<Geometry>::type> type;
 501 };
 502
 503 template <typename AnyTag, typename Point, typename Geometry>
 504 struct default_strategy<point_tag, AnyTag, point_tag, areal_tag, spherical_tag, spherical_tag, Point, Geometry>
 505 {
 506     typedef strategy::within::winding<Point, typename geometry::point_type<Geometry>::type> type;
 507 };
 508
 509 // TODO: use linear_tag and pointlike_tag the same way how areal_tag is used
 510
 511 template <typename Point, typename Geometry, typename AnyTag>
 512 struct default_strategy<point_tag, AnyTag, point_tag, AnyTag, cartesian_tag, cartesian_tag, Point, Geometry>
 513 {
 514     typedef strategy::within::winding<Point, typename geometry::point_type<Geometry>::type> type;
 515 };
 516
 517 template <typename Point, typename Geometry, typename AnyTag>
 518 struct default_strategy<point_tag, AnyTag, point_tag, AnyTag, spherical_tag, spherical_tag, Point, Geometry>
 519 {
 520     typedef strategy::within::winding<Point, typename geometry::point_type<Geometry>::type> type;
 521 };
 522
 523 }}} // namespace strategy::covered_by::services
 524 #endif
 525
 526
 527 }} // namespace boost::geometry
 528
 529
 530 #endif // BOOST_GEOMETRY_STRATEGY_AGNOSTIC_POINT_IN_POLY_WINDING_HPP