From 48b3bd9afe0ffccf06b2c8dceef1b636834ba876 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 26 May 2020 15:49:03 +0200 Subject: Implement wink2 inverse by generic inversion of forward method - Move the generic method initiated from adams_ws2 to a pj_generic_inverse_2d() method - Use it in adams_ws2 - Use it in wink2 Fixes https://github.com/qgis/QGIS/issues/35512 --- src/Makefile.am | 1 + src/generic_inverse.cpp | 114 ++++++++++++++++++++++++++++++++++++++++++++++ src/lib_proj.cmake | 1 + src/proj_internal.h | 2 + src/projections/adams.cpp | 80 +------------------------------- src/projections/wink2.cpp | 14 +++++- 6 files changed, 131 insertions(+), 81 deletions(-) create mode 100644 src/generic_inverse.cpp (limited to 'src') diff --git a/src/Makefile.am b/src/Makefile.am index 3f5b2968..bbd1e49e 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -212,6 +212,7 @@ libproj_la_SOURCES = \ tsfn.cpp units.cpp ctx.cpp log.cpp zpoly1.cpp rtodms.cpp \ release.cpp gauss.cpp \ fileapi.cpp \ + generic_inverse.cpp \ \ datums.cpp datum_set.cpp transform.cpp \ utils.cpp \ diff --git a/src/generic_inverse.cpp b/src/generic_inverse.cpp new file mode 100644 index 00000000..a15cfae1 --- /dev/null +++ b/src/generic_inverse.cpp @@ -0,0 +1,114 @@ +/****************************************************************************** + * + * Project: PROJ + * Purpose: Generic method to compute inverse projection from forward method + * Author: Even Rouault + * + ****************************************************************************** + * Copyright (c) 2018, Even Rouault + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + ****************************************************************************/ + +#include "proj_internal.h" + +#include +#include + +/** Compute (lam, phi) corresponding to input (xy.x, xy.y) for projection P. + * + * Uses Newton-Raphson method, extended to 2D variables, that is using + * inversion of the Jacobian 2D matrix of partial derivatives. The derivatives + * are estimated numerically from the P->fwd method evaluated at close points. + * + * Note: thresholds used have been verified to work with adams_ws2 and wink2 + * + * Starts with initial guess provided by user in lpInitial + */ +PJ_LP pj_generic_inverse_2d(PJ_XY xy, PJ *P, PJ_LP lpInitial) { + PJ_LP lp = lpInitial; + double deriv_lam_X = 0; + double deriv_lam_Y = 0; + double deriv_phi_X = 0; + double deriv_phi_Y = 0; + for (int i = 0; i < 15; i++) { + PJ_XY xyApprox = P->fwd(lp, P); + const double deltaX = xyApprox.x - xy.x; + const double deltaY = xyApprox.y - xy.y; + if (fabs(deltaX) < 1e-10 && fabs(deltaY) < 1e-10) { + return lp; + } + + if (fabs(deltaX) > 1e-6 || fabs(deltaY) > 1e-6) { + // Compute Jacobian matrix (only if we aren't close to the final + // result to speed things a bit) + PJ_LP lp2; + PJ_XY xy2; + const double dLam = lp.lam > 0 ? -1e-6 : 1e-6; + lp2.lam = lp.lam + dLam; + lp2.phi = lp.phi; + xy2 = P->fwd(lp2, P); + const double deriv_X_lam = (xy2.x - xyApprox.x) / dLam; + const double deriv_Y_lam = (xy2.y - xyApprox.y) / dLam; + + const double dPhi = lp.phi > 0 ? -1e-6 : 1e-6; + lp2.lam = lp.lam; + lp2.phi = lp.phi + dPhi; + xy2 = P->fwd(lp2, P); + const double deriv_X_phi = (xy2.x - xyApprox.x) / dPhi; + const double deriv_Y_phi = (xy2.y - xyApprox.y) / dPhi; + + // Inverse of Jacobian matrix + const double det = + deriv_X_lam * deriv_Y_phi - deriv_X_phi * deriv_Y_lam; + if (det != 0) { + deriv_lam_X = deriv_Y_phi / det; + deriv_lam_Y = -deriv_X_phi / det; + deriv_phi_X = -deriv_Y_lam / det; + deriv_phi_Y = deriv_X_lam / det; + } + } + + if (xy.x != 0) { + // Limit the amplitude of correction to avoid overshoots due to + // bad initial guess + const double delta_lam = std::max( + std::min(deltaX * deriv_lam_X + deltaY * deriv_lam_Y, 0.3), + -0.3); + lp.lam -= delta_lam; + if (lp.lam < -M_PI) + lp.lam = -M_PI; + else if (lp.lam > M_PI) + lp.lam = M_PI; + } + + if (xy.y != 0) { + const double delta_phi = std::max( + std::min(deltaX * deriv_phi_X + deltaY * deriv_phi_Y, 0.3), + -0.3); + lp.phi -= delta_phi; + if (lp.phi < -M_HALFPI) + lp.phi = -M_HALFPI; + else if (lp.phi > M_HALFPI) + lp.phi = M_HALFPI; + } + } + pj_ctx_set_errno(P->ctx, PJD_ERR_NON_CONVERGENT); + return lp; +} diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake index 4eca0256..abc9cc4b 100644 --- a/src/lib_proj.cmake +++ b/src/lib_proj.cmake @@ -220,6 +220,7 @@ set(SRC_LIBPROJ_CORE fileapi.cpp fwd.cpp gauss.cpp + generic_inverse.cpp geodesic.c init.cpp initcache.cpp diff --git a/src/proj_internal.h b/src/proj_internal.h index 618133bb..738e8ac0 100644 --- a/src/proj_internal.h +++ b/src/proj_internal.h @@ -902,6 +902,8 @@ const PJ_UNITS *pj_list_angular_units(); void pj_clear_hgridshift_knowngrids_cache(); void pj_clear_vgridshift_knowngrids_cache(); +PJ_LP pj_generic_inverse_2d(PJ_XY xy, PJ *P, PJ_LP lpInitial); + /* classic public API */ #include "proj_api.h" diff --git a/src/projections/adams.cpp b/src/projections/adams.cpp index e704c687..4f7d1a03 100644 --- a/src/projections/adams.cpp +++ b/src/projections/adams.cpp @@ -182,84 +182,6 @@ static PJ_XY adams_forward(PJ_LP lp, PJ *P) { return xy; } -static PJ_LP inverse_with_newton_raphson(PJ_XY xy, - PJ *P, - PJ_LP lp, // initial guess - PJ_XY (*fwd)(PJ_LP, PJ *)) -{ - double deriv_lam_X = 0; - double deriv_lam_Y = 0; - double deriv_phi_X = 0; - double deriv_phi_Y = 0; - for(int i = 0; i < 15; i++ ) - { - PJ_XY xyApprox = fwd(lp, P); - const double deltaX = xyApprox.x - xy.x; - const double deltaY = xyApprox.y - xy.y; - if( fabs(deltaX) < 1e-10 && fabs(deltaY) < 1e-10 ) - { - return lp; - } - - if( fabs(deltaX) > 1e-6 || fabs(deltaY) > 1e-6 ) - { - // Compute Jacobian matrix (only if we aren't close to the final - // result to speed things a bit) - PJ_LP lp2; - PJ_XY xy2; - const double dLam = lp.lam > 0 ? -1e-6 : 1e-6; - lp2.lam = lp.lam + dLam; - lp2.phi = lp.phi; - xy2 = fwd(lp2, P); - const double deriv_X_lam = (xy2.x - xyApprox.x) / dLam; - const double deriv_Y_lam = (xy2.y - xyApprox.y) / dLam; - - const double dPhi = lp.phi > 0 ? -1e-6 : 1e-6; - lp2.lam = lp.lam; - lp2.phi = lp.phi + dPhi; - xy2 = fwd(lp2, P); - const double deriv_X_phi = (xy2.x - xyApprox.x) / dPhi; - const double deriv_Y_phi = (xy2.y - xyApprox.y) / dPhi; - - // Inverse of Jacobian matrix - const double det = deriv_X_lam * deriv_Y_phi - deriv_X_phi * deriv_Y_lam; - if( det != 0 ) - { - deriv_lam_X = deriv_Y_phi / det; - deriv_lam_Y = -deriv_X_phi / det; - deriv_phi_X = -deriv_Y_lam / det; - deriv_phi_Y = deriv_X_lam / det; - } - } - - if( xy.x != 0 ) - { - // Limit the amplitude of correction to avoid overshoots due to - // bad initial guess - const double delta_lam = std::max(std::min( - deltaX * deriv_lam_X + deltaY * deriv_lam_Y, 0.3), -0.3); - lp.lam -= delta_lam; - if( lp.lam < -M_PI ) - lp.lam = -M_PI; - else if( lp.lam > M_PI ) - lp.lam = M_PI; - } - - if( xy.y != 0 ) - { - const double delta_phi = std::max(std::min( - deltaX * deriv_phi_X + deltaY * deriv_phi_Y, 0.3), -0.3); - lp.phi -= delta_phi; - if( lp.phi < -M_HALFPI ) - lp.phi = -M_HALFPI ; - else if( lp.phi > M_HALFPI ) - lp.phi = M_HALFPI; - } - } - pj_ctx_set_errno( P->ctx, PJD_ERR_NON_CONVERGENT ); - return lp; -} - static PJ_LP adams_inverse(PJ_XY xy, PJ *P) { PJ_LP lp; @@ -277,7 +199,7 @@ static PJ_LP adams_inverse(PJ_XY xy, PJ *P) lp.lam = fabs(lp.phi) >= M_HALFPI ? 0 : std::max(std::min(xy.x / 2.62205760 / cos(lp.phi), 1.0), -1.0) * M_PI; - return inverse_with_newton_raphson(xy, P, lp, adams_forward); + return pj_generic_inverse_2d(xy, P, lp); } diff --git a/src/projections/wink2.cpp b/src/projections/wink2.cpp index 7fd08649..d457f842 100644 --- a/src/projections/wink2.cpp +++ b/src/projections/wink2.cpp @@ -6,7 +6,7 @@ #include "proj.h" #include "proj_internal.h" -PROJ_HEAD(wink2, "Winkel II") "\n\tPCyl, Sph, no inv\n\tlat_1="; +PROJ_HEAD(wink2, "Winkel II") "\n\tPCyl, Sph\n\tlat_1="; namespace { // anonymous namespace struct pj_opaque { @@ -41,6 +41,16 @@ static PJ_XY wink2_s_forward (PJ_LP lp, PJ *P) { /* Spheroidal, forwar return xy; } +static PJ_LP wink2_s_inverse(PJ_XY xy, PJ *P) +{ + PJ_LP lpInit; + + lpInit.phi = xy.y; + lpInit.lam = xy.x; + + return pj_generic_inverse_2d(xy, P, lpInit); +} + PJ *PROJECTION(wink2) { struct pj_opaque *Q = static_cast(pj_calloc (1, sizeof (struct pj_opaque))); @@ -50,8 +60,8 @@ PJ *PROJECTION(wink2) { static_cast(P->opaque)->cosphi1 = cos(pj_param(P->ctx, P->params, "rlat_1").f); P->es = 0.; - P->inv = nullptr; P->fwd = wink2_s_forward; + P->inv = wink2_s_inverse; return P; } -- cgit v1.2.3