boost/geometry/algorithms/detail/direction_code.hpp

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