summaryrefslogtreecommitdiff
path: root/boost/geometry/strategies
diff options
context:
space:
mode:
authorDongHun Kwak <dh0128.kwak@samsung.com>2019-12-05 15:11:01 +0900
committerDongHun Kwak <dh0128.kwak@samsung.com>2019-12-05 15:11:01 +0900
commit3fdc3e5ee96dca5b11d1694975a65200787eab86 (patch)
tree5c1733853892b8397d67706fa453a9bd978d2102 /boost/geometry/strategies
parent88e602c57797660ebe0f9e15dbd64c1ff16dead3 (diff)
downloadboost-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')
-rw-r--r--boost/geometry/strategies/agnostic/buffer_distance_symmetric.hpp2
-rw-r--r--boost/geometry/strategies/agnostic/point_in_point.hpp2
-rw-r--r--boost/geometry/strategies/agnostic/point_in_poly_winding.hpp549
-rw-r--r--boost/geometry/strategies/cartesian/intersection.hpp5
-rw-r--r--boost/geometry/strategies/cartesian/point_in_poly_winding.hpp296
-rw-r--r--boost/geometry/strategies/cartesian/side_by_triangle.hpp16
-rw-r--r--boost/geometry/strategies/compare.hpp248
-rw-r--r--boost/geometry/strategies/geographic/area.hpp36
-rw-r--r--boost/geometry/strategies/geographic/distance.hpp41
-rw-r--r--boost/geometry/strategies/geographic/distance_cross_track.hpp641
-rw-r--r--boost/geometry/strategies/geographic/intersection.hpp21
-rw-r--r--boost/geometry/strategies/geographic/parameters.hpp51
-rw-r--r--boost/geometry/strategies/geographic/point_in_poly_winding.hpp80
-rw-r--r--boost/geometry/strategies/spherical/area.hpp11
-rw-r--r--boost/geometry/strategies/spherical/compare.hpp321
-rw-r--r--boost/geometry/strategies/spherical/compare_circular.hpp152
-rw-r--r--boost/geometry/strategies/spherical/distance_cross_track.hpp61
-rw-r--r--boost/geometry/strategies/spherical/intersection.hpp5
-rw-r--r--boost/geometry/strategies/spherical/point_in_poly_winding.hpp581
-rw-r--r--boost/geometry/strategies/spherical/side_by_cross_track.hpp23
-rw-r--r--boost/geometry/strategies/strategies.hpp6
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>