diff options
Diffstat (limited to 'boost/geometry/formulas/gnomonic_intersection.hpp')
-rw-r--r-- | boost/geometry/formulas/gnomonic_intersection.hpp | 148 |
1 files changed, 148 insertions, 0 deletions
diff --git a/boost/geometry/formulas/gnomonic_intersection.hpp b/boost/geometry/formulas/gnomonic_intersection.hpp new file mode 100644 index 0000000000..7e1b7bcdab --- /dev/null +++ b/boost/geometry/formulas/gnomonic_intersection.hpp @@ -0,0 +1,148 @@ +// Boost.Geometry + +// Copyright (c) 2016 Oracle and/or its affiliates. + +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_GEOMETRY_FORMULAS_GNOMONIC_INTERSECTION_HPP +#define BOOST_GEOMETRY_FORMULAS_GNOMONIC_INTERSECTION_HPP + +#include <boost/geometry/core/access.hpp> +#include <boost/geometry/core/cs.hpp> + +#include <boost/geometry/arithmetic/cross_product.hpp> +#include <boost/geometry/formulas/gnomonic_spheroid.hpp> +#include <boost/geometry/geometries/point.hpp> +#include <boost/geometry/util/math.hpp> + + +namespace boost { namespace geometry { namespace formula +{ + +/*! +\brief The intersection of two geodesics using spheroidal gnomonic projection + as proposed by Karney. +\author See + - Charles F.F Karney, Algorithms for geodesics, 2011 + https://arxiv.org/pdf/1109.4448.pdf + - GeographicLib forum thread: Intersection between two geodesic lines + https://sourceforge.net/p/geographiclib/discussion/1026621/thread/21aaff9f/ +*/ +template +< + typename CT, + template <typename, bool, bool, bool, bool, bool> class Inverse, + template <typename, bool, bool, bool, bool> class Direct +> +class gnomonic_intersection +{ +public: + template <typename T1, typename T2, typename Spheroid> + static inline bool 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) + { + CT const lon_a1 = lona1; + CT const lat_a1 = lata1; + CT const lon_a2 = lona2; + CT const lat_a2 = lata2; + CT const lon_b1 = lonb1; + CT const lat_b1 = latb1; + CT const lon_b2 = lonb2; + CT const lat_b2 = latb2; + + return apply(lon_a1, lat_a1, lon_a2, lat_a2, lon_b1, lat_b1, lon_b2, lat_b2, lon, lat, spheroid); + } + + template <typename Spheroid> + static inline bool 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) + { + typedef gnomonic_spheroid<CT, Inverse, Direct> gnom_t; + + lon = (lona1 + lona2 + lonb1 + lonb2) / 4; + lat = (lata1 + lata2 + latb1 + latb2) / 4; + // TODO: consider normalizing lon + + for (int i = 0; i < 10; ++i) + { + CT xa1, ya1, xa2, ya2; + CT xb1, yb1, xb2, yb2; + CT x, y; + double lat1, lon1; + + bool ok = gnom_t::forward(lon, lat, lona1, lata1, xa1, ya1, spheroid) + && gnom_t::forward(lon, lat, lona2, lata2, xa2, ya2, spheroid) + && gnom_t::forward(lon, lat, lonb1, latb1, xb1, yb1, spheroid) + && gnom_t::forward(lon, lat, lonb2, latb2, xb2, yb2, spheroid) + && intersect(xa1, ya1, xa2, ya2, xb1, yb1, xb2, yb2, x, y) + && gnom_t::inverse(lon, lat, x, y, lon1, lat1, spheroid); + + if (! ok) + { + return false; + } + + if (math::equals(lat1, lat) && math::equals(lon1, lon)) + { + break; + } + + lat = lat1; + lon = lon1; + } + + // NOTE: true is also returned if the number of iterations is too great + // which means that the accuracy of the result is low + return true; + } + +private: + static inline bool 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) + { + typedef model::point<CT, 3, cs::cartesian> v3d_t; + + CT const c0 = 0; + CT const c1 = 1; + + v3d_t const va1(xa1, ya1, c1); + v3d_t const va2(xa2, ya2, c1); + v3d_t const vb1(xb1, yb1, c1); + v3d_t const vb2(xb2, yb2, c1); + + v3d_t const la = cross_product(va1, va2); + v3d_t const lb = cross_product(vb1, vb2); + v3d_t const p = cross_product(la, lb); + + CT const z = get<2>(p); + + if (math::equals(z, c0)) + { + // degenerated or collinear segments + return false; + } + + x = get<0>(p) / z; + y = get<1>(p) / z; + + return true; + } +}; + +}}} // namespace boost::geometry::formula + + +#endif // BOOST_GEOMETRY_FORMULAS_GNOMONIC_INTERSECTION_HPP |