From a48401d94a8b4dd10590631c1948a392ab834dda Mon Sep 17 00:00:00 2001 From: Adam Wulkiewicz Date: Tue, 9 Nov 2021 21:58:10 +0100 Subject: [PATCH] [envelope] Simplify the pole in ring check with side_of_pole strategy. --- .../strategy/geographic/envelope_range.hpp | 12 +- .../strategy/spherical/envelope_range.hpp | 114 +++++++++++++++--- .../envelope_expand/envelope_on_spheroid.cpp | 24 ++-- 3 files changed, 119 insertions(+), 31 deletions(-) diff --git a/include/boost/geometry/strategy/geographic/envelope_range.hpp b/include/boost/geometry/strategy/geographic/envelope_range.hpp index a3cbbe1fa..39739107f 100644 --- a/include/boost/geometry/strategy/geographic/envelope_range.hpp +++ b/include/boost/geometry/strategy/geographic/envelope_range.hpp @@ -14,8 +14,8 @@ #include #include -// TEMP - get rid of this dependency -#include +// Get rid of this dependency? +#include namespace boost { namespace geometry { @@ -96,11 +96,11 @@ public: < FormulaPolicy, Spheroid, CalculationType >(m_spheroid), - within::geographic_winding + within::detail::spherical_winding_base < - void, void, - FormulaPolicy, Spheroid, CalculationType - >(m_spheroid)); + envelope::detail::side_of_pole, + CalculationType + >()); } Spheroid model() const diff --git a/include/boost/geometry/strategy/spherical/envelope_range.hpp b/include/boost/geometry/strategy/spherical/envelope_range.hpp index d0c9d3b89..48b5d79c8 100644 --- a/include/boost/geometry/strategy/spherical/envelope_range.hpp +++ b/include/boost/geometry/strategy/spherical/envelope_range.hpp @@ -18,8 +18,7 @@ #include #include -// TEMP - get rid of these dependencies -#include +// Get rid of this dependency? #include @@ -73,8 +72,87 @@ inline void spheroidal_linestring(Range const& range, Box& mbr, } -template -inline bool pole_within(bool north_pole, Ring const& ring, WithinStrategy const& within_strategy) +// Side of pole WRT segment which doesn't contain it. +template +struct side_of_pole +{ + typedef spherical_tag cs_tag; + + template + static inline int apply(P const& p1, P const& p2, P const& p) + { + using calc_t = typename promote_floating_point + < + typename select_calculation_type_alt + < + CalculationType, P + >::type + >::type; + + using units_t = typename geometry::detail::cs_angular_units

::type; + using constants = math::detail::constants_on_spheroid; + + calc_t const c0 = 0; + calc_t const pi = constants::half_period(); + + calc_t const lon1 = get<0>(p1); + calc_t const lat1 = get<1>(p1); + calc_t const lon2 = get<0>(p2); + calc_t const lat2 = get<1>(p2); + calc_t const lat = get<1>(p); + + calc_t const s_lon_diff = math::longitude_distance_signed(lon1, lon2); + bool const s_vertical = math::equals(s_lon_diff, c0) + || math::equals(s_lon_diff, pi); + // Side of vertical segment is 0 for both poles. + if (s_vertical) + { + return 0; + } + + // This strategy shouldn't be called in this case but just in case + // check if segment starts at a pole + if (math::equals(lat, lat1) || math::equals(lat, lat2)) + { + return 0; + } + + // -1 is rhs + // 1 is lhs + if (lat >= c0) // north pole + { + return s_lon_diff < c0 ? -1 : 1; + } + else // south pole + { + return s_lon_diff > c0 ? -1 : 1; + } + } +}; + + +template +inline int point_in_range(Point const& point, Range const& range, Strategy const& strategy) +{ + typename Strategy::state_type state; + + auto it = boost::begin(range); + auto const end = boost::end(range); + for (auto previous = it++ ; it != end ; ++previous, ++it ) + { + if ( ! strategy.apply(point, *previous, *it, state) ) + { + break; + } + } + + return strategy.result(state); +} + + +template +inline bool pole_within(bool north_pole, Ring const& ring, + PoleWithinStrategy const& pole_within_strategy) { if (boost::size(ring) < core_detail::closure::minimum_ring_size < @@ -99,14 +177,18 @@ inline bool pole_within(bool north_pole, Ring const& ring, WithinStrategy const& geometry::set<1>(point, constants_t::min_latitude()); } geometry::detail::closed_clockwise_view view(ring); - return geometry::detail::within::point_in_range(point, view, within_strategy) > 0; + return point_in_range(point, view, pole_within_strategy) > 0; } -template +template +< + typename Range, typename Box, + typename EnvelopeStrategy, typename ExpandStrategy, typename PoleWithinStrategy +> inline void spheroidal_ring(Range const& range, Box& mbr, EnvelopeStrategy const& envelope_strategy, ExpandStrategy const& expand_strategy, - WithinStrategy const& within_strategy) + PoleWithinStrategy const& pole_within_strategy) { geometry::detail::closed_view closed_range(range); @@ -121,9 +203,9 @@ inline void spheroidal_ring(Range const& range, Box& mbr, coord_t const lon_max = geometry::get<1, 0>(mbr); // If box covers the whole longitude range it is possible that the ring contains // one of the poles. - // TODO: Technically it is possible that a reversed ring may cover more than + // Technically it is possible that a reversed ring may cover more than // half of the globe and mbr of it's linear ring may be small and not cover the - // longitude range. + // longitude range. We currently don't support such rings. if (lon_max - lon_min >= two_pi) { coord_t const lat_n_pole = constants_t::max_latitude(); @@ -140,20 +222,16 @@ inline void spheroidal_ring(Range const& range, Box& mbr, lat_max = lat_n_pole; } - // TODO - implement something simpler than within strategy because here - // we know that neither min nor max is a pole so there is no segment which - // contains a pole, no endpoint, no vertex at pole, there are no antipodal - // points. So many special cases can be ignored. if (lat_max < lat_n_pole) { - if (pole_within(true, range, within_strategy)) + if (pole_within(true, range, pole_within_strategy)) { lat_max = lat_n_pole; } } if (lat_min > lat_s_pole) { - if (pole_within(false, range, within_strategy)) + if (pole_within(false, range, pole_within_strategy)) { lat_min = lat_s_pole; } @@ -190,7 +268,11 @@ public: detail::spheroidal_ring(range, mbr, envelope::spherical_segment(), expand::spherical_segment(), - within::spherical_winding()); + within::detail::spherical_winding_base + < + envelope::detail::side_of_pole, + CalculationType + >()); } }; diff --git a/test/algorithms/envelope_expand/envelope_on_spheroid.cpp b/test/algorithms/envelope_expand/envelope_on_spheroid.cpp index e270fdba4..92d73e361 100644 --- a/test/algorithms/envelope_expand/envelope_on_spheroid.cpp +++ b/test/algorithms/envelope_expand/envelope_on_spheroid.cpp @@ -2736,14 +2736,20 @@ BOOST_AUTO_TEST_CASE( envelope_cw_ring ) tester::apply("r28cw", from_wkt("POLYGON((),(0 -80,-90 -80,-180 -80,90 -80,0 -80))"), -180, -90, 180, -80); - - typedef bg::model::polygon G2; - typedef test_envelope_on_sphere_or_spheroid tester2; - tester2::apply("r27ccw", - from_wkt("POLYGON((),(0 80,-90 80,-180 80,90 80,0 80))"), - -180, 80, 180, 90); - tester2::apply("r28ccw", - from_wkt("POLYGON((),(0 -80,90 -80,-180 -80,-90 -80,0 -80))"), - -180, -90, 180, -80); } +BOOST_AUTO_TEST_CASE(envelope_ccw_ring) +{ + typedef bg::cs::spherical_equatorial coordinate_system_type; + typedef bg::model::point point_type; + typedef bg::model::polygon G; + typedef bg::model::box B; + typedef test_envelope_on_sphere_or_spheroid tester; + + tester::apply("r27ccw", + from_wkt("POLYGON((),(0 80,-90 80,-180 80,90 80,0 80))"), + -180, 80, 180, 90); + tester::apply("r28ccw", + from_wkt("POLYGON((),(0 -80,90 -80,-180 -80,-90 -80,0 -80))"), + -180, -90, 180, -80); +} \ No newline at end of file