1 /******************************************************************************
2 *
3 * Project: PROJ
4 * Purpose: Generic method to compute inverse projection from forward method
5 * Author: Even Rouault <even dot rouault at spatialys dot com>
6 *
7 ******************************************************************************
8 * Copyright (c) 2018, Even Rouault <even dot rouault at spatialys dot com>
9 *
10 * Permission is hereby granted, free of charge, to any person obtaining a
11 * copy of this software and associated documentation files (the "Software"),
12 * to deal in the Software without restriction, including without limitation
13 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
14 * and/or sell copies of the Software, and to permit persons to whom the
15 * Software is furnished to do so, subject to the following conditions:
16 *
17 * The above copyright notice and this permission notice shall be included
18 * in all copies or substantial portions of the Software.
19 *
20 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
21 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
23 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
25 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
26 * DEALINGS IN THE SOFTWARE.
27 ****************************************************************************/
28
29 #include "proj_internal.h"
30
31 #include <algorithm>
32 #include <cmath>
33
34 /** Compute (lam, phi) corresponding to input (xy.x, xy.y) for projection P.
35 *
36 * Uses Newton-Raphson method, extended to 2D variables, that is using
37 * inversion of the Jacobian 2D matrix of partial derivatives. The derivatives
38 * are estimated numerically from the P->fwd method evaluated at close points.
39 *
40 * Note: thresholds used have been verified to work with adams_ws2 and wink2
41 *
42 * Starts with initial guess provided by user in lpInitial
43 */
pj_generic_inverse_2d(PJ_XY xy,PJ * P,PJ_LP lpInitial)44 PJ_LP pj_generic_inverse_2d(PJ_XY xy, PJ *P, PJ_LP lpInitial) {
45 PJ_LP lp = lpInitial;
46 double deriv_lam_X = 0;
47 double deriv_lam_Y = 0;
48 double deriv_phi_X = 0;
49 double deriv_phi_Y = 0;
50 for (int i = 0; i < 15; i++) {
51 PJ_XY xyApprox = P->fwd(lp, P);
52 const double deltaX = xyApprox.x - xy.x;
53 const double deltaY = xyApprox.y - xy.y;
54 if (fabs(deltaX) < 1e-10 && fabs(deltaY) < 1e-10) {
55 return lp;
56 }
57
58 if (fabs(deltaX) > 1e-6 || fabs(deltaY) > 1e-6) {
59 // Compute Jacobian matrix (only if we aren't close to the final
60 // result to speed things a bit)
61 PJ_LP lp2;
62 PJ_XY xy2;
63 const double dLam = lp.lam > 0 ? -1e-6 : 1e-6;
64 lp2.lam = lp.lam + dLam;
65 lp2.phi = lp.phi;
66 xy2 = P->fwd(lp2, P);
67 const double deriv_X_lam = (xy2.x - xyApprox.x) / dLam;
68 const double deriv_Y_lam = (xy2.y - xyApprox.y) / dLam;
69
70 const double dPhi = lp.phi > 0 ? -1e-6 : 1e-6;
71 lp2.lam = lp.lam;
72 lp2.phi = lp.phi + dPhi;
73 xy2 = P->fwd(lp2, P);
74 const double deriv_X_phi = (xy2.x - xyApprox.x) / dPhi;
75 const double deriv_Y_phi = (xy2.y - xyApprox.y) / dPhi;
76
77 // Inverse of Jacobian matrix
78 const double det =
79 deriv_X_lam * deriv_Y_phi - deriv_X_phi * deriv_Y_lam;
80 if (det != 0) {
81 deriv_lam_X = deriv_Y_phi / det;
82 deriv_lam_Y = -deriv_X_phi / det;
83 deriv_phi_X = -deriv_Y_lam / det;
84 deriv_phi_Y = deriv_X_lam / det;
85 }
86 }
87
88 if (xy.x != 0) {
89 // Limit the amplitude of correction to avoid overshoots due to
90 // bad initial guess
91 const double delta_lam = std::max(
92 std::min(deltaX * deriv_lam_X + deltaY * deriv_lam_Y, 0.3),
93 -0.3);
94 lp.lam -= delta_lam;
95 if (lp.lam < -M_PI)
96 lp.lam = -M_PI;
97 else if (lp.lam > M_PI)
98 lp.lam = M_PI;
99 }
100
101 if (xy.y != 0) {
102 const double delta_phi = std::max(
103 std::min(deltaX * deriv_phi_X + deltaY * deriv_phi_Y, 0.3),
104 -0.3);
105 lp.phi -= delta_phi;
106 if (lp.phi < -M_HALFPI)
107 lp.phi = -M_HALFPI;
108 else if (lp.phi > M_HALFPI)
109 lp.phi = M_HALFPI;
110 }
111 }
112 pj_ctx_set_errno(P->ctx, PJD_ERR_NON_CONVERGENT);
113 return lp;
114 }
115