algorithms/detail/direction_code.hpp

0001 // Boost.Geometry (aka GGL, Generic Geometry Library)
0002
0003 // Copyright (c) 2015-2020 Barend Gehrels, Amsterdam, the Netherlands.
0004
0005 // This file was modified by Oracle on 2015-2020.
0006 // Modifications copyright (c) 2015-2020 Oracle and/or its affiliates.
0007
0008 // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
0009 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0010
0011 // Use, modification and distribution is subject to the Boost Software License,
0012 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0013 // http://www.boost.org/LICENSE_1_0.txt)
0014
0015 #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_DIRECTION_CODE_HPP
0016 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_DIRECTION_CODE_HPP
0017
0018
0019 #include <type_traits>
0020
0021 #include <boost/geometry/core/access.hpp>
0022 #include <boost/geometry/core/static_assert.hpp>
0023 #include <boost/geometry/arithmetic/infinite_line_functions.hpp>
0024 #include <boost/geometry/algorithms/detail/make/make.hpp>
0025 #include <boost/geometry/util/math.hpp>
0026 #include <boost/geometry/util/select_coordinate_type.hpp>
0027 #include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>
0028
0029
0030 namespace boost { namespace geometry
0031 {
0032
0033
0034 #ifndef DOXYGEN_NO_DETAIL
0035 namespace detail
0036 {
0037
0038 template <typename CSTag>
0039 struct direction_code_impl
0040 {
0041     BOOST_GEOMETRY_STATIC_ASSERT_FALSE(
0042         "Not implemented for this coordinate system.",
0043         CSTag);
0044 };
0045
0046 template <>
0047 struct direction_code_impl<cartesian_tag>
0048 {
0049     template <typename PointSegmentA, typename PointSegmentB, typename Point2>
0050     static inline int apply(PointSegmentA const& segment_a, PointSegmentB const& segment_b,
0051                             Point2 const& point)
0052     {
0053         using calc_t = typename geometry::select_coordinate_type
0054             <
0055                 PointSegmentA, PointSegmentB, Point2
0056             >::type;
0057
0058         using line_type = model::infinite_line<calc_t>;
0059
0060         // Situation and construction of perpendicular line
0061         //
0062         //     P1     a--------------->b   P2
0063         //                             |
0064         //                             |
0065         //                             v
0066         //
0067         // P1 is located right of the (directional) perpendicular line
0068         // and therefore gets a negative side_value, and returns -1.
0069         // P2 is to the left of the perpendicular line and returns 1.
0070         // If the specified point is located on top of b, it returns 0.
0071
0072         line_type const line
0073             = detail::make::make_perpendicular_line<calc_t>(segment_a,
0074                 segment_b, segment_b);
0075
0076         if (arithmetic::is_degenerate(line))
0077         {
0078             return 0;
0079         }
0080
0081         calc_t const sv = arithmetic::side_value(line, point);
0082         static calc_t const zero = 0;
0083         return sv == zero ? 0 : sv > zero ? 1 : -1;
0084     }
0085 };
0086
0087 template <>
0088 struct direction_code_impl<spherical_equatorial_tag>
0089 {
0090     template <typename PointSegmentA, typename PointSegmentB, typename Point2>
0091     static inline int apply(PointSegmentA const& segment_a, PointSegmentB const& segment_b,
0092                             Point2 const& p)
0093     {
0094         {
0095             using units_sa_t =  typename cs_angular_units<PointSegmentA>::type;
0096             using units_sb_t =  typename cs_angular_units<PointSegmentB>::type;
0097             using units_p_t = typename cs_angular_units<Point2>::type;
0098             BOOST_GEOMETRY_STATIC_ASSERT(
0099                 (std::is_same<units_sa_t, units_sb_t>::value),
0100                 "Not implemented for different units.",
0101                 units_sa_t, units_sb_t);
0102             BOOST_GEOMETRY_STATIC_ASSERT(
0103                 (std::is_same<units_sa_t, units_p_t>::value),
0104                 "Not implemented for different units.",
0105                 units_sa_t, units_p_t);
0106         }
0107
0108         using coor_sa_t = typename coordinate_type<PointSegmentA>::type;
0109         using coor_sb_t = typename coordinate_type<PointSegmentB>::type;
0110         using coor_p_t = typename coordinate_type<Point2>::type;
0111
0112         // Declare unit type (equal for all types) and calc type (coerced to most precise)
0113         using units_t = typename cs_angular_units<Point2>::type;
0114         using calc_t = typename geometry::select_coordinate_type
0115             <
0116                 PointSegmentA, PointSegmentB, Point2
0117             >::type;
0118         using constants_sa_t = math::detail::constants_on_spheroid<coor_sa_t, units_t>;
0119         using constants_sb_t = math::detail::constants_on_spheroid<coor_sb_t, units_t>;
0120         using constants_p_t = math::detail::constants_on_spheroid<coor_p_t, units_t>;
0121
0122         static coor_sa_t const pi_half_sa = constants_sa_t::max_latitude();
0123         static coor_sb_t const pi_half_sb = constants_sb_t::max_latitude();
0124         static coor_p_t const pi_half_p = constants_p_t::max_latitude();
0125         static calc_t const c0 = 0;
0126
0127         coor_sa_t const a0 = geometry::get<0>(segment_a);
0128         coor_sa_t const a1 = geometry::get<1>(segment_a);
0129         coor_sb_t const b0 = geometry::get<0>(segment_b);
0130         coor_sb_t const b1 = geometry::get<1>(segment_b);
0131         coor_p_t const p0 = geometry::get<0>(p);
0132         coor_p_t const p1 = geometry::get<1>(p);
0133
0134         if ( (math::equals(b0, a0) && math::equals(b1, a1))
0135           || (math::equals(b0, p0) && math::equals(b1, p1)) )
0136         {
0137             return 0;
0138         }
0139
0140         bool const is_a_pole = math::equals(pi_half_sa, math::abs(a1));
0141         bool const is_b_pole = math::equals(pi_half_sb, math::abs(b1));
0142         bool const is_p_pole = math::equals(pi_half_p, math::abs(p1));
0143
0144         if ( is_b_pole && ((is_a_pole && math::sign(b1) == math::sign(a1))
0145                         || (is_p_pole && math::sign(b1) == math::sign(p1))) )
0146         {
0147             return 0;
0148         }
0149
0150         // NOTE: as opposed to the implementation for cartesian CS
0151         // here point b is the origin
0152
0153         calc_t const dlon1 = math::longitude_distance_signed<units_t, calc_t>(b0, a0);
0154         calc_t const dlon2 = math::longitude_distance_signed<units_t, calc_t>(b0, p0);
0155
0156         bool is_antilon1 = false, is_antilon2 = false;
0157         calc_t const dlat1 = latitude_distance_signed<units_t, calc_t>(b1, a1, dlon1, is_antilon1);
0158         calc_t const dlat2 = latitude_distance_signed<units_t, calc_t>(b1, p1, dlon2, is_antilon2);
0159
0160         calc_t const mx = is_a_pole || is_b_pole || is_p_pole
0161             ? c0
0162             : (std::min)(is_antilon1 ? c0 : math::abs(dlon1),
0163                          is_antilon2 ? c0 : math::abs(dlon2));
0164         calc_t const my = (std::min)(math::abs(dlat1),
0165                                      math::abs(dlat2));
0166
0167         int s1 = 0, s2 = 0;
0168         if (mx >= my)
0169         {
0170             s1 = dlon1 > 0 ? 1 : -1;
0171             s2 = dlon2 > 0 ? 1 : -1;
0172         }
0173         else
0174         {
0175             s1 = dlat1 > 0 ? 1 : -1;
0176             s2 = dlat2 > 0 ? 1 : -1;
0177         }
0178
0179         return s1 == s2 ? -1 : 1;
0180     }
0181
0182     template <typename Units, typename T>
0183     static inline T latitude_distance_signed(T const& lat1, T const& lat2, T const& lon_ds, bool & is_antilon)
0184     {
0185         using constants = math::detail::constants_on_spheroid<T, Units>;
0186         static T const pi = constants::half_period();
0187         static T const c0 = 0;
0188
0189         T res = lat2 - lat1;
0190
0191         is_antilon = math::equals(math::abs(lon_ds), pi);
0192         if (is_antilon)
0193         {
0194             res = lat2 + lat1;
0195             if (res >= c0)
0196                 res = pi - res;
0197             else
0198                 res = -pi - res;
0199         }
0200
0201         return res;
0202     }
0203 };
0204
0205 template <>
0206 struct direction_code_impl<spherical_polar_tag>
0207 {
0208     template <typename PointSegmentA, typename PointSegmentB, typename Point2>
0209     static inline int apply(PointSegmentA segment_a, PointSegmentB segment_b,
0210                             Point2 p)
0211     {
0212         using constants_sa_t = math::detail::constants_on_spheroid
0213             <
0214                 typename coordinate_type<PointSegmentA>::type,
0215                 typename cs_angular_units<PointSegmentA>::type
0216             >;
0217         using constants_p_t = math::detail::constants_on_spheroid
0218             <
0219                 typename coordinate_type<Point2>::type,
0220                 typename cs_angular_units<Point2>::type
0221             >;
0222
0223         geometry::set<1>(segment_a,
0224             constants_sa_t::max_latitude() - geometry::get<1>(segment_a));
0225         geometry::set<1>(segment_b,
0226             constants_sa_t::max_latitude() - geometry::get<1>(segment_b));
0227         geometry::set<1>(p,
0228             constants_p_t::max_latitude() - geometry::get<1>(p));
0229
0230         return direction_code_impl
0231                 <
0232                     spherical_equatorial_tag
0233                 >::apply(segment_a, segment_b, p);
0234     }
0235 };
0236
0237 // if spherical_tag is passed then pick cs_tag based on PointSegmentA type
0238 // with spherical_equatorial_tag as the default
0239 template <>
0240 struct direction_code_impl<spherical_tag>
0241 {
0242     template <typename PointSegmentA, typename PointSegmentB, typename Point2>
0243     static inline int apply(PointSegmentA segment_a, PointSegmentB segment_b,
0244                             Point2 p)
0245     {
0246         return direction_code_impl
0247             <
0248                 std::conditional_t
0249                     <
0250                         std::is_same
0251                             <
0252                                 typename geometry::cs_tag<PointSegmentA>::type,
0253                                 spherical_polar_tag
0254                             >::value,
0255                         spherical_polar_tag,
0256                         spherical_equatorial_tag
0257                     >
0258             >::apply(segment_a, segment_b, p);
0259     }
0260 };
0261
0262 template <>
0263 struct direction_code_impl<geographic_tag>
0264     : direction_code_impl<spherical_equatorial_tag>
0265 {};
0266
0267 // Gives sense of direction for point p, collinear w.r.t. segment (a,b)
0268 // Returns -1 if p goes backward w.r.t (a,b), so goes from b in direction of a
0269 // Returns 1 if p goes forward, so extends (a,b)
0270 // Returns 0 if p is equal with b, or if (a,b) is degenerate
0271 // Note that it does not do any collinearity test, that should be done before
0272 // In some cases the "segment" consists of different source points, and therefore
0273 // their types might differ.
0274 template <typename CSTag, typename PointSegmentA, typename PointSegmentB, typename Point2>
0275 inline int direction_code(PointSegmentA const& segment_a, PointSegmentB const& segment_b,
0276                           Point2 const& p)
0277 {
0278     return direction_code_impl<CSTag>::apply(segment_a, segment_b, p);
0279 }
0280
0281
0282 } // namespace detail
0283 #endif //DOXYGEN_NO_DETAIL
0284
0285
0286 }} // namespace boost::geometry
0287
0288 #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_DIRECTION_CODE_HPP