diff options
Diffstat (limited to 'boost/geometry/strategies/agnostic/point_in_poly_winding.hpp')
| -rw-r--r-- | boost/geometry/strategies/agnostic/point_in_poly_winding.hpp | 330 |
1 files changed, 216 insertions, 114 deletions
diff --git a/boost/geometry/strategies/agnostic/point_in_poly_winding.hpp b/boost/geometry/strategies/agnostic/point_in_poly_winding.hpp index 641533fc6a..9e2ec2c4ff 100644 --- a/boost/geometry/strategies/agnostic/point_in_poly_winding.hpp +++ b/boost/geometry/strategies/agnostic/point_in_poly_winding.hpp @@ -3,8 +3,9 @@ // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. // Copyright (c) 2013 Adam Wulkiewicz, Lodz, Poland. -// This file was modified by Oracle on 2013, 2014. -// Modifications copyright (c) 2013, 2014 Oracle and/or its affiliates. +// This file was modified by Oracle on 2013, 2014, 2016. +// Modifications copyright (c) 2013-2016 Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands. @@ -13,8 +14,6 @@ // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at // http://www.boost.org/LICENSE_1_0.txt) -// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle - #ifndef BOOST_GEOMETRY_STRATEGY_AGNOSTIC_POINT_IN_POLY_WINDING_HPP #define BOOST_GEOMETRY_STRATEGY_AGNOSTIC_POINT_IN_POLY_WINDING_HPP @@ -35,15 +34,34 @@ namespace boost { namespace geometry namespace strategy { namespace within { +// 1 deg or pi/180 rad +template <typename Point, + typename CalculationType = typename coordinate_type<Point>::type> +struct winding_small_angle +{ + typedef typename coordinate_system<Point>::type cs_t; + typedef math::detail::constants_on_spheroid + < + CalculationType, + typename cs_t::units + > constants; + + static inline CalculationType apply() + { + return constants::half_period() / CalculationType(180); + } +}; + // Fix for https://svn.boost.org/trac/boost/ticket/9628 -// For floating point coordinates, the <1> coordinate of a point is compared +// For floating point coordinates, the <D> coordinate of a point is compared // with the segment's points using some EPS. If the coordinates are "equal" // the sides are calculated. Therefore we can treat a segment as a long areal // geometry having some width. There is a small ~triangular area somewhere // between the segment's effective area and a segment's line used in sides // calculation where the segment is on the one side of the line but on the // other side of a segment (due to the width). +// Below picture assuming D = 1, if D = 0 horiz<->vert, E<->N, RIGHT<->UP. // For the s1 of a segment going NE the real side is RIGHT but the point may // be detected as LEFT, like this: // RIGHT @@ -54,6 +72,9 @@ namespace strategy { namespace within // _____7 // _____/ // _____/ +// In the code below actually D = 0, so segments are nearly-vertical +// Called when the point is on the same level as one of the segment's points +// but the point is not aligned with a vertical segment template <typename CSTag> struct winding_side_equal { @@ -62,122 +83,228 @@ struct winding_side_equal CSTag >::type strategy_side_type; - template <size_t D, typename Point, typename PointOfSegment> + template <typename Point, typename PointOfSegment> static inline int apply(Point const& point, PointOfSegment const& se, int count) { - // Create a vertical segment intersecting the original segment's endpoint - // equal to the point, with the derived direction (UP/DOWN). - // Set only the 2 first coordinates, the other ones are ignored + typedef typename coordinate_type<PointOfSegment>::type scoord_t; + typedef typename coordinate_system<PointOfSegment>::type::units units_t; + + if (math::equals(get<1>(point), get<1>(se))) + return 0; + + // Create a horizontal segment intersecting the original segment's endpoint + // equal to the point, with the derived direction (E/W). PointOfSegment ss1, ss2; - set<1-D>(ss1, get<1-D>(se)); - set<1-D>(ss2, get<1-D>(se)); - if (count > 0) // UP + set<1>(ss1, get<1>(se)); + set<0>(ss1, get<0>(se)); + set<1>(ss2, get<1>(se)); + scoord_t ss20 = get<0>(se); + if (count > 0) { - set<D>(ss1, 0); - set<D>(ss2, 1); + ss20 += winding_small_angle<PointOfSegment>::apply(); } - else // DOWN + else { - set<D>(ss1, 1); - set<D>(ss2, 0); + ss20 -= winding_small_angle<PointOfSegment>::apply(); } + math::normalize_longitude<units_t>(ss20); + set<0>(ss2, ss20); + // Check the side using this vertical segment return strategy_side_type::apply(ss1, ss2, point); } }; - // The optimization for cartesian template <> struct winding_side_equal<cartesian_tag> { - template <size_t D, typename Point, typename PointOfSegment> + template <typename Point, typename PointOfSegment> static inline int apply(Point const& point, PointOfSegment const& se, int count) { - return math::equals(get<1-D>(point), get<1-D>(se)) ? + // NOTE: for D=0 the signs would be reversed + return math::equals(get<1>(point), get<1>(se)) ? 0 : - get<1-D>(point) < get<1-D>(se) ? + get<1>(point) < get<1>(se) ? // assuming count is equal to 1 or -1 - count : // ( count > 0 ? 1 : -1) : - -count; // ( count > 0 ? -1 : 1) ; + -count : // ( count > 0 ? -1 : 1) : + count; // ( count > 0 ? 1 : -1) ; } }; -template <typename CSTag> -struct winding_side_between +template <typename Point, + typename CalculationType, + typename CSTag = typename cs_tag<Point>::type> +struct winding_check_touch { - typedef typename strategy::side::services::default_strategy - < - CSTag - >::type strategy_side_type; + typedef CalculationType calc_t; + typedef typename coordinate_system<Point>::type::units units_t; + typedef math::detail::constants_on_spheroid<CalculationType, units_t> constants; - template <size_t D, typename Point, typename PointOfSegment> + template <typename PointOfSegment, typename State> static inline int apply(Point const& point, - PointOfSegment const& s1, PointOfSegment const& s2, - int count) + PointOfSegment const& seg1, + PointOfSegment const& seg2, + State& state, + bool& eq1, + bool& eq2) { - // Create a vertical segment intersecting the original segment's endpoint - // equal to the point, with the derived direction (UP/DOWN). - // Set only the 2 first coordinates, the other ones are ignored - PointOfSegment ss1, ss2; - set<1-D>(ss1, get<1-D>(s1)); - set<1-D>(ss2, get<1-D>(s1)); - - if (count > 0) // UP + calc_t const pi = constants::half_period(); + calc_t const pi2 = pi / calc_t(2); + + calc_t const px = get<0>(point); + calc_t const s1x = get<0>(seg1); + calc_t const s2x = get<0>(seg2); + calc_t const py = get<1>(point); + calc_t const s1y = get<1>(seg1); + calc_t const s2y = get<1>(seg2); + + // NOTE: lat in {-90, 90} and arbitrary lon + // it doesn't matter what lon it is if it's a pole + // so e.g. if one of the segment endpoints is a pole + // then only the other lon matters + + bool eq1_strict = math::equals(s1x, px); + bool eq2_strict = math::equals(s2x, px); + + eq1 = eq1_strict // lon strictly equal to s1 + || math::equals(s1y, pi2) || math::equals(s1y, -pi2); // s1 is pole + eq2 = eq2_strict // lon strictly equal to s2 + || math::equals(s2y, pi2) || math::equals(s2y, -pi2); // s2 is pole + + // segment overlapping pole + calc_t s1x_anti = s1x + constants::half_period(); + math::normalize_longitude<units_t, calc_t>(s1x_anti); + bool antipodal = math::equals(s2x, s1x_anti); + if (antipodal) { - set<D>(ss1, 0); - set<D>(ss2, 1); + eq1 = eq2 = eq1 || eq2; + + // segment overlapping pole and point is pole + if (math::equals(py, pi2) || math::equals(py, -pi2)) + { + eq1 = eq2 = true; + } } - else // DOWN + + // Both equal p -> segment vertical + // The only thing which has to be done is check if point is ON segment + if (eq1 && eq2) { - set<D>(ss1, 1); - set<D>(ss2, 0); + // segment endpoints on the same sides of the globe + if (! antipodal + // p's lat between segment endpoints' lats + ? (s1y <= py && s2y >= py) || (s2y <= py && s1y >= py) + // going through north or south pole? + : (pi - s1y - s2y <= pi + ? (eq1_strict && s1y <= py) || (eq2_strict && s2y <= py) // north + || math::equals(py, pi2) // point on north pole + : (eq1_strict && s1y >= py) || (eq2_strict && s2y >= py)) // south + || math::equals(py, -pi2) // point on south pole + ) + { + state.m_touches = true; + } + return true; } + return false; + } +}; +// The optimization for cartesian +template <typename Point, typename CalculationType> +struct winding_check_touch<Point, CalculationType, cartesian_tag> +{ + typedef CalculationType calc_t; - int const seg_side = strategy_side_type::apply(ss1, ss2, s2); + template <typename PointOfSegment, typename State> + static inline bool apply(Point const& point, + PointOfSegment const& seg1, + PointOfSegment const& seg2, + State& state, + bool& eq1, + bool& eq2) + { + calc_t const px = get<0>(point); + calc_t const s1x = get<0>(seg1); + calc_t const s2x = get<0>(seg2); - if (seg_side != 0) // segment not vertical + eq1 = math::equals(s1x, px); + eq2 = math::equals(s2x, px); + + // Both equal p -> segment vertical + // The only thing which has to be done is check if point is ON segment + if (eq1 && eq2) { - if (strategy_side_type::apply(ss1, ss2, point) == -seg_side) // point on the opposite side than s2 - { - return -seg_side; - } - else + calc_t const py = get<1>(point); + calc_t const s1y = get<1>(seg1); + calc_t const s2y = get<1>(seg2); + if ((s1y <= py && s2y >= py) || (s2y <= py && s1y >= py)) { - set<1-D>(ss1, get<1-D>(s2)); - set<1-D>(ss2, get<1-D>(s2)); - - if (strategy_side_type::apply(ss1, ss2, point) == seg_side) // point behind s2 - { - return seg_side; - } + state.m_touches = true; } + return true; } - - // segment is vertical or point is between p1 and p2 - return strategy_side_type::apply(s1, s2, point); + return false; } }; -// The specialization for cartesian -template <> -struct winding_side_between<cartesian_tag> + +// Called if point is not aligned with a vertical segment +template <typename Point, + typename CalculationType, + typename CSTag = typename cs_tag<Point>::type> +struct winding_calculate_count { - typedef strategy::side::services::default_strategy - < - cartesian_tag - >::type strategy_side_type; + typedef CalculationType calc_t; + typedef typename coordinate_system<Point>::type::units units_t; - template <size_t D, typename Point, typename PointOfSegment> - static inline int apply(Point const& point, - PointOfSegment const& s1, PointOfSegment const& s2, - int /*count*/) + static inline bool greater(calc_t const& l, calc_t const& r) + { + calc_t diff = l - r; + math::normalize_longitude<units_t, calc_t>(diff); + return diff > calc_t(0); + } + + static inline int apply(calc_t const& p, + calc_t const& s1, calc_t const& s2, + bool eq1, bool eq2) + { + // Probably could be optimized by avoiding normalization for some comparisons + // e.g. s1 > p could be calculated from p > s1 + + // If both segment endpoints were poles below checks wouldn't be enough + // but this means that either both are the same or that they are N/S poles + // and therefore the segment is not valid. + // If needed (eq1 && eq2 ? 0) could be returned + + return + eq1 ? (greater(s2, p) ? 1 : -1) // Point on level s1, E/W depending on s2 + : eq2 ? (greater(s1, p) ? -1 : 1) // idem + : greater(p, s1) && greater(s2, p) ? 2 // Point between s1 -> s2 --> E + : greater(p, s2) && greater(s1, p) ? -2 // Point between s2 -> s1 --> W + : 0; + } +}; +// The optimization for cartesian +template <typename Point, typename CalculationType> +struct winding_calculate_count<Point, CalculationType, cartesian_tag> +{ + typedef CalculationType calc_t; + + static inline int apply(calc_t const& p, + calc_t const& s1, calc_t const& s2, + bool eq1, bool eq2) { - return strategy_side_type::apply(s1, s2, point); + return + eq1 ? (s2 > p ? 1 : -1) // Point on level s1, E/W depending on s2 + : eq2 ? (s1 > p ? -1 : 1) // idem + : s1 < p && s2 > p ? 2 // Point between s1 -> s2 --> E + : s2 < p && s1 > p ? -2 // Point between s2 -> s1 --> W + : 0; } }; @@ -234,6 +361,9 @@ class winding public : friend class winding; + template <typename P, typename CT, typename CST> + friend struct winding_check_touch; + inline counter() : m_count(0) , m_touches(false) @@ -242,48 +372,21 @@ class winding }; - template <size_t D> - static inline int check_touch(Point const& point, - PointOfSegment const& seg1, PointOfSegment const& seg2, - counter& state) - { - calculation_type const p = get<D>(point); - calculation_type const s1 = get<D>(seg1); - calculation_type const s2 = get<D>(seg2); - if ((s1 <= p && s2 >= p) || (s2 <= p && s1 >= p)) - { - state.m_touches = true; - } - return 0; - } - - - template <size_t D> static inline int check_segment(Point const& point, PointOfSegment const& seg1, PointOfSegment const& seg2, counter& state, bool& eq1, bool& eq2) { - calculation_type const p = get<D>(point); - calculation_type const s1 = get<D>(seg1); - calculation_type const s2 = get<D>(seg2); - - // Check if one of segment endpoints is at same level of point - eq1 = math::equals(s1, p); - eq2 = math::equals(s2, p); - - if (eq1 && eq2) + if (winding_check_touch<Point, calculation_type> + ::apply(point, seg1, seg2, state, eq1, eq2)) { - // Both equal p -> segment is horizontal (or vertical for D=0) - // The only thing which has to be done is check if point is ON segment - return check_touch<1 - D>(point, seg1, seg2, state); + return 0; } - return - eq1 ? (s2 > p ? 1 : -1) // Point on level s1, UP/DOWN depending on s2 - : eq2 ? (s1 > p ? -1 : 1) // idem - : s1 < p && s2 > p ? 2 // Point between s1 -> s2 --> UP - : s2 < p && s1 > p ? -2 // Point between s2 -> s1 --> DOWN - : 0; + calculation_type const p = get<0>(point); + calculation_type const s1 = get<0>(seg1); + calculation_type const s2 = get<0>(seg2); + return winding_calculate_count<Point, calculation_type> + ::apply(p, s1, s2, eq1, eq2); } @@ -304,19 +407,18 @@ public : bool eq2 = false; boost::ignore_unused(eq2); - int count = check_segment<1>(point, s1, s2, state, eq1, eq2); + int count = check_segment(point, s1, s2, state, eq1, eq2); if (count != 0) { int side = 0; if (count == 1 || count == -1) { - side = winding_side_equal<cs_t> - ::template apply<1>(point, eq1 ? s1 : s2, count); + side = winding_side_equal<cs_t>::apply(point, eq1 ? s1 : s2, count); } else // count == 2 || count == -2 { - side = winding_side_between<cs_t> - ::template apply<1>(point, s1, s2, count); + // 1 left, -1 right + side = strategy_side_type::apply(s1, s2, point); } if (side == 0) |