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

2 files changed, 13 insertions, 81 deletions

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<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque)));
@@ -50,8 +60,8 @@ PJ *PROJECTION(wink2) {
 
     static_cast<struct pj_opaque*>(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;
 }