Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: Generic method to compute inverse coordinate transformation from
5 : * forward method
6 : * Author: Even Rouault <even dot rouault at spatialys dot com>
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2023, Even Rouault <even dot rouault at spatialys dot com>
10 : *
11 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include <algorithm>
31 : #include <cmath>
32 :
33 : #include "gdalgenericinverse.h"
34 : #include <cstdio>
35 :
36 : /** Compute (xOut, yOut) corresponding to input (xIn, yIn) using
37 : * the provided forward transformation to emulate the reverse direction.
38 : *
39 : * Uses Newton-Raphson method, extended to 2D variables, that is using
40 : * inversion of the Jacobian 2D matrix of partial derivatives. The derivatives
41 : * are estimated numerically from the forward method evaluated at close points.
42 : *
43 : * Starts with initial guess provided by user in (guessedXOut, guessedYOut)
44 : *
45 : * It iterates at most for 15 iterations or as soon as we get below the specified
46 : * tolerance (on input coordinates)
47 : */
48 70066 : bool GDALGenericInverse2D(double xIn, double yIn, double guessedXOut,
49 : double guessedYOut,
50 : GDALForwardCoordTransformer pfnForwardTranformer,
51 : void *pfnForwardTranformerUserData, double &xOut,
52 : double &yOut,
53 : bool computeJacobianMatrixOnlyAtFirstIter,
54 : double toleranceOnInputCoordinates,
55 : double toleranceOnOutputCoordinates)
56 : {
57 70066 : const double dfAbsValOut = std::max(fabs(guessedXOut), fabs(guessedYOut));
58 70066 : const double dfEps = dfAbsValOut > 0 ? dfAbsValOut * 1e-6 : 1e-6;
59 70066 : if (toleranceOnInputCoordinates == 0)
60 : {
61 70066 : const double dfAbsValIn = std::max(fabs(xIn), fabs(yIn));
62 70066 : toleranceOnInputCoordinates =
63 70066 : dfAbsValIn > 0 ? dfAbsValIn * 1e-12 : 1e-12;
64 : }
65 70066 : xOut = guessedXOut;
66 70066 : yOut = guessedYOut;
67 70066 : double deriv_lam_X = 0;
68 70066 : double deriv_lam_Y = 0;
69 70066 : double deriv_phi_X = 0;
70 70066 : double deriv_phi_Y = 0;
71 333611 : for (int i = 0; i < 15; i++)
72 : {
73 : double xApprox;
74 : double yApprox;
75 333585 : if (!pfnForwardTranformer(xOut, yOut, xApprox, yApprox,
76 : pfnForwardTranformerUserData))
77 70040 : return false;
78 333585 : const double deltaX = xApprox - xIn;
79 333585 : const double deltaY = yApprox - yIn;
80 333585 : if (fabs(deltaX) < toleranceOnInputCoordinates &&
81 83854 : fabs(deltaY) < toleranceOnInputCoordinates)
82 : {
83 70040 : return true;
84 : }
85 :
86 263545 : if (i == 0 || !computeJacobianMatrixOnlyAtFirstIter)
87 : {
88 : // Compute Jacobian matrix
89 63865 : double xTmp = xOut + dfEps;
90 63865 : double yTmp = yOut;
91 : double xTmpOut;
92 : double yTmpOut;
93 63865 : if (!pfnForwardTranformer(xTmp, yTmp, xTmpOut, yTmpOut,
94 : pfnForwardTranformerUserData))
95 0 : return false;
96 63865 : const double deriv_X_lam = (xTmpOut - xApprox) / dfEps;
97 63865 : const double deriv_Y_lam = (yTmpOut - yApprox) / dfEps;
98 :
99 63865 : xTmp = xOut;
100 63865 : yTmp = yOut + dfEps;
101 63865 : if (!pfnForwardTranformer(xTmp, yTmp, xTmpOut, yTmpOut,
102 : pfnForwardTranformerUserData))
103 0 : return false;
104 63865 : const double deriv_X_phi = (xTmpOut - xApprox) / dfEps;
105 63865 : const double deriv_Y_phi = (yTmpOut - yApprox) / dfEps;
106 :
107 : // Inverse of Jacobian matrix
108 63865 : const double det =
109 63865 : deriv_X_lam * deriv_Y_phi - deriv_X_phi * deriv_Y_lam;
110 63865 : if (det != 0)
111 : {
112 63865 : deriv_lam_X = deriv_Y_phi / det;
113 63865 : deriv_lam_Y = -deriv_X_phi / det;
114 63865 : deriv_phi_X = -deriv_Y_lam / det;
115 63865 : deriv_phi_Y = deriv_X_lam / det;
116 : }
117 : else
118 : {
119 0 : return false;
120 : }
121 : }
122 :
123 263545 : const double xOutDelta = deltaX * deriv_lam_X + deltaY * deriv_lam_Y;
124 263545 : const double yOutDelta = deltaX * deriv_phi_X + deltaY * deriv_phi_Y;
125 263545 : xOut -= xOutDelta;
126 263545 : yOut -= yOutDelta;
127 263545 : if (toleranceOnOutputCoordinates > 0 &&
128 0 : fabs(xOutDelta) < toleranceOnOutputCoordinates &&
129 0 : fabs(yOutDelta) < toleranceOnOutputCoordinates)
130 : {
131 0 : return true;
132 : }
133 : }
134 26 : return false;
135 : }
|
