diff options
author | DongHun Kwak <dh0128.kwak@samsung.com> | 2019-12-05 15:11:01 +0900 |
---|---|---|
committer | DongHun Kwak <dh0128.kwak@samsung.com> | 2019-12-05 15:11:01 +0900 |
commit | 3fdc3e5ee96dca5b11d1694975a65200787eab86 (patch) | |
tree | 5c1733853892b8397d67706fa453a9bd978d2102 /boost/geometry/strategies | |
parent | 88e602c57797660ebe0f9e15dbd64c1ff16dead3 (diff) | |
download | boost-3fdc3e5ee96dca5b11d1694975a65200787eab86.tar.gz boost-3fdc3e5ee96dca5b11d1694975a65200787eab86.tar.bz2 boost-3fdc3e5ee96dca5b11d1694975a65200787eab86.zip |
Imported Upstream version 1.66.0upstream/1.66.0
Diffstat (limited to 'boost/geometry/strategies')
21 files changed, 2313 insertions, 835 deletions
diff --git a/boost/geometry/strategies/agnostic/buffer_distance_symmetric.hpp b/boost/geometry/strategies/agnostic/buffer_distance_symmetric.hpp index 73bd21ac73..fe973e3f38 100644 --- a/boost/geometry/strategies/agnostic/buffer_distance_symmetric.hpp +++ b/boost/geometry/strategies/agnostic/buffer_distance_symmetric.hpp @@ -68,7 +68,7 @@ public : return negative() ? -1 : 1; } - //! Returns true if distance is negative + //! Returns true if distance is negative (aka deflate) inline bool negative() const { return m_distance < 0; diff --git a/boost/geometry/strategies/agnostic/point_in_point.hpp b/boost/geometry/strategies/agnostic/point_in_point.hpp index 1a9274149a..d4692766c5 100644 --- a/boost/geometry/strategies/agnostic/point_in_point.hpp +++ b/boost/geometry/strategies/agnostic/point_in_point.hpp @@ -32,7 +32,7 @@ struct point_in_point { static inline bool apply(Point1 const& point1, Point2 const& point2) { - return detail::equals::equals_point_point(point1, point2); + return geometry::detail::equals::equals_point_point(point1, point2); } }; 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 diff --git a/boost/geometry/strategies/cartesian/intersection.hpp b/boost/geometry/strategies/cartesian/intersection.hpp index 50e903885b..233bb50b64 100644 --- a/boost/geometry/strategies/cartesian/intersection.hpp +++ b/boost/geometry/strategies/cartesian/intersection.hpp @@ -33,10 +33,10 @@ #include <boost/geometry/util/promote_integral.hpp> #include <boost/geometry/util/select_calculation_type.hpp> -#include <boost/geometry/strategies/agnostic/point_in_poly_winding.hpp> #include <boost/geometry/strategies/cartesian/area_surveyor.hpp> #include <boost/geometry/strategies/cartesian/distance_pythagoras.hpp> #include <boost/geometry/strategies/cartesian/envelope_segment.hpp> +#include <boost/geometry/strategies/cartesian/point_in_poly_winding.hpp> #include <boost/geometry/strategies/cartesian/side_by_triangle.hpp> #include <boost/geometry/strategies/covered_by.hpp> #include <boost/geometry/strategies/intersection.hpp> @@ -81,11 +81,10 @@ struct cartesian_segments template <typename Geometry1, typename Geometry2> struct point_in_geometry_strategy { - typedef strategy::within::winding + typedef strategy::within::cartesian_winding < typename point_type<Geometry1>::type, typename point_type<Geometry2>::type, - side_strategy_type, CalculationType > type; }; diff --git a/boost/geometry/strategies/cartesian/point_in_poly_winding.hpp b/boost/geometry/strategies/cartesian/point_in_poly_winding.hpp new file mode 100644 index 0000000000..c41bc9b83d --- /dev/null +++ b/boost/geometry/strategies/cartesian/point_in_poly_winding.hpp @@ -0,0 +1,296 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// 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, 2016, 2017. +// Modifications copyright (c) 2013-2017 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. + +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_GEOMETRY_STRATEGY_CARTESIAN_POINT_IN_POLY_WINDING_HPP +#define BOOST_GEOMETRY_STRATEGY_CARTESIAN_POINT_IN_POLY_WINDING_HPP + + +#include <boost/geometry/core/tags.hpp> + +#include <boost/geometry/util/math.hpp> +#include <boost/geometry/util/select_calculation_type.hpp> + +#include <boost/geometry/strategies/cartesian/side_by_triangle.hpp> +#include <boost/geometry/strategies/covered_by.hpp> +#include <boost/geometry/strategies/within.hpp> + + +namespace boost { namespace geometry +{ + +namespace strategy { namespace within +{ + + +/*! +\brief Within detection using winding rule in cartesian coordinate system. +\ingroup strategies +\tparam Point \tparam_point +\tparam PointOfSegment \tparam_segment_point +\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 + +\qbk{ +[heading See also] +[link geometry.reference.algorithms.within.within_3_with_strategy within (with strategy)] +} + */ +template +< + typename Point, + typename PointOfSegment = Point, + typename CalculationType = void +> +class cartesian_winding +{ + typedef side::side_by_triangle<CalculationType> side_strategy_type; + + typedef typename select_calculation_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 cartesian_winding; + + inline counter() + : m_count(0) + , m_touches(false) + {} + + }; + +public: + typedef typename side_strategy_type::envelope_strategy_type envelope_strategy_type; + + static inline envelope_strategy_type get_envelope_strategy() + { + return side_strategy_type::get_envelope_strategy(); + } + + typedef typename side_strategy_type::disjoint_strategy_type disjoint_strategy_type; + + static inline disjoint_strategy_type get_disjoint_strategy() + { + return side_strategy_type::get_disjoint_strategy(); + } + + // Typedefs and static methods to fulfill the concept + typedef Point point_type; + typedef PointOfSegment segment_point_type; + typedef counter state_type; + + static inline bool apply(Point const& point, + PointOfSegment const& s1, PointOfSegment const& s2, + counter& state) + { + bool eq1 = false; + bool eq2 = false; + + int count = check_segment(point, s1, s2, state, eq1, eq2); + if (count != 0) + { + int side = 0; + if (count == 1 || count == -1) + { + side = side_equal(point, eq1 ? s1 : s2, count); + } + else // count == 2 || count == -2 + { + // 1 left, -1 right + side = side_strategy_type::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 left, 2 for right (or -1/1) + // Side positive thus means RIGHT and LEFTSIDE or LEFT 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: + static inline int check_segment(Point const& point, + PointOfSegment const& seg1, + PointOfSegment const& seg2, + counter& state, + bool& eq1, bool& eq2) + { + if (check_touch(point, seg1, seg2, state, eq1, eq2)) + { + return 0; + } + + return calculate_count(point, seg1, seg2, eq1, eq2); + } + + static inline bool check_touch(Point const& point, + PointOfSegment const& seg1, + PointOfSegment const& seg2, + counter& state, + bool& eq1, bool& eq2) + { + calculation_type const px = get<0>(point); + calculation_type const s1x = get<0>(seg1); + calculation_type 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) + { + calculation_type const py = get<1>(point); + calculation_type const s1y = get<1>(seg1); + calculation_type const s2y = get<1>(seg2); + if ((s1y <= py && s2y >= py) || (s2y <= py && s1y >= py)) + { + state.m_touches = true; + } + return true; + } + return false; + } + + static inline int calculate_count(Point const& point, + PointOfSegment const& seg1, + PointOfSegment const& seg2, + bool eq1, bool eq2) + { + calculation_type const p = get<0>(point); + calculation_type const s1 = get<0>(seg1); + calculation_type const s2 = get<0>(seg2); + + 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; + } + + static inline int side_equal(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) ; + } +}; + + +#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 cartesian_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 cartesian_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::cartesian_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::cartesian_winding + < + typename geometry::point_type<PointLike>::type, + typename geometry::point_type<Geometry>::type + > type; +}; + +}}} // namespace strategy::covered_by::services +#endif + + +}} // namespace boost::geometry + + +#endif // BOOST_GEOMETRY_STRATEGY_CARTESIAN_POINT_IN_POLY_WINDING_HPP diff --git a/boost/geometry/strategies/cartesian/side_by_triangle.hpp b/boost/geometry/strategies/cartesian/side_by_triangle.hpp index 4d1d97520f..91bbee7bb2 100644 --- a/boost/geometry/strategies/cartesian/side_by_triangle.hpp +++ b/boost/geometry/strategies/cartesian/side_by_triangle.hpp @@ -4,10 +4,11 @@ // Copyright (c) 2008-2015 Bruno Lalande, Paris, France. // Copyright (c) 2009-2015 Mateusz Loskot, London, UK. -// This file was modified by Oracle on 2015. -// Modifications copyright (c) 2015, Oracle and/or its affiliates. +// This file was modified by Oracle on 2015, 2017. +// Modifications copyright (c) 2015-2017, Oracle and/or its affiliates. // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle +// 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. @@ -29,9 +30,9 @@ #include <boost/geometry/strategies/cartesian/disjoint_segment_box.hpp> #include <boost/geometry/strategies/cartesian/envelope_segment.hpp> +#include <boost/geometry/strategies/compare.hpp> #include <boost/geometry/strategies/side.hpp> -#include <boost/geometry/algorithms/detail/relate/less.hpp> #include <boost/geometry/algorithms/detail/equals/point_point.hpp> @@ -182,10 +183,11 @@ public : // arguments, we cyclically permute them so that the first // argument is always the lexicographically smallest point. - geometry::detail::relate::less less; - if (less(p, p1)) + typedef compare::cartesian<compare::less> less; + + if (less::apply(p, p1)) { - if (less(p, p2)) + if (less::apply(p, p2)) { // p is the lexicographically smallest return side_value<CoordinateType, PromotedType>(p, p1, p2, epsp); @@ -197,7 +199,7 @@ public : } } - if (less(p1, p2)) + if (less::apply(p1, p2)) { // p1 is the lexicographically smallest return side_value<CoordinateType, PromotedType>(p1, p2, p, epsp); diff --git a/boost/geometry/strategies/compare.hpp b/boost/geometry/strategies/compare.hpp index 2958319229..b196b75ec2 100644 --- a/boost/geometry/strategies/compare.hpp +++ b/boost/geometry/strategies/compare.hpp @@ -4,6 +4,11 @@ // Copyright (c) 2008-2012 Bruno Lalande, Paris, France. // Copyright (c) 2009-2012 Mateusz Loskot, London, UK. +// This file was modified by Oracle on 2017. +// Modifications copyright (c) 2017, 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. @@ -15,155 +20,204 @@ #ifndef BOOST_GEOMETRY_STRATEGIES_COMPARE_HPP #define BOOST_GEOMETRY_STRATEGIES_COMPARE_HPP + #include <cstddef> #include <functional> -#include <boost/mpl/if.hpp> +#include <boost/mpl/assert.hpp> +#include <boost/mpl/min.hpp> +#include <boost/geometry/core/access.hpp> #include <boost/geometry/core/cs.hpp> #include <boost/geometry/core/coordinate_type.hpp> +#include <boost/geometry/core/coordinate_dimension.hpp> -#include <boost/geometry/strategies/tags.hpp> +#include <boost/geometry/util/math.hpp> namespace boost { namespace geometry { -/*! - \brief Traits class binding a comparing strategy to a coordinate system - \ingroup util - \tparam Tag tag of coordinate system of point-type - \tparam Direction direction to compare on: 1 for less (-> ascending order) - and -1 for greater (-> descending order) - \tparam Point point-type - \tparam CoordinateSystem coordinate sytem of point - \tparam Dimension: the dimension to compare on -*/ -template -< - typename Tag, - int Direction, - typename Point, - typename CoordinateSystem, - std::size_t Dimension -> -struct strategy_compare +namespace strategy { namespace compare { - typedef strategy::not_implemented type; + + +struct less +{ + template <typename T1, typename T2> + static inline bool apply(T1 const& l, T2 const& r) + { + return l < r; + } }; +struct greater +{ + template <typename T1, typename T2> + static inline bool apply(T1 const& l, T2 const& r) + { + return l > r; + } +}; + +struct equal_to +{ + template <typename T1, typename T2> + static inline bool apply(T1 const& , T2 const& ) + { + return false; + } +}; + + +#ifndef DOXYGEN_NO_DETAIL +namespace detail +{ -#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS -// For compare we add defaults specializations, -// because they defaultly redirect to std::less / greater / equal_to template < - typename Tag, - typename Point, - typename CoordinateSystem, - std::size_t Dimension + typename ComparePolicy, + std::size_t Dimension, + std::size_t DimensionCount > -struct strategy_compare<Tag, 1, Point, CoordinateSystem, Dimension> +struct compare_loop { - typedef std::less<typename coordinate_type<Point>::type> type; + template <typename Point1, typename Point2> + static inline bool apply(Point1 const& left, Point2 const& right) + { + typename geometry::coordinate_type<Point1>::type const& + cleft = geometry::get<Dimension>(left); + typename geometry::coordinate_type<Point2>::type const& + cright = geometry::get<Dimension>(right); + + if (math::equals(cleft, cright)) + { + return compare_loop + < + ComparePolicy, + Dimension + 1, DimensionCount + >::apply(left, right); + } + else + { + return ComparePolicy::apply(cleft, cright); + } + } }; - template < - typename Tag, - typename Point, - typename CoordinateSystem, - std::size_t Dimension + typename ComparePolicy, + std::size_t DimensionCount > -struct strategy_compare<Tag, -1, Point, CoordinateSystem, Dimension> +struct compare_loop<ComparePolicy, DimensionCount, DimensionCount> { - typedef std::greater<typename coordinate_type<Point>::type> type; + template <typename Point1, typename Point2> + static inline bool apply(Point1 const& , Point2 const& ) + { + // On coming here, points are equal. + // Return false for less/greater. + return false; + } }; - template < - typename Tag, - typename Point, - typename CoordinateSystem, - std::size_t Dimension + std::size_t DimensionCount > -struct strategy_compare<Tag, 0, Point, CoordinateSystem, Dimension> +struct compare_loop<strategy::compare::equal_to, DimensionCount, DimensionCount> { - typedef std::equal_to<typename coordinate_type<Point>::type> type; + template <typename Point1, typename Point2> + static inline bool apply(Point1 const& , Point2 const& ) + { + // On coming here, points are equal. + // Return true for equal_to. + return true; + } }; - -#endif +} // namespace detail +#endif // DOXYGEN_NO_DETAIL -namespace strategy { namespace compare +template +< + typename ComparePolicy, + int Dimension = -1 +> +struct cartesian { + template <typename Point1, typename Point2> + static inline bool apply(Point1 const& left, Point2 const& right) + { + return compare::detail::compare_loop + < + ComparePolicy, Dimension, Dimension + 1 + >::apply(left, right); + } +}; - -/*! - \brief Default strategy, indicates the default strategy for comparisons - \details The default strategy for comparisons defer in most cases - to std::less (for ascending) and std::greater (for descending). - However, if a spherical coordinate system is used, and comparison - is done on longitude, it will take another strategy handling circular -*/ -struct default_strategy {}; - - -#ifndef DOXYGEN_NO_DETAIL -namespace detail +#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS +template +< + typename ComparePolicy +> +struct cartesian<ComparePolicy, -1> { + template <typename Point1, typename Point2> + static inline bool apply(Point1 const& left, Point2 const& right) + { + return compare::detail::compare_loop + < + ComparePolicy, + 0, + boost::mpl::min + < + geometry::dimension<Point1>, + geometry::dimension<Point2> + >::type::value + >::apply(left, right); + } +}; +#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS -template <typename Type> -struct is_default : boost::false_type -{}; - - -template <> -struct is_default<default_strategy> : boost::true_type -{}; +namespace services +{ -/*! - \brief Meta-function to select strategy - \details If "default_strategy" is specified, it will take the - traits-registered class for the specified coordinate system. - If another strategy is explicitly specified, it takes that one. -*/ template < - typename Strategy, - int Direction, - typename Point, - std::size_t Dimension + typename ComparePolicy, + typename Point1, + typename Point2 = Point1, + int Dimension = -1, + typename CSTag1 = typename cs_tag<Point1>::type, + typename CSTag2 = typename cs_tag<Point2>::type > -struct select_strategy +struct default_strategy { - typedef typename - boost::mpl::if_ - < - is_default<Strategy>, - typename strategy_compare - < - typename cs_tag<Point>::type, - Direction, - Point, - typename coordinate_system<Point>::type, - Dimension - >::type, - Strategy - >::type type; + BOOST_MPL_ASSERT_MSG + ( + false, + NOT_IMPLEMENTED_FOR_THESE_TYPES, + (types<CSTag1, CSTag2>) + ); }; -} // namespace detail -#endif // DOXYGEN_NO_DETAIL + +template <typename ComparePolicy, typename Point1, typename Point2, int Dimension> +struct default_strategy<ComparePolicy, Point1, Point2, Dimension, cartesian_tag, cartesian_tag> +{ + typedef compare::cartesian<ComparePolicy, Dimension> type; +}; + + +} // namespace services -}} // namespace strategy::compare +}} // namespace strategy compare }} // namespace boost::geometry diff --git a/boost/geometry/strategies/geographic/area.hpp b/boost/geometry/strategies/geographic/area.hpp index 44dc2e6945..40d6c8243e 100644 --- a/boost/geometry/strategies/geographic/area.hpp +++ b/boost/geometry/strategies/geographic/area.hpp @@ -15,12 +15,11 @@ #include <boost/geometry/core/srs.hpp> #include <boost/geometry/formulas/area_formulas.hpp> -#include <boost/geometry/formulas/flattening.hpp> +#include <boost/geometry/formulas/authalic_radius_sqr.hpp> +#include <boost/geometry/formulas/eccentricity_sqr.hpp> #include <boost/geometry/strategies/geographic/parameters.hpp> -#include <boost/math/special_functions/atanh.hpp> - namespace boost { namespace geometry { @@ -88,31 +87,16 @@ protected : inline spheroid_constants(Spheroid const& spheroid) : m_spheroid(spheroid) , m_a2(math::sqr(get_radius<0>(spheroid))) - , m_e2(formula::flattening<CT>(spheroid) - * (CT(2.0) - CT(formula::flattening<CT>(spheroid)))) + , m_e2(formula::eccentricity_sqr<CT>(spheroid)) , m_ep2(m_e2 / (CT(1.0) - m_e2)) , m_ep(math::sqrt(m_ep2)) - , m_c2(authalic_radius(spheroid, m_a2, m_e2)) + , m_c2(formula_dispatch::authalic_radius_sqr + < + CT, Spheroid, srs_spheroid_tag + >::apply(m_a2, m_e2)) {} }; - static inline CT authalic_radius(Spheroid const& sph, CT const& a2, CT const& e2) - { - CT const c0 = 0; - - if (math::equals(e2, c0)) - { - return a2; - } - - CT const sqrt_e2 = math::sqrt(e2); - CT const c2 = 2; - - return (a2 / c2) + - ((math::sqr(get_radius<2>(sph)) * boost::math::atanh(sqrt_e2)) - / (c2 * sqrt_e2)); - } - struct area_sums { CT m_excess_sum; @@ -190,8 +174,10 @@ public : state.m_correction_sum += result.ellipsoidal_term; // Keep track whenever a segment crosses the prime meridian - geometry::formula::area_formulas<CT> - ::crosses_prime_meridian(p1, p2, state); + if (area_formulas::crosses_prime_meridian(p1, p2)) + { + state.m_crosses_prime_meridian++; + } } } diff --git a/boost/geometry/strategies/geographic/distance.hpp b/boost/geometry/strategies/geographic/distance.hpp index d3656f449c..98ee46a369 100644 --- a/boost/geometry/strategies/geographic/distance.hpp +++ b/boost/geometry/strategies/geographic/distance.hpp @@ -5,6 +5,7 @@ // This file was modified by Oracle on 2014-2017. // Modifications copyright (c) 2014-2017 Oracle and/or its affiliates. +// Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle // Use, modification and distribution is subject to the Boost Software License, @@ -22,12 +23,14 @@ #include <boost/geometry/core/srs.hpp> #include <boost/geometry/formulas/andoyer_inverse.hpp> +#include <boost/geometry/formulas/elliptic_arc_length.hpp> #include <boost/geometry/formulas/flattening.hpp> #include <boost/geometry/strategies/distance.hpp> #include <boost/geometry/strategies/geographic/parameters.hpp> #include <boost/geometry/util/math.hpp> +#include <boost/geometry/util/normalize_spheroidal_coordinates.hpp> #include <boost/geometry/util/promote_floating_point.hpp> #include <boost/geometry/util/select_calculation_type.hpp> @@ -70,17 +73,41 @@ public : : m_spheroid(spheroid) {} + template <typename CT> + static inline CT apply(CT lon1, CT lat1, CT lon2, CT lat2, + Spheroid const& spheroid) + { + typedef typename formula::elliptic_arc_length + < + CT, strategy::default_order<FormulaPolicy>::value + > elliptic_arc_length; + + typename elliptic_arc_length::result res = + elliptic_arc_length::apply(lon1, lat1, lon2, lat2, spheroid); + + if (res.meridian) + { + return res.distance; + } + + return FormulaPolicy::template inverse + < + CT, true, false, false, false, false + >::apply(lon1, lat1, lon2, lat2, spheroid).distance; + } + template <typename Point1, typename Point2> inline typename calculation_type<Point1, Point2>::type apply(Point1 const& point1, Point2 const& point2) const { - return FormulaPolicy::template inverse - < - typename calculation_type<Point1, Point2>::type, - true, false, false, false, false - >::apply(get_as_radian<0>(point1), get_as_radian<1>(point1), - get_as_radian<0>(point2), get_as_radian<1>(point2), - m_spheroid).distance; + typedef typename calculation_type<Point1, Point2>::type CT; + + CT lon1 = get_as_radian<0>(point1); + CT lat1 = get_as_radian<1>(point1); + CT lon2 = get_as_radian<0>(point2); + CT lat2 = get_as_radian<1>(point2); + + return apply(lon1, lat1, lon2, lat2, m_spheroid); } inline Spheroid const& model() const diff --git a/boost/geometry/strategies/geographic/distance_cross_track.hpp b/boost/geometry/strategies/geographic/distance_cross_track.hpp new file mode 100644 index 0000000000..799208a96d --- /dev/null +++ b/boost/geometry/strategies/geographic/distance_cross_track.hpp @@ -0,0 +1,641 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2016-2017, Oracle and/or its affiliates. + +// Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle + +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_DISTANCE_CROSS_TRACK_HPP +#define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_DISTANCE_CROSS_TRACK_HPP + +#include <algorithm> + +#include <boost/config.hpp> +#include <boost/concept_check.hpp> +#include <boost/mpl/if.hpp> +#include <boost/type_traits/is_void.hpp> + +#include <boost/geometry/core/cs.hpp> +#include <boost/geometry/core/access.hpp> +#include <boost/geometry/core/radian_access.hpp> +#include <boost/geometry/core/tags.hpp> + +#include <boost/geometry/strategies/distance.hpp> +#include <boost/geometry/strategies/concepts/distance_concept.hpp> +#include <boost/geometry/strategies/spherical/distance_haversine.hpp> +#include <boost/geometry/strategies/geographic/azimuth.hpp> +#include <boost/geometry/strategies/geographic/parameters.hpp> + +#include <boost/geometry/formulas/vincenty_direct.hpp> + +#include <boost/geometry/util/math.hpp> +#include <boost/geometry/util/promote_floating_point.hpp> +#include <boost/geometry/util/select_calculation_type.hpp> +#include <boost/geometry/util/normalize_spheroidal_coordinates.hpp> + +#include <boost/geometry/formulas/result_direct.hpp> +#include <boost/geometry/formulas/mean_radius.hpp> + +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK +# include <boost/geometry/io/dsv/write.hpp> +#endif + +#ifndef BOOST_GEOMETRY_DETAIL_POINT_SEGMENT_DISTANCE_MAX_STEPS +#define BOOST_GEOMETRY_DETAIL_POINT_SEGMENT_DISTANCE_MAX_STEPS 100 +#endif + +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK +#include <iostream> +#endif + +namespace boost { namespace geometry +{ + +namespace strategy { namespace distance +{ + +/*! +\brief Strategy functor for distance point to segment calculation on ellipsoid + Algorithm uses direct and inverse geodesic problems as subroutines. + The algorithm approximates the distance by an iterative Newton method. +\ingroup strategies +\details Class which calculates the distance of a point to a segment, for points +on the ellipsoid +\see C.F.F.Karney - Geodesics on an ellipsoid of revolution, + https://arxiv.org/abs/1102.1215 +\tparam FormulaPolicy underlying point-point distance strategy +\tparam Spheroid is the spheroidal model used +\tparam CalculationType \tparam_calculation +\tparam EnableClosestPoint computes the closest point on segment if true +*/ +template +< + typename FormulaPolicy = strategy::andoyer, + typename Spheroid = srs::spheroid<double>, + typename CalculationType = void, + bool EnableClosestPoint = false +> +class geographic_cross_track +{ +public : + template <typename Point, typename PointOfSegment> + struct return_type + : promote_floating_point + < + typename select_calculation_type + < + Point, + PointOfSegment, + CalculationType + >::type + > + {}; + + struct distance_strategy + { + typedef geographic<FormulaPolicy, Spheroid, CalculationType> type; + }; + + inline typename distance_strategy::type get_distance_strategy() const + { + typedef typename distance_strategy::type distance_type; + return distance_type(m_spheroid); + } + + explicit geographic_cross_track(Spheroid const& spheroid = Spheroid()) + : m_spheroid(spheroid) + {} + + template <typename Point, typename PointOfSegment> + inline typename return_type<Point, PointOfSegment>::type + apply(Point const& p, PointOfSegment const& sp1, PointOfSegment const& sp2) const + { + typedef typename coordinate_system<Point>::type::units units_type; + + return (apply<units_type>(get<0>(sp1), get<1>(sp1), + get<0>(sp2), get<1>(sp2), + get<0>(p), get<1>(p), + m_spheroid)).distance; + } + +private : + + template <typename CT> + struct result_distance_point_segment + { + result_distance_point_segment() + : distance(0) + , closest_point_lon(0) + , closest_point_lat(0) + {} + + CT distance; + CT closest_point_lon; + CT closest_point_lat; + }; + + template <typename CT> + result_distance_point_segment<CT> + static inline non_iterative_case(CT lon, CT lat, CT distance) + { + result_distance_point_segment<CT> result; + result.distance = distance; + + if (EnableClosestPoint) + { + result.closest_point_lon = lon; + result.closest_point_lat = lat; + } + return result; + } + + template <typename CT> + result_distance_point_segment<CT> + static inline non_iterative_case(CT lon1, CT lat1, //p1 + CT lon2, CT lat2, //p2 + Spheroid const& spheroid) + { + CT distance = geometry::strategy::distance::geographic<FormulaPolicy, Spheroid, CT> + ::apply(lon1, lat1, lon2, lat2, spheroid); + + return non_iterative_case(lon1, lat1, distance); + } + + template <typename CT> + CT static inline normalize(CT g4) + { + CT const pi = math::pi<CT>(); + if (g4 < 0 && g4 < -pi)//close to -270 + { + return g4 + 1.5 * pi; + } + else if (g4 > 0 && g4 > pi)//close to 270 + { + return - g4 + 1.5 * pi; + } + else if (g4 < 0 && g4 > -pi)//close to -90 + { + return -g4 - pi/2; + } + return g4 - pi/2; + } + + template <typename Units, typename CT> + result_distance_point_segment<CT> + static inline apply(CT lon1, CT lat1, //p1 + CT lon2, CT lat2, //p2 + CT lon3, CT lat3, //query point p3 + Spheroid const& spheroid) + { + typedef typename FormulaPolicy::template inverse<CT, true, false, false, true, true> + inverse_distance_quantities_type; + typedef typename FormulaPolicy::template inverse<CT, false, true, false, false, false> + inverse_azimuth_type; + typedef typename FormulaPolicy::template inverse<CT, false, true, true, false, false> + inverse_azimuth_reverse_type; + typedef typename FormulaPolicy::template direct<CT, true, false, false, false> + direct_distance_type; + + CT const earth_radius = geometry::formula::mean_radius<CT>(spheroid); + + result_distance_point_segment<CT> result; + + // Constants + CT const f = geometry::formula::flattening<CT>(spheroid); + CT const pi = math::pi<CT>(); + CT const half_pi = pi / CT(2); + CT const c0 = CT(0); + + // Convert to radians + lon1 = math::as_radian<Units>(lon1); + lat1 = math::as_radian<Units>(lat1); + lon2 = math::as_radian<Units>(lon2); + lat2 = math::as_radian<Units>(lat2); + lon3 = math::as_radian<Units>(lon3); + lat3 = math::as_radian<Units>(lat3); + + if (lon1 > lon2) + { + std::swap(lon1, lon2); + std::swap(lat1, lat2); + } + + //segment on equator + //Note: antipodal points on equator does not define segment on equator + //but pass by the pole + CT diff = geometry::math::longitude_distance_signed<geometry::radian>(lon1, lon2); + if (math::equals(lat1, c0) && math::equals(lat2, c0) + && !math::equals(math::abs(diff), pi)) + { +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK + std::cout << "Equatorial segment" << std::endl; +#endif + if (lon3 <= lon1) + { + return non_iterative_case(lon1, lat1, lon3, lat3, spheroid); + } + if (lon3 >= lon2) + { + return non_iterative_case(lon2, lat2, lon3, lat3, spheroid); + } + return non_iterative_case(lon3, lat1, lon3, lat3, spheroid); + } + + if (math::equals(math::abs(diff), pi)) + { +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK + std::cout << "Meridian segment" << std::endl; +#endif + result_distance_point_segment<CT> d1 = apply<geometry::radian>(lon1, lat1, lon1, half_pi, lon3, lat3, spheroid); + result_distance_point_segment<CT> d2 = apply<geometry::radian>(lon2, lat2, lon2, half_pi, lon3, lat3, spheroid); + if (d1.distance < d2.distance) + { + return d1; + } + else + { + return d2; + } + } + + CT d1 = geometry::strategy::distance::geographic<FormulaPolicy, Spheroid, CT> + ::apply(lon1, lat1, lon3, lat3, spheroid); + + CT d3 = geometry::strategy::distance::geographic<FormulaPolicy, Spheroid, CT> + ::apply(lon1, lat1, lon2, lat2, spheroid); + + if (geometry::math::equals(d3, c0)) + { +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK + std::cout << "Degenerate segment" << std::endl; + std::cout << "distance between points=" << d1 << std::endl; +#endif + return non_iterative_case(lon1, lat2, d1); + } + + CT d2 = geometry::strategy::distance::geographic<FormulaPolicy, Spheroid, CT> + ::apply(lon2, lat2, lon3, lat3, spheroid); + + // Compute a12 (GEO) + geometry::formula::result_inverse<CT> res12 = + inverse_azimuth_reverse_type::apply(lon1, lat1, lon2, lat2, spheroid); + CT a12 = res12.azimuth; + CT a13 = inverse_azimuth_type::apply(lon1, lat1, lon3, lat3, spheroid).azimuth; + + CT a312 = a13 - a12; + + if (geometry::math::equals(a312, c0)) + { +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK + std::cout << "point on segment" << std::endl; +#endif + return non_iterative_case(lon3, lat3, c0); + } + + CT projection1 = cos( a312 ) * d1 / d3; + +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK + std::cout << "segment=(" << lon1 * math::r2d<CT>(); + std::cout << "," << lat1 * math::r2d<CT>(); + std::cout << "),(" << lon2 * math::r2d<CT>(); + std::cout << "," << lat2 * math::r2d<CT>(); + std::cout << ")\np=(" << lon3 * math::r2d<CT>(); + std::cout << "," << lat3 * math::r2d<CT>(); + std::cout << ")\na1=" << a12 * math::r2d<CT>() << std::endl; + std::cout << "a13=" << a13 * math::r2d<CT>() << std::endl; + std::cout << "a312=" << a312 * math::r2d<CT>() << std::endl; + std::cout << "cos(a312)=" << cos(a312) << std::endl; +#endif + if (projection1 < 0.0) + { +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK + std::cout << "projection closer to p1" << std::endl; +#endif + // projection of p3 on geodesic spanned by segment (p1,p2) fall + // outside of segment on the side of p1 + return non_iterative_case(lon1, lat1, lon3, lat3, spheroid); + } + + CT a21 = res12.reverse_azimuth - pi; + CT a23 = inverse_azimuth_type::apply(lon2, lat2, lon3, lat3, spheroid).azimuth; + + CT a321 = a23 - a21; + +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK + std::cout << "a21=" << a21 * math::r2d<CT>() << std::endl; + std::cout << "a23=" << a23 * math::r2d<CT>() << std::endl; + std::cout << "a321=" << a321 * math::r2d<CT>() << std::endl; + std::cout << "cos(a321)=" << cos(a321) << std::endl; +#endif + CT projection2 = cos( a321 ) * d2 / d3; + + if (projection2 < 0.0) + { +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK + std::cout << "projection closer to p2" << std::endl; +#endif + // projection of p3 on geodesic spanned by segment (p1,p2) fall + // outside of segment on the side of p2 + return non_iterative_case(lon2, lat2, lon3, lat3, spheroid); + } + + // Guess s14 (SPHERICAL) + typedef geometry::model::point + < + CT, 2, + geometry::cs::spherical_equatorial<geometry::radian> + > point; + + CT bet1 = atan((1 - f) * tan(lon1)); + CT bet2 = atan((1 - f) * tan(lon2)); + CT bet3 = atan((1 - f) * tan(lon3)); + point p1 = point(bet1, lat1); + point p2 = point(bet2, lat2); + point p3 = point(bet3, lat3); + + geometry::strategy::distance::cross_track<CT> cross_track(earth_radius); + CT s34 = cross_track.apply(p3, p1, p2); + + geometry::strategy::distance::haversine<CT> str(earth_radius); + CT s13 = str.apply(p1, p3); + CT s14 = acos( cos(s13/earth_radius) / cos(s34/earth_radius) ) * earth_radius; + +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK + std::cout << "s34=" << s34 << std::endl; + std::cout << "s13=" << s13 << std::endl; + std::cout << "s14=" << s14 << std::endl; + std::cout << "===============" << std::endl; +#endif + + // Update s14 (using Newton method) + CT prev_distance = 0; + geometry::formula::result_direct<CT> res14; + geometry::formula::result_inverse<CT> res34; + + int counter = 0; // robustness + CT g4; + CT delta_g4; + + do{ + prev_distance = res34.distance; + + // Solve the direct problem to find p4 (GEO) + res14 = direct_distance_type::apply(lon1, lat1, s14, a12, spheroid); + + // Solve an inverse problem to find g4 + // g4 is the angle between segment (p1,p2) and segment (p3,p4) that meet on p4 (GEO) + + CT a4 = inverse_azimuth_type::apply(res14.lon2, res14.lat2, + lon2, lat2, spheroid).azimuth; + res34 = inverse_distance_quantities_type::apply(res14.lon2, res14.lat2, + lon3, lat3, spheroid); + g4 = res34.azimuth - a4; + + delta_g4 = normalize(g4); + + CT M43 = res34.geodesic_scale; // cos(s14/earth_radius) is the spherical limit + CT m34 = res34.reduced_length; + CT der = (M43 / m34) * sin(g4); + s14 = s14 - delta_g4 / der; + +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK + std::cout << "p4=" << res14.lon2 * math::r2d<CT>() << + "," << res14.lat2 * math::r2d<CT>() << std::endl; + std::cout << "delta_g4=" << delta_g4 << std::endl; + std::cout << "g4=" << g4 * math::r2d<CT>() << std::endl; + std::cout << "der=" << der << std::endl; + std::cout << "M43=" << M43 << std::endl; + std::cout << "spherical limit=" << cos(s14/earth_radius) << std::endl; + std::cout << "m34=" << m34 << std::endl; + std::cout << "new_s14=" << s14 << std::endl; + std::cout << std::setprecision(16) << "dist =" << res34.distance << std::endl; + std::cout << "---------end of step " << counter << std::endl<< std::endl; +#endif + result.distance = prev_distance; + +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK + if (g4 == half_pi) + { + std::cout << "Stop msg: g4 == half_pi" << std::endl; + } + if (res34.distance >= prev_distance && prev_distance != 0) + { + std::cout << "Stop msg: res34.distance >= prev_distance" << std::endl; + } + if (delta_g4 == 0) + { + std::cout << "Stop msg: delta_g4 == 0" << std::endl; + } + if (counter == 19) + { + std::cout << "Stop msg: counter" << std::endl; + } +#endif + + } while (g4 != half_pi + && (prev_distance > res34.distance || prev_distance == 0) + && delta_g4 != 0 + && ++counter < BOOST_GEOMETRY_DETAIL_POINT_SEGMENT_DISTANCE_MAX_STEPS ) ; + +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK + std::cout << "distance=" << res34.distance << std::endl; + + point p4(res14.lon2, res14.lat2); + CT s34_sph = str.apply(p4, p3); + + std::cout << "s34(sph) =" << s34_sph << std::endl; + std::cout << "s34(geo) =" + << inverse_distance_quantities_type::apply(get<0>(p4), get<1>(p4), lon3, lat3, spheroid).distance + << ", p4=(" << get<0>(p4) * math::r2d<double>() << "," + << get<1>(p4) * math::r2d<double>() << ")" + << std::endl; + + CT s31 = inverse_distance_quantities_type::apply(lon3, lat3, lon1, lat1, spheroid).distance; + CT s32 = inverse_distance_quantities_type::apply(lon3, lat3, lon2, lat2, spheroid).distance; + + CT a4 = inverse_azimuth_type::apply(get<0>(p4), get<1>(p4), lon2, lat2, spheroid).azimuth; + geometry::formula::result_direct<CT> res4 = direct_distance_type::apply(get<0>(p4), get<1>(p4), .04, a4, spheroid); + CT p4_plus = inverse_distance_quantities_type::apply(res4.lon2, res4.lat2, lon3, lat3, spheroid).distance; + + geometry::formula::result_direct<CT> res1 = direct_distance_type::apply(lon1, lat1, s14-.04, a12, spheroid); + CT p4_minus = inverse_distance_quantities_type::apply(res1.lon2, res1.lat2, lon3, lat3, spheroid).distance; + + std::cout << "s31=" << s31 << "\ns32=" << s32 + << "\np4_plus=" << p4_plus << ", p4=(" << res4.lon2 * math::r2d<double>() << "," << res4.lat2 * math::r2d<double>() << ")" + << "\np4_minus=" << p4_minus << ", p4=(" << res1.lon2 * math::r2d<double>() << "," << res1.lat2 * math::r2d<double>() << ")" + << std::endl; + + if (res34.distance <= p4_plus && res34.distance <= p4_minus) + { + std::cout << "Closest point computed" << std::endl; + } + else + { + std::cout << "There is a closer point nearby" << std::endl; + } +#endif + + return result; + } + + Spheroid m_spheroid; +}; + + + +#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS +namespace services +{ + +//tags +template <typename FormulaPolicy> +struct tag<geographic_cross_track<FormulaPolicy> > +{ + typedef strategy_tag_distance_point_segment type; +}; + +template +< + typename FormulaPolicy, + typename Spheroid +> +struct tag<geographic_cross_track<FormulaPolicy, Spheroid> > +{ + typedef strategy_tag_distance_point_segment type; +}; + +template +< + typename FormulaPolicy, + typename Spheroid, + typename CalculationType +> +struct tag<geographic_cross_track<FormulaPolicy, Spheroid, CalculationType> > +{ + typedef strategy_tag_distance_point_segment type; +}; + + +//return types +template <typename FormulaPolicy, typename P, typename PS> +struct return_type<geographic_cross_track<FormulaPolicy>, P, PS> + : geographic_cross_track<FormulaPolicy>::template return_type<P, PS> +{}; + +template +< + typename FormulaPolicy, + typename Spheroid, + typename P, + typename PS +> +struct return_type<geographic_cross_track<FormulaPolicy, Spheroid>, P, PS> + : geographic_cross_track<FormulaPolicy, Spheroid>::template return_type<P, PS> +{}; + +template +< + typename FormulaPolicy, + typename Spheroid, + typename CalculationType, + typename P, + typename PS +> +struct return_type<geographic_cross_track<FormulaPolicy, Spheroid, CalculationType>, P, PS> + : geographic_cross_track<FormulaPolicy, Spheroid, CalculationType>::template return_type<P, PS> +{}; + +//comparable types +template +< + typename FormulaPolicy, + typename Spheroid, + typename CalculationType +> +struct comparable_type<geographic_cross_track<FormulaPolicy, Spheroid, CalculationType> > +{ + typedef geographic_cross_track + < + FormulaPolicy, Spheroid, CalculationType + > type; +}; + +template +< + typename FormulaPolicy, + typename Spheroid, + typename CalculationType +> +struct get_comparable<geographic_cross_track<FormulaPolicy, Spheroid, CalculationType> > +{ + typedef typename comparable_type + < + geographic_cross_track<FormulaPolicy, Spheroid, CalculationType> + >::type comparable_type; +public : + static inline comparable_type + apply(geographic_cross_track<FormulaPolicy, Spheroid, CalculationType> const& ) + { + return comparable_type(); + } +}; + + +template +< + typename FormulaPolicy, + typename P, + typename PS +> +struct result_from_distance<geographic_cross_track<FormulaPolicy>, P, PS> +{ +private : + typedef typename geographic_cross_track + < + FormulaPolicy + >::template return_type<P, PS>::type return_type; +public : + template <typename T> + static inline return_type + apply(geographic_cross_track<FormulaPolicy> const& , T const& distance) + { + return distance; + } +}; + + +template <typename Point, typename PointOfSegment> +struct default_strategy + < + point_tag, segment_tag, Point, PointOfSegment, + geographic_tag, geographic_tag + > +{ + typedef geographic_cross_track<> type; +}; + + +template <typename PointOfSegment, typename Point> +struct default_strategy + < + segment_tag, point_tag, PointOfSegment, Point, + geographic_tag, geographic_tag + > +{ + typedef typename default_strategy + < + point_tag, segment_tag, Point, PointOfSegment, + geographic_tag, geographic_tag + >::type type; +}; + +} // namespace services +#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + +}} // namespace strategy::distance + +}} // namespace boost::geometry +#endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_DISTANCE_CROSS_TRACK_HPP diff --git a/boost/geometry/strategies/geographic/intersection.hpp b/boost/geometry/strategies/geographic/intersection.hpp index 59a40f281d..e91659d40e 100644 --- a/boost/geometry/strategies/geographic/intersection.hpp +++ b/boost/geometry/strategies/geographic/intersection.hpp @@ -26,6 +26,7 @@ #include <boost/geometry/formulas/andoyer_inverse.hpp> #include <boost/geometry/formulas/sjoberg_intersection.hpp> #include <boost/geometry/formulas/spherical.hpp> +#include <boost/geometry/formulas/unit_spheroid.hpp> #include <boost/geometry/geometries/concepts/point_concept.hpp> #include <boost/geometry/geometries/concepts/segment_concept.hpp> @@ -36,6 +37,7 @@ #include <boost/geometry/strategies/geographic/distance.hpp> #include <boost/geometry/strategies/geographic/envelope_segment.hpp> #include <boost/geometry/strategies/geographic/parameters.hpp> +#include <boost/geometry/strategies/geographic/point_in_poly_winding.hpp> #include <boost/geometry/strategies/geographic/side.hpp> #include <boost/geometry/strategies/intersection.hpp> #include <boost/geometry/strategies/intersection_result.hpp> @@ -78,11 +80,12 @@ struct geographic_segments template <typename Geometry1, typename Geometry2> struct point_in_geometry_strategy { - typedef strategy::within::winding + typedef strategy::within::geographic_winding < typename point_type<Geometry1>::type, typename point_type<Geometry2>::type, - side_strategy_type, + FormulaPolicy, + Spheroid, CalculationType > type; }; @@ -95,7 +98,7 @@ struct geographic_segments < Geometry1, Geometry2 >::type strategy_type; - return strategy_type(get_side_strategy()); + return strategy_type(m_spheroid); } template <typename Geometry> @@ -289,10 +292,12 @@ private: typedef typename select_calculation_type <Segment1, Segment2, CalculationType>::type calc_t; + typedef srs::spheroid<calc_t> spheroid_type; + static const calc_t c0 = 0; // normalized spheroid - srs::spheroid<calc_t> spheroid = normalized_spheroid<calc_t>(m_spheroid); + spheroid_type spheroid = formula::unit_spheroid<spheroid_type>(m_spheroid); // TODO: check only 2 first coordinates here? using geometry::detail::equals::equals_point_point; @@ -958,14 +963,6 @@ private: ip_flag; } - template <typename CalcT, typename SpheroidT> - static inline srs::spheroid<CalcT> normalized_spheroid(SpheroidT const& spheroid) - { - return srs::spheroid<CalcT>(CalcT(1), - CalcT(get_radius<2>(spheroid)) // b/a - / CalcT(get_radius<0>(spheroid))); - } - private: Spheroid m_spheroid; }; diff --git a/boost/geometry/strategies/geographic/parameters.hpp b/boost/geometry/strategies/geographic/parameters.hpp index 5638db50fa..4896e42e8d 100644 --- a/boost/geometry/strategies/geographic/parameters.hpp +++ b/boost/geometry/strategies/geographic/parameters.hpp @@ -12,7 +12,9 @@ #include <boost/geometry/formulas/andoyer_inverse.hpp> +#include <boost/geometry/formulas/thomas_direct.hpp> #include <boost/geometry/formulas/thomas_inverse.hpp> +#include <boost/geometry/formulas/vincenty_direct.hpp> #include <boost/geometry/formulas/vincenty_inverse.hpp> #include <boost/mpl/assert.hpp> @@ -24,6 +26,23 @@ namespace boost { namespace geometry { namespace strategy struct andoyer { + //TODO: this should be replaced by an andoyer direct formula + template + < + typename CT, + bool EnableCoordinates = true, + bool EnableReverseAzimuth = false, + bool EnableReducedLength = false, + bool EnableGeodesicScale = false + > + struct direct + : formula::thomas_direct + < + CT, EnableCoordinates, EnableReverseAzimuth, + EnableReducedLength, EnableGeodesicScale + > + {}; + template < typename CT, @@ -48,6 +67,22 @@ struct thomas template < typename CT, + bool EnableCoordinates = true, + bool EnableReverseAzimuth = false, + bool EnableReducedLength = false, + bool EnableGeodesicScale = false + > + struct direct + : formula::thomas_direct + < + CT, EnableCoordinates, EnableReverseAzimuth, + EnableReducedLength, EnableGeodesicScale + > + {}; + + template + < + typename CT, bool EnableDistance, bool EnableAzimuth, bool EnableReverseAzimuth = false, @@ -69,6 +104,22 @@ struct vincenty template < typename CT, + bool EnableCoordinates = true, + bool EnableReverseAzimuth = false, + bool EnableReducedLength = false, + bool EnableGeodesicScale = false + > + struct direct + : formula::vincenty_direct + < + CT, EnableCoordinates, EnableReverseAzimuth, + EnableReducedLength, EnableGeodesicScale + > + {}; + + template + < + typename CT, bool EnableDistance, bool EnableAzimuth, bool EnableReverseAzimuth = false, diff --git a/boost/geometry/strategies/geographic/point_in_poly_winding.hpp b/boost/geometry/strategies/geographic/point_in_poly_winding.hpp new file mode 100644 index 0000000000..95a1961474 --- /dev/null +++ b/boost/geometry/strategies/geographic/point_in_poly_winding.hpp @@ -0,0 +1,80 @@ +// Boost.Geometry + +// Copyright (c) 2017 Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_GEOMETRY_STRATEGY_GEOGRAPHIC_POINT_IN_POLY_WINDING_HPP +#define BOOST_GEOMETRY_STRATEGY_GEOGRAPHIC_POINT_IN_POLY_WINDING_HPP + + +#include <boost/geometry/strategies/geographic/side.hpp> +#include <boost/geometry/strategies/spherical/point_in_poly_winding.hpp> + + +namespace boost { namespace geometry +{ + +namespace strategy { namespace within +{ + + +/*! +\brief Within detection using winding rule in geographic coordinate system. +\ingroup strategies +\tparam Point \tparam_point +\tparam PointOfSegment \tparam_segment_point +\tparam FormulaPolicy Geodesic formula policy +\tparam Spheroid Spheroid model +\tparam CalculationType \tparam_calculation + +\qbk{ +[heading See also] +[link geometry.reference.algorithms.within.within_3_with_strategy within (with strategy)] +} + */ +template +< + typename Point, + typename PointOfSegment = Point, + typename FormulaPolicy = strategy::andoyer, + typename Spheroid = srs::spheroid<double>, + typename CalculationType = void +> +class geographic_winding + : public within::detail::spherical_winding_base + < + Point, + PointOfSegment, + side::geographic<FormulaPolicy, Spheroid, CalculationType>, + CalculationType + > +{ + typedef within::detail::spherical_winding_base + < + Point, + PointOfSegment, + side::geographic<FormulaPolicy, Spheroid, CalculationType>, + CalculationType + > base_t; + +public: + geographic_winding() + {} + + explicit geographic_winding(Spheroid const& model) + : base_t(model) + {} +}; + + +}} // namespace strategy::within + + +}} // namespace boost::geometry + + +#endif // BOOST_GEOMETRY_STRATEGY_GEOGRAPHIC_POINT_IN_POLY_WINDING_HPP diff --git a/boost/geometry/strategies/spherical/area.hpp b/boost/geometry/strategies/spherical/area.hpp index 206b734548..1ab61b644f 100644 --- a/boost/geometry/strategies/spherical/area.hpp +++ b/boost/geometry/strategies/spherical/area.hpp @@ -127,14 +127,15 @@ public : { if (! geometry::math::equals(get<0>(p1), get<0>(p2))) { + typedef geometry::formula::area_formulas<CT> area_formulas; - state.m_sum += geometry::formula::area_formulas - <CT>::template spherical<LongSegment>(p1, p2); + state.m_sum += area_formulas::template spherical<LongSegment>(p1, p2); // Keep track whenever a segment crosses the prime meridian - geometry::formula::area_formulas - <CT>::crosses_prime_meridian(p1, p2, state); - + if (area_formulas::crosses_prime_meridian(p1, p2)) + { + state.m_crosses_prime_meridian++; + } } } diff --git a/boost/geometry/strategies/spherical/compare.hpp b/boost/geometry/strategies/spherical/compare.hpp new file mode 100644 index 0000000000..26163f7406 --- /dev/null +++ b/boost/geometry/strategies/spherical/compare.hpp @@ -0,0 +1,321 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2017. +// Modifications copyright (c) 2017, Oracle and/or its affiliates. + +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_GEOMETRY_STRATEGIES_SPHERICAL_COMPARE_HPP +#define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_COMPARE_HPP + + +#include <boost/mpl/if.hpp> +#include <boost/mpl/min.hpp> + +#include <boost/type_traits/is_same.hpp> + +#include <boost/geometry/core/access.hpp> +#include <boost/geometry/core/coordinate_dimension.hpp> +#include <boost/geometry/core/coordinate_system.hpp> +#include <boost/geometry/core/coordinate_type.hpp> +#include <boost/geometry/core/cs.hpp> +#include <boost/geometry/core/tags.hpp> + +#include <boost/geometry/strategies/compare.hpp> + +#include <boost/geometry/util/math.hpp> +#include <boost/geometry/util/normalize_spheroidal_coordinates.hpp> + + +namespace boost { namespace geometry +{ + + +namespace strategy { namespace compare +{ + + +#ifndef DOXYGEN_NO_DETAIL +namespace detail +{ + +template <std::size_t I, typename P> +static inline typename geometry::coordinate_type<P>::type +get(P const& p, boost::true_type /*same units*/) +{ + return geometry::get<I>(p); +} + +template <std::size_t I, typename P> +static inline typename geometry::coordinate_type<P>::type +get(P const& p, boost::false_type /*different units*/) +{ + return geometry::get_as_radian<I>(p); +} + +template +< + typename ComparePolicy, + typename Point1, + typename Point2, + std::size_t DimensionCount +> +struct spherical_latitude +{ + typedef typename geometry::coordinate_type<Point1>::type coordinate1_type; + typedef typename geometry::coordinate_system<Point1>::type::units units1_type; + typedef typename geometry::coordinate_type<Point2>::type coordinate2_type; + typedef typename geometry::coordinate_system<Point2>::type::units units2_type; + typedef typename boost::is_same<units1_type, units2_type>::type same_units_type; + + template <typename T1, typename T2> + static inline bool apply(Point1 const& left, Point2 const& right, + T1 const& l1, T2 const& r1) + { + // latitudes equal + if (math::equals(l1, r1)) + { + return compare::detail::compare_loop + < + ComparePolicy, 2, DimensionCount + >::apply(left, right); + } + else + { + return ComparePolicy::apply(l1, r1); + } + } + + static inline bool apply(Point1 const& left, Point2 const& right) + { + coordinate1_type const& l1 = compare::detail::get<1>(left, same_units_type()); + coordinate2_type const& r1 = compare::detail::get<1>(right, same_units_type()); + + return apply(left, right, l1, r1); + } +}; + +template +< + typename ComparePolicy, + typename Point1, + typename Point2 +> +struct spherical_latitude<ComparePolicy, Point1, Point2, 1> +{ + template <typename T1, typename T2> + static inline bool apply(Point1 const& left, Point2 const& right, + T1 const& , T2 const& ) + { + return apply(left, right); + } + + static inline bool apply(Point1 const& left, Point2 const& right) + { + return compare::detail::compare_loop + < + ComparePolicy, 1, 1 + >::apply(left, right); + } +}; + +template +< + typename ComparePolicy, + typename Point1, + typename Point2, + std::size_t DimensionCount +> +struct spherical_longitude +{ + typedef typename geometry::coordinate_type<Point1>::type coordinate1_type; + typedef typename geometry::coordinate_system<Point1>::type::units units1_type; + typedef typename geometry::coordinate_type<Point2>::type coordinate2_type; + typedef typename geometry::coordinate_system<Point2>::type::units units2_type; + typedef typename boost::is_same<units1_type, units2_type>::type same_units_type; + typedef typename boost::mpl::if_<same_units_type, units1_type, geometry::radian>::type units_type; + + static const bool is_equatorial = ! boost::is_same + < + typename geometry::cs_tag<Point1>::type, + geometry::spherical_polar_tag + >::value; + + static inline bool are_both_at_antimeridian(coordinate1_type const& l0, + coordinate2_type const& r0, + bool & is_left_at, + bool & is_right_at) + { + is_left_at = math::is_longitude_antimeridian<units_type>(l0); + is_right_at = math::is_longitude_antimeridian<units_type>(r0); + return is_left_at && is_right_at; + } + + static inline bool apply(Point1 const& left, Point2 const& right) + { + // if units are different the coordinates are in radians + coordinate1_type const& l0 = compare::detail::get<0>(left, same_units_type()); + coordinate2_type const& r0 = compare::detail::get<0>(right, same_units_type()); + coordinate1_type const& l1 = compare::detail::get<1>(left, same_units_type()); + coordinate2_type const& r1 = compare::detail::get<1>(right, same_units_type()); + + bool is_left_at_antimeridian = false; + bool is_right_at_antimeridian = false; + + // longitudes equal + if (math::equals(l0, r0) + // both at antimeridian + || are_both_at_antimeridian(l0, r0, is_left_at_antimeridian, is_right_at_antimeridian) + // both at pole + || (math::equals(l1, r1) + && math::is_latitude_pole<units_type, is_equatorial>(l1))) + { + return spherical_latitude + < + ComparePolicy, Point1, Point2, DimensionCount + >::apply(left, right, l1, r1); + } + // if left is at antimeridian and right is not at antimeridian + // then left is greater than right + else if (is_left_at_antimeridian) + { + // less/equal_to -> false, greater -> true + return ComparePolicy::apply(1, 0); + } + // if right is at antimeridian and left is not at antimeridian + // then left is lesser than right + else if (is_right_at_antimeridian) + { + // less -> true, equal_to/greater -> false + return ComparePolicy::apply(0, 1); + } + else + { + return ComparePolicy::apply(l0, r0); + } + } +}; + + +} // namespace detail +#endif // DOXYGEN_NO_DETAIL + + +/*! +\brief Compare strategy for spherical coordinates +\ingroup strategies +\tparam Point point-type +\tparam Dimension dimension +*/ +template +< + typename ComparePolicy, + int Dimension = -1 +> +struct spherical + : cartesian<ComparePolicy, Dimension> +{}; + +#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS +// all dimensions starting from longitude +template <typename ComparePolicy> +struct spherical<ComparePolicy, -1> +{ + template <typename Point1, typename Point2> + static inline bool apply(Point1 const& left, Point2 const& right) + { + return compare::detail::spherical_longitude + < + ComparePolicy, + Point1, + Point2, + boost::mpl::min + < + geometry::dimension<Point1>, + geometry::dimension<Point2> + >::type::value + >::apply(left, right); + } +}; + +// only longitudes (and latitudes to check poles) +template <typename ComparePolicy> +struct spherical<ComparePolicy, 0> +{ + template <typename Point1, typename Point2> + static inline bool apply(Point1 const& left, Point2 const& right) + { + return compare::detail::spherical_longitude + < + ComparePolicy, Point1, Point2, 1 + >::apply(left, right); + } +}; + +// only latitudes +template <typename ComparePolicy> +struct spherical<ComparePolicy, 1> +{ + template <typename Point1, typename Point2> + static inline bool apply(Point1 const& left, Point2 const& right) + { + return compare::detail::spherical_latitude + < + ComparePolicy, Point1, Point2, 2 + >::apply(left, right); + } +}; + +#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + + +namespace services +{ + + +template <typename ComparePolicy, typename Point1, typename Point2, int Dimension> +struct default_strategy + < + ComparePolicy, Point1, Point2, Dimension, + spherical_polar_tag, spherical_polar_tag + > +{ + typedef compare::spherical<ComparePolicy, Dimension> type; +}; + +template <typename ComparePolicy, typename Point1, typename Point2, int Dimension> +struct default_strategy + < + ComparePolicy, Point1, Point2, Dimension, + spherical_equatorial_tag, spherical_equatorial_tag + > +{ + typedef compare::spherical<ComparePolicy, Dimension> type; +}; + +template <typename ComparePolicy, typename Point1, typename Point2, int Dimension> +struct default_strategy + < + ComparePolicy, Point1, Point2, Dimension, + geographic_tag, geographic_tag + > +{ + typedef compare::spherical<ComparePolicy, Dimension> type; +}; + + +} // namespace services + + +}} // namespace strategy::compare + + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_COMPARE_HPP diff --git a/boost/geometry/strategies/spherical/compare_circular.hpp b/boost/geometry/strategies/spherical/compare_circular.hpp deleted file mode 100644 index 2f890dfd87..0000000000 --- a/boost/geometry/strategies/spherical/compare_circular.hpp +++ /dev/null @@ -1,152 +0,0 @@ -// Boost.Geometry (aka GGL, Generic Geometry Library) - -// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. - -// Use, modification and distribution is subject to the Boost Software License, -// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at -// http://www.boost.org/LICENSE_1_0.txt) - -#ifndef BOOST_GEOMETRY_STRATEGIES_SPHERICAL_COMPARE_SPHERICAL_HPP -#define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_COMPARE_SPHERICAL_HPP - -#include <boost/math/constants/constants.hpp> - -#include <boost/geometry/core/cs.hpp> -#include <boost/geometry/core/tags.hpp> -#include <boost/geometry/strategies/compare.hpp> -#include <boost/geometry/util/math.hpp> - - -namespace boost { namespace geometry -{ - - -namespace strategy { namespace compare -{ - - -#ifndef DOXYGEN_NO_DETAIL -namespace detail -{ - -template <typename Units> -struct shift -{ -}; - -template <> -struct shift<degree> -{ - static inline double full() { return 360.0; } - static inline double half() { return 180.0; } -}; - -template <> -struct shift<radian> -{ - static inline double full() { return 2.0 * boost::math::constants::pi<double>(); } - static inline double half() { return boost::math::constants::pi<double>(); } -}; - -} // namespace detail -#endif - -/*! -\brief Compare (in one direction) strategy for spherical coordinates -\ingroup strategies -\tparam Point point-type -\tparam Dimension dimension -*/ -template <typename CoordinateType, typename Units, typename Compare> -struct circular_comparator -{ - static inline CoordinateType put_in_range(CoordinateType const& c, - double min_border, double max_border) - { - CoordinateType value = c; - while (value < min_border) - { - value += detail::shift<Units>::full(); - } - while (value > max_border) - { - value -= detail::shift<Units>::full(); - } - return value; - } - - inline bool operator()(CoordinateType const& c1, CoordinateType const& c2) const - { - Compare compare; - - // Check situation that one of them is e.g. std::numeric_limits. - static const double full = detail::shift<Units>::full(); - double mx = 10.0 * full; - if (c1 < -mx || c1 > mx || c2 < -mx || c2 > mx) - { - // do normal comparison, using circular is not useful - return compare(c1, c2); - } - - static const double half = full / 2.0; - CoordinateType v1 = put_in_range(c1, -half, half); - CoordinateType v2 = put_in_range(c2, -half, half); - - // Two coordinates on a circle are - // at max <= half a circle away from each other. - // So if it is more, shift origin. - CoordinateType diff = geometry::math::abs(v1 - v2); - if (diff > half) - { - v1 = put_in_range(v1, 0, full); - v2 = put_in_range(v2, 0, full); - } - - return compare(v1, v2); - } -}; - -}} // namespace strategy::compare - -#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS - -// Specialize for the longitude (dim 0) -template -< - typename Point, - template<typename> class CoordinateSystem, - typename Units -> -struct strategy_compare<spherical_polar_tag, 1, Point, CoordinateSystem<Units>, 0> -{ - typedef typename coordinate_type<Point>::type coordinate_type; - typedef strategy::compare::circular_comparator - < - coordinate_type, - Units, - std::less<coordinate_type> - > type; -}; - -template -< - typename Point, - template<typename> class CoordinateSystem, - typename Units -> -struct strategy_compare<spherical_polar_tag, -1, Point, CoordinateSystem<Units>, 0> -{ - typedef typename coordinate_type<Point>::type coordinate_type; - typedef strategy::compare::circular_comparator - < - coordinate_type, - Units, - std::greater<coordinate_type> - > type; -}; - -#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS - -}} // namespace boost::geometry - -#endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_COMPARE_SPHERICAL_HPP diff --git a/boost/geometry/strategies/spherical/distance_cross_track.hpp b/boost/geometry/strategies/spherical/distance_cross_track.hpp index 7daafa4a16..82a2fb3d6f 100644 --- a/boost/geometry/strategies/spherical/distance_cross_track.hpp +++ b/boost/geometry/strategies/spherical/distance_cross_track.hpp @@ -2,9 +2,10 @@ // Copyright (c) 2007-2014 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2014. -// Modifications copyright (c) 2014, Oracle and/or its affiliates. +// This file was modified by Oracle on 2014-2017. +// Modifications copyright (c) 2014-2017, Oracle and/or its affiliates. +// Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle // Use, modification and distribution is subject to the Boost Software License, @@ -26,8 +27,6 @@ #include <boost/geometry/core/radian_access.hpp> #include <boost/geometry/core/tags.hpp> -#include <boost/geometry/algorithms/detail/course.hpp> - #include <boost/geometry/strategies/distance.hpp> #include <boost/geometry/strategies/concepts/distance_concept.hpp> #include <boost/geometry/strategies/spherical/distance_haversine.hpp> @@ -373,19 +372,6 @@ public : typedef typename return_type<Point, PointOfSegment>::type return_type; -#ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK - std::cout << "Course " << dsv(sp1) << " to " << dsv(p) << " " - << crs_AD * geometry::math::r2d<return_type>() << std::endl; - std::cout << "Course " << dsv(sp1) << " to " << dsv(sp2) << " " - << crs_AB * geometry::math::r2d<return_type>() << std::endl; - std::cout << "Course " << dsv(sp2) << " to " << dsv(p) << " " - << crs_BD * geometry::math::r2d << std::endl; - std::cout << "Projection AD-AB " << projection1 << " : " - << d_crs1 * geometry::math::r2d<return_type>() << std::endl; - std::cout << "Projection BD-BA " << projection2 << " : " - << d_crs2 * geometry::math::r2d<return_type>() << std::endl; -#endif - // http://williams.best.vwh.net/avform.htm#XTE return_type d1 = m_strategy.apply(sp1, p); return_type d3 = m_strategy.apply(sp1, sp2); @@ -398,10 +384,25 @@ public : return_type d2 = m_strategy.apply(sp2, p); - return_type crs_AD = geometry::detail::course<return_type>(sp1, p); - return_type crs_AB = geometry::detail::course<return_type>(sp1, sp2); - return_type crs_BA = crs_AB - geometry::math::pi<return_type>(); - return_type crs_BD = geometry::detail::course<return_type>(sp2, p); + return_type lon1 = geometry::get_as_radian<0>(sp1); + return_type lat1 = geometry::get_as_radian<1>(sp1); + return_type lon2 = geometry::get_as_radian<0>(sp2); + return_type lat2 = geometry::get_as_radian<1>(sp2); + return_type lon = geometry::get_as_radian<0>(p); + return_type lat = geometry::get_as_radian<1>(p); + + return_type crs_AD = geometry::formula::spherical_azimuth<return_type, false> + (lon1, lat1, lon, lat).azimuth; + + geometry::formula::result_spherical<return_type> result = + geometry::formula::spherical_azimuth<return_type, true> + (lon1, lat1, lon2, lat2); + return_type crs_AB = result.azimuth; + return_type crs_BA = result.reverse_azimuth - geometry::math::pi<return_type>(); + + return_type crs_BD = geometry::formula::spherical_azimuth<return_type, false> + (lon2, lat2, lon, lat).azimuth; + return_type d_crs1 = crs_AD - crs_AB; return_type d_crs2 = crs_BD - crs_BA; @@ -409,6 +410,24 @@ public : return_type projection1 = cos( d_crs1 ) * d1 / d3; return_type projection2 = cos( d_crs2 ) * d2 / d3; +#ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK + std::cout << "Course " << dsv(sp1) << " to " << dsv(p) << " " + << crs_AD * geometry::math::r2d<return_type>() << std::endl; + std::cout << "Course " << dsv(sp1) << " to " << dsv(sp2) << " " + << crs_AB * geometry::math::r2d<return_type>() << std::endl; + std::cout << "Course " << dsv(sp2) << " to " << dsv(sp1) << " " + << crs_BA * geometry::math::r2d<return_type>() << std::endl; + std::cout << "Course " << dsv(sp2) << " to " << dsv(p) << " " + << crs_BD * geometry::math::r2d<return_type>() << std::endl; + std::cout << "Projection AD-AB " << projection1 << " : " + << d_crs1 * geometry::math::r2d<return_type>() << std::endl; + std::cout << "Projection BD-BA " << projection2 << " : " + << d_crs2 * geometry::math::r2d<return_type>() << std::endl; + std::cout << " d1: " << (d1 ) + << " d2: " << (d2 ) + << std::endl; +#endif + if (projection1 > 0.0 && projection2 > 0.0) { #ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK diff --git a/boost/geometry/strategies/spherical/intersection.hpp b/boost/geometry/strategies/spherical/intersection.hpp index 44b1cc62bf..b5ae878b85 100644 --- a/boost/geometry/strategies/spherical/intersection.hpp +++ b/boost/geometry/strategies/spherical/intersection.hpp @@ -33,7 +33,6 @@ #include <boost/geometry/policies/robustness/segment_ratio.hpp> -#include <boost/geometry/strategies/agnostic/point_in_poly_winding.hpp> #include <boost/geometry/strategies/covered_by.hpp> #include <boost/geometry/strategies/intersection.hpp> #include <boost/geometry/strategies/intersection_result.hpp> @@ -42,6 +41,7 @@ #include <boost/geometry/strategies/spherical/area.hpp> #include <boost/geometry/strategies/spherical/distance_haversine.hpp> #include <boost/geometry/strategies/spherical/envelope_segment.hpp> +#include <boost/geometry/strategies/spherical/point_in_poly_winding.hpp> #include <boost/geometry/strategies/spherical/ssf.hpp> #include <boost/geometry/strategies/within.hpp> @@ -94,11 +94,10 @@ struct ecef_segments template <typename Geometry1, typename Geometry2> struct point_in_geometry_strategy { - typedef strategy::within::winding + typedef strategy::within::spherical_winding < typename point_type<Geometry1>::type, typename point_type<Geometry2>::type, - side_strategy_type, CalculationType > type; }; diff --git a/boost/geometry/strategies/spherical/point_in_poly_winding.hpp b/boost/geometry/strategies/spherical/point_in_poly_winding.hpp new file mode 100644 index 0000000000..0f1a901d13 --- /dev/null +++ b/boost/geometry/strategies/spherical/point_in_poly_winding.hpp @@ -0,0 +1,581 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// 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, 2016, 2017. +// Modifications copyright (c) 2013-2017 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. + +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_GEOMETRY_STRATEGY_SPHERICAL_POINT_IN_POLY_WINDING_HPP +#define BOOST_GEOMETRY_STRATEGY_SPHERICAL_POINT_IN_POLY_WINDING_HPP + + +#include <boost/geometry/core/access.hpp> +#include <boost/geometry/core/coordinate_system.hpp> +#include <boost/geometry/core/cs.hpp> +#include <boost/geometry/core/tags.hpp> + +#include <boost/geometry/util/math.hpp> +#include <boost/geometry/util/select_calculation_type.hpp> +#include <boost/geometry/util/normalize_spheroidal_coordinates.hpp> + +#include <boost/geometry/strategies/covered_by.hpp> +#include <boost/geometry/strategies/side.hpp> +#include <boost/geometry/strategies/spherical/ssf.hpp> +#include <boost/geometry/strategies/within.hpp> + + +namespace boost { namespace geometry +{ + +namespace strategy { namespace within +{ + + +#ifndef DOXYGEN_NO_DETAIL +namespace detail +{ + +template +< + typename Point, + typename PointOfSegment = Point, + typename SideStrategy = typename strategy::side::services::default_strategy + < + typename cs_tag<Point>::type + >::type, + typename CalculationType = void +> +class spherical_winding_base +{ + typedef typename select_calculation_type + < + Point, + PointOfSegment, + CalculationType + >::type calculation_type; + + typedef typename coordinate_system<Point>::type::units units_t; + typedef math::detail::constants_on_spheroid<calculation_type, units_t> constants; + + /*! subclass to keep state */ + class counter + { + int m_count; + //int m_count_n; + int m_count_s; + int m_raw_count; + int m_raw_count_anti; + bool m_touches; + + inline int code() const + { + if (m_touches) + { + return 0; + } + + if (m_raw_count != 0 && m_raw_count_anti != 0) + { + if (m_raw_count > 0) // right, wrap around south pole + { + return (m_count + m_count_s) == 0 ? -1 : 1; + } + else // left, wrap around north pole + { + //return (m_count + m_count_n) == 0 ? -1 : 1; + // m_count_n is 0 + return m_count == 0 ? -1 : 1; + } + } + + return m_count == 0 ? -1 : 1; + } + + public : + friend class spherical_winding_base; + + inline counter() + : m_count(0) + //, m_count_n(0) + , m_count_s(0) + , m_raw_count(0) + , m_raw_count_anti(0) + , m_touches(false) + {} + + }; + + struct count_info + { + explicit count_info(int c = 0, bool ia = false) + : count(c) + , is_anti(ia) + {} + + int count; + bool is_anti; + }; + +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(); + } + + typedef typename SideStrategy::disjoint_strategy_type disjoint_strategy_type; + + inline disjoint_strategy_type get_disjoint_strategy() const + { + return m_side_strategy.get_disjoint_strategy(); + } + + spherical_winding_base() + {} + + template <typename Model> + explicit spherical_winding_base(Model const& model) + : m_side_strategy(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 + { + bool eq1 = false; + bool eq2 = false; + bool s_antipodal = false; + + count_info ci = check_segment(point, s1, s2, state, eq1, eq2, s_antipodal); + if (ci.count != 0) + { + if (! ci.is_anti) + { + int side = 0; + if (ci.count == 1 || ci.count == -1) + { + side = side_equal(point, eq1 ? s1 : s2, ci, s1, s2); + } + else // count == 2 || count == -2 + { + if (! s_antipodal) + { + // 1 left, -1 right + side = m_side_strategy.apply(s1, s2, point); + } + else + { + calculation_type const pi = constants::half_period(); + calculation_type const s1_lat = get<1>(s1); + calculation_type const s2_lat = get<1>(s2); + + side = math::sign(ci.count) + * (pi - s1_lat - s2_lat <= pi // segment goes through north pole + ? -1 // going right all points will be on right side + : 1); // going right all points will be on left side + } + } + + 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 left, 2 for right (or -1/1) + // Side positive thus means RIGHT and LEFTSIDE or LEFT and RIGHTSIDE + // See accompagnying figure (TODO) + if (side * ci.count > 0) + { + state.m_count += ci.count; + } + + state.m_raw_count += ci.count; + } + else + { + // Count negated because the segment is on the other side of the globe + // so it is reversed to match this side of the globe + + // Assuming geometry wraps around north pole, for segments on the other side of the globe + // the point will always be RIGHT+RIGHTSIDE or LEFT+LEFTSIDE, so side*-count always < 0 + //state.m_count_n -= 0; + + // Assuming geometry wraps around south pole, for segments on the other side of the globe + // the point will always be RIGHT+LEFTSIDE or LEFT+RIGHTSIDE, so side*-count always > 0 + state.m_count_s -= ci.count; + + state.m_raw_count_anti -= ci.count; + } + } + return ! state.m_touches; + } + + static inline int result(counter const& state) + { + return state.code(); + } + +private: + + static inline count_info check_segment(Point const& point, + PointOfSegment const& seg1, + PointOfSegment const& seg2, + counter& state, + bool& eq1, bool& eq2, bool& s_antipodal) + { + if (check_touch(point, seg1, seg2, state, eq1, eq2, s_antipodal)) + { + return count_info(0, false); + } + + return calculate_count(point, seg1, seg2, eq1, eq2, s_antipodal); + } + + static inline int check_touch(Point const& point, + PointOfSegment const& seg1, + PointOfSegment const& seg2, + counter& state, + bool& eq1, + bool& eq2, + bool& s_antipodal) + { + calculation_type const c0 = 0; + calculation_type const c2 = 2; + calculation_type const pi = constants::half_period(); + calculation_type const half_pi = pi / c2; + + calculation_type const p_lon = get<0>(point); + calculation_type const s1_lon = get<0>(seg1); + calculation_type const s2_lon = get<0>(seg2); + calculation_type const p_lat = get<1>(point); + calculation_type const s1_lat = get<1>(seg1); + calculation_type const s2_lat = 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 = longitudes_equal(s1_lon, p_lon); + bool eq2_strict = longitudes_equal(s2_lon, p_lon); + bool eq1_anti = false; + bool eq2_anti = false; + + calculation_type const anti_p_lon = p_lon + (p_lon <= c0 ? pi : -pi); + + eq1 = eq1_strict // lon strictly equal to s1 + || (eq1_anti = longitudes_equal(s1_lon, anti_p_lon)) // anti-lon strictly equal to s1 + || math::equals(math::abs(s1_lat), half_pi); // s1 is pole + eq2 = eq2_strict // lon strictly equal to s2 + || (eq2_anti = longitudes_equal(s2_lon, anti_p_lon)) // anti-lon strictly equal to s2 + || math::equals(math::abs(s2_lat), half_pi); // s2 is pole + + // segment overlapping pole + calculation_type const s_lon_diff = math::longitude_distance_signed<units_t>(s1_lon, s2_lon); + s_antipodal = math::equals(s_lon_diff, pi); + if (s_antipodal) + { + eq1 = eq2 = eq1 || eq2; + + // segment overlapping pole and point is pole + if (math::equals(math::abs(p_lat), half_pi)) + { + 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 (! s_antipodal) + { + // p's lat between segment endpoints' lats + if ( (s1_lat <= p_lat && s2_lat >= p_lat) || (s2_lat <= p_lat && s1_lat >= p_lat) ) + { + if (!eq1_anti || !eq2_anti) + { + state.m_touches = true; + } + } + } + else + { + // going through north or south pole? + if (pi - s1_lat - s2_lat <= pi) + { + if ( (eq1_strict && s1_lat <= p_lat) || (eq2_strict && s2_lat <= p_lat) // north + || math::equals(p_lat, half_pi) ) // point on north pole + { + state.m_touches = true; + } + else if (! eq1_strict && ! eq2_strict && math::equals(p_lat, -half_pi) ) // point on south pole + { + return false; + } + } + else // south pole + { + if ( (eq1_strict && s1_lat >= p_lat) || (eq2_strict && s2_lat >= p_lat) // south + || math::equals(p_lat, -half_pi) ) // point on south pole + { + state.m_touches = true; + } + else if (! eq1_strict && ! eq2_strict && math::equals(p_lat, half_pi) ) // point on north pole + { + return false; + } + } + } + + return true; + } + + return false; + } + + // Called if point is not aligned with a vertical segment + static inline count_info calculate_count(Point const& point, + PointOfSegment const& seg1, + PointOfSegment const& seg2, + bool eq1, bool eq2, bool s_antipodal) + { + // 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 + + calculation_type const c0 = 0; + calculation_type const pi = constants::half_period(); + + calculation_type const p = get<0>(point); + calculation_type const s1 = get<0>(seg1); + calculation_type const s2 = get<0>(seg2); + + calculation_type const s1_p = math::longitude_distance_signed<units_t>(s1, p); + + if (s_antipodal) + { + return count_info(s1_p < c0 ? -2 : 2, false); // choose W/E + } + + calculation_type const s1_s2 = math::longitude_distance_signed<units_t>(s1, s2); + + if (eq1 || eq2) // Point on level s1 or s2 + { + return count_info(s1_s2 < c0 ? -1 : 1, // choose W/E + longitudes_equal(p + pi, (eq1 ? s1 : s2))); + } + + // Point between s1 and s2 + if ( math::sign(s1_p) == math::sign(s1_s2) + && math::abs(s1_p) < math::abs(s1_s2) ) + { + return count_info(s1_s2 < c0 ? -2 : 2, false); // choose W/E + } + + calculation_type const s1_p_anti = math::longitude_distance_signed<units_t>(s1, p + pi); + + // Anti-Point between s1 and s2 + if ( math::sign(s1_p_anti) == math::sign(s1_s2) + && math::abs(s1_p_anti) < math::abs(s1_s2) ) + { + return count_info(s1_s2 < c0 ? -2 : 2, true); // choose W/E + } + + return count_info(0, false); + } + + + // 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 + inline int side_equal(Point const& point, + PointOfSegment const& se, + count_info const& ci, + PointOfSegment const& s1, PointOfSegment const& s2) const + { + 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 (ci.count > 0) + { + ss20 += small_angle(); + } + else + { + ss20 -= small_angle(); + } + math::normalize_longitude<units_t>(ss20); + set<0>(ss2, ss20); + + // Check the side using this vertical segment + return m_side_strategy.apply(ss1, ss2, point); + } + + // 1 deg or pi/180 rad + static inline calculation_type small_angle() + { + return constants::half_period() / calculation_type(180); + }; + + static inline bool longitudes_equal(calculation_type const& lon1, calculation_type const& lon2) + { + return math::equals( + math::longitude_distance_signed<units_t>(lon1, lon2), + calculation_type(0)); + } + + SideStrategy m_side_strategy; +}; + + +} // namespace detail +#endif // DOXYGEN_NO_DETAIL + + +/*! +\brief Within detection using winding rule in spherical coordinate system. +\ingroup strategies +\tparam Point \tparam_point +\tparam PointOfSegment \tparam_segment_point +\tparam CalculationType \tparam_calculation + +\qbk{ +[heading See also] +[link geometry.reference.algorithms.within.within_3_with_strategy within (with strategy)] +} + */ +template +< + typename Point, + typename PointOfSegment = Point, + typename CalculationType = void +> +class spherical_winding + : public within::detail::spherical_winding_base + < + Point, + PointOfSegment, + side::spherical_side_formula<CalculationType>, + CalculationType + > +{}; + + +#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, spherical_tag, spherical_tag> +{ + typedef within::detail::spherical_winding_base + < + 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::detail::spherical_winding_base + < + 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, spherical_tag, spherical_tag> +{ + typedef within::detail::spherical_winding_base + < + 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::detail::spherical_winding_base + < + typename geometry::point_type<PointLike>::type, + typename geometry::point_type<Geometry>::type + > type; +}; + +}}} // namespace strategy::covered_by::services +#endif + + +}} // namespace boost::geometry + + +#endif // BOOST_GEOMETRY_STRATEGY_SPHERICAL_POINT_IN_POLY_WINDING_HPP diff --git a/boost/geometry/strategies/spherical/side_by_cross_track.hpp b/boost/geometry/strategies/spherical/side_by_cross_track.hpp index 3f7be05557..2b195e96f0 100644 --- a/boost/geometry/strategies/spherical/side_by_cross_track.hpp +++ b/boost/geometry/strategies/spherical/side_by_cross_track.hpp @@ -2,9 +2,10 @@ // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2014. -// Modifications copyright (c) 2014, Oracle and/or its affiliates. +// This file was modified by Oracle on 2014-2017. +// Modifications copyright (c) 2014-2017, Oracle and/or its affiliates. +// Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle // Use, modification and distribution is subject to the Boost Software License, @@ -18,7 +19,7 @@ #include <boost/geometry/core/access.hpp> #include <boost/geometry/core/radian_access.hpp> -#include <boost/geometry/algorithms/detail/course.hpp> +#include <boost/geometry/formulas/spherical.hpp> #include <boost/geometry/util/math.hpp> #include <boost/geometry/util/promote_floating_point.hpp> @@ -59,8 +60,20 @@ public : >::type calc_t; calc_t d1 = 0.001; // m_strategy.apply(sp1, p); - calc_t crs_AD = geometry::detail::course<calc_t>(p1, p); - calc_t crs_AB = geometry::detail::course<calc_t>(p1, p2); + + calc_t lon1 = geometry::get_as_radian<0>(p1); + calc_t lat1 = geometry::get_as_radian<1>(p1); + calc_t lon2 = geometry::get_as_radian<0>(p2); + calc_t lat2 = geometry::get_as_radian<1>(p2); + calc_t lon = geometry::get_as_radian<0>(p); + calc_t lat = geometry::get_as_radian<1>(p); + + calc_t crs_AD = geometry::formula::spherical_azimuth<calc_t, false> + (lon1, lat1, lon, lat).azimuth; + + calc_t crs_AB = geometry::formula::spherical_azimuth<calc_t, false> + (lon1, lat1, lon2, lat2).azimuth; + calc_t XTD = asin(sin(d1) * sin(crs_AD - crs_AB)); return math::equals(XTD, 0) ? 0 : XTD < 0 ? 1 : -1; diff --git a/boost/geometry/strategies/strategies.hpp b/boost/geometry/strategies/strategies.hpp index 27a025c7c4..e7f3604abf 100644 --- a/boost/geometry/strategies/strategies.hpp +++ b/boost/geometry/strategies/strategies.hpp @@ -65,6 +65,7 @@ #include <boost/geometry/strategies/cartesian/point_in_box.hpp> #include <boost/geometry/strategies/cartesian/point_in_poly_franklin.hpp> #include <boost/geometry/strategies/cartesian/point_in_poly_crossings_multiply.hpp> +#include <boost/geometry/strategies/cartesian/point_in_poly_winding.hpp> #include <boost/geometry/strategies/cartesian/side_by_triangle.hpp> #include <boost/geometry/strategies/spherical/area.hpp> @@ -73,9 +74,10 @@ #include <boost/geometry/strategies/spherical/distance_haversine.hpp> #include <boost/geometry/strategies/spherical/distance_cross_track.hpp> #include <boost/geometry/strategies/spherical/distance_cross_track_point_box.hpp> -#include <boost/geometry/strategies/spherical/compare_circular.hpp> +#include <boost/geometry/strategies/spherical/compare.hpp> #include <boost/geometry/strategies/spherical/envelope_segment.hpp> #include <boost/geometry/strategies/spherical/intersection.hpp> +#include <boost/geometry/strategies/spherical/point_in_poly_winding.hpp> #include <boost/geometry/strategies/spherical/ssf.hpp> #include <boost/geometry/strategies/geographic/area.hpp> @@ -83,11 +85,13 @@ #include <boost/geometry/strategies/geographic/disjoint_segment_box.hpp> #include <boost/geometry/strategies/geographic/distance.hpp> #include <boost/geometry/strategies/geographic/distance_andoyer.hpp> +#include <boost/geometry/strategies/geographic/distance_cross_track.hpp> #include <boost/geometry/strategies/geographic/distance_thomas.hpp> #include <boost/geometry/strategies/geographic/distance_vincenty.hpp> #include <boost/geometry/strategies/geographic/envelope_segment.hpp> #include <boost/geometry/strategies/geographic/intersection.hpp> //#include <boost/geometry/strategies/geographic/intersection_elliptic.hpp> +#include <boost/geometry/strategies/geographic/point_in_poly_winding.hpp> #include <boost/geometry/strategies/geographic/side.hpp> #include <boost/geometry/strategies/geographic/side_andoyer.hpp> #include <boost/geometry/strategies/geographic/side_thomas.hpp> |