1 // Boost.Geometry 2 3 // Copyright (c) 2016 Oracle and/or its affiliates. 4 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_GNOMONIC_INTERSECTION_HPP 12 #define BOOST_GEOMETRY_FORMULAS_GNOMONIC_INTERSECTION_HPP 13 14 #include <boost/geometry/core/access.hpp> 15 #include <boost/geometry/core/cs.hpp> 16 17 #include <boost/geometry/arithmetic/cross_product.hpp> 18 #include <boost/geometry/formulas/gnomonic_spheroid.hpp> 19 #include <boost/geometry/geometries/point.hpp> 20 #include <boost/geometry/util/math.hpp> 21 22 23 namespace boost { namespace geometry { namespace formula 24 { 25 26 /*! 27 \brief The intersection of two geodesics using spheroidal gnomonic projection 28 as proposed by Karney. 29 \author See 30 - Charles F.F Karney, Algorithms for geodesics, 2011 31 https://arxiv.org/pdf/1109.4448.pdf 32 - GeographicLib forum thread: Intersection between two geodesic lines 33 https://sourceforge.net/p/geographiclib/discussion/1026621/thread/21aaff9f/ 34 */ 35 template 36 < 37 typename CT, 38 template <typename, bool, bool, bool, bool, bool> class Inverse, 39 template <typename, bool, bool, bool, bool> class Direct 40 > 41 class gnomonic_intersection 42 { 43 public: 44 template <typename T1, typename T2, typename Spheroid> apply(T1 const & lona1,T1 const & lata1,T1 const & lona2,T1 const & lata2,T2 const & lonb1,T2 const & latb1,T2 const & lonb2,T2 const & latb2,CT & lon,CT & lat,Spheroid const & spheroid)45 static inline bool apply(T1 const& lona1, T1 const& lata1, 46 T1 const& lona2, T1 const& lata2, 47 T2 const& lonb1, T2 const& latb1, 48 T2 const& lonb2, T2 const& latb2, 49 CT & lon, CT & lat, 50 Spheroid const& spheroid) 51 { 52 CT const lon_a1 = lona1; 53 CT const lat_a1 = lata1; 54 CT const lon_a2 = lona2; 55 CT const lat_a2 = lata2; 56 CT const lon_b1 = lonb1; 57 CT const lat_b1 = latb1; 58 CT const lon_b2 = lonb2; 59 CT const lat_b2 = latb2; 60 61 return apply(lon_a1, lat_a1, lon_a2, lat_a2, lon_b1, lat_b1, lon_b2, lat_b2, lon, lat, spheroid); 62 } 63 64 template <typename Spheroid> apply(CT const & lona1,CT const & lata1,CT const & lona2,CT const & lata2,CT const & lonb1,CT const & latb1,CT const & lonb2,CT const & latb2,CT & lon,CT & lat,Spheroid const & spheroid)65 static inline bool apply(CT const& lona1, CT const& lata1, 66 CT const& lona2, CT const& lata2, 67 CT const& lonb1, CT const& latb1, 68 CT const& lonb2, CT const& latb2, 69 CT & lon, CT & lat, 70 Spheroid const& spheroid) 71 { 72 typedef gnomonic_spheroid<CT, Inverse, Direct> gnom_t; 73 74 lon = (lona1 + lona2 + lonb1 + lonb2) / 4; 75 lat = (lata1 + lata2 + latb1 + latb2) / 4; 76 // TODO: consider normalizing lon 77 78 for (int i = 0; i < 10; ++i) 79 { 80 CT xa1, ya1, xa2, ya2; 81 CT xb1, yb1, xb2, yb2; 82 CT x, y; 83 double lat1, lon1; 84 85 bool ok = gnom_t::forward(lon, lat, lona1, lata1, xa1, ya1, spheroid) 86 && gnom_t::forward(lon, lat, lona2, lata2, xa2, ya2, spheroid) 87 && gnom_t::forward(lon, lat, lonb1, latb1, xb1, yb1, spheroid) 88 && gnom_t::forward(lon, lat, lonb2, latb2, xb2, yb2, spheroid) 89 && intersect(xa1, ya1, xa2, ya2, xb1, yb1, xb2, yb2, x, y) 90 && gnom_t::inverse(lon, lat, x, y, lon1, lat1, spheroid); 91 92 if (! ok) 93 { 94 return false; 95 } 96 97 if (math::equals(lat1, lat) && math::equals(lon1, lon)) 98 { 99 break; 100 } 101 102 lat = lat1; 103 lon = lon1; 104 } 105 106 // NOTE: true is also returned if the number of iterations is too great 107 // which means that the accuracy of the result is low 108 return true; 109 } 110 111 private: intersect(CT const & xa1,CT const & ya1,CT const & xa2,CT const & ya2,CT const & xb1,CT const & yb1,CT const & xb2,CT const & yb2,CT & x,CT & y)112 static inline bool intersect(CT const& xa1, CT const& ya1, CT const& xa2, CT const& ya2, 113 CT const& xb1, CT const& yb1, CT const& xb2, CT const& yb2, 114 CT & x, CT & y) 115 { 116 typedef model::point<CT, 3, cs::cartesian> v3d_t; 117 118 CT const c0 = 0; 119 CT const c1 = 1; 120 121 v3d_t const va1(xa1, ya1, c1); 122 v3d_t const va2(xa2, ya2, c1); 123 v3d_t const vb1(xb1, yb1, c1); 124 v3d_t const vb2(xb2, yb2, c1); 125 126 v3d_t const la = cross_product(va1, va2); 127 v3d_t const lb = cross_product(vb1, vb2); 128 v3d_t const p = cross_product(la, lb); 129 130 CT const z = get<2>(p); 131 132 if (math::equals(z, c0)) 133 { 134 // degenerated or collinear segments 135 return false; 136 } 137 138 x = get<0>(p) / z; 139 y = get<1>(p) / z; 140 141 return true; 142 } 143 }; 144 145 }}} // namespace boost::geometry::formula 146 147 148 #endif // BOOST_GEOMETRY_FORMULAS_GNOMONIC_INTERSECTION_HPP 149