ceph/src/boost/boost/geometry/formulas/spherical.hpp

   1 // Boost.Geometry
   2
   3 // Copyright (c) 2016, Oracle and/or its affiliates.
   4 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
   5 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
   6
   7 // Use, modification and distribution is subject to the Boost Software License,
   8 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
   9 // http://www.boost.org/LICENSE_1_0.txt)
  10
  11 #ifndef BOOST_GEOMETRY_FORMULAS_SPHERICAL_HPP
  12 #define BOOST_GEOMETRY_FORMULAS_SPHERICAL_HPP
  13
  14 #include <boost/geometry/core/coordinate_system.hpp>
  15 #include <boost/geometry/core/coordinate_type.hpp>
  16 #include <boost/geometry/core/access.hpp>
  17 #include <boost/geometry/core/radian_access.hpp>
  18
  19 //#include <boost/geometry/arithmetic/arithmetic.hpp>
  20 #include <boost/geometry/arithmetic/cross_product.hpp>
  21 #include <boost/geometry/arithmetic/dot_product.hpp>
  22
  23 #include <boost/geometry/util/math.hpp>
  24 #include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>
  25 #include <boost/geometry/util/select_coordinate_type.hpp>
  26
  27 namespace boost { namespace geometry {
  28
  29 namespace formula {
  30
  31 template <typename T>
  32 struct result_spherical
  33 {
  34     result_spherical()
  35         : azimuth(0)
  36         , reverse_azimuth(0)
  37     {}
  38
  39     T azimuth;
  40     T reverse_azimuth;
  41 };
  42
  43 template <typename T>
  44 static inline void sph_to_cart3d(T const& lon, T const& lat, T & x, T & y, T & z)
  45 {
  46     T const cos_lat = cos(lat);
  47     x = cos_lat * cos(lon);
  48     y = cos_lat * sin(lon);
  49     z = sin(lat);
  50 }
  51
  52 template <typename Point3d, typename PointSph>
  53 static inline Point3d sph_to_cart3d(PointSph const& point_sph)
  54 {
  55     typedef typename coordinate_type<Point3d>::type calc_t;
  56
  57     calc_t const lon = get_as_radian<0>(point_sph);
  58     calc_t const lat = get_as_radian<1>(point_sph);
  59     calc_t x, y, z;
  60     sph_to_cart3d(lon, lat, x, y, z);
  61
  62     Point3d res;
  63     set<0>(res, x);
  64     set<1>(res, y);
  65     set<2>(res, z);
  66
  67     return res;
  68 }
  69
  70 template <typename T>
  71 static inline void cart3d_to_sph(T const& x, T const& y, T const& z, T & lon, T & lat)
  72 {
  73     lon = atan2(y, x);
  74     lat = asin(z);
  75 }
  76
  77 template <typename PointSph, typename Point3d>
  78 static inline PointSph cart3d_to_sph(Point3d const& point_3d)
  79 {
  80     typedef typename coordinate_type<PointSph>::type coord_t;
  81     typedef typename coordinate_type<Point3d>::type calc_t;
  82
  83     calc_t const x = get<0>(point_3d);
  84     calc_t const y = get<1>(point_3d);
  85     calc_t const z = get<2>(point_3d);
  86     calc_t lonr, latr;
  87     cart3d_to_sph(x, y, z, lonr, latr);
  88
  89     PointSph res;
  90     set_from_radian<0>(res, lonr);
  91     set_from_radian<1>(res, latr);
  92
  93     coord_t lon = get<0>(res);
  94     coord_t lat = get<1>(res);
  95
  96     math::normalize_spheroidal_coordinates
  97         <
  98             typename coordinate_system<PointSph>::type::units,
  99             coord_t
 100         >(lon, lat);
 101
 102     set<0>(res, lon);
 103     set<1>(res, lat);
 104
 105     return res;
 106 }
 107
 108 // -1 right
 109 // 1 left
 110 // 0 on
 111 template <typename Point3d1, typename Point3d2>
 112 static inline int sph_side_value(Point3d1 const& norm, Point3d2 const& pt)
 113 {
 114     typedef typename select_coordinate_type<Point3d1, Point3d2>::type calc_t;
 115     calc_t c0 = 0;
 116     calc_t d = dot_product(norm, pt);
 117     return math::equals(d, c0) ? 0
 118         : d > c0 ? 1
 119         : -1; // d < 0
 120 }
 121
 122 template <typename CT, bool ReverseAzimuth, typename T1, typename T2>
 123 static inline result_spherical<CT> spherical_azimuth(T1 const& lon1,
 124                                                      T1 const& lat1,
 125                                                      T2 const& lon2,
 126                                                      T2 const& lat2)
 127 {
 128     typedef result_spherical<CT> result_type;
 129     result_type result;
 130
 131     // http://williams.best.vwh.net/avform.htm#Crs
 132     // https://en.wikipedia.org/wiki/Great-circle_navigation
 133     CT dlon = lon2 - lon1;
 134
 135     // An optimization which should kick in often for Boxes
 136     //if ( math::equals(dlon, ReturnType(0)) )
 137     //if ( get<0>(p1) == get<0>(p2) )
 138     //{
 139     //    return - sin(get_as_radian<1>(p1)) * cos_p2lat);
 140     //}
 141
 142     CT const cos_dlon = cos(dlon);
 143     CT const sin_dlon = sin(dlon);
 144     CT const cos_lat1 = cos(lat1);
 145     CT const cos_lat2 = cos(lat2);
 146     CT const sin_lat1 = sin(lat1);
 147     CT const sin_lat2 = sin(lat2);
 148
 149     {
 150         // "An alternative formula, not requiring the pre-computation of d"
 151         // In the formula below dlon is used as "d"
 152         CT const y = sin_dlon * cos_lat2;
 153         CT const x = cos_lat1 * sin_lat2 - sin_lat1 * cos_lat2 * cos_dlon;
 154         result.azimuth = atan2(y, x);
 155     }
 156
 157     if (ReverseAzimuth)
 158     {
 159         CT const y = sin_dlon * cos_lat1;
 160         CT const x = sin_lat2 * cos_lat1 * cos_dlon - cos_lat2 * sin_lat1;
 161         result.reverse_azimuth = atan2(y, x);
 162     }
 163
 164     return result;
 165 }
 166
 167 template <typename ReturnType, typename T1, typename T2>
 168 inline ReturnType spherical_azimuth(T1 const& lon1, T1 const& lat1,
 169                                     T2 const& lon2, T2 const& lat2)
 170 {
 171     return spherical_azimuth<ReturnType, false>(lon1, lat1, lon2, lat2).azimuth;
 172 }
 173
 174 template <typename T>
 175 inline T spherical_azimuth(T const& lon1, T const& lat1, T const& lon2, T const& lat2)
 176 {
 177     return spherical_azimuth<T, false>(lon1, lat1, lon2, lat2).azimuth;
 178 }
 179
 180 template <typename T>
 181 inline int azimuth_side_value(T const& azi_a1_p, T const& azi_a1_a2)
 182 {
 183     T const pi = math::pi<T>();
 184     T const two_pi = math::two_pi<T>();
 185
 186     // instead of the formula from XTD
 187     //calc_t a_diff = asin(sin(azi_a1_p - azi_a1_a2));
 188
 189     T a_diff = azi_a1_p - azi_a1_a2;
 190     // normalize, angle in [-pi, pi]
 191     while (a_diff > pi)
 192         a_diff -= two_pi;
 193     while (a_diff < -pi)
 194         a_diff += two_pi;
 195
 196     // NOTE: in general it shouldn't be required to support the pi/-pi case
 197     // because in non-cartesian systems it makes sense to check the side
 198     // only "between" the endpoints.
 199     // However currently the winding strategy calls the side strategy
 200     // for vertical segments to check if the point is "between the endpoints.
 201     // This could be avoided since the side strategy is not required for that
 202     // because meridian is the shortest path. So a difference of
 203     // longitudes would be sufficient (of course normalized to [-pi, pi]).
 204
 205     // NOTE: with the above said, the pi/-pi check is temporary
 206     // however in case if this was required
 207     // the geodesics on ellipsoid aren't "symmetrical"
 208     // therefore instead of comparing a_diff to pi and -pi
 209     // one should probably use inverse azimuths and compare
 210     // the difference to 0 as well
 211
 212     // positive azimuth is on the right side
 213     return math::equals(a_diff, 0)
 214         || math::equals(a_diff, pi)
 215         || math::equals(a_diff, -pi) ? 0
 216         : a_diff > 0 ? -1 // right
 217         : 1; // left
 218 }
 219
 220 } // namespace formula
 221
 222 }} // namespace boost::geometry
 223
 224 #endif // BOOST_GEOMETRY_FORMULAS_SPHERICAL_HPP