[formulas] Add ellipsoidal gnomonic projection.

The projection was proposed by C.F.F. Karney in "Algorithms for geodesics". On ellipsoid of revolution (spheroid). The projection takes geodesic inverse and direct formulas used internally as template parameters.
2026-02-01 20:42:10 +00:00 · 2016-07-14 20:56:07 +02:00
parent a31669b958
commit 1d8938d53f
1 changed files with 120 additions and 0 deletions
--- a/include/boost/geometry/formulas/gnomonic_spheroid.hpp
+++ b/include/boost/geometry/formulas/gnomonic_spheroid.hpp
@@ -0,0 +1,120 @@
+// 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_SPHEROID_HPP
+#define BOOST_GEOMETRY_FORMULAS_GNOMONIC_SPHEROID_HPP
+
+
+#include <boost/geometry/core/radius.hpp>
+
+#include <boost/geometry/algorithms/detail/flattening.hpp>
+
+#include <boost/geometry/util/condition.hpp>
+#include <boost/geometry/util/math.hpp>
+
+#include <boost/geometry/formulas/andoyer_inverse.hpp>
+#include <boost/geometry/formulas/thomas_inverse.hpp>
+#include <boost/geometry/formulas/vincenty_direct.hpp>
+#include <boost/geometry/formulas/vincenty_inverse.hpp>
+
+
+namespace boost { namespace geometry { namespace formula
+{
+
+/*!
+\brief Gnomonic projection on spheroid (ellipsoid of revolution).
+\author See
+- Charles F.F Karney, Algorithms for geodesics, 2011
+https://arxiv.org/pdf/1109.4448.pdf
+*/
+template <
+    typename CT,
+    template <typename, bool, bool, bool, bool ,bool> class Inverse,
+    template <typename, bool, bool, bool, bool> class Direct
+>
+class gnomonic_spheroid
+{
+    typedef Inverse<CT, false, true, true, true, true> inverse_type;
+    typedef typename inverse_type::result_type inverse_result;
+
+    typedef Direct<CT, false, false, true, true> direct_quantities_type;
+    typedef Direct<CT, true, false, false, false> direct_coordinates_type;
+    typedef typename direct_coordinates_type::result_type direct_result;
+
+public:
+    template <typename Spheroid>
+    static inline bool forward(CT const& lon0, CT const& lat0,
+                               CT const& lon, CT const& lat,
+                               CT & x, CT & y,
+                               Spheroid const& spheroid)
+    {
+        inverse_result i_res = inverse_type::apply(lat0, lon0, lat, lon, spheroid);
+
+        if (math::smaller_or_equals(i_res.geodesic_scale, CT(0)))
+        {
+            return false;
+        }
+
+        CT rho = i_res.reduced_length / i_res.geodesic_scale;
+        x = sin(i_res.azimuth) * rho;
+        y = cos(i_res.azimuth) * rho;
+        
+        return true;
+    }
+
+    template <typename Spheroid>
+    static inline bool inverse(CT const& lon0, CT const& lat0,
+                               CT const& x, CT const& y,
+                               CT & lon, CT & lat,
+                               Spheroid const& spheroid)
+    {
+        CT const a = get_radius<0>(spheroid);
+        CT const b = get_radius<2>(spheroid);
+        CT const f = detail::flattening<CT>(spheroid);
+        CT const ds_threshold = a * std::numeric_limits<CT>::epsilon(); // TODO: 0 for non-fundamental type
+        
+        CT const azimuth = atan2(x, y);
+        CT const rho = math::sqrt(math::sqr(x) + math::sqr(y)); // use hypot?
+        CT distance = a * atan(rho / a);
+
+        bool found = false;
+        for (int i = 0 ; i < 10 ; ++i)
+        {
+            direct_result d_res = direct_quantities_type::apply(lon0, lat0, distance, azimuth, spheroid);
+            CT const& m = d_res.reduced_length;
+            CT const& M = d_res.geodesic_scale;
+            
+            CT const drho = m / M - rho; // rho = m / M
+            CT const ds = drho * math::sqr(M); // drho/ds = 1/M^2
+            distance -= ds;
+
+            // ds_threshold may be 0
+            if (math::abs(ds) <= ds_threshold)
+            {
+                found = true;
+                break;
+            }
+        }
+
+        if (found)
+        {
+            direct_result d_res = direct_coordinates_type::apply(lon0, lat0, distance, azimuth, spheroid);
+            lon = d_res.lon2;
+            lat = d_res.lat2;
+        }
+
+        return found;
+    }
+};
+
+}}} // namespace boost::geometry::formula
+
+
+#endif // BOOST_GEOMETRY_FORMULAS_GNOMONIC_SPHEROID_HPP