3 // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
5 // This file was modified by Oracle on 2014, 2015.
6 // Modifications copyright (c) 2014-2015 Oracle and/or its affiliates.
8 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
10 // Use, modification and distribution is subject to the Boost Software License,
11 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
12 // http://www.boost.org/LICENSE_1_0.txt)
14 #ifndef BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_SIDE_DETAIL_HPP
15 #define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_SIDE_DETAIL_HPP
17 #include <boost/geometry/core/cs.hpp>
18 #include <boost/geometry/core/access.hpp>
19 #include <boost/geometry/core/radian_access.hpp>
20 #include <boost/geometry/core/radius.hpp>
22 #include <boost/geometry/util/math.hpp>
23 #include <boost/geometry/util/promote_floating_point.hpp>
24 #include <boost/geometry/util/select_calculation_type.hpp>
26 #include <boost/geometry/strategies/side.hpp>
27 //#include <boost/geometry/strategies/concepts/side_concept.hpp>
30 namespace boost { namespace geometry
34 namespace strategy { namespace side
37 #ifndef DOXYGEN_NO_DETAIL
42 \brief Check at which side of a segment a point lies
43 left of segment (> 0), right of segment (< 0), on segment (0)
45 \tparam InverseFormula Geodesic inverse solution formula.
46 \tparam Model Reference model of coordinate system.
47 \tparam CalculationType \tparam_calculation
49 template <template<typename, bool, bool> class InverseFormula,
51 typename CalculationType = void>
55 by_azimuth(Model const& model = Model())
59 template <typename P1, typename P2, typename P>
60 inline int apply(P1 const& p1, P2 const& p2, P const& p)
62 typedef typename promote_floating_point
64 typename select_calculation_type_alt
71 typedef InverseFormula<calc_t, false, true> inverse_formula;
73 calc_t a1p = azimuth<calc_t, inverse_formula>(p1, p, m_model);
74 calc_t a12 = azimuth<calc_t, inverse_formula>(p1, p2, m_model);
76 calc_t const pi = math::pi<calc_t>();
78 // instead of the formula from XTD
79 //calc_t a_diff = asin(sin(a1p - a12));
81 calc_t a_diff = a1p - a12;
82 // normalize, angle in [-pi, pi]
84 a_diff -= calc_t(2) * pi;
85 while ( a_diff < -pi )
86 a_diff += calc_t(2) * pi;
88 // NOTE: in general it shouldn't be required to support the pi/-pi case
89 // because in non-cartesian systems it makes sense to check the side
90 // only "between" the endpoints.
91 // However currently the winding strategy calls the side strategy
92 // for vertical segments to check if the point is "between the endpoints.
93 // This could be avoided since the side strategy is not required for that
94 // because meridian is the shortest path. So a difference of
95 // longitudes would be sufficient (of course normalized to [-pi, pi]).
97 // NOTE: with the above said, the pi/-pi check is temporary
98 // however in case if this was required
99 // the geodesics on ellipsoid aren't "symmetrical"
100 // therefore instead of comparing a_diff to pi and -pi
101 // one should probably use inverse azimuths and compare
102 // the difference to 0 as well
104 // positive azimuth is on the right side
105 return math::equals(a_diff, 0)
106 || math::equals(a_diff, pi)
107 || math::equals(a_diff, -pi) ? 0
108 : a_diff > 0 ? -1 // right
113 template <typename ResultType,
114 typename InverseFormulaType,
118 static inline ResultType azimuth(Point1 const& point1, Point2 const& point2, ModelT const& model)
120 return InverseFormulaType::apply(get_as_radian<0>(point1),
121 get_as_radian<1>(point1),
122 get_as_radian<0>(point2),
123 get_as_radian<1>(point2),
131 #endif // DOXYGEN_NO_DETAIL
133 }} // namespace strategy::side
136 }} // namespace boost::geometry
139 #endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_SIDE_DETAIL_HPP