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             :  * SPDX-License-Identifier: MIT
      12             :  ****************************************************************************/
      13             : 
      14             : #include <algorithm>
      15             : #include <cmath>
      16             : 
      17             : #include "gdalgenericinverse.h"
      18             : #include <cstdio>
      19             : 
      20             : /** Compute (xOut, yOut) corresponding to input (xIn, yIn) using
      21             :  * the provided forward transformation to emulate the reverse direction.
      22             :  *
      23             :  * Uses Newton-Raphson method, extended to 2D variables, that is using
      24             :  * inversion of the Jacobian 2D matrix of partial derivatives. The derivatives
      25             :  * are estimated numerically from the forward method evaluated at close points.
      26             :  *
      27             :  * Starts with initial guess provided by user in (guessedXOut, guessedYOut)
      28             :  *
      29             :  * It iterates at most for 15 iterations or as soon as we get below the specified
      30             :  * tolerance (on input coordinates)
      31             :  */
      32       70066 : bool GDALGenericInverse2D(double xIn, double yIn, double guessedXOut,
      33             :                           double guessedYOut,
      34             :                           GDALForwardCoordTransformer pfnForwardTranformer,
      35             :                           void *pfnForwardTranformerUserData, double &xOut,
      36             :                           double &yOut,
      37             :                           bool computeJacobianMatrixOnlyAtFirstIter,
      38             :                           double toleranceOnInputCoordinates,
      39             :                           double toleranceOnOutputCoordinates)
      40             : {
      41       70066 :     const double dfAbsValOut = std::max(fabs(guessedXOut), fabs(guessedYOut));
      42       70066 :     const double dfEps = dfAbsValOut > 0 ? dfAbsValOut * 1e-6 : 1e-6;
      43       70066 :     if (toleranceOnInputCoordinates == 0)
      44             :     {
      45       70066 :         const double dfAbsValIn = std::max(fabs(xIn), fabs(yIn));
      46       70066 :         toleranceOnInputCoordinates =
      47       70066 :             dfAbsValIn > 0 ? dfAbsValIn * 1e-12 : 1e-12;
      48             :     }
      49       70066 :     xOut = guessedXOut;
      50       70066 :     yOut = guessedYOut;
      51       70066 :     double deriv_lam_X = 0;
      52       70066 :     double deriv_lam_Y = 0;
      53       70066 :     double deriv_phi_X = 0;
      54       70066 :     double deriv_phi_Y = 0;
      55      333611 :     for (int i = 0; i < 15; i++)
      56             :     {
      57             :         double xApprox;
      58             :         double yApprox;
      59      333585 :         if (!pfnForwardTranformer(xOut, yOut, xApprox, yApprox,
      60             :                                   pfnForwardTranformerUserData))
      61       70040 :             return false;
      62      333585 :         const double deltaX = xApprox - xIn;
      63      333585 :         const double deltaY = yApprox - yIn;
      64      333585 :         if (fabs(deltaX) < toleranceOnInputCoordinates &&
      65       83854 :             fabs(deltaY) < toleranceOnInputCoordinates)
      66             :         {
      67       70040 :             return true;
      68             :         }
      69             : 
      70      263545 :         if (i == 0 || !computeJacobianMatrixOnlyAtFirstIter)
      71             :         {
      72             :             // Compute Jacobian matrix
      73       63865 :             double xTmp = xOut + dfEps;
      74       63865 :             double yTmp = yOut;
      75             :             double xTmpOut;
      76             :             double yTmpOut;
      77       63865 :             if (!pfnForwardTranformer(xTmp, yTmp, xTmpOut, yTmpOut,
      78             :                                       pfnForwardTranformerUserData))
      79           0 :                 return false;
      80       63865 :             const double deriv_X_lam = (xTmpOut - xApprox) / dfEps;
      81       63865 :             const double deriv_Y_lam = (yTmpOut - yApprox) / dfEps;
      82             : 
      83       63865 :             xTmp = xOut;
      84       63865 :             yTmp = yOut + dfEps;
      85       63865 :             if (!pfnForwardTranformer(xTmp, yTmp, xTmpOut, yTmpOut,
      86             :                                       pfnForwardTranformerUserData))
      87           0 :                 return false;
      88       63865 :             const double deriv_X_phi = (xTmpOut - xApprox) / dfEps;
      89       63865 :             const double deriv_Y_phi = (yTmpOut - yApprox) / dfEps;
      90             : 
      91             :             // Inverse of Jacobian matrix
      92       63865 :             const double det =
      93       63865 :                 deriv_X_lam * deriv_Y_phi - deriv_X_phi * deriv_Y_lam;
      94       63865 :             if (det != 0)
      95             :             {
      96       63865 :                 deriv_lam_X = deriv_Y_phi / det;
      97       63865 :                 deriv_lam_Y = -deriv_X_phi / det;
      98       63865 :                 deriv_phi_X = -deriv_Y_lam / det;
      99       63865 :                 deriv_phi_Y = deriv_X_lam / det;
     100             :             }
     101             :             else
     102             :             {
     103           0 :                 return false;
     104             :             }
     105             :         }
     106             : 
     107      263545 :         const double xOutDelta = deltaX * deriv_lam_X + deltaY * deriv_lam_Y;
     108      263545 :         const double yOutDelta = deltaX * deriv_phi_X + deltaY * deriv_phi_Y;
     109      263545 :         xOut -= xOutDelta;
     110      263545 :         yOut -= yOutDelta;
     111      263545 :         if (toleranceOnOutputCoordinates > 0 &&
     112           0 :             fabs(xOutDelta) < toleranceOnOutputCoordinates &&
     113           0 :             fabs(yOutDelta) < toleranceOnOutputCoordinates)
     114             :         {
     115           0 :             return true;
     116             :         }
     117             :     }
     118          26 :     return false;
     119             : }