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 | 549 |
1 files changed, 54 insertions, 495 deletions
diff --git a/boost/geometry/strategies/agnostic/point_in_poly_winding.hpp b/boost/geometry/strategies/agnostic/point_in_poly_winding.hpp index 0a797ac0f0..774294b57e 100644 --- a/boost/geometry/strategies/agnostic/point_in_poly_winding.hpp +++ b/boost/geometry/strategies/agnostic/point_in_poly_winding.hpp @@ -18,14 +18,15 @@ #define BOOST_GEOMETRY_STRATEGY_AGNOSTIC_POINT_IN_POLY_WINDING_HPP -#include <boost/core/ignore_unused.hpp> +#include <boost/mpl/assert.hpp> -#include <boost/geometry/util/math.hpp> -#include <boost/geometry/util/select_calculation_type.hpp> +#include <boost/geometry/core/cs.hpp> +#include <boost/geometry/core/tag_cast.hpp> +#include <boost/geometry/core/tags.hpp> +#include <boost/geometry/strategies/cartesian/point_in_poly_winding.hpp> #include <boost/geometry/strategies/side.hpp> -#include <boost/geometry/strategies/covered_by.hpp> -#include <boost/geometry/strategies/within.hpp> +#include <boost/geometry/strategies/spherical/point_in_poly_winding.hpp> namespace boost { namespace geometry @@ -34,292 +35,63 @@ 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 <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 -// ___-----> -// ^ O Pt __ __ -// EPS __ __ -// v__ __ BUT DETECTED AS LEFT OF THIS LINE -// _____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 +#ifndef DOXYGEN_NO_DETAIL +namespace detail { - typedef typename strategy::side::services::default_strategy - < - CSTag - >::type strategy_side_type; - - template <typename Point, typename PointOfSegment> - static inline int apply(Point const& point, - PointOfSegment const& se, - int count) - { - 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>(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) - { - ss20 += winding_small_angle<PointOfSegment>::apply(); - } - else - { - 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 +< + typename Point, + typename PointOfSegment, + typename CalculationType, + typename CSTag = typename tag_cast + < + typename cs_tag<Point>::type, + spherical_tag + >::type +> +struct winding_base_type { - template <typename Point, typename PointOfSegment> - static inline int apply(Point const& point, - PointOfSegment const& se, - int count) - { - // NOTE: for D=0 the signs would be reversed - return math::equals(get<1>(point), get<1>(se)) ? - 0 : - get<1>(point) < get<1>(se) ? - // assuming count is equal to 1 or -1 - -count : // ( count > 0 ? -1 : 1) : - count; // ( count > 0 ? 1 : -1) ; - } + BOOST_MPL_ASSERT_MSG(false, + NOT_IMPLEMENTED_FOR_THIS_COORDINATE_SYSTEM, + (CSTag)); }; - -template <typename Point, - typename CalculationType, - typename CSTag = typename cs_tag<Point>::type> -struct winding_check_touch +template <typename Point, typename PointOfSegment, typename CalculationType> +struct winding_base_type<Point, PointOfSegment, CalculationType, cartesian_tag> { - typedef CalculationType calc_t; - typedef typename coordinate_system<Point>::type::units units_t; - typedef math::detail::constants_on_spheroid<CalculationType, units_t> constants; - - template <typename PointOfSegment, typename State> - static inline int apply(Point const& point, - PointOfSegment const& seg1, - PointOfSegment const& seg2, - State& state, - bool& eq1, - bool& eq2) - { - 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) - { - eq1 = eq2 = eq1 || eq2; - - // segment overlapping pole and point is pole - if (math::equals(py, pi2) || math::equals(py, -pi2)) - { - eq1 = eq2 = true; - } - } - - // Both equal p -> segment vertical - // The only thing which has to be done is check if point is ON segment - if (eq1 && eq2) - { - // 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; - - 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); - - 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) - { - 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)) - { - state.m_touches = true; - } - return true; - } - return false; - } + typedef within::cartesian_winding<Point, PointOfSegment, CalculationType> type; }; - -// 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 +template <typename Point, typename PointOfSegment, typename CalculationType> +struct winding_base_type<Point, PointOfSegment, CalculationType, spherical_tag> { - typedef CalculationType calc_t; - typedef typename coordinate_system<Point>::type::units units_t; - - 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 + typedef within::detail::spherical_winding_base + < + Point, + PointOfSegment, + typename strategy::side::services::default_strategy + < + typename cs_tag<Point>::type + >::type, + CalculationType + > type; +}; - // 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 - 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; - } -}; +} // namespace detail +#endif // DOXYGEN_NO_DETAIL /*! -\brief Within detection using winding rule +\brief Within detection using winding rule. Side strategy used internally is + choosen based on Point's coordinate system. \ingroup strategies \tparam Point \tparam_point \tparam PointOfSegment \tparam_segment_point -\tparam SideStrategy Side strategy \tparam CalculationType \tparam_calculation -\author Barend Gehrels -\note The implementation is inspired by terralib http://www.terralib.org (LGPL) -\note but totally revised afterwards, especially for cases on segments -\note Only dependant on "side", -> agnostic, suitable for spherical/latlong \qbk{ [heading See also] @@ -330,245 +102,32 @@ template < typename Point, typename PointOfSegment = Point, - typename SideStrategy = typename strategy::side::services::default_strategy - < - typename cs_tag<Point>::type - >::type, typename CalculationType = void > class winding + : public within::detail::winding_base_type + < + Point, PointOfSegment, CalculationType + >::type { - typedef typename select_calculation_type + typedef typename within::detail::winding_base_type < - Point, - PointOfSegment, - CalculationType - >::type calculation_type; - - /*! subclass to keep state */ - class counter - { - int m_count; - bool m_touches; - - inline int code() const - { - return m_touches ? 0 : m_count == 0 ? -1 : 1; - } - - public : - friend class winding; - - template <typename P, typename CT, typename CST> - friend struct winding_check_touch; - - inline counter() - : m_count(0) - , m_touches(false) - {} - - }; - - static inline int check_segment(Point const& point, - PointOfSegment const& seg1, PointOfSegment const& seg2, - counter& state, bool& eq1, bool& eq2) - { - if (winding_check_touch<Point, calculation_type> - ::apply(point, seg1, seg2, state, eq1, eq2)) - { - return 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); - } - + Point, PointOfSegment, CalculationType + >::type base_t; public: - typedef typename SideStrategy::envelope_strategy_type envelope_strategy_type; - - inline envelope_strategy_type get_envelope_strategy() const - { - return m_side_strategy.get_envelope_strategy(); - } + winding() {} - typedef typename SideStrategy::disjoint_strategy_type disjoint_strategy_type; - - inline disjoint_strategy_type get_disjoint_strategy() const - { - return m_side_strategy.get_disjoint_strategy(); - } - - winding() - {} - - explicit winding(SideStrategy const& side_strategy) - : m_side_strategy(side_strategy) + template <typename Model> + explicit winding(Model const& model) + : base_t(model) {} - - // Typedefs and static methods to fulfill the concept - typedef Point point_type; - typedef PointOfSegment segment_point_type; - typedef counter state_type; - - inline bool apply(Point const& point, - PointOfSegment const& s1, PointOfSegment const& s2, - counter& state) const - { - typedef typename cs_tag<Point>::type cs_t; - - bool eq1 = false; - bool eq2 = false; - boost::ignore_unused(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>::apply(point, eq1 ? s1 : s2, count); - } - else // count == 2 || count == -2 - { - // 1 left, -1 right - side = m_side_strategy.apply(s1, s2, point); - } - - if (side == 0) - { - // Point is lying on segment - state.m_touches = true; - state.m_count = 0; - return false; - } - - // Side is NEG for right, POS for left. - // The count is -2 for down, 2 for up (or -1/1) - // Side positive thus means UP and LEFTSIDE or DOWN and RIGHTSIDE - // See accompagnying figure (TODO) - if (side * count > 0) - { - state.m_count += count; - } - } - return ! state.m_touches; - } - - static inline int result(counter const& state) - { - return state.code(); - } - -private: - SideStrategy m_side_strategy; -}; - - -#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS - -namespace services -{ - -template <typename PointLike, typename Geometry, typename AnyTag1, typename AnyTag2> -struct default_strategy<PointLike, Geometry, AnyTag1, AnyTag2, pointlike_tag, polygonal_tag, cartesian_tag, cartesian_tag> -{ - typedef winding - < - typename geometry::point_type<PointLike>::type, - typename geometry::point_type<Geometry>::type - > type; }; -template <typename PointLike, typename Geometry, typename AnyTag1, typename AnyTag2> -struct default_strategy<PointLike, Geometry, AnyTag1, AnyTag2, pointlike_tag, polygonal_tag, spherical_tag, spherical_tag> -{ - typedef winding - < - typename geometry::point_type<PointLike>::type, - typename geometry::point_type<Geometry>::type - > type; -}; - -template <typename PointLike, typename Geometry, typename AnyTag1, typename AnyTag2> -struct default_strategy<PointLike, Geometry, AnyTag1, AnyTag2, pointlike_tag, linear_tag, cartesian_tag, cartesian_tag> -{ - typedef winding - < - typename geometry::point_type<PointLike>::type, - typename geometry::point_type<Geometry>::type - > type; -}; - -template <typename PointLike, typename Geometry, typename AnyTag1, typename AnyTag2> -struct default_strategy<PointLike, Geometry, AnyTag1, AnyTag2, pointlike_tag, linear_tag, spherical_tag, spherical_tag> -{ - typedef winding - < - typename geometry::point_type<PointLike>::type, - typename geometry::point_type<Geometry>::type - > type; -}; - -} // namespace services - -#endif - }} // namespace strategy::within -#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS -namespace strategy { namespace covered_by { namespace services -{ - -template <typename PointLike, typename Geometry, typename AnyTag1, typename AnyTag2> -struct default_strategy<PointLike, Geometry, AnyTag1, AnyTag2, pointlike_tag, polygonal_tag, cartesian_tag, cartesian_tag> -{ - typedef within::winding - < - typename geometry::point_type<PointLike>::type, - typename geometry::point_type<Geometry>::type - > type; -}; - -template <typename PointLike, typename Geometry, typename AnyTag1, typename AnyTag2> -struct default_strategy<PointLike, Geometry, AnyTag1, AnyTag2, pointlike_tag, polygonal_tag, spherical_tag, spherical_tag> -{ - typedef within::winding - < - typename geometry::point_type<PointLike>::type, - typename geometry::point_type<Geometry>::type - > type; -}; - -template <typename PointLike, typename Geometry, typename AnyTag1, typename AnyTag2> -struct default_strategy<PointLike, Geometry, AnyTag1, AnyTag2, pointlike_tag, linear_tag, cartesian_tag, cartesian_tag> -{ - typedef within::winding - < - typename geometry::point_type<PointLike>::type, - typename geometry::point_type<Geometry>::type - > type; -}; - -template <typename PointLike, typename Geometry, typename AnyTag1, typename AnyTag2> -struct default_strategy<PointLike, Geometry, AnyTag1, AnyTag2, pointlike_tag, linear_tag, spherical_tag, spherical_tag> -{ - typedef within::winding - < - typename geometry::point_type<PointLike>::type, - typename geometry::point_type<Geometry>::type - > type; -}; - -}}} // namespace strategy::covered_by::services -#endif - - }} // namespace boost::geometry |