geometry/formulas/gnomonic_intersection.hpp

0001 // Boost.Geometry
0002
0003 // Copyright (c) 2023 Adam Wulkiewicz, Lodz, Poland.
0004
0005 // Copyright (c) 2016 Oracle and/or its affiliates.
0006 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0007
0008 // Use, modification and distribution is subject to the Boost Software License,
0009 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0010 // http://www.boost.org/LICENSE_1_0.txt)
0011
0012 #ifndef BOOST_GEOMETRY_FORMULAS_GNOMONIC_INTERSECTION_HPP
0013 #define BOOST_GEOMETRY_FORMULAS_GNOMONIC_INTERSECTION_HPP
0014
0015 #include <boost/geometry/core/access.hpp>
0016 #include <boost/geometry/core/cs.hpp>
0017
0018 #include <boost/geometry/arithmetic/cross_product.hpp>
0019 #include <boost/geometry/formulas/gnomonic_spheroid.hpp>
0020 #include <boost/geometry/geometries/point.hpp>
0021 #include <boost/geometry/util/math.hpp>
0022
0023
0024 namespace boost { namespace geometry { namespace formula
0025 {
0026
0027 /*!
0028 \brief The intersection of two geodesics using spheroidal gnomonic projection
0029        as proposed by Karney.
0030 \author See
0031     - Charles F.F Karney, Algorithms for geodesics, 2011
0032       https://arxiv.org/pdf/1109.4448.pdf
0033     - GeographicLib forum thread: Intersection between two geodesic lines
0034       https://sourceforge.net/p/geographiclib/discussion/1026621/thread/21aaff9f/
0035 */
0036 template
0037 <
0038     typename CT,
0039     template <typename, bool, bool, bool, bool, bool> class Inverse,
0040     template <typename, bool, bool, bool, bool> class Direct
0041 >
0042 class gnomonic_intersection
0043 {
0044 public:
0045     template <typename T1, typename T2, typename Spheroid>
0046     static inline bool apply(T1 const& lona1, T1 const& lata1,
0047                              T1 const& lona2, T1 const& lata2,
0048                              T2 const& lonb1, T2 const& latb1,
0049                              T2 const& lonb2, T2 const& latb2,
0050                              CT & lon, CT & lat,
0051                              Spheroid const& spheroid)
0052     {
0053         CT const lon_a1 = lona1;
0054         CT const lat_a1 = lata1;
0055         CT const lon_a2 = lona2;
0056         CT const lat_a2 = lata2;
0057         CT const lon_b1 = lonb1;
0058         CT const lat_b1 = latb1;
0059         CT const lon_b2 = lonb2;
0060         CT const lat_b2 = latb2;
0061
0062         return apply(lon_a1, lat_a1, lon_a2, lat_a2, lon_b1, lat_b1, lon_b2, lat_b2, lon, lat, spheroid);
0063     }
0064
0065     template <typename Spheroid>
0066     static inline bool apply(CT const& lona1, CT const& lata1,
0067                              CT const& lona2, CT const& lata2,
0068                              CT const& lonb1, CT const& latb1,
0069                              CT const& lonb2, CT const& latb2,
0070                              CT & lon, CT & lat,
0071                              Spheroid const& spheroid)
0072     {
0073         typedef gnomonic_spheroid<CT, Inverse, Direct> gnom_t;
0074
0075         lon = (lona1 + lona2 + lonb1 + lonb2) / 4;
0076         lat = (lata1 + lata2 + latb1 + latb2) / 4;
0077         // TODO: consider normalizing lon
0078
0079         for (int i = 0; i < 10; ++i)
0080         {
0081             CT xa1, ya1, xa2, ya2;
0082             CT xb1, yb1, xb2, yb2;
0083             CT x, y;
0084             CT lon1, lat1;
0085
0086             bool ok = gnom_t::forward(lon, lat, lona1, lata1, xa1, ya1, spheroid)
0087                    && gnom_t::forward(lon, lat, lona2, lata2, xa2, ya2, spheroid)
0088                    && gnom_t::forward(lon, lat, lonb1, latb1, xb1, yb1, spheroid)
0089                    && gnom_t::forward(lon, lat, lonb2, latb2, xb2, yb2, spheroid)
0090                    && intersect(xa1, ya1, xa2, ya2, xb1, yb1, xb2, yb2, x, y)
0091                    && gnom_t::inverse(lon, lat, x, y, lon1, lat1, spheroid);
0092
0093             if (! ok)
0094             {
0095                 return false;
0096             }
0097
0098             if (math::equals(lat1, lat) && math::equals(lon1, lon))
0099             {
0100                 break;
0101             }
0102
0103             lat = lat1;
0104             lon = lon1;
0105         }
0106
0107         // NOTE: true is also returned if the number of iterations is too great
0108         //       which means that the accuracy of the result is low
0109         return true;
0110     }
0111
0112 private:
0113     static inline bool intersect(CT const& xa1, CT const& ya1, CT const& xa2, CT const& ya2,
0114                                  CT const& xb1, CT const& yb1, CT const& xb2, CT const& yb2,
0115                                  CT & x, CT & y)
0116     {
0117         typedef model::point<CT, 3, cs::cartesian> v3d_t;
0118
0119         CT const c0 = 0;
0120         CT const c1 = 1;
0121
0122         v3d_t const va1(xa1, ya1, c1);
0123         v3d_t const va2(xa2, ya2, c1);
0124         v3d_t const vb1(xb1, yb1, c1);
0125         v3d_t const vb2(xb2, yb2, c1);
0126
0127         v3d_t const la = cross_product(va1, va2);
0128         v3d_t const lb = cross_product(vb1, vb2);
0129         v3d_t const p = cross_product(la, lb);
0130
0131         CT const z = get<2>(p);
0132
0133         if (math::equals(z, c0))
0134         {
0135             // degenerated or collinear segments
0136             return false;
0137         }
0138
0139         x = get<0>(p) / z;
0140         y = get<1>(p) / z;
0141
0142         return true;
0143     }
0144 };
0145
0146 }}} // namespace boost::geometry::formula
0147
0148
0149 #endif // BOOST_GEOMETRY_FORMULAS_GNOMONIC_INTERSECTION_HPP