ceph/src/boost/boost/geometry/algorithms/detail/direction_code.hpp

   1 // Boost.Geometry (aka GGL, Generic Geometry Library)
   2
   3 // Copyright (c) 2015-2020 Barend Gehrels, Amsterdam, the Netherlands.
   4
   5 // This file was modified by Oracle on 2015-2020.
   6 // Modifications copyright (c) 2015-2020 Oracle and/or its affiliates.
   7
   8 // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
   9 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
  10
  11 // Use, modification and distribution is subject to the Boost Software License,
  12 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  13 // http://www.boost.org/LICENSE_1_0.txt)
  14
  15 #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_DIRECTION_CODE_HPP
  16 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_DIRECTION_CODE_HPP
  17
  18
  19 #include <type_traits>
  20
  21 #include <boost/geometry/core/access.hpp>
  22 #include <boost/geometry/core/static_assert.hpp>
  23 #include <boost/geometry/arithmetic/infinite_line_functions.hpp>
  24 #include <boost/geometry/algorithms/detail/make/make.hpp>
  25 #include <boost/geometry/util/math.hpp>
  26 #include <boost/geometry/util/select_coordinate_type.hpp>
  27 #include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>
  28
  29
  30 namespace boost { namespace geometry
  31 {
  32
  33
  34 #ifndef DOXYGEN_NO_DETAIL
  35 namespace detail
  36 {
  37
  38 template <typename CSTag>
  39 struct direction_code_impl
  40 {
  41     BOOST_GEOMETRY_STATIC_ASSERT_FALSE(
  42         "Not implemented for this coordinate system.",
  43         CSTag);
  44 };
  45
  46 template <>
  47 struct direction_code_impl<cartesian_tag>
  48 {
  49     template <typename Point1, typename Point2>
  50     static inline int apply(Point1 const& segment_a, Point1 const& segment_b,
  51                             Point2 const& point)
  52     {
  53         typedef typename geometry::select_coordinate_type
  54             <
  55                 Point1, Point2
  56             >::type calc_t;
  57
  58         typedef model::infinite_line<calc_t> line_type;
  59
  60         // Situation and construction of perpendicular line
  61         //
  62         //     P1     a--------------->b   P2
  63         //                             |
  64         //                             |
  65         //                             v
  66         //
  67         // P1 is located right of the (directional) perpendicular line
  68         // and therefore gets a negative side_value, and returns -1.
  69         // P2 is to the left of the perpendicular line and returns 1.
  70         // If the specified point is located on top of b, it returns 0.
  71
  72         line_type const line
  73             = detail::make::make_perpendicular_line<calc_t>(segment_a,
  74                 segment_b, segment_b);
  75
  76         if (arithmetic::is_degenerate(line))
  77         {
  78             return 0;
  79         }
  80
  81         calc_t const sv = arithmetic::side_value(line, point);
  82         return sv == 0 ? 0 : sv > 0 ? 1 : -1;
  83     }
  84 };
  85
  86 template <>
  87 struct direction_code_impl<spherical_equatorial_tag>
  88 {
  89     template <typename Point1, typename Point2>
  90     static inline int apply(Point1 const& segment_a, Point1 const& segment_b,
  91                             Point2 const& p)
  92     {
  93         typedef typename coordinate_type<Point1>::type coord1_t;
  94         typedef typename coordinate_type<Point2>::type coord2_t;
  95         typedef typename cs_angular_units<Point1>::type units_t;
  96         typedef typename cs_angular_units<Point2>::type units2_t;
  97         BOOST_GEOMETRY_STATIC_ASSERT(
  98             (std::is_same<units_t, units2_t>::value),
  99             "Not implemented for different units.",
 100             units_t, units2_t);
 101
 102         typedef typename geometry::select_coordinate_type <Point1, Point2>::type calc_t;
 103         typedef math::detail::constants_on_spheroid<coord1_t, units_t> constants1;
 104         typedef math::detail::constants_on_spheroid<coord2_t, units_t> constants2;
 105         static coord1_t const pi_half1 = constants1::max_latitude();
 106         static coord2_t const pi_half2 = constants2::max_latitude();
 107         static calc_t const c0 = 0;
 108
 109         coord1_t const a0 = geometry::get<0>(segment_a);
 110         coord1_t const a1 = geometry::get<1>(segment_a);
 111         coord1_t const b0 = geometry::get<0>(segment_b);
 112         coord1_t const b1 = geometry::get<1>(segment_b);
 113         coord2_t const p0 = geometry::get<0>(p);
 114         coord2_t const p1 = geometry::get<1>(p);
 115
 116         if ( (math::equals(b0, a0) && math::equals(b1, a1))
 117           || (math::equals(b0, p0) && math::equals(b1, p1)) )
 118         {
 119             return 0;
 120         }
 121
 122         bool const is_a_pole = math::equals(pi_half1, math::abs(a1));
 123         bool const is_b_pole = math::equals(pi_half1, math::abs(b1));
 124         bool const is_p_pole = math::equals(pi_half2, math::abs(p1));
 125
 126         if ( is_b_pole && ((is_a_pole && math::sign(b1) == math::sign(a1))
 127                         || (is_p_pole && math::sign(b1) == math::sign(p1))) )
 128         {
 129             return 0;
 130         }
 131
 132         // NOTE: as opposed to the implementation for cartesian CS
 133         // here point b is the origin
 134
 135         calc_t const dlon1 = math::longitude_distance_signed<units_t, calc_t>(b0, a0);
 136         calc_t const dlon2 = math::longitude_distance_signed<units_t, calc_t>(b0, p0);
 137
 138         bool is_antilon1 = false, is_antilon2 = false;
 139         calc_t const dlat1 = latitude_distance_signed<units_t, calc_t>(b1, a1, dlon1, is_antilon1);
 140         calc_t const dlat2 = latitude_distance_signed<units_t, calc_t>(b1, p1, dlon2, is_antilon2);
 141
 142         calc_t mx = is_a_pole || is_b_pole || is_p_pole ?
 143                     c0 :
 144                     (std::min)(is_antilon1 ? c0 : math::abs(dlon1),
 145                                is_antilon2 ? c0 : math::abs(dlon2));
 146         calc_t my = (std::min)(math::abs(dlat1),
 147                                math::abs(dlat2));
 148
 149         int s1 = 0, s2 = 0;
 150         if (mx >= my)
 151         {
 152             s1 = dlon1 > 0 ? 1 : -1;
 153             s2 = dlon2 > 0 ? 1 : -1;
 154         }
 155         else
 156         {
 157             s1 = dlat1 > 0 ? 1 : -1;
 158             s2 = dlat2 > 0 ? 1 : -1;
 159         }
 160
 161         return s1 == s2 ? -1 : 1;
 162     }
 163
 164     template <typename Units, typename T>
 165     static inline T latitude_distance_signed(T const& lat1, T const& lat2, T const& lon_ds, bool & is_antilon)
 166     {
 167         typedef math::detail::constants_on_spheroid<T, Units> constants;
 168         static T const pi = constants::half_period();
 169         static T const c0 = 0;
 170
 171         T res = lat2 - lat1;
 172
 173         is_antilon = math::equals(math::abs(lon_ds), pi);
 174         if (is_antilon)
 175         {
 176             res = lat2 + lat1;
 177             if (res >= c0)
 178                 res = pi - res;
 179             else
 180                 res = -pi - res;
 181         }
 182
 183         return res;
 184     }
 185 };
 186
 187 template <>
 188 struct direction_code_impl<spherical_polar_tag>
 189 {
 190     template <typename Point1, typename Point2>
 191     static inline int apply(Point1 segment_a, Point1 segment_b,
 192                             Point2 p)
 193     {
 194         typedef math::detail::constants_on_spheroid
 195             <
 196                 typename coordinate_type<Point1>::type,
 197                 typename cs_angular_units<Point1>::type
 198             > constants1;
 199         typedef math::detail::constants_on_spheroid
 200             <
 201                 typename coordinate_type<Point2>::type,
 202                 typename cs_angular_units<Point2>::type
 203             > constants2;
 204
 205         geometry::set<1>(segment_a,
 206             constants1::max_latitude() - geometry::get<1>(segment_a));
 207         geometry::set<1>(segment_b,
 208             constants1::max_latitude() - geometry::get<1>(segment_b));
 209         geometry::set<1>(p,
 210             constants2::max_latitude() - geometry::get<1>(p));
 211
 212         return direction_code_impl
 213                 <
 214                     spherical_equatorial_tag
 215                 >::apply(segment_a, segment_b, p);
 216     }
 217 };
 218
 219 // if spherical_tag is passed then pick cs_tag based on Point1 type
 220 // with spherical_equatorial_tag as the default
 221 template <>
 222 struct direction_code_impl<spherical_tag>
 223 {
 224     template <typename Point1, typename Point2>
 225     static inline int apply(Point1 segment_a, Point1 segment_b,
 226                             Point2 p)
 227     {
 228         return direction_code_impl
 229             <
 230                 std::conditional_t
 231                     <
 232                         std::is_same
 233                             <
 234                                 typename geometry::cs_tag<Point1>::type,
 235                                 spherical_polar_tag
 236                             >::value,
 237                         spherical_polar_tag,
 238                         spherical_equatorial_tag
 239                     >
 240             >::apply(segment_a, segment_b, p);
 241     }
 242 };
 243
 244 template <>
 245 struct direction_code_impl<geographic_tag>
 246     : direction_code_impl<spherical_equatorial_tag>
 247 {};
 248
 249 // Gives sense of direction for point p, collinear w.r.t. segment (a,b)
 250 // Returns -1 if p goes backward w.r.t (a,b), so goes from b in direction of a
 251 // Returns 1 if p goes forward, so extends (a,b)
 252 // Returns 0 if p is equal with b, or if (a,b) is degenerate
 253 // Note that it does not do any collinearity test, that should be done before
 254 template <typename CSTag, typename Point1, typename Point2>
 255 inline int direction_code(Point1 const& segment_a, Point1 const& segment_b,
 256                           Point2 const& p)
 257 {
 258     return direction_code_impl<CSTag>::apply(segment_a, segment_b, p);
 259 }
 260
 261
 262 } // namespace detail
 263 #endif //DOXYGEN_NO_DETAIL
 264
 265
 266 }} // namespace boost::geometry
 267
 268 #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_DIRECTION_CODE_HPP