projections/impl/pj_generic_inverse.hpp

0001 // Boost.Geometry
0002
0003 // Copyright (c) 2022 Adam Wulkiewicz, Lodz, Poland.
0004
0005 // Copyright (c) 2022, Oracle and/or its affiliates.
0006 // Contributed and/or modified by Vissarion Fysikopoulos, 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 // This file is converted from PROJ, https://github.com/OSGeo/PROJ
0013
0014 // Last updated version of proj: 9.0.0
0015
0016 // Original copyright notice:
0017
0018 /******************************************************************************
0019  *
0020  * Project:  PROJ
0021  * Purpose:  Generic method to compute inverse projection from forward method
0022  * Author:   Even Rouault <even dot rouault at spatialys dot com>
0023  *
0024  ******************************************************************************
0025  * Copyright (c) 2018, Even Rouault <even dot rouault at spatialys dot com>
0026  *
0027  * Permission is hereby granted, free of charge, to any person obtaining a
0028  * copy of this software and associated documentation files (the "Software"),
0029  * to deal in the Software without restriction, including without limitation
0030  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
0031  * and/or sell copies of the Software, and to permit persons to whom the
0032  * Software is furnished to do so, subject to the following conditions:
0033  *
0034  * The above copyright notice and this permission notice shall be included
0035  * in all copies or substantial portions of the Software.
0036  *
0037  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
0038  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
0039  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
0040  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
0041  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
0042  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
0043  * DEALINGS IN THE SOFTWARE.
0044  ****************************************************************************/
0045
0046 #ifndef BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_GENERIC_INVERSE_HPP
0047 #define BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_GENERIC_INVERSE_HPP
0048
0049 #include <algorithm>
0050 #include <cmath>
0051
0052 #include <boost/geometry/util/math.hpp>
0053
0054 /** Compute (lam, phi) corresponding to input (xy_x, xy_y) for projection P.
0055  *
0056  * Uses Newton-Raphson method, extended to 2D variables, that is using
0057  * inversion of the Jacobian 2D matrix of partial derivatives. The derivatives
0058  * are estimated numerically from the P->fwd method evaluated at close points.
0059  *
0060  * Note: thresholds used have been verified to work with adams_ws2 and wink2
0061  *
0062  * Starts with initial guess provided by user in lpInitial
0063  */
0064
0065
0066 namespace boost { namespace geometry { namespace projections
0067 {
0068
0069 namespace detail
0070 {
0071
0072 template <typename T, typename Parameters, typename Projection>
0073 void pj_generic_inverse_2d(T const& xy_x,
0074                            T const& xy_y,
0075                            Parameters const& par,
0076                            Projection const& proj,
0077                            T& lp_lat,
0078                            T& lp_lon)
0079 {
0080     T deriv_lam_X = 0;
0081     T deriv_lam_Y = 0;
0082     T deriv_phi_X = 0;
0083     T deriv_phi_Y = 0;
0084
0085     for (int i = 0; i < 15; i++)
0086     {
0087         T xyApprox_x;
0088         T xyApprox_y;
0089         proj->fwd(par, lp_lon, lp_lat, xyApprox_x, xyApprox_y);
0090         T const deltaX = xyApprox_x - xy_x;
0091         T const deltaY = xyApprox_y - xy_y;
0092         if (fabs(deltaX) < 1e-10 && fabs(deltaY) < 1e-10) return;
0093
0094         if (fabs(deltaX) > 1e-6 || fabs(deltaY) > 1e-6)
0095         {
0096             // Compute Jacobian matrix (only if we aren't close to the final
0097             // result to speed things a bit)
0098             T lp2_lat;
0099             T lp2_lon;
0100             T xy2_x;
0101             T xy2_y;
0102             T const dLam = lp_lat > 0 ? -1e-6 : 1e-6;
0103             lp2_lat = lp_lat + dLam;
0104             lp2_lon = lp_lon;
0105             proj->fwd(par, lp2_lon, lp2_lat, xy2_x, xy2_y);
0106             //xy2 = P->fwd(lp2, P);
0107             T const deriv_X_lam = (xy2_x - xyApprox_x) / dLam;
0108             T const deriv_Y_lam = (xy2_y - xyApprox_y) / dLam;
0109
0110             T const dPhi = lp_lon > 0 ? -1e-6 : 1e-6;
0111             lp2_lat = lp_lat;
0112             lp2_lon = lp_lon + dPhi;
0113             proj->fwd(par, lp2_lon, lp2_lat, xy2_x, xy2_y);
0114             //xy2 = P->fwd(lp2, P);
0115             T const deriv_X_phi = (xy2_x - xyApprox_x) / dPhi;
0116             T const deriv_Y_phi = (xy2_y - xyApprox_y) / dPhi;
0117
0118             // Inverse of Jacobian matrix
0119             T const det = deriv_X_lam * deriv_Y_phi - deriv_X_phi * deriv_Y_lam;
0120             if (det != 0)
0121             {
0122                 deriv_lam_X = deriv_Y_phi / det;
0123                 deriv_lam_Y = -deriv_X_phi / det;
0124                 deriv_phi_X = -deriv_Y_lam / det;
0125                 deriv_phi_Y = deriv_X_lam / det;
0126             }
0127         }
0128
0129         if (xy_x != 0)
0130         {
0131             // Limit the amplitude of correction to avoid overshoots due to
0132             // bad initial guess
0133             T const delta_lam = (std::max)(
0134                 (std::min)(deltaX * deriv_lam_X + deltaY * deriv_lam_Y, 0.3),
0135                 -0.3);
0136             lp_lat -= delta_lam;
0137             if (lp_lat < -math::pi<T>())
0138                 lp_lat = -math::pi<T>();
0139             else if (lp_lat > math::pi<T>())
0140                 lp_lat = math::pi<T>();
0141         }
0142
0143         if (xy_y != 0)
0144         {
0145             T const delta_phi = (std::max)(
0146                 (std::min)(deltaX * deriv_phi_X + deltaY * deriv_phi_Y, 0.3),
0147                 -0.3);
0148             lp_lon -= delta_phi;
0149             static T const half_pi = math::half_pi<T>();
0150             if (lp_lon < -half_pi)
0151                 lp_lon = -half_pi;
0152             else if (lp_lon > half_pi)
0153                 lp_lon = half_pi;
0154         }
0155     }
0156     //pj_ctx_set_errno(P->ctx, PJD_ERR_NON_CONVERGENT);
0157 }
0158
0159 } // namespace detail
0160
0161 }}} // namespace boost::geometry::projections
0162
0163 #endif // BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_GENERIC_INVERSE_HPP