LCOV - code coverage report
Current view: top level - alg - gdaltransformer.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 1671 1891 88.4 %
Date: 2025-07-09 17:50:03 Functions: 62 66 93.9 %

          Line data    Source code
       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  Mapinfo Image Warper
       4             :  * Purpose:  Implementation of one or more GDALTrasformerFunc types, including
       5             :  *           the GenImgProj (general image reprojector) transformer.
       6             :  * Author:   Frank Warmerdam, warmerdam@pobox.com
       7             :  *
       8             :  ******************************************************************************
       9             :  * Copyright (c) 2002, i3 - information integration and imaging
      10             :  *                          Fort Collin, CO
      11             :  * Copyright (c) 2008-2013, Even Rouault <even dot rouault at spatialys.com>
      12             :  * Copyright (c) 2021, CLS
      13             :  *
      14             :  * SPDX-License-Identifier: MIT
      15             :  ****************************************************************************/
      16             : 
      17             : #include "cpl_port.h"
      18             : #include "gdal_alg.h"
      19             : #include "gdal_alg_priv.h"
      20             : 
      21             : #include <climits>
      22             : #include <cmath>
      23             : #include <cstddef>
      24             : #include <cstdlib>
      25             : #include <cstring>
      26             : 
      27             : #include <algorithm>
      28             : #include <limits>
      29             : #include <utility>
      30             : 
      31             : #include "cpl_conv.h"
      32             : #include "cpl_error.h"
      33             : #include "cpl_list.h"
      34             : #include "cpl_minixml.h"
      35             : #include "cpl_multiproc.h"
      36             : #include "cpl_string.h"
      37             : #include "cpl_vsi.h"
      38             : #include "gdal.h"
      39             : #include "gdal_priv.h"
      40             : #include "ogr_core.h"
      41             : #include "ogr_spatialref.h"
      42             : #include "ogr_srs_api.h"
      43             : 
      44             : CPL_C_START
      45             : void *GDALDeserializeGCPTransformer(CPLXMLNode *psTree);
      46             : void *GDALDeserializeTPSTransformer(CPLXMLNode *psTree);
      47             : void *GDALDeserializeGeoLocTransformer(CPLXMLNode *psTree);
      48             : void *GDALDeserializeRPCTransformer(CPLXMLNode *psTree);
      49             : void *GDALDeserializeHomographyTransformer(CPLXMLNode *psTree);
      50             : CPL_C_END
      51             : 
      52             : static CPLXMLNode *GDALSerializeReprojectionTransformer(void *pTransformArg);
      53             : static void *GDALDeserializeReprojectionTransformer(CPLXMLNode *psTree);
      54             : 
      55             : static CPLXMLNode *GDALSerializeGenImgProjTransformer(void *pTransformArg);
      56             : static void *GDALDeserializeGenImgProjTransformer(CPLXMLNode *psTree);
      57             : 
      58             : static void *GDALCreateApproxTransformer2(GDALTransformerFunc pfnRawTransformer,
      59             :                                           void *pRawTransformerArg,
      60             :                                           double dfMaxErrorForward,
      61             :                                           double dfMaxErrorReverse);
      62             : 
      63             : /************************************************************************/
      64             : /*                            GDALIsTransformer()                       */
      65             : /************************************************************************/
      66             : 
      67       13861 : bool GDALIsTransformer(void *hTransformerArg, const char *pszClassName)
      68             : {
      69       13861 :     if (!hTransformerArg)
      70         558 :         return false;
      71             :     // All transformers should have a GDALTransformerInfo member as their first members
      72       13303 :     GDALTransformerInfo *psInfo =
      73             :         static_cast<GDALTransformerInfo *>(hTransformerArg);
      74       13303 :     return memcmp(psInfo->abySignature, GDAL_GTI2_SIGNATURE,
      75       25712 :                   strlen(GDAL_GTI2_SIGNATURE)) == 0 &&
      76       25712 :            strcmp(psInfo->pszClassName, pszClassName) == 0;
      77             : }
      78             : 
      79             : /************************************************************************/
      80             : /*                          GDALTransformFunc                           */
      81             : /*                                                                      */
      82             : /*      Documentation for GDALTransformFunc typedef.                    */
      83             : /************************************************************************/
      84             : 
      85             : /*!
      86             : 
      87             : \typedef typedef int (*GDALTransformerFunc)( void *pTransformerArg, int
      88             : bDstToSrc, int nPointCount, double *x, double *y, double *z, int *panSuccess );
      89             : 
      90             : Generic signature for spatial point transformers.
      91             : 
      92             : This function signature is used for a variety of functions that accept
      93             : passed in functions used to transform point locations between two coordinate
      94             : spaces.
      95             : 
      96             : The GDALCreateGenImgProjTransformer(), GDALCreateReprojectionTransformerEx(),
      97             : GDALCreateGCPTransformer() and GDALCreateApproxTransformer() functions can
      98             : be used to prepare argument data for some built-in transformers.  As well,
      99             : applications can implement their own transformers to the following signature.
     100             : 
     101             : \code
     102             : typedef int
     103             : (*GDALTransformerFunc)( void *pTransformerArg,
     104             :                         int bDstToSrc, int nPointCount,
     105             :                         double *x, double *y, double *z, int *panSuccess );
     106             : \endcode
     107             : 
     108             : @param pTransformerArg application supplied callback data used by the
     109             : transformer.
     110             : 
     111             : @param bDstToSrc if TRUE the transformation will be from the destination
     112             : coordinate space to the source coordinate system, otherwise the transformation
     113             : will be from the source coordinate system to the destination coordinate system.
     114             : 
     115             : @param nPointCount number of points in the x, y and z arrays.
     116             : 
     117             : @param[in,out] x input X coordinates.  Results returned in same array.
     118             : 
     119             : @param[in,out] y input Y coordinates.  Results returned in same array.
     120             : 
     121             : @param[in,out] z input Z coordinates.  Results returned in same array.
     122             : 
     123             : @param[out] panSuccess array of ints in which success (TRUE) or failure (FALSE)
     124             : flags are returned for the translation of each point. Must not be NULL.
     125             : 
     126             : @return TRUE if all points have been successfully transformed (changed in 3.11,
     127             : previously was TRUE if some points have been successfully transformed)
     128             : 
     129             : */
     130             : 
     131             : /************************************************************************/
     132             : /*                      GDALSuggestedWarpOutput()                       */
     133             : /************************************************************************/
     134             : 
     135             : /**
     136             :  * Suggest output file size.
     137             :  *
     138             :  * This function is used to suggest the size, and georeferenced extents
     139             :  * appropriate given the indicated transformation and input file.  It walks
     140             :  * the edges of the input file (approximately 20 sample points along each
     141             :  * edge) transforming into output coordinates in order to get an extents box.
     142             :  *
     143             :  * Then a resolution is computed with the intent that the length of the
     144             :  * distance from the top left corner of the output imagery to the bottom right
     145             :  * corner would represent the same number of pixels as in the source image.
     146             :  * Note that if the image is somewhat rotated the diagonal taken isn't of the
     147             :  * whole output bounding rectangle, but instead of the locations where the
     148             :  * top/left and bottom/right corners transform.  The output pixel size is
     149             :  * always square.  This is intended to approximately preserve the resolution
     150             :  * of the input data in the output file.
     151             :  *
     152             :  * The values returned in padfGeoTransformOut, pnPixels and pnLines are
     153             :  * the suggested number of pixels and lines for the output file, and the
     154             :  * geotransform relating those pixels to the output georeferenced coordinates.
     155             :  *
     156             :  * The trickiest part of using the function is ensuring that the
     157             :  * transformer created is from source file pixel/line coordinates to
     158             :  * output file georeferenced coordinates.  This can be accomplished with
     159             :  * GDALCreateGenImgProjTransformer() by passing a NULL for the hDstDS.
     160             :  *
     161             :  * @param hSrcDS the input image (it is assumed the whole input image is
     162             :  * being transformed).
     163             :  * @param pfnTransformer the transformer function.
     164             :  * @param pTransformArg the callback data for the transformer function.
     165             :  * @param padfGeoTransformOut the array of six doubles in which the suggested
     166             :  * geotransform is returned.
     167             :  * @param pnPixels int in which the suggest pixel width of output is returned.
     168             :  * @param pnLines int in which the suggest pixel height of output is returned.
     169             :  *
     170             :  * @return CE_None if successful or CE_Failure otherwise.
     171             :  */
     172             : 
     173          41 : CPLErr CPL_STDCALL GDALSuggestedWarpOutput(GDALDatasetH hSrcDS,
     174             :                                            GDALTransformerFunc pfnTransformer,
     175             :                                            void *pTransformArg,
     176             :                                            double *padfGeoTransformOut,
     177             :                                            int *pnPixels, int *pnLines)
     178             : 
     179             : {
     180          41 :     VALIDATE_POINTER1(hSrcDS, "GDALSuggestedWarpOutput", CE_Failure);
     181             : 
     182          41 :     double adfExtent[4] = {};
     183             : 
     184          41 :     return GDALSuggestedWarpOutput2(hSrcDS, pfnTransformer, pTransformArg,
     185             :                                     padfGeoTransformOut, pnPixels, pnLines,
     186          41 :                                     adfExtent, 0);
     187             : }
     188             : 
     189         580 : static bool GDALSuggestedWarpOutput2_MustAdjustForRightBorder(
     190             :     GDALTransformerFunc pfnTransformer, void *pTransformArg, double *padfExtent,
     191             :     int /* nPixels*/, int nLines, double dfPixelSizeX, double dfPixelSizeY)
     192             : {
     193         580 :     double adfX[21] = {};
     194         580 :     double adfY[21] = {};
     195             : 
     196         580 :     const double dfMaxXOut = padfExtent[2];
     197         580 :     const double dfMaxYOut = padfExtent[3];
     198             : 
     199             :     // Take 20 steps.
     200         580 :     int nSamplePoints = 0;
     201       12760 :     for (double dfRatio = 0.0; dfRatio <= 1.01; dfRatio += 0.05)
     202             :     {
     203             :         // Ensure we end exactly at the end.
     204       12180 :         if (dfRatio > 0.99)
     205         580 :             dfRatio = 1.0;
     206             : 
     207             :         // Along right.
     208       12180 :         adfX[nSamplePoints] = dfMaxXOut;
     209       12180 :         adfY[nSamplePoints] = dfMaxYOut - dfPixelSizeY * dfRatio * nLines;
     210       12180 :         nSamplePoints++;
     211             :     }
     212         580 :     double adfZ[21] = {};
     213             : 
     214         580 :     int abSuccess[21] = {};
     215             : 
     216         580 :     pfnTransformer(pTransformArg, TRUE, nSamplePoints, adfX, adfY, adfZ,
     217             :                    abSuccess);
     218             : 
     219         580 :     int abSuccess2[21] = {};
     220             : 
     221         580 :     pfnTransformer(pTransformArg, FALSE, nSamplePoints, adfX, adfY, adfZ,
     222             :                    abSuccess2);
     223             : 
     224         580 :     nSamplePoints = 0;
     225         580 :     int nBadCount = 0;
     226       12760 :     for (double dfRatio = 0.0; dfRatio <= 1.01; dfRatio += 0.05)
     227             :     {
     228       12180 :         const double expected_x = dfMaxXOut;
     229       12180 :         const double expected_y = dfMaxYOut - dfPixelSizeY * dfRatio * nLines;
     230       12180 :         if (!abSuccess[nSamplePoints] || !abSuccess2[nSamplePoints] ||
     231        9368 :             fabs(adfX[nSamplePoints] - expected_x) > dfPixelSizeX ||
     232        8078 :             fabs(adfY[nSamplePoints] - expected_y) > dfPixelSizeY)
     233             :         {
     234        4102 :             nBadCount++;
     235             :         }
     236       12180 :         nSamplePoints++;
     237             :     }
     238             : 
     239         580 :     return nBadCount == nSamplePoints;
     240             : }
     241             : 
     242         487 : static bool GDALSuggestedWarpOutput2_MustAdjustForBottomBorder(
     243             :     GDALTransformerFunc pfnTransformer, void *pTransformArg, double *padfExtent,
     244             :     int nPixels, int /* nLines */, double dfPixelSizeX, double dfPixelSizeY)
     245             : {
     246         487 :     double adfX[21] = {};
     247         487 :     double adfY[21] = {};
     248             : 
     249         487 :     const double dfMinXOut = padfExtent[0];
     250         487 :     const double dfMinYOut = padfExtent[1];
     251             : 
     252             :     // Take 20 steps.
     253         487 :     int nSamplePoints = 0;
     254       10714 :     for (double dfRatio = 0.0; dfRatio <= 1.01; dfRatio += 0.05)
     255             :     {
     256             :         // Ensure we end exactly at the end.
     257       10227 :         if (dfRatio > 0.99)
     258         487 :             dfRatio = 1.0;
     259             : 
     260             :         // Along right.
     261       10227 :         adfX[nSamplePoints] = dfMinXOut + dfPixelSizeX * dfRatio * nPixels;
     262       10227 :         adfY[nSamplePoints] = dfMinYOut;
     263       10227 :         nSamplePoints++;
     264             :     }
     265         487 :     double adfZ[21] = {};
     266             : 
     267         487 :     int abSuccess[21] = {};
     268             : 
     269         487 :     pfnTransformer(pTransformArg, TRUE, nSamplePoints, adfX, adfY, adfZ,
     270             :                    abSuccess);
     271             : 
     272         487 :     int abSuccess2[21] = {};
     273             : 
     274         487 :     pfnTransformer(pTransformArg, FALSE, nSamplePoints, adfX, adfY, adfZ,
     275             :                    abSuccess2);
     276             : 
     277         487 :     nSamplePoints = 0;
     278         487 :     int nBadCount = 0;
     279       10714 :     for (double dfRatio = 0.0; dfRatio <= 1.01; dfRatio += 0.05)
     280             :     {
     281       10227 :         const double expected_x = dfMinXOut + dfPixelSizeX * dfRatio * nPixels;
     282       10227 :         const double expected_y = dfMinYOut;
     283       10227 :         if (!abSuccess[nSamplePoints] || !abSuccess2[nSamplePoints] ||
     284        8299 :             fabs(adfX[nSamplePoints] - expected_x) > dfPixelSizeX ||
     285        8010 :             fabs(adfY[nSamplePoints] - expected_y) > dfPixelSizeY)
     286             :         {
     287        2217 :             nBadCount++;
     288             :         }
     289       10227 :         nSamplePoints++;
     290             :     }
     291             : 
     292         487 :     return nBadCount == nSamplePoints;
     293             : }
     294             : 
     295             : /************************************************************************/
     296             : /*                      GDALSuggestedWarpOutput2()                      */
     297             : /************************************************************************/
     298             : 
     299             : /**
     300             :  * Suggest output file size.
     301             :  *
     302             :  * This function is used to suggest the size, and georeferenced extents
     303             :  * appropriate given the indicated transformation and input file.  It walks
     304             :  * the edges of the input file (approximately 20 sample points along each
     305             :  * edge) transforming into output coordinates in order to get an extents box.
     306             :  *
     307             :  * Then a resolution is computed with the intent that the length of the
     308             :  * distance from the top left corner of the output imagery to the bottom right
     309             :  * corner would represent the same number of pixels as in the source image.
     310             :  * Note that if the image is somewhat rotated the diagonal taken isn't of the
     311             :  * whole output bounding rectangle, but instead of the locations where the
     312             :  * top/left and bottom/right corners transform.  The output pixel size is
     313             :  * always square.  This is intended to approximately preserve the resolution
     314             :  * of the input data in the output file.
     315             :  *
     316             :  * The values returned in padfGeoTransformOut, pnPixels and pnLines are
     317             :  * the suggested number of pixels and lines for the output file, and the
     318             :  * geotransform relating those pixels to the output georeferenced coordinates.
     319             :  *
     320             :  * The trickiest part of using the function is ensuring that the
     321             :  * transformer created is from source file pixel/line coordinates to
     322             :  * output file georeferenced coordinates.  This can be accomplished with
     323             :  * GDALCreateGenImgProjTransformer() by passing a NULL for the hDstDS.
     324             :  *
     325             :  * @param hSrcDS the input image (it is assumed the whole input image is
     326             :  * being transformed).
     327             :  * @param pfnTransformer the transformer function.
     328             :  * @param pTransformArg the callback data for the transformer function.
     329             :  * @param padfGeoTransformOut the array of six doubles in which the suggested
     330             :  * geotransform is returned.
     331             :  * @param pnPixels int in which the suggest pixel width of output is returned.
     332             :  * @param pnLines int in which the suggest pixel height of output is returned.
     333             :  * @param padfExtent Four entry array to return extents as (xmin, ymin, xmax,
     334             :  * ymax).
     335             :  * @param nOptions Options flags. Zero or GDAL_SWO_ROUND_UP_SIZE  to ask *pnPixels
     336             :  * and *pnLines to be rounded up instead of being rounded to the closes integer, or
     337             :  * GDAL_SWO_FORCE_SQUARE_PIXEL to indicate that the generated pixel size is a square.
     338             :  *
     339             :  * @return CE_None if successful or CE_Failure otherwise.
     340             :  */
     341             : 
     342        1042 : CPLErr CPL_STDCALL GDALSuggestedWarpOutput2(GDALDatasetH hSrcDS,
     343             :                                             GDALTransformerFunc pfnTransformer,
     344             :                                             void *pTransformArg,
     345             :                                             double *padfGeoTransformOut,
     346             :                                             int *pnPixels, int *pnLines,
     347             :                                             double *padfExtent, int nOptions)
     348             : {
     349        1042 :     VALIDATE_POINTER1(hSrcDS, "GDALSuggestedWarpOutput2", CE_Failure);
     350             : 
     351             :     const bool bIsGDALGenImgProjTransform{
     352        2084 :         pTransformArg &&
     353        1042 :         GDALIsTransformer(pTransformArg, GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME)};
     354             : 
     355             :     /* -------------------------------------------------------------------- */
     356             :     /*      Setup sample points all around the edge of the input raster.    */
     357             :     /* -------------------------------------------------------------------- */
     358        1042 :     if (bIsGDALGenImgProjTransform)
     359             :     {
     360             :         // In case CHECK_WITH_INVERT_PROJ has been modified.
     361        1042 :         GDALRefreshGenImgProjTransformer(pTransformArg);
     362             :     }
     363           0 :     else if (GDALIsTransformer(pTransformArg,
     364             :                                GDAL_APPROX_TRANSFORMER_CLASS_NAME))
     365             :     {
     366             :         // In case CHECK_WITH_INVERT_PROJ has been modified.
     367           0 :         GDALRefreshApproxTransformer(pTransformArg);
     368             :     }
     369             : 
     370        1042 :     const int nInXSize = GDALGetRasterXSize(hSrcDS);
     371        1042 :     const int nInYSize = GDALGetRasterYSize(hSrcDS);
     372             : 
     373             :     /* ------------------------------------------------------------- */
     374             :     /* Special case for warping on the same (or null) CRS.           */
     375             :     /* ------------------------------------------------------------- */
     376        1042 :     if ((!nOptions || (nOptions & GDAL_SWO_FORCE_SQUARE_PIXEL) == 0) &&
     377        1041 :         pTransformArg && bIsGDALGenImgProjTransform)
     378             :     {
     379        1041 :         const GDALGenImgProjTransformInfo *psInfo =
     380             :             static_cast<const GDALGenImgProjTransformInfo *>(pTransformArg);
     381             : 
     382        1041 :         if (!psInfo->sSrcParams.pTransformer &&
     383         968 :             !psInfo->bHasCustomTransformationPipeline &&
     384         964 :             !psInfo->sDstParams.pTransformer &&
     385         964 :             psInfo->sSrcParams.adfGeoTransform[2] == 0 &&
     386         964 :             psInfo->sSrcParams.adfGeoTransform[4] == 0 &&
     387         964 :             psInfo->sDstParams.adfGeoTransform[0] == 0 &&
     388         952 :             psInfo->sDstParams.adfGeoTransform[1] == 1 &&
     389         952 :             psInfo->sDstParams.adfGeoTransform[2] == 0 &&
     390         952 :             psInfo->sDstParams.adfGeoTransform[3] == 0 &&
     391         952 :             psInfo->sDstParams.adfGeoTransform[4] == 0 &&
     392         952 :             psInfo->sDstParams.adfGeoTransform[5] == 1)
     393             :         {
     394         952 :             const OGRSpatialReference *poSourceCRS = nullptr;
     395         952 :             const OGRSpatialReference *poTargetCRS = nullptr;
     396             : 
     397         952 :             if (psInfo->pReprojectArg)
     398             :             {
     399         680 :                 const GDALReprojectionTransformInfo *psRTI =
     400             :                     static_cast<const GDALReprojectionTransformInfo *>(
     401             :                         psInfo->pReprojectArg);
     402         680 :                 poSourceCRS = psRTI->poForwardTransform->GetSourceCS();
     403         680 :                 poTargetCRS = psRTI->poForwardTransform->GetTargetCS();
     404             :             }
     405             : 
     406        1632 :             if ((!poSourceCRS && !poTargetCRS) ||
     407         680 :                 (poSourceCRS && poTargetCRS &&
     408         680 :                  poSourceCRS->IsSame(poTargetCRS)))
     409             :             {
     410             : 
     411         620 :                 const bool bNorthUp{psInfo->sSrcParams.adfGeoTransform[5] <
     412             :                                     0.0};
     413             : 
     414         620 :                 memcpy(padfGeoTransformOut, psInfo->sSrcParams.adfGeoTransform,
     415             :                        sizeof(double) * 6);
     416             : 
     417         620 :                 if (!bNorthUp)
     418             :                 {
     419          62 :                     padfGeoTransformOut[3] = padfGeoTransformOut[3] +
     420          62 :                                              nInYSize * padfGeoTransformOut[5];
     421          62 :                     padfGeoTransformOut[5] = -padfGeoTransformOut[5];
     422             :                 }
     423             : 
     424         620 :                 *pnPixels = nInXSize;
     425         620 :                 *pnLines = nInYSize;
     426             : 
     427             :                 // Calculate extent from hSrcDS
     428         620 :                 if (padfExtent)
     429             :                 {
     430         620 :                     padfExtent[0] = psInfo->sSrcParams.adfGeoTransform[0];
     431         620 :                     padfExtent[1] =
     432         620 :                         psInfo->sSrcParams.adfGeoTransform[3] +
     433         620 :                         nInYSize * psInfo->sSrcParams.adfGeoTransform[5];
     434         620 :                     padfExtent[2] =
     435         620 :                         psInfo->sSrcParams.adfGeoTransform[0] +
     436         620 :                         nInXSize * psInfo->sSrcParams.adfGeoTransform[1];
     437         620 :                     padfExtent[3] = psInfo->sSrcParams.adfGeoTransform[3];
     438         620 :                     if (!bNorthUp)
     439             :                     {
     440          62 :                         std::swap(padfExtent[1], padfExtent[3]);
     441             :                     }
     442             :                 }
     443         620 :                 return CE_None;
     444             :             }
     445             :         }
     446             :     }
     447             : 
     448         422 :     const int N_PIXELSTEP = 50;
     449             :     int nSteps = static_cast<int>(
     450         422 :         static_cast<double>(std::min(nInYSize, nInXSize)) / N_PIXELSTEP + 0.5);
     451         422 :     if (nSteps < 20)
     452         397 :         nSteps = 20;
     453          25 :     else if (nSteps > 100)
     454          10 :         nSteps = 100;
     455             : 
     456             :     // TODO(rouault): How is this goto retry supposed to work?  Added in r20537.
     457             :     // Does redoing the same malloc multiple times work?  If it is needed, can
     458             :     // it be converted to a tigher while loop around the MALLOC3s and free?  Is
     459             :     // the point to try with the full requested steps.  Then, if there is not
     460             :     // enough memory, back off and try with just 20 steps?
     461         422 : retry:
     462         422 :     int nStepsPlusOne = nSteps + 1;
     463         422 :     int nSampleMax = nStepsPlusOne * nStepsPlusOne;
     464             : 
     465         422 :     double dfStep = 1.0 / nSteps;
     466         422 :     double *padfY = nullptr;
     467         422 :     double *padfZ = nullptr;
     468         422 :     double *padfYRevert = nullptr;
     469         422 :     double *padfZRevert = nullptr;
     470             : 
     471             :     int *pabSuccess = static_cast<int *>(
     472         422 :         VSI_MALLOC3_VERBOSE(sizeof(int), nStepsPlusOne, nStepsPlusOne));
     473             :     double *padfX = static_cast<double *>(
     474         422 :         VSI_MALLOC3_VERBOSE(sizeof(double) * 3, nStepsPlusOne, nStepsPlusOne));
     475             :     double *padfXRevert = static_cast<double *>(
     476         422 :         VSI_MALLOC3_VERBOSE(sizeof(double) * 3, nStepsPlusOne, nStepsPlusOne));
     477         422 :     if (pabSuccess == nullptr || padfX == nullptr || padfXRevert == nullptr)
     478             :     {
     479           0 :         CPLFree(padfX);
     480           0 :         CPLFree(padfXRevert);
     481           0 :         CPLFree(pabSuccess);
     482           0 :         if (nSteps > 20)
     483             :         {
     484           0 :             nSteps = 20;
     485           0 :             goto retry;
     486             :         }
     487           0 :         return CE_Failure;
     488             :     }
     489             : 
     490         422 :     padfY = padfX + nSampleMax;
     491         422 :     padfZ = padfX + nSampleMax * 2;
     492         422 :     padfYRevert = padfXRevert + nSampleMax;
     493         422 :     padfZRevert = padfXRevert + nSampleMax * 2;
     494             : 
     495             :     // Take N_STEPS steps.
     496       10187 :     for (int iStep = 0; iStep <= nSteps; iStep++)
     497             :     {
     498        9765 :         double dfRatio = (iStep == nSteps) ? 1.0 : iStep * dfStep;
     499        9765 :         int iStep2 = iStep;
     500             : 
     501             :         // Along top.
     502        9765 :         padfX[iStep2] = dfRatio * nInXSize;
     503        9765 :         padfY[iStep2] = 0.0;
     504        9765 :         padfZ[iStep2] = 0.0;
     505             : 
     506             :         // Along bottom.
     507        9765 :         iStep2 += nStepsPlusOne;
     508        9765 :         padfX[iStep2] = dfRatio * nInXSize;
     509        9765 :         padfY[iStep2] = nInYSize;
     510        9765 :         padfZ[iStep2] = 0.0;
     511             : 
     512             :         // Along left.
     513        9765 :         iStep2 += nStepsPlusOne;
     514        9765 :         padfX[iStep2] = 0.0;
     515        9765 :         padfY[iStep2] = dfRatio * nInYSize;
     516        9765 :         padfZ[iStep2] = 0.0;
     517             : 
     518             :         // Along right.
     519        9765 :         iStep2 += nStepsPlusOne;
     520        9765 :         padfX[iStep2] = nInXSize;
     521        9765 :         padfY[iStep2] = dfRatio * nInYSize;
     522        9765 :         padfZ[iStep2] = 0.0;
     523             :     }
     524             : 
     525         422 :     int nSamplePoints = 4 * nStepsPlusOne;
     526             : 
     527         422 :     memset(pabSuccess, 1, sizeof(int) * nSampleMax);
     528             : 
     529             :     /* -------------------------------------------------------------------- */
     530             :     /*      Transform them to the output coordinate system.                 */
     531             :     /* -------------------------------------------------------------------- */
     532             :     {
     533         844 :         CPLTurnFailureIntoWarningBackuper oErrorsToWarnings{};
     534         422 :         pfnTransformer(pTransformArg, FALSE, nSamplePoints, padfX, padfY, padfZ,
     535             :                        pabSuccess);
     536             :     }
     537         422 :     constexpr int SIGN_FINAL_UNINIT = -2;
     538         422 :     constexpr int SIGN_FINAL_INVALID = 0;
     539         422 :     int iSignDiscontinuity = SIGN_FINAL_UNINIT;
     540         422 :     int nFailedCount = 0;
     541         422 :     const int iSignArray[2] = {-1, 1};
     542       39482 :     for (int i = 0; i < nSamplePoints; i++)
     543             :     {
     544       39060 :         if (pabSuccess[i])
     545             :         {
     546             :             // Fix for https://trac.osgeo.org/gdal/ticket/7243
     547             :             // where echo "-2050000.000 2050000.000" |
     548             :             //              gdaltransform -s_srs EPSG:3411 -t_srs EPSG:4326
     549             :             // gives "-180 63.691332898492"
     550             :             // but we would rather like 180
     551       34999 :             if (iSignDiscontinuity == 1 || iSignDiscontinuity == -1)
     552             :             {
     553       10631 :                 if (!((iSignDiscontinuity * padfX[i] > 0 &&
     554       10563 :                        iSignDiscontinuity * padfX[i] <= 180.0) ||
     555          69 :                       (fabs(padfX[i] - iSignDiscontinuity * -180.0) < 1e-8)))
     556             :                 {
     557          46 :                     iSignDiscontinuity = SIGN_FINAL_INVALID;
     558             :                 }
     559             :             }
     560       24368 :             else if (iSignDiscontinuity == SIGN_FINAL_UNINIT)
     561             :             {
     562         987 :                 for (const auto &iSign : iSignArray)
     563             :                 {
     564         727 :                     if ((iSign * padfX[i] > 0 && iSign * padfX[i] <= 180.0) ||
     565         567 :                         (fabs(padfX[i] - iSign * -180.0) < 1e-8))
     566             :                     {
     567         160 :                         iSignDiscontinuity = iSign;
     568         160 :                         break;
     569             :                     }
     570             :                 }
     571         420 :                 if (iSignDiscontinuity == SIGN_FINAL_UNINIT)
     572             :                 {
     573         260 :                     iSignDiscontinuity = SIGN_FINAL_INVALID;
     574             :                 }
     575             :             }
     576             :         }
     577             :         else
     578             :         {
     579        4061 :             nFailedCount++;
     580             :         }
     581             :     }
     582             : 
     583         422 :     if (iSignDiscontinuity == 1 || iSignDiscontinuity == -1)
     584             :     {
     585       11018 :         for (int i = 0; i < nSamplePoints; i++)
     586             :         {
     587       10904 :             if (pabSuccess[i])
     588             :             {
     589       10287 :                 if (fabs(padfX[i] - iSignDiscontinuity * -180.0) < 1e-8)
     590             :                 {
     591           2 :                     double axTemp[2] = {iSignDiscontinuity * -180.0,
     592           2 :                                         iSignDiscontinuity * 180.0};
     593           2 :                     double ayTemp[2] = {padfY[i], padfY[i]};
     594           2 :                     double azTemp[2] = {padfZ[i], padfZ[i]};
     595           2 :                     int abSuccess[2] = {FALSE, FALSE};
     596           4 :                     CPLTurnFailureIntoWarningBackuper oErrorsToWarnings{};
     597           2 :                     if (pfnTransformer(pTransformArg, TRUE, 2, axTemp, ayTemp,
     598           2 :                                        azTemp, abSuccess) &&
     599           4 :                         fabs(axTemp[0] - axTemp[1]) < 1e-8 &&
     600           2 :                         fabs(ayTemp[0] - ayTemp[1]) < 1e-8)
     601             :                     {
     602           2 :                         padfX[i] = iSignDiscontinuity * 180.0;
     603             :                     }
     604             :                 }
     605             :             }
     606             :         }
     607             :     }
     608             : 
     609             :     /* -------------------------------------------------------------------- */
     610             :     /*      Check if the computed target coordinates are revertable.        */
     611             :     /*      If not, try the detailed grid sampling.                         */
     612             :     /* -------------------------------------------------------------------- */
     613         422 :     if (nFailedCount)
     614             :     {
     615          86 :         CPLDebug("WARP", "At least one point failed after direct transform");
     616             :     }
     617             :     else
     618             :     {
     619         336 :         memcpy(padfXRevert, padfX, nSamplePoints * sizeof(double));
     620         336 :         memcpy(padfYRevert, padfY, nSamplePoints * sizeof(double));
     621         336 :         memcpy(padfZRevert, padfZ, nSamplePoints * sizeof(double));
     622             :         {
     623         672 :             CPLTurnFailureIntoWarningBackuper oErrorsToWarnings{};
     624         336 :             pfnTransformer(pTransformArg, TRUE, nSamplePoints, padfXRevert,
     625             :                            padfYRevert, padfZRevert, pabSuccess);
     626             :         }
     627             : 
     628       30054 :         for (int i = 0; nFailedCount == 0 && i < nSamplePoints; i++)
     629             :         {
     630       29734 :             if (!pabSuccess[i])
     631             :             {
     632          16 :                 nFailedCount++;
     633          16 :                 break;
     634             :             }
     635             : 
     636       29718 :             double dfRatio = (i % nStepsPlusOne) * dfStep;
     637       29718 :             if (dfRatio > 0.99)
     638        1266 :                 dfRatio = 1.0;
     639             : 
     640       29718 :             double dfExpectedX = 0.0;
     641       29718 :             double dfExpectedY = 0.0;
     642       29718 :             if (i < nStepsPlusOne)
     643             :             {
     644        7563 :                 dfExpectedX = dfRatio * nInXSize;
     645             :             }
     646       22155 :             else if (i < 2 * nStepsPlusOne)
     647             :             {
     648        7513 :                 dfExpectedX = dfRatio * nInXSize;
     649        7513 :                 dfExpectedY = nInYSize;
     650             :             }
     651       14642 :             else if (i < 3 * nStepsPlusOne)
     652             :             {
     653        7383 :                 dfExpectedY = dfRatio * nInYSize;
     654             :             }
     655             :             else
     656             :             {
     657        7259 :                 dfExpectedX = nInXSize;
     658        7259 :                 dfExpectedY = dfRatio * nInYSize;
     659             :             }
     660             : 
     661       29718 :             if (fabs(padfXRevert[i] - dfExpectedX) >
     662       29718 :                     nInXSize / static_cast<double>(nSteps) ||
     663       29706 :                 fabs(padfYRevert[i] - dfExpectedY) >
     664       29706 :                     nInYSize / static_cast<double>(nSteps))
     665          12 :                 nFailedCount++;
     666             :         }
     667         336 :         if (nFailedCount != 0)
     668          28 :             CPLDebug("WARP",
     669             :                      "At least one point failed after revert transform");
     670             :     }
     671             : 
     672             :     /* -------------------------------------------------------------------- */
     673             :     /*      If any of the edge points failed to transform, we need to       */
     674             :     /*      build a fairly detailed internal grid of points instead to      */
     675             :     /*      help identify the area that is transformable.                   */
     676             :     /* -------------------------------------------------------------------- */
     677         422 :     if (nFailedCount)
     678             :     {
     679         114 :         nSamplePoints = 0;
     680             : 
     681             :         // Take N_STEPS steps.
     682        2662 :         for (int iStep = 0; iStep <= nSteps; iStep++)
     683             :         {
     684        2548 :             double dfRatio = (iStep == nSteps) ? 1.0 : iStep * dfStep;
     685             : 
     686       69582 :             for (int iStep2 = 0; iStep2 <= nSteps; iStep2++)
     687             :             {
     688       67034 :                 const double dfRatio2 =
     689       67034 :                     iStep2 == nSteps ? 1.0 : iStep2 * dfStep;
     690             : 
     691             :                 // From top to bottom, from left to right.
     692       67034 :                 padfX[nSamplePoints] = dfRatio2 * nInXSize;
     693       67034 :                 padfY[nSamplePoints] = dfRatio * nInYSize;
     694       67034 :                 padfZ[nSamplePoints] = 0.0;
     695       67034 :                 nSamplePoints++;
     696             :             }
     697             :         }
     698             : 
     699         114 :         CPLAssert(nSamplePoints == nSampleMax);
     700             : 
     701             :         {
     702         228 :             CPLTurnFailureIntoWarningBackuper oErrorsToWarnings{};
     703         114 :             pfnTransformer(pTransformArg, FALSE, nSamplePoints, padfX, padfY,
     704             :                            padfZ, pabSuccess);
     705             :         }
     706             :     }
     707             : 
     708             :     /* -------------------------------------------------------------------- */
     709             :     /*      Collect the bounds, ignoring any failed points.                 */
     710             :     /* -------------------------------------------------------------------- */
     711         422 :     double dfMinXOut = 0.0;
     712         422 :     double dfMinYOut = 0.0;
     713         422 :     double dfMaxXOut = 0.0;
     714         422 :     double dfMaxYOut = 0.0;
     715         422 :     bool bGotInitialPoint = false;
     716             : 
     717         422 :     nFailedCount = 0;
     718       96324 :     for (int i = 0; i < nSamplePoints; i++)
     719             :     {
     720       95902 :         int x_i = 0;
     721       95902 :         int y_i = 0;
     722             : 
     723       95902 :         if (nSamplePoints == nSampleMax)
     724             :         {
     725       67034 :             x_i = i % nStepsPlusOne;
     726       67034 :             y_i = i / nStepsPlusOne;
     727             :         }
     728             :         else
     729             :         {
     730       28868 :             if (i < 2 * nStepsPlusOne)
     731             :             {
     732       14434 :                 x_i = i % nStepsPlusOne;
     733       14434 :                 y_i = (i < nStepsPlusOne) ? 0 : nSteps;
     734             :             }
     735             :         }
     736             : 
     737       95902 :         if (x_i > 0 && (pabSuccess[i - 1] || pabSuccess[i]))
     738             :         {
     739       65625 :             double x_out_before = padfX[i - 1];
     740       65625 :             double x_out_after = padfX[i];
     741       65625 :             int nIter = 0;
     742       65625 :             double x_in_before =
     743       65625 :                 static_cast<double>(x_i - 1) * nInXSize / nSteps;
     744       65625 :             double x_in_after = static_cast<double>(x_i) * nInXSize / nSteps;
     745       65625 :             int invalid_before = !(pabSuccess[i - 1]);
     746       65625 :             int invalid_after = !(pabSuccess[i]);
     747             : 
     748             :             // Detect discontinuity in target coordinates when the target x
     749             :             // coordinates change sign. This may be a false positive when the
     750             :             // target tx is around 0 Dichotomic search to reduce the interval
     751             :             // to near the discontinuity and get a better out extent.
     752       82886 :             while ((invalid_before || invalid_after ||
     753      173439 :                     x_out_before * x_out_after < 0.0) &&
     754             :                    nIter < 16)
     755             :             {
     756       24928 :                 double x = (x_in_before + x_in_after) / 2.0;
     757       24928 :                 double y = static_cast<double>(y_i) * nInYSize / nSteps;
     758       24928 :                 double z = 0.0;
     759       24928 :                 int bSuccess = TRUE;
     760       24928 :                 if (pfnTransformer(pTransformArg, FALSE, 1, &x, &y, &z,
     761       43046 :                                    &bSuccess) &&
     762       18118 :                     bSuccess)
     763             :                 {
     764       18118 :                     if (bGotInitialPoint)
     765             :                     {
     766       18098 :                         dfMinXOut = std::min(dfMinXOut, x);
     767       18098 :                         dfMinYOut = std::min(dfMinYOut, y);
     768       18098 :                         dfMaxXOut = std::max(dfMaxXOut, x);
     769       18098 :                         dfMaxYOut = std::max(dfMaxYOut, y);
     770             :                     }
     771             :                     else
     772             :                     {
     773          20 :                         bGotInitialPoint = true;
     774          20 :                         dfMinXOut = x;
     775          20 :                         dfMaxXOut = x;
     776          20 :                         dfMinYOut = y;
     777          20 :                         dfMaxYOut = y;
     778             :                     }
     779             : 
     780       18118 :                     if (invalid_before || x_out_before * x < 0)
     781             :                     {
     782        9857 :                         invalid_after = FALSE;
     783        9857 :                         x_in_after = (x_in_before + x_in_after) / 2.0;
     784        9857 :                         x_out_after = x;
     785             :                     }
     786             :                     else
     787             :                     {
     788        8261 :                         invalid_before = FALSE;
     789        8261 :                         x_out_before = x;
     790        8261 :                         x_in_before = (x_in_before + x_in_after) / 2.0;
     791             :                     }
     792             :                 }
     793             :                 else
     794             :                 {
     795        6810 :                     if (invalid_before)
     796             :                     {
     797        3402 :                         x_in_before = (x_in_before + x_in_after) / 2.0;
     798             :                     }
     799        3408 :                     else if (invalid_after)
     800             :                     {
     801        3408 :                         x_in_after = (x_in_before + x_in_after) / 2.0;
     802             :                     }
     803             :                     else
     804             :                     {
     805           0 :                         break;
     806             :                     }
     807             :                 }
     808       24928 :                 nIter++;
     809             :             }
     810             :         }
     811             : 
     812       95902 :         if (!pabSuccess[i])
     813             :         {
     814       13718 :             nFailedCount++;
     815       13718 :             continue;
     816             :         }
     817             : 
     818       82184 :         if (bGotInitialPoint)
     819             :         {
     820       81784 :             dfMinXOut = std::min(dfMinXOut, padfX[i]);
     821       81784 :             dfMinYOut = std::min(dfMinYOut, padfY[i]);
     822       81784 :             dfMaxXOut = std::max(dfMaxXOut, padfX[i]);
     823       81784 :             dfMaxYOut = std::max(dfMaxYOut, padfY[i]);
     824             :         }
     825             :         else
     826             :         {
     827         400 :             bGotInitialPoint = true;
     828         400 :             dfMinXOut = padfX[i];
     829         400 :             dfMaxXOut = padfX[i];
     830         400 :             dfMinYOut = padfY[i];
     831         400 :             dfMaxYOut = padfY[i];
     832             :         }
     833             :     }
     834             : 
     835         422 :     if (nFailedCount > nSamplePoints - 10)
     836             :     {
     837           5 :         CPLError(CE_Failure, CPLE_AppDefined,
     838             :                  "Too many points (%d out of %d) failed to transform, "
     839             :                  "unable to compute output bounds.",
     840             :                  nFailedCount, nSamplePoints);
     841             : 
     842           5 :         CPLFree(padfX);
     843           5 :         CPLFree(padfXRevert);
     844           5 :         CPLFree(pabSuccess);
     845             : 
     846           5 :         return CE_Failure;
     847             :     }
     848             : 
     849         417 :     if (nFailedCount)
     850          82 :         CPLDebug("GDAL",
     851             :                  "GDALSuggestedWarpOutput(): %d out of %d points failed to "
     852             :                  "transform.",
     853             :                  nFailedCount, nSamplePoints);
     854             : 
     855         417 :     bool bIsGeographicCoordsDeg = false;
     856         417 :     if (bIsGDALGenImgProjTransform)
     857             :     {
     858         417 :         const GDALGenImgProjTransformInfo *pGIPTI =
     859             :             static_cast<const GDALGenImgProjTransformInfo *>(pTransformArg);
     860         417 :         if (pGIPTI->sSrcParams.pTransformer == GDALGeoLocTransform &&
     861          28 :             pGIPTI->sDstParams.pTransformer == nullptr &&
     862          28 :             pGIPTI->sDstParams.adfGeoTransform[0] == 0 &&
     863          26 :             pGIPTI->sDstParams.adfGeoTransform[1] == 1 &&
     864          26 :             pGIPTI->sDstParams.adfGeoTransform[2] == 0 &&
     865          26 :             pGIPTI->sDstParams.adfGeoTransform[3] == 0 &&
     866          26 :             pGIPTI->sDstParams.adfGeoTransform[4] == 0 &&
     867          26 :             pGIPTI->sDstParams.adfGeoTransform[5] == 1)
     868             :         {
     869             :             /* --------------------------------------------------------------------
     870             :              */
     871             :             /*      Special case for geolocation array, to quickly find the
     872             :              * bounds. */
     873             :             /* --------------------------------------------------------------------
     874             :              */
     875          26 :             const GDALGeoLocTransformInfo *pGLTI =
     876             :                 static_cast<const GDALGeoLocTransformInfo *>(
     877             :                     pGIPTI->sSrcParams.pTransformArg);
     878             : 
     879          26 :             if (pGIPTI->pReproject == nullptr)
     880             :             {
     881             :                 const char *pszGLSRS =
     882          26 :                     CSLFetchNameValue(pGLTI->papszGeolocationInfo, "SRS");
     883          26 :                 if (pszGLSRS == nullptr)
     884             :                 {
     885           4 :                     bIsGeographicCoordsDeg = true;
     886             :                 }
     887             :                 else
     888             :                 {
     889          44 :                     OGRSpatialReference oSRS;
     890          22 :                     if (oSRS.SetFromUserInput(pszGLSRS) == OGRERR_NONE &&
     891          42 :                         oSRS.IsGeographic() &&
     892          20 :                         std::fabs(oSRS.GetAngularUnits() -
     893          20 :                                   CPLAtof(SRS_UA_DEGREE_CONV)) < 1e-9)
     894             :                     {
     895          20 :                         bIsGeographicCoordsDeg = true;
     896             :                     }
     897             :                 }
     898             :             }
     899             : 
     900         208 :             for (const auto &xy :
     901          26 :                  {std::pair<double, double>(pGLTI->dfMinX, pGLTI->dfYAtMinX),
     902          26 :                   std::pair<double, double>(pGLTI->dfXAtMinY, pGLTI->dfMinY),
     903          26 :                   std::pair<double, double>(pGLTI->dfMaxX, pGLTI->dfYAtMaxX),
     904         130 :                   std::pair<double, double>(pGLTI->dfXAtMaxY, pGLTI->dfMaxY)})
     905             :             {
     906         104 :                 double x = xy.first;
     907         104 :                 double y = xy.second;
     908         104 :                 if (pGLTI->bSwapXY)
     909             :                 {
     910           4 :                     std::swap(x, y);
     911             :                 }
     912         104 :                 double xOut = std::numeric_limits<double>::quiet_NaN();
     913         104 :                 double yOut = std::numeric_limits<double>::quiet_NaN();
     914         104 :                 if (pGIPTI->pReproject == nullptr ||
     915           0 :                     pGIPTI->pReproject(pGIPTI->pReprojectArg, false, 1, &x, &y,
     916             :                                        nullptr, nullptr))
     917             :                 {
     918         104 :                     xOut = x;
     919         104 :                     yOut = y;
     920             :                 }
     921         104 :                 dfMinXOut = std::min(dfMinXOut, xOut);
     922         104 :                 dfMinYOut = std::min(dfMinYOut, yOut);
     923         104 :                 dfMaxXOut = std::max(dfMaxXOut, xOut);
     924         104 :                 dfMaxYOut = std::max(dfMaxYOut, yOut);
     925          26 :             }
     926             :         }
     927         391 :         else if (pGIPTI->sSrcParams.pTransformer == nullptr &&
     928         348 :                  pGIPTI->sDstParams.pTransformer == nullptr &&
     929         348 :                  pGIPTI->pReproject == GDALReprojectionTransform &&
     930         337 :                  pGIPTI->sDstParams.adfGeoTransform[0] == 0 &&
     931         335 :                  pGIPTI->sDstParams.adfGeoTransform[1] == 1 &&
     932         335 :                  pGIPTI->sDstParams.adfGeoTransform[2] == 0 &&
     933         335 :                  pGIPTI->sDstParams.adfGeoTransform[3] == 0 &&
     934         335 :                  pGIPTI->sDstParams.adfGeoTransform[4] == 0 &&
     935         335 :                  pGIPTI->sDstParams.adfGeoTransform[5] == 1)
     936             :         {
     937             :             /* ------------------------------------------------------------- */
     938             :             /* Special case for warping using source geotransform and        */
     939             :             /* reprojection to deal with the poles.                          */
     940             :             /* ------------------------------------------------------------- */
     941         335 :             const GDALReprojectionTransformInfo *psRTI =
     942             :                 static_cast<const GDALReprojectionTransformInfo *>(
     943             :                     pGIPTI->pReprojectArg);
     944             :             const OGRSpatialReference *poSourceCRS =
     945         335 :                 psRTI->poForwardTransform->GetSourceCS();
     946             :             const OGRSpatialReference *poTargetCRS =
     947         335 :                 psRTI->poForwardTransform->GetTargetCS();
     948         669 :             if (poTargetCRS != nullptr &&
     949         334 :                 psRTI->poReverseTransform != nullptr &&
     950         334 :                 poTargetCRS->IsGeographic() &&
     951         116 :                 fabs(poTargetCRS->GetAngularUnits() -
     952         785 :                      CPLAtof(SRS_UA_DEGREE_CONV)) < 1e-9 &&
     953         116 :                 (!poSourceCRS || !poSourceCRS->IsGeographic()))
     954             :             {
     955          89 :                 bIsGeographicCoordsDeg = true;
     956             : 
     957          89 :                 std::unique_ptr<CPLConfigOptionSetter> poSetter;
     958          89 :                 if (pGIPTI->bCheckWithInvertPROJ)
     959             :                 {
     960             :                     // CHECK_WITH_INVERT_PROJ=YES prevent reliable
     961             :                     // transformation of poles.
     962           4 :                     poSetter = std::make_unique<CPLConfigOptionSetter>(
     963           4 :                         "CHECK_WITH_INVERT_PROJ", "NO", false);
     964           4 :                     GDALRefreshGenImgProjTransformer(pTransformArg);
     965             :                     // GDALRefreshGenImgProjTransformer() has invalidated psRTI
     966           4 :                     psRTI = static_cast<const GDALReprojectionTransformInfo *>(
     967             :                         pGIPTI->pReprojectArg);
     968             :                 }
     969             : 
     970         267 :                 for (const auto &sign : iSignArray)
     971             :                 {
     972         178 :                     double X = 0.0;
     973         178 :                     const double Yinit = 90.0 * sign;
     974         178 :                     double Y = Yinit;
     975         178 :                     if (psRTI->poReverseTransform->Transform(1, &X, &Y))
     976             :                     {
     977         122 :                         const auto invGT =
     978             :                             pGIPTI->sSrcParams.adfInvGeoTransform;
     979         122 :                         const double x = invGT[0] + X * invGT[1] + Y * invGT[2];
     980         122 :                         const double y = invGT[3] + X * invGT[4] + Y * invGT[5];
     981         122 :                         constexpr double EPSILON = 1e-5;
     982         122 :                         if (x >= -EPSILON && x <= nInXSize + EPSILON &&
     983          26 :                             y >= -EPSILON && y <= nInYSize + EPSILON)
     984             :                         {
     985           6 :                             if (psRTI->poForwardTransform->Transform(1, &X,
     986          12 :                                                                      &Y) &&
     987           6 :                                 fabs(Y - Yinit) <= 1e-6)
     988             :                             {
     989           6 :                                 bool bMinXMaxXSet = false;
     990           6 :                                 if (poSourceCRS)
     991             :                                 {
     992             :                                     const char *pszProjection =
     993           6 :                                         poSourceCRS->GetAttrValue("PROJECTION");
     994           6 :                                     if (pszProjection &&
     995           6 :                                         EQUAL(pszProjection,
     996             :                                               SRS_PT_ORTHOGRAPHIC))
     997             :                                     {
     998             :                                         const double dfLon0 =
     999           4 :                                             poSourceCRS->GetNormProjParm(
    1000             :                                                 SRS_PP_CENTRAL_MERIDIAN, 0.0);
    1001           4 :                                         dfMinXOut = dfLon0 - 90;
    1002           4 :                                         dfMaxXOut = dfLon0 + 90;
    1003           4 :                                         bMinXMaxXSet = true;
    1004             :                                     }
    1005             :                                 }
    1006           6 :                                 if (!bMinXMaxXSet)
    1007             :                                 {
    1008           2 :                                     dfMinXOut = -180;
    1009           2 :                                     dfMaxXOut = 180;
    1010             :                                 }
    1011           6 :                                 if (sign < 0)
    1012           2 :                                     dfMinYOut = Yinit;
    1013             :                                 else
    1014           4 :                                     dfMaxYOut = Yinit;
    1015             :                             }
    1016             :                         }
    1017             :                     }
    1018             :                 }
    1019             : 
    1020          89 :                 if (poSetter)
    1021             :                 {
    1022           4 :                     poSetter.reset();
    1023           4 :                     GDALRefreshGenImgProjTransformer(pTransformArg);
    1024           4 :                     pGIPTI = static_cast<const GDALGenImgProjTransformInfo *>(
    1025             :                         pTransformArg);
    1026           4 :                     psRTI = static_cast<const GDALReprojectionTransformInfo *>(
    1027             :                         pGIPTI->pReprojectArg);
    1028           4 :                     poSourceCRS = psRTI->poForwardTransform->GetSourceCS();
    1029           4 :                     poTargetCRS = psRTI->poForwardTransform->GetTargetCS();
    1030             :                 }
    1031             :             }
    1032             : 
    1033             :             // Use TransformBounds() to handle more particular cases
    1034         335 :             if (poSourceCRS != nullptr && poTargetCRS != nullptr &&
    1035         334 :                 pGIPTI->sSrcParams.adfGeoTransform[1] != 0 &&
    1036         334 :                 pGIPTI->sSrcParams.adfGeoTransform[2] == 0 &&
    1037         334 :                 pGIPTI->sSrcParams.adfGeoTransform[4] == 0 &&
    1038         334 :                 pGIPTI->sSrcParams.adfGeoTransform[5] != 0)
    1039             :             {
    1040         334 :                 const double dfULX = pGIPTI->sSrcParams.adfGeoTransform[0];
    1041         334 :                 const double dfULY = pGIPTI->sSrcParams.adfGeoTransform[3];
    1042         334 :                 const double dfLRX =
    1043         334 :                     dfULX + pGIPTI->sSrcParams.adfGeoTransform[1] * nInXSize;
    1044         334 :                 const double dfLRY =
    1045         334 :                     dfULY + pGIPTI->sSrcParams.adfGeoTransform[5] * nInYSize;
    1046         334 :                 const double dfMinSrcX = std::min(dfULX, dfLRX);
    1047         334 :                 const double dfMinSrcY = std::min(dfULY, dfLRY);
    1048         334 :                 const double dfMaxSrcX = std::max(dfULX, dfLRX);
    1049         334 :                 const double dfMaxSrcY = std::max(dfULY, dfLRY);
    1050         334 :                 double dfTmpMinXOut = std::numeric_limits<double>::max();
    1051         334 :                 double dfTmpMinYOut = std::numeric_limits<double>::max();
    1052         334 :                 double dfTmpMaxXOut = std::numeric_limits<double>::min();
    1053         334 :                 double dfTmpMaxYOut = std::numeric_limits<double>::min();
    1054         668 :                 if (psRTI->poForwardTransform->TransformBounds(
    1055             :                         dfMinSrcX, dfMinSrcY, dfMaxSrcX, dfMaxSrcY,
    1056             :                         &dfTmpMinXOut, &dfTmpMinYOut, &dfTmpMaxXOut,
    1057             :                         &dfTmpMaxYOut,
    1058         334 :                         2))  // minimum number of points as we already have a
    1059             :                              // logic above to sample
    1060             :                 {
    1061         327 :                     dfMinXOut = std::min(dfMinXOut, dfTmpMinXOut);
    1062         327 :                     dfMinYOut = std::min(dfMinYOut, dfTmpMinYOut);
    1063         327 :                     dfMaxXOut = std::max(dfMaxXOut, dfTmpMaxXOut);
    1064         327 :                     dfMaxYOut = std::max(dfMaxYOut, dfTmpMaxYOut);
    1065             :                 }
    1066             :             }
    1067             :         }
    1068             :     }
    1069             : 
    1070             :     /* -------------------------------------------------------------------- */
    1071             :     /*      Compute the distance in "georeferenced" units from the top      */
    1072             :     /*      corner of the transformed input image to the bottom left        */
    1073             :     /*      corner of the transformed input.  Use this distance to          */
    1074             :     /*      compute an approximate pixel size in the output                 */
    1075             :     /*      georeferenced coordinates.                                      */
    1076             :     /* -------------------------------------------------------------------- */
    1077         417 :     double dfDiagonalDist = 0.0;
    1078         417 :     double dfDeltaX = 0.0;
    1079         417 :     double dfDeltaY = 0.0;
    1080             : 
    1081         417 :     if (pabSuccess[0] && pabSuccess[nSamplePoints - 1])
    1082             :     {
    1083         340 :         dfDeltaX = padfX[nSamplePoints - 1] - padfX[0];
    1084         340 :         dfDeltaY = padfY[nSamplePoints - 1] - padfY[0];
    1085             :         // In some cases this can result in 0 values. See #5980
    1086             :         // Fallback to safer method in that case.
    1087             :     }
    1088         417 :     if (dfDeltaX == 0.0 || dfDeltaY == 0.0)
    1089             :     {
    1090          83 :         dfDeltaX = dfMaxXOut - dfMinXOut;
    1091          83 :         dfDeltaY = dfMaxYOut - dfMinYOut;
    1092             :     }
    1093             : 
    1094         417 :     dfDiagonalDist = sqrt(dfDeltaX * dfDeltaX + dfDeltaY * dfDeltaY);
    1095             : 
    1096             :     /* -------------------------------------------------------------------- */
    1097             :     /*      Compute a pixel size from this.                                 */
    1098             :     /* -------------------------------------------------------------------- */
    1099             :     const double dfPixelSize =
    1100         417 :         dfDiagonalDist / sqrt(static_cast<double>(nInXSize) * nInXSize +
    1101         417 :                               static_cast<double>(nInYSize) * nInYSize);
    1102             : 
    1103         417 :     const double dfPixels = (dfMaxXOut - dfMinXOut) / dfPixelSize;
    1104         417 :     const double dfLines = (dfMaxYOut - dfMinYOut) / dfPixelSize;
    1105             : 
    1106         417 :     const int knIntMaxMinusOne = std::numeric_limits<int>::max() - 1;
    1107         417 :     if (dfPixels > knIntMaxMinusOne || dfLines > knIntMaxMinusOne)
    1108             :     {
    1109           0 :         CPLError(CE_Failure, CPLE_AppDefined,
    1110             :                  "Computed dimensions are too big : %.0f x %.0f",
    1111             :                  dfPixels + 0.5, dfLines + 0.5);
    1112             : 
    1113           0 :         CPLFree(padfX);
    1114           0 :         CPLFree(padfXRevert);
    1115           0 :         CPLFree(pabSuccess);
    1116             : 
    1117           0 :         return CE_Failure;
    1118             :     }
    1119             : 
    1120         417 :     if ((nOptions & GDAL_SWO_ROUND_UP_SIZE) != 0)
    1121             :     {
    1122           8 :         constexpr double EPS = 1e-5;
    1123           8 :         *pnPixels = static_cast<int>(std::ceil(dfPixels - EPS));
    1124           8 :         *pnLines = static_cast<int>(std::ceil(dfLines - EPS));
    1125             :     }
    1126             :     else
    1127             :     {
    1128         409 :         *pnPixels = static_cast<int>(dfPixels + 0.5);
    1129         409 :         *pnLines = static_cast<int>(dfLines + 0.5);
    1130             :     }
    1131             : 
    1132         417 :     double dfPixelSizeX = dfPixelSize;
    1133         417 :     double dfPixelSizeY = dfPixelSize;
    1134             : 
    1135         417 :     const double adfRatioArray[] = {0.000, 0.001, 0.010, 0.100, 1.000};
    1136             : 
    1137             :     /* -------------------------------------------------------------------- */
    1138             :     /*      Check that the right border is not completely out of source     */
    1139             :     /*      image. If so, adjust the x pixel size a bit in the hope it will */
    1140             :     /*      fit.                                                            */
    1141             :     /* -------------------------------------------------------------------- */
    1142         588 :     for (const auto &dfRatio : adfRatioArray)
    1143             :     {
    1144         580 :         const double dfTryPixelSizeX =
    1145         580 :             dfPixelSizeX - dfPixelSizeX * dfRatio / *pnPixels;
    1146         580 :         double adfExtent[4] = {dfMinXOut, dfMaxYOut - (*pnLines) * dfPixelSizeY,
    1147         580 :                                dfMinXOut + (*pnPixels) * dfTryPixelSizeX,
    1148         580 :                                dfMaxYOut};
    1149         580 :         if (!GDALSuggestedWarpOutput2_MustAdjustForRightBorder(
    1150             :                 pfnTransformer, pTransformArg, adfExtent, *pnPixels, *pnLines,
    1151             :                 dfTryPixelSizeX, dfPixelSizeY))
    1152             :         {
    1153         409 :             dfPixelSizeX = dfTryPixelSizeX;
    1154         409 :             break;
    1155             :         }
    1156             :     }
    1157             : 
    1158             :     /* -------------------------------------------------------------------- */
    1159             :     /*      Check that the bottom border is not completely out of source    */
    1160             :     /*      image. If so, adjust the y pixel size a bit in the hope it will */
    1161             :     /*      fit.                                                            */
    1162             :     /* -------------------------------------------------------------------- */
    1163         497 :     for (const auto &dfRatio : adfRatioArray)
    1164             :     {
    1165         487 :         const double dfTryPixelSizeY =
    1166         487 :             dfPixelSizeY - dfPixelSizeY * dfRatio / *pnLines;
    1167             :         double adfExtent[4] = {
    1168         487 :             dfMinXOut, dfMaxYOut - (*pnLines) * dfTryPixelSizeY,
    1169         487 :             dfMinXOut + (*pnPixels) * dfPixelSizeX, dfMaxYOut};
    1170         487 :         if (!GDALSuggestedWarpOutput2_MustAdjustForBottomBorder(
    1171             :                 pfnTransformer, pTransformArg, adfExtent, *pnPixels, *pnLines,
    1172             :                 dfPixelSizeX, dfTryPixelSizeY))
    1173             :         {
    1174         407 :             dfPixelSizeY = dfTryPixelSizeY;
    1175         407 :             break;
    1176             :         }
    1177             :     }
    1178             : 
    1179             :     /* -------------------------------------------------------------------- */
    1180             :     /*      Recompute some bounds so that all return values are consistent  */
    1181             :     /* -------------------------------------------------------------------- */
    1182         417 :     double dfMaxXOutNew = dfMinXOut + (*pnPixels) * dfPixelSizeX;
    1183         417 :     if (bIsGeographicCoordsDeg &&
    1184         113 :         ((dfMaxXOut <= 180 && dfMaxXOutNew > 180) || dfMaxXOut == 180))
    1185             :     {
    1186           3 :         dfMaxXOut = 180;
    1187           3 :         dfPixelSizeX = (dfMaxXOut - dfMinXOut) / *pnPixels;
    1188             :     }
    1189             :     else
    1190             :     {
    1191         414 :         dfMaxXOut = dfMaxXOutNew;
    1192             :     }
    1193             : 
    1194         417 :     double dfMinYOutNew = dfMaxYOut - (*pnLines) * dfPixelSizeY;
    1195         417 :     if (bIsGeographicCoordsDeg && dfMinYOut >= -90 && dfMinYOutNew < -90)
    1196             :     {
    1197           0 :         dfMinYOut = -90;
    1198           0 :         dfPixelSizeY = (dfMaxYOut - dfMinYOut) / *pnLines;
    1199             :     }
    1200             :     else
    1201             :     {
    1202         417 :         dfMinYOut = dfMinYOutNew;
    1203             :     }
    1204             : 
    1205             :     /* -------------------------------------------------------------------- */
    1206             :     /*      Return raw extents.                                             */
    1207             :     /* -------------------------------------------------------------------- */
    1208         417 :     padfExtent[0] = dfMinXOut;
    1209         417 :     padfExtent[1] = dfMinYOut;
    1210         417 :     padfExtent[2] = dfMaxXOut;
    1211         417 :     padfExtent[3] = dfMaxYOut;
    1212             : 
    1213             :     /* -------------------------------------------------------------------- */
    1214             :     /*      Set the output geotransform.                                    */
    1215             :     /* -------------------------------------------------------------------- */
    1216         417 :     padfGeoTransformOut[0] = dfMinXOut;
    1217         417 :     padfGeoTransformOut[1] = dfPixelSizeX;
    1218         417 :     padfGeoTransformOut[2] = 0.0;
    1219         417 :     padfGeoTransformOut[3] = dfMaxYOut;
    1220         417 :     padfGeoTransformOut[4] = 0.0;
    1221         417 :     padfGeoTransformOut[5] = -dfPixelSizeY;
    1222             : 
    1223         417 :     CPLFree(padfX);
    1224         417 :     CPLFree(padfXRevert);
    1225         417 :     CPLFree(pabSuccess);
    1226             : 
    1227         417 :     return CE_None;
    1228             : }
    1229             : 
    1230             : /************************************************************************/
    1231             : /*                    GetCurrentCheckWithInvertPROJ()                   */
    1232             : /************************************************************************/
    1233             : 
    1234        3073 : static bool GetCurrentCheckWithInvertPROJ()
    1235             : {
    1236        3073 :     return CPLTestBool(CPLGetConfigOption("CHECK_WITH_INVERT_PROJ", "NO"));
    1237             : }
    1238             : 
    1239             : /************************************************************************/
    1240             : /*               GDALCreateGenImgProjTransformerInternal()              */
    1241             : /************************************************************************/
    1242             : 
    1243             : static void *GDALCreateSimilarGenImgProjTransformer(void *hTransformArg,
    1244             :                                                     double dfRatioX,
    1245             :                                                     double dfRatioY);
    1246             : 
    1247        1900 : static GDALGenImgProjTransformInfo *GDALCreateGenImgProjTransformerInternal()
    1248             : {
    1249             :     /* -------------------------------------------------------------------- */
    1250             :     /*      Initialize the transform info.                                  */
    1251             :     /* -------------------------------------------------------------------- */
    1252             :     GDALGenImgProjTransformInfo *psInfo =
    1253             :         static_cast<GDALGenImgProjTransformInfo *>(
    1254        1900 :             CPLCalloc(sizeof(GDALGenImgProjTransformInfo), 1));
    1255             : 
    1256        1900 :     memcpy(psInfo->sTI.abySignature, GDAL_GTI2_SIGNATURE,
    1257             :            strlen(GDAL_GTI2_SIGNATURE));
    1258        1900 :     psInfo->sTI.pszClassName = GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME;
    1259        1900 :     psInfo->sTI.pfnTransform = GDALGenImgProjTransform;
    1260        1900 :     psInfo->sTI.pfnCleanup = GDALDestroyGenImgProjTransformer;
    1261        1900 :     psInfo->sTI.pfnSerialize = GDALSerializeGenImgProjTransformer;
    1262        1900 :     psInfo->sTI.pfnCreateSimilar = GDALCreateSimilarGenImgProjTransformer;
    1263             : 
    1264        1900 :     psInfo->bCheckWithInvertPROJ = GetCurrentCheckWithInvertPROJ();
    1265        1900 :     psInfo->bHasCustomTransformationPipeline = false;
    1266             : 
    1267        1900 :     return psInfo;
    1268             : }
    1269             : 
    1270             : /************************************************************************/
    1271             : /*                GDALCreateSimilarGenImgProjTransformer()              */
    1272             : /************************************************************************/
    1273             : 
    1274          40 : static void *GDALCreateSimilarGenImgProjTransformer(void *hTransformArg,
    1275             :                                                     double dfRatioX,
    1276             :                                                     double dfRatioY)
    1277             : {
    1278          40 :     VALIDATE_POINTER1(hTransformArg, "GDALCreateSimilarGenImgProjTransformer",
    1279             :                       nullptr);
    1280             : 
    1281          40 :     GDALGenImgProjTransformInfo *psInfo =
    1282             :         static_cast<GDALGenImgProjTransformInfo *>(hTransformArg);
    1283             : 
    1284             :     GDALGenImgProjTransformInfo *psClonedInfo =
    1285          40 :         GDALCreateGenImgProjTransformerInternal();
    1286             : 
    1287          40 :     memcpy(psClonedInfo, psInfo, sizeof(GDALGenImgProjTransformInfo));
    1288             : 
    1289          40 :     psClonedInfo->bCheckWithInvertPROJ = GetCurrentCheckWithInvertPROJ();
    1290             : 
    1291          40 :     if (psClonedInfo->sSrcParams.pTransformArg)
    1292           7 :         psClonedInfo->sSrcParams.pTransformArg = GDALCreateSimilarTransformer(
    1293             :             psInfo->sSrcParams.pTransformArg, dfRatioX, dfRatioY);
    1294          33 :     else if (dfRatioX != 1.0 || dfRatioY != 1.0)
    1295             :     {
    1296          10 :         if (psClonedInfo->sSrcParams.adfGeoTransform[2] == 0.0 &&
    1297          10 :             psClonedInfo->sSrcParams.adfGeoTransform[4] == 0.0)
    1298             :         {
    1299          10 :             psClonedInfo->sSrcParams.adfGeoTransform[1] *= dfRatioX;
    1300          10 :             psClonedInfo->sSrcParams.adfGeoTransform[5] *= dfRatioY;
    1301             :         }
    1302             :         else
    1303             :         {
    1304             :             // If the x and y ratios are not equal, then we cannot really
    1305             :             // compute a geotransform.
    1306           0 :             psClonedInfo->sSrcParams.adfGeoTransform[1] *= dfRatioX;
    1307           0 :             psClonedInfo->sSrcParams.adfGeoTransform[2] *= dfRatioX;
    1308           0 :             psClonedInfo->sSrcParams.adfGeoTransform[4] *= dfRatioX;
    1309           0 :             psClonedInfo->sSrcParams.adfGeoTransform[5] *= dfRatioX;
    1310             :         }
    1311          10 :         if (!GDALInvGeoTransform(psClonedInfo->sSrcParams.adfGeoTransform,
    1312          10 :                                  psClonedInfo->sSrcParams.adfInvGeoTransform))
    1313             :         {
    1314           0 :             CPLError(CE_Failure, CPLE_AppDefined, "Cannot invert geotransform");
    1315           0 :             GDALDestroyGenImgProjTransformer(psClonedInfo);
    1316           0 :             return nullptr;
    1317             :         }
    1318             :     }
    1319             : 
    1320          40 :     if (psClonedInfo->pReprojectArg)
    1321          19 :         psClonedInfo->pReprojectArg =
    1322          19 :             GDALCloneTransformer(psInfo->pReprojectArg);
    1323             : 
    1324          40 :     if (psClonedInfo->sDstParams.pTransformArg)
    1325           0 :         psClonedInfo->sDstParams.pTransformArg =
    1326           0 :             GDALCloneTransformer(psInfo->sDstParams.pTransformArg);
    1327             : 
    1328          40 :     return psClonedInfo;
    1329             : }
    1330             : 
    1331             : /************************************************************************/
    1332             : /*                  GDALCreateGenImgProjTransformer()                   */
    1333             : /************************************************************************/
    1334             : 
    1335             : /**
    1336             :  * Create image to image transformer.
    1337             :  *
    1338             :  * This function creates a transformation object that maps from pixel/line
    1339             :  * coordinates on one image to pixel/line coordinates on another image.  The
    1340             :  * images may potentially be georeferenced in different coordinate systems,
    1341             :  * and may used GCPs to map between their pixel/line coordinates and
    1342             :  * georeferenced coordinates (as opposed to the default assumption that their
    1343             :  * geotransform should be used).
    1344             :  *
    1345             :  * This transformer potentially performs three concatenated transformations.
    1346             :  *
    1347             :  * The first stage is from source image pixel/line coordinates to source
    1348             :  * image georeferenced coordinates, and may be done using the geotransform,
    1349             :  * or if not defined using a polynomial model derived from GCPs.  If GCPs
    1350             :  * are used this stage is accomplished using GDALGCPTransform().
    1351             :  *
    1352             :  * The second stage is to change projections from the source coordinate system
    1353             :  * to the destination coordinate system, assuming they differ.  This is
    1354             :  * accomplished internally using GDALReprojectionTransform().
    1355             :  *
    1356             :  * The third stage is converting from destination image georeferenced
    1357             :  * coordinates to destination image coordinates.  This is done using the
    1358             :  * destination image geotransform, or if not available, using a polynomial
    1359             :  * model derived from GCPs. If GCPs are used this stage is accomplished using
    1360             :  * GDALGCPTransform().  This stage is skipped if hDstDS is NULL when the
    1361             :  * transformation is created.
    1362             :  *
    1363             :  * @param hSrcDS source dataset, or NULL.
    1364             :  * @param pszSrcWKT the coordinate system for the source dataset.  If NULL,
    1365             :  * it will be read from the dataset itself.
    1366             :  * @param hDstDS destination dataset (or NULL).
    1367             :  * @param pszDstWKT the coordinate system for the destination dataset.  If
    1368             :  * NULL, and hDstDS not NULL, it will be read from the destination dataset.
    1369             :  * @param bGCPUseOK TRUE if GCPs should be used if the geotransform is not
    1370             :  * available on the source dataset (not destination).
    1371             :  * @param dfGCPErrorThreshold ignored/deprecated.
    1372             :  * @param nOrder the maximum order to use for GCP derived polynomials if
    1373             :  * possible.  Use 0 to autoselect, or -1 for thin plate splines.
    1374             :  *
    1375             :  * @return handle suitable for use GDALGenImgProjTransform(), and to be
    1376             :  * deallocated with GDALDestroyGenImgProjTransformer().
    1377             :  */
    1378             : 
    1379          62 : void *GDALCreateGenImgProjTransformer(GDALDatasetH hSrcDS,
    1380             :                                       const char *pszSrcWKT,
    1381             :                                       GDALDatasetH hDstDS,
    1382             :                                       const char *pszDstWKT, int bGCPUseOK,
    1383             :                                       CPL_UNUSED double dfGCPErrorThreshold,
    1384             :                                       int nOrder)
    1385             : {
    1386          62 :     char **papszOptions = nullptr;
    1387             : 
    1388          62 :     if (pszSrcWKT != nullptr)
    1389           3 :         papszOptions = CSLSetNameValue(papszOptions, "SRC_SRS", pszSrcWKT);
    1390          62 :     if (pszDstWKT != nullptr)
    1391           3 :         papszOptions = CSLSetNameValue(papszOptions, "DST_SRS", pszDstWKT);
    1392          62 :     if (!bGCPUseOK)
    1393           0 :         papszOptions = CSLSetNameValue(papszOptions, "GCPS_OK", "FALSE");
    1394          62 :     if (nOrder != 0)
    1395           0 :         papszOptions = CSLSetNameValue(papszOptions, "MAX_GCP_ORDER",
    1396           0 :                                        CPLString().Printf("%d", nOrder));
    1397             : 
    1398          62 :     void *pRet = GDALCreateGenImgProjTransformer2(hSrcDS, hDstDS, papszOptions);
    1399          62 :     CSLDestroy(papszOptions);
    1400             : 
    1401          62 :     return pRet;
    1402             : }
    1403             : 
    1404             : /************************************************************************/
    1405             : /*                          InsertCenterLong()                          */
    1406             : /*                                                                      */
    1407             : /*      Insert a CENTER_LONG Extension entry on a GEOGCS to indicate    */
    1408             : /*      the center longitude of the dataset for wrapping purposes.      */
    1409             : /************************************************************************/
    1410             : 
    1411         863 : static void InsertCenterLong(GDALDatasetH hDS, OGRSpatialReference *poSRS,
    1412             :                              CPLStringList &aosOptions)
    1413             : 
    1414             : {
    1415        1488 :     if (!poSRS->IsGeographic() || std::fabs(poSRS->GetAngularUnits() -
    1416         625 :                                             CPLAtof(SRS_UA_DEGREE_CONV)) > 1e-9)
    1417             :     {
    1418         239 :         return;
    1419             :     }
    1420             : 
    1421         624 :     if (poSRS->GetExtension(nullptr, "CENTER_LONG"))
    1422           0 :         return;
    1423             : 
    1424             :     /* -------------------------------------------------------------------- */
    1425             :     /*      For now we only do this if we have a geotransform since         */
    1426             :     /*      other forms require a bunch of extra work.                      */
    1427             :     /* -------------------------------------------------------------------- */
    1428         624 :     double adfGeoTransform[6] = {};
    1429             : 
    1430         624 :     if (GDALGetGeoTransform(hDS, adfGeoTransform) != CE_None)
    1431           0 :         return;
    1432             : 
    1433             :     /* -------------------------------------------------------------------- */
    1434             :     /*      Compute min/max longitude based on testing the four corners.    */
    1435             :     /* -------------------------------------------------------------------- */
    1436         624 :     const int nXSize = GDALGetRasterXSize(hDS);
    1437         624 :     const int nYSize = GDALGetRasterYSize(hDS);
    1438             : 
    1439             :     const double dfMinLong =
    1440        1248 :         std::min(std::min(adfGeoTransform[0] + 0 * adfGeoTransform[1] +
    1441         624 :                               0 * adfGeoTransform[2],
    1442        1248 :                           adfGeoTransform[0] + nXSize * adfGeoTransform[1] +
    1443         624 :                               0 * adfGeoTransform[2]),
    1444        1248 :                  std::min(adfGeoTransform[0] + 0 * adfGeoTransform[1] +
    1445         624 :                               nYSize * adfGeoTransform[2],
    1446        1248 :                           adfGeoTransform[0] + nXSize * adfGeoTransform[1] +
    1447         624 :                               nYSize * adfGeoTransform[2]));
    1448             :     const double dfMaxLong =
    1449        1248 :         std::max(std::max(adfGeoTransform[0] + 0 * adfGeoTransform[1] +
    1450         624 :                               0 * adfGeoTransform[2],
    1451        1248 :                           adfGeoTransform[0] + nXSize * adfGeoTransform[1] +
    1452         624 :                               0 * adfGeoTransform[2]),
    1453        1248 :                  std::max(adfGeoTransform[0] + 0 * adfGeoTransform[1] +
    1454         624 :                               nYSize * adfGeoTransform[2],
    1455        1248 :                           adfGeoTransform[0] + nXSize * adfGeoTransform[1] +
    1456         624 :                               nYSize * adfGeoTransform[2]));
    1457             : 
    1458             :     const double dfEpsilon =
    1459         624 :         std::max(std::fabs(adfGeoTransform[1]), std::fabs(adfGeoTransform[2]));
    1460             :     // If the raster covers more than 360 degree (allow an extra pixel),
    1461             :     // give up
    1462         624 :     constexpr double RELATIVE_EPSILON = 0.05;  // for numeric precision issues
    1463         624 :     if (dfMaxLong - dfMinLong > 360.0 + dfEpsilon * (1 + RELATIVE_EPSILON))
    1464           0 :         return;
    1465             : 
    1466             :     /* -------------------------------------------------------------------- */
    1467             :     /*      Insert center long.                                             */
    1468             :     /* -------------------------------------------------------------------- */
    1469         624 :     const double dfCenterLong = (dfMaxLong + dfMinLong) / 2.0;
    1470         624 :     aosOptions.SetNameValue("CENTER_LONG", CPLSPrintf("%g", dfCenterLong));
    1471             : }
    1472             : 
    1473             : /************************************************************************/
    1474             : /*                      GDALComputeAreaOfInterest()                     */
    1475             : /************************************************************************/
    1476             : 
    1477        1193 : bool GDALComputeAreaOfInterest(const OGRSpatialReference *poSRS,
    1478             :                                double adfGT[6], int nXSize, int nYSize,
    1479             :                                double &dfWestLongitudeDeg,
    1480             :                                double &dfSouthLatitudeDeg,
    1481             :                                double &dfEastLongitudeDeg,
    1482             :                                double &dfNorthLatitudeDeg)
    1483             : {
    1484        1193 :     bool ret = false;
    1485             : 
    1486        1193 :     if (!poSRS)
    1487           0 :         return false;
    1488             : 
    1489        1193 :     OGRSpatialReference oSrcSRSHoriz(*poSRS);
    1490        1193 :     if (oSrcSRSHoriz.IsCompound())
    1491             :     {
    1492          17 :         oSrcSRSHoriz.StripVertical();
    1493             :     }
    1494             : 
    1495        1193 :     OGRSpatialReference *poGeog = oSrcSRSHoriz.CloneGeogCS();
    1496        1193 :     if (poGeog)
    1497             :     {
    1498        1193 :         poGeog->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
    1499        1193 :         poGeog->SetAngularUnits(SRS_UA_DEGREE, CPLAtof(SRS_UA_DEGREE_CONV));
    1500             : 
    1501        1193 :         auto poCT = OGRCreateCoordinateTransformation(&oSrcSRSHoriz, poGeog);
    1502        1193 :         if (poCT)
    1503             :         {
    1504        1193 :             poCT->SetEmitErrors(false);
    1505             : 
    1506             :             double x[4], y[4];
    1507        1193 :             x[0] = adfGT[0];
    1508        1193 :             y[0] = adfGT[3];
    1509        1193 :             x[1] = adfGT[0] + nXSize * adfGT[1];
    1510        1193 :             y[1] = adfGT[3];
    1511        1193 :             x[2] = adfGT[0];
    1512        1193 :             y[2] = adfGT[3] + nYSize * adfGT[5];
    1513        1193 :             x[3] = x[1];
    1514        1193 :             y[3] = y[2];
    1515        1193 :             int validity[4] = {false, false, false, false};
    1516        1193 :             poCT->Transform(4, x, y, nullptr, validity);
    1517        1193 :             dfWestLongitudeDeg = std::numeric_limits<double>::max();
    1518        1193 :             dfSouthLatitudeDeg = std::numeric_limits<double>::max();
    1519        1193 :             dfEastLongitudeDeg = -std::numeric_limits<double>::max();
    1520        1193 :             dfNorthLatitudeDeg = -std::numeric_limits<double>::max();
    1521        5965 :             for (int i = 0; i < 4; i++)
    1522             :             {
    1523        4772 :                 if (validity[i])
    1524             :                 {
    1525        4756 :                     ret = true;
    1526        4756 :                     dfWestLongitudeDeg = std::min(dfWestLongitudeDeg, x[i]);
    1527        4756 :                     dfSouthLatitudeDeg = std::min(dfSouthLatitudeDeg, y[i]);
    1528        4756 :                     dfEastLongitudeDeg = std::max(dfEastLongitudeDeg, x[i]);
    1529        4756 :                     dfNorthLatitudeDeg = std::max(dfNorthLatitudeDeg, y[i]);
    1530             :                 }
    1531             :             }
    1532        1193 :             if (validity[0] && validity[1] && x[0] > x[1])
    1533             :             {
    1534           3 :                 dfWestLongitudeDeg = x[0];
    1535           3 :                 dfEastLongitudeDeg = x[1];
    1536             :             }
    1537        1193 :             if (ret && std::fabs(dfWestLongitudeDeg) <= 180 &&
    1538        1188 :                 std::fabs(dfEastLongitudeDeg) <= 180 &&
    1539        1184 :                 std::fabs(dfSouthLatitudeDeg) <= 90 &&
    1540        1180 :                 std::fabs(dfNorthLatitudeDeg) <= 90)
    1541             :             {
    1542        1180 :                 CPLDebug("GDAL", "Computing area of interest: %g, %g, %g, %g",
    1543             :                          dfWestLongitudeDeg, dfSouthLatitudeDeg,
    1544             :                          dfEastLongitudeDeg, dfNorthLatitudeDeg);
    1545             :             }
    1546             :             else
    1547             :             {
    1548          13 :                 CPLDebug("GDAL", "Could not compute area of interest");
    1549          13 :                 dfWestLongitudeDeg = 0;
    1550          13 :                 dfSouthLatitudeDeg = 0;
    1551          13 :                 dfEastLongitudeDeg = 0;
    1552          13 :                 dfNorthLatitudeDeg = 0;
    1553             :             }
    1554        1193 :             OGRCoordinateTransformation::DestroyCT(poCT);
    1555             :         }
    1556             : 
    1557        1193 :         delete poGeog;
    1558             :     }
    1559             : 
    1560        1193 :     return ret;
    1561             : }
    1562             : 
    1563           5 : bool GDALComputeAreaOfInterest(const OGRSpatialReference *poSRS, double dfX1,
    1564             :                                double dfY1, double dfX2, double dfY2,
    1565             :                                double &dfWestLongitudeDeg,
    1566             :                                double &dfSouthLatitudeDeg,
    1567             :                                double &dfEastLongitudeDeg,
    1568             :                                double &dfNorthLatitudeDeg)
    1569             : {
    1570           5 :     bool ret = false;
    1571             : 
    1572           5 :     if (!poSRS)
    1573           0 :         return false;
    1574             : 
    1575           5 :     OGRSpatialReference oSrcSRSHoriz(*poSRS);
    1576           5 :     if (oSrcSRSHoriz.IsCompound())
    1577             :     {
    1578           0 :         oSrcSRSHoriz.StripVertical();
    1579             :     }
    1580             : 
    1581           5 :     OGRSpatialReference *poGeog = oSrcSRSHoriz.CloneGeogCS();
    1582           5 :     if (poGeog)
    1583             :     {
    1584           5 :         poGeog->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
    1585             : 
    1586           5 :         auto poCT = OGRCreateCoordinateTransformation(&oSrcSRSHoriz, poGeog);
    1587           5 :         if (poCT)
    1588             :         {
    1589             :             double x[4], y[4];
    1590           5 :             x[0] = dfX1;
    1591           5 :             y[0] = dfY1;
    1592           5 :             x[1] = dfX2;
    1593           5 :             y[1] = dfY1;
    1594           5 :             x[2] = dfX1;
    1595           5 :             y[2] = dfY2;
    1596           5 :             x[3] = dfX2;
    1597           5 :             y[3] = dfY2;
    1598           5 :             int validity[4] = {false, false, false, false};
    1599           5 :             poCT->Transform(4, x, y, nullptr, validity);
    1600           5 :             dfWestLongitudeDeg = std::numeric_limits<double>::max();
    1601           5 :             dfSouthLatitudeDeg = std::numeric_limits<double>::max();
    1602           5 :             dfEastLongitudeDeg = -std::numeric_limits<double>::max();
    1603           5 :             dfNorthLatitudeDeg = -std::numeric_limits<double>::max();
    1604          25 :             for (int i = 0; i < 4; i++)
    1605             :             {
    1606          20 :                 if (validity[i])
    1607             :                 {
    1608          20 :                     ret = true;
    1609          20 :                     dfWestLongitudeDeg = std::min(dfWestLongitudeDeg, x[i]);
    1610          20 :                     dfSouthLatitudeDeg = std::min(dfSouthLatitudeDeg, y[i]);
    1611          20 :                     dfEastLongitudeDeg = std::max(dfEastLongitudeDeg, x[i]);
    1612          20 :                     dfNorthLatitudeDeg = std::max(dfNorthLatitudeDeg, y[i]);
    1613             :                 }
    1614             :             }
    1615           5 :             if (validity[0] && validity[1] && (dfX1 - dfX2) * (x[0] - x[1]) < 0)
    1616             :             {
    1617           0 :                 dfWestLongitudeDeg = x[0];
    1618           0 :                 dfEastLongitudeDeg = x[1];
    1619             :             }
    1620           5 :             if (ret)
    1621             :             {
    1622           5 :                 CPLDebug("GDAL", "Computing area of interest: %g, %g, %g, %g",
    1623             :                          dfWestLongitudeDeg, dfSouthLatitudeDeg,
    1624             :                          dfEastLongitudeDeg, dfNorthLatitudeDeg);
    1625             :             }
    1626             :             else
    1627             :             {
    1628           0 :                 CPLDebug("GDAL", "Could not compute area of interest");
    1629           0 :                 dfWestLongitudeDeg = 0;
    1630           0 :                 dfSouthLatitudeDeg = 0;
    1631           0 :                 dfEastLongitudeDeg = 0;
    1632           0 :                 dfNorthLatitudeDeg = 0;
    1633             :             }
    1634           5 :             delete poCT;
    1635             :         }
    1636             : 
    1637           5 :         delete poGeog;
    1638             :     }
    1639             : 
    1640           5 :     return ret;
    1641             : }
    1642             : 
    1643             : /************************************************************************/
    1644             : /*                    GDALGCPAntimeridianUnwrap()                       */
    1645             : /************************************************************************/
    1646             : 
    1647             : /* Deal with discontinuties of dfGCPX longitudes around the anti-meridian.
    1648             :  * Cf https://github.com/OSGeo/gdal/issues/8371
    1649             :  */
    1650          41 : static void GDALGCPAntimeridianUnwrap(int nGCPCount, GDAL_GCP *pasGCPList,
    1651             :                                       const OGRSpatialReference &oSRS,
    1652             :                                       CSLConstList papszOptions)
    1653             : {
    1654             :     const char *pszGCPAntimeridianUnwrap =
    1655          41 :         CSLFetchNameValueDef(papszOptions, "GCP_ANTIMERIDIAN_UNWRAP", "AUTO");
    1656         122 :     const bool bForced = EQUAL(pszGCPAntimeridianUnwrap, "YES") ||
    1657          40 :                          EQUAL(pszGCPAntimeridianUnwrap, "ON") ||
    1658         121 :                          EQUAL(pszGCPAntimeridianUnwrap, "TRUE") ||
    1659          40 :                          EQUAL(pszGCPAntimeridianUnwrap, "1");
    1660          47 :     if (bForced || (!oSRS.IsEmpty() && oSRS.IsGeographic() &&
    1661           6 :                     fabs(oSRS.GetAngularUnits(nullptr) -
    1662           6 :                          CPLAtof(SRS_UA_DEGREE_CONV)) < 1e-8 &&
    1663           6 :                     EQUAL(pszGCPAntimeridianUnwrap, "AUTO")))
    1664             :     {
    1665           6 :         if (!bForced)
    1666             :         {
    1667             :             // Proceed to unwrapping only if the longitudes are within
    1668             :             // [-180, -170] or [170, 180]
    1669         425 :             for (int i = 0; i < nGCPCount; ++i)
    1670             :             {
    1671         423 :                 const double dfLongAbs = fabs(pasGCPList[i].dfGCPX);
    1672         423 :                 if (dfLongAbs > 180 || dfLongAbs < 170)
    1673             :                 {
    1674           3 :                     return;
    1675             :                 }
    1676             :             }
    1677             :         }
    1678             : 
    1679           3 :         bool bDone = false;
    1680         633 :         for (int i = 0; i < nGCPCount; ++i)
    1681             :         {
    1682         630 :             if (pasGCPList[i].dfGCPX < 0)
    1683             :             {
    1684          48 :                 if (!bDone)
    1685             :                 {
    1686           3 :                     bDone = true;
    1687           3 :                     CPLDebug("WARP", "GCP longitude unwrapping");
    1688             :                 }
    1689          48 :                 pasGCPList[i].dfGCPX += 360;
    1690             :             }
    1691             :         }
    1692             :     }
    1693             : }
    1694             : 
    1695             : /************************************************************************/
    1696             : /*              GDALGetGenImgProjTranformerOptionList()                 */
    1697             : /************************************************************************/
    1698             : 
    1699             : /** Return a XML string describing options accepted by
    1700             :  * GDALCreateGenImgProjTransformer2().
    1701             :  *
    1702             :  * @since 3.11
    1703             :  */
    1704        1646 : const char *GDALGetGenImgProjTranformerOptionList(void)
    1705             : {
    1706             :     return "<OptionList>"
    1707             :            "<Option name='SRC_SRS' type='string' description='WKT SRS, or any "
    1708             :            "string recognized by OGRSpatialReference::SetFromUserInput(), to "
    1709             :            "be used as an override for CRS of input dataset'/>"
    1710             :            "<Option name='DST_SRS' type='string' description='WKT SRS, or any "
    1711             :            "string recognized by OGRSpatialReference::SetFromUserInput(), to "
    1712             :            "be used as an override for CRS of output dataset'/>"
    1713             :            "<Option name='PROMOTE_TO_3D' type='boolean' description='"
    1714             :            "Whether to promote SRC_SRS / DST_SRS to 3D.' "
    1715             :            "default='NO'/>"
    1716             :            "<Option name='COORDINATE_OPERATION' type='string' description='"
    1717             :            "Coordinate operation, as a PROJ or WKT string, used as an override "
    1718             :            "over the normally computed pipeline. The pipeline must take into "
    1719             :            "account the axis order of the source and target SRS.'/>"
    1720             :            "<Option name='ALLOW_BALLPARK' type='boolean' description='"
    1721             :            "Whether ballpark coordinate operations are allowed.' "
    1722             :            "default='YES'/>"
    1723             :            "<Option name='ONLY_BEST' type='string-select' "
    1724             :            "description='"
    1725             :            "By default (at least in the PROJ 9.x series), PROJ may use "
    1726             :            "coordinate operations that are not the \"best\" if resources "
    1727             :            "(typically grids) needed to use them are missing. It will then "
    1728             :            "fallback to other coordinate operations that have a lesser "
    1729             :            "accuracy, for example using Helmert transformations, or in the "
    1730             :            "absence of such operations, to ones with potential very rough "
    1731             :            " accuracy, using \"ballpark\" transformations (see "
    1732             :            "https://proj.org/glossary.html). "
    1733             :            "When calling this method with YES, PROJ will only consider the "
    1734             :            "\"best\" operation, and error out (at Transform() time) if they "
    1735             :            "cannot be used. This method may be used together with "
    1736             :            "ALLOW_BALLPARK=NO to only allow best operations that have a known "
    1737             :            "accuracy. Note that this method has no effect on PROJ versions "
    1738             :            "before 9.2. The default value for this option can be also set with "
    1739             :            "the PROJ_ONLY_BEST_DEFAULT environment variable, or with the "
    1740             :            "\"only_best_default\" setting of proj.ini. Setting "
    1741             :            "ONLY_BEST=YES/NO overrides such default value' default='AUTO'>"
    1742             :            "  <Value>AUTO</Value>"
    1743             :            "  <Value>YES</Value>"
    1744             :            "  <Value>NO</Value>"
    1745             :            "</Option>"
    1746             :            "<Option name='COORDINATE_EPOCH' type='float' description='"
    1747             :            "Coordinate epoch, expressed as a decimal year. Useful for "
    1748             :            "time-dependent coordinate operations.'/>"
    1749             :            "<Option name='SRC_COORDINATE_EPOCH' type='float' description='"
    1750             :            "Coordinate epoch of source CRS, expressed as a decimal year. "
    1751             :            "Useful for time-dependent coordinate operations.'/>"
    1752             :            "<Option name='DST_COORDINATE_EPOCH' type='float' description='"
    1753             :            "Coordinate epoch of target CRS, expressed as a decimal year. "
    1754             :            "Useful for time-dependent coordinate operations.'/>"
    1755             :            "<Option name='GCPS_OK' type='boolean' description='"
    1756             :            "Allow use of GCPs.' default='YES'/>"
    1757             :            "<Option name='REFINE_MINIMUM_GCPS' type='int' description='"
    1758             :            "The minimum amount of GCPs that should be available after the "
    1759             :            "refinement'/>"
    1760             :            "<Option name='REFINE_TOLERANCE' type='float' description='"
    1761             :            "The tolerance that specifies when a GCP will be eliminated.'/>"
    1762             :            "<Option name='MAX_GCP_ORDER' type='int' description='"
    1763             :            "The maximum order to use for GCP derived polynomials if possible. "
    1764             :            "The default is to autoselect based on the number of GCPs. A value "
    1765             :            "of -1 triggers use of Thin Plate Spline instead of polynomials.'/>"
    1766             :            "<Option name='GCP_ANTIMERIDIAN_UNWRAP' type='string-select' "
    1767             :            "description='"
    1768             :            "Whether to \"unwrap\" longitudes of ground control points that "
    1769             :            "span the antimeridian. For datasets with GCPs in "
    1770             :            "longitude/latitude coordinate space spanning the antimeridian, "
    1771             :            "longitudes will have a discontinuity on +/- 180 deg, and will "
    1772             :            "result in a subset of the GCPs with longitude in the [-180,-170] "
    1773             :            "range and another subset in [170, 180]. By default (AUTO), that "
    1774             :            "situation will be detected and longitudes in [-180,-170] will be "
    1775             :            "shifted to [180, 190] to get a continuous set. This option can be "
    1776             :            "set to YES to force that behavior (useful if no SRS information is "
    1777             :            "available), or to NO to disable it.' default='AUTO'>"
    1778             :            "  <Value>AUTO</Value>"
    1779             :            "  <Value>YES</Value>"
    1780             :            "  <Value>NO</Value>"
    1781             :            "</Option>"
    1782             :            "<Option name='SRC_METHOD' alias='METHOD' type='string-select' "
    1783             :            "description='"
    1784             :            "Force only one geolocation method to be considered on the source "
    1785             :            "dataset. Will be used for pixel/line to georef transformation on "
    1786             :            "the source dataset. NO_GEOTRANSFORM can be used to specify the "
    1787             :            "identity geotransform (ungeoreferenced image)'>"
    1788             :            "  <Value>GEOTRANSFORM</Value>"
    1789             :            "  <Value>GCP_POLYNOMIAL</Value>"
    1790             :            "  <Value>GCP_TPS</Value>"
    1791             :            "  <Value>GCP_HOMOGRAPHY</Value>"
    1792             :            "  <Value>GEOLOC_ARRAY</Value>"
    1793             :            "  <Value>RPC</Value>"
    1794             :            "  <Value>NO_GEOTRANSFORM</Value>"
    1795             :            "</Option>"
    1796             :            "<Option name='DST_METHOD' type='string-select' description='"
    1797             :            "Force only one geolocation method to be considered on the target "
    1798             :            "dataset. Will be used for pixel/line to georef transformation on "
    1799             :            "the targe dataset. NO_GEOTRANSFORM can be used to specify the "
    1800             :            "identity geotransform (ungeoreferenced image)'>"
    1801             :            "  <Value>GEOTRANSFORM</Value>"
    1802             :            "  <Value>GCP_POLYNOMIAL</Value>"
    1803             :            "  <Value>GCP_TPS</Value>"
    1804             :            "  <Value>GCP_HOMOGRAPHY</Value>"
    1805             :            "  <Value>GEOLOC_ARRAY</Value>"
    1806             :            "  <Value>RPC</Value>"
    1807             :            "  <Value>NO_GEOTRANSFORM</Value>"
    1808             :            "</Option>"
    1809             :            "<Option name='RPC_HEIGHT' type='float' description='"
    1810             :            "A fixed height to be used with RPC calculations. If RPC_HEIGHT and "
    1811             :            "RPC_DEM are not specified but that the RPC metadata domain contains"
    1812             :            " a HEIGHT_DEFAULT item (for example, the DIMAP driver may fill it),"
    1813             :            "this value will be used as the RPC_HEIGHT. Otherwise, if none of "
    1814             :            "RPC_HEIGHT and RPC_DEM are specified as transformer options and "
    1815             :            "if HEIGHT_DEFAULT is no available, a height of 0 will be used.'/>"
    1816             :            "<Option name='RPC_DEM' type='string' description='"
    1817             :            "Name of a GDAL dataset (a DEM file typically) used to extract "
    1818             :            "elevation offsets from. In this situation the Z passed into the "
    1819             :            "transformation function is assumed to be height above ground. "
    1820             :            "This option should be used in replacement of RPC_HEIGHT to provide "
    1821             :            "a way of defining a non uniform ground for the target scene.'/>"
    1822             :            "<Option name='RPC_HEIGHT_SCALE' type='float' description='"
    1823             :            "Factor used to multiply heights above ground. Useful when "
    1824             :            "elevation offsets of the DEM are not expressed in meters.'/>"
    1825             :            "<Option name='RPC_DEMINTERPOLATION' type='string-select' "
    1826             :            "description='DEM interpolation method' default='BILINEAR'>"
    1827             :            "  <Value>NEAR</Value>"
    1828             :            "  <Value>BILINEAR</Value>"
    1829             :            "  <Value>CUBIC</Value>"
    1830             :            "</Option>"
    1831             :            "<Option name='RPC_DEM_MISSING_VALUE' type='float' description='"
    1832             :            "Value of DEM height that must be used in case the DEM has nodata "
    1833             :            "value at the sampling point, or if its extent does not cover the "
    1834             :            "requested coordinate. When not specified, missing values will "
    1835             :            "cause a failed transform.'/>"
    1836             :            "<Option name='RPC_DEM_SRS' type='string' description='"
    1837             :            "WKT SRS, or any string recognized by "
    1838             :            "OGRSpatialReference::SetFromUserInput(), to be used as an "
    1839             :            "override for DEM SRS. Useful if DEM SRS does not have an explicit "
    1840             :            "vertical component.'/>"
    1841             :            "<Option name='RPC_DEM_APPLY_VDATUM_SHIFT' type='boolean' "
    1842             :            "description='"
    1843             :            "Whether the vertical component of a compound SRS for the DEM "
    1844             :            "should be used (when it is present). This is useful so as to "
    1845             :            "be able to transform the raw values from the DEM expressed with "
    1846             :            "respect to a geoid to the heights with respect to the WGS84 "
    1847             :            "ellipsoid. When this is enabled, the GTIFF_REPORT_COMPD_CS "
    1848             :            "configuration option will be also set temporarily so as to get "
    1849             :            "the vertical information from GeoTIFF files.' default='YES'/>"
    1850             :            "<Option name='RPC_PIXEL_ERROR_THRESHOLD' type='float' description='"
    1851             :            "Overrides the dfPixErrThreshold parameter, i.e. the error "
    1852             :            "(measured in pixels) allowed in the iterative solution of "
    1853             :            "pixel/line to lat/long computations (the other way is always "
    1854             :            "exact given the equations).'/>"
    1855             :            "<Option name='RPC_MAX_ITERATIONS' type='int' description='"
    1856             :            "Maximum number of iterations allowed in the iterative solution of "
    1857             :            "pixel/line to lat/long computations. Default value is 10 in the "
    1858             :            "absence of a DEM, or 20 if there is a DEM.'/>"
    1859             :            "<Option name='RPC_FOOTPRINT' type='string' description='"
    1860             :            "WKT or GeoJSON polygon (in long / lat coordinate space) with a "
    1861             :            "validity footprint for the RPC. Any coordinate transformation that "
    1862             :            "goes from or arrive outside this footprint will be considered "
    1863             :            "invalid. This* is useful in situations where the RPC values become "
    1864             :            "highly unstable outside of the area on which they have been "
    1865             :            "computed for, potentially leading to undesirable \"echoes\" / "
    1866             :            "false positives. This requires GDAL to be built against GEOS..'/>"
    1867             :            "<Option name='RPC_MAX_ITERATIONS' type='int' description='"
    1868             :            "Maximum number of iterations allowed in the iterative solution of "
    1869             :            "pixel/line to lat/long computations. Default value is 10 in the "
    1870             :            "absence of a DEM, or 20 if there is a DEM.'/>"
    1871             :            "<Option name='INSERT_CENTER_LONG' type='boolean' description='"
    1872             :            "May be set to FALSE to disable setting up a CENTER_LONG value on "
    1873             :            "the coordinate system to rewrap things around the center of the "
    1874             :            "image.' default='YES'/>"
    1875             :            "<Option name='SRC_APPROX_ERROR_IN_SRS_UNIT' type='float' "
    1876             :            "description='"
    1877             :            "Use an approximate transformer for the source transformer. Must be "
    1878             :            "defined together with SRC_APPROX_ERROR_IN_PIXEL to be taken into "
    1879             :            "account.'/>"
    1880             :            "<Option name='SRC_APPROX_ERROR_IN_PIXEL' type='float' "
    1881             :            "description='"
    1882             :            "Use an approximate transformer for the source transformer. Must be "
    1883             :            "defined together with SRC_APPROX_ERROR_IN_SRS_UNIT to be taken "
    1884             :            "into "
    1885             :            "account.'/>"
    1886             :            "<Option name='DST_APPROX_ERROR_IN_SRS_UNIT' type='float' "
    1887             :            "description='"
    1888             :            "Use an approximate transformer for the target transformer. Must be "
    1889             :            "defined together with DST_APPROX_ERROR_IN_PIXEL to be taken into "
    1890             :            "account.'/>"
    1891             :            "<Option name='DST_APPROX_ERROR_IN_PIXEL' type='float' "
    1892             :            "description='"
    1893             :            "Use an approximate transformer for the target transformer. Must be "
    1894             :            "defined together with DST_APPROX_ERROR_IN_SRS_UNIT to be taken "
    1895             :            "into "
    1896             :            "account.'/>"
    1897             :            "<Option name='REPROJECTION_APPROX_ERROR_IN_SRC_SRS_UNIT' "
    1898             :            "type='float' "
    1899             :            "description='"
    1900             :            "Use an approximate transformer for the coordinate reprojection. "
    1901             :            "Must be used together with "
    1902             :            "REPROJECTION_APPROX_ERROR_IN_DST_SRS_UNIT to be taken into "
    1903             :            "account.'/>"
    1904             :            "<Option name='REPROJECTION_APPROX_ERROR_IN_DST_SRS_UNIT' "
    1905             :            "type='float' "
    1906             :            "description='"
    1907             :            "Use an approximate transformer for the coordinate reprojection. "
    1908             :            "Must be used together with "
    1909             :            "REPROJECTION_APPROX_ERROR_IN_SRC_SRS_UNIT to be taken into "
    1910             :            "account.'/>"
    1911             :            "<Option name='AREA_OF_INTEREST' type='string' "
    1912             :            "description='"
    1913             :            "Area of interest, as "
    1914             :            "west_lon_deg,south_lat_deg,east_lon_deg,north_lat_deg, used to "
    1915             :            "compute the best coordinate operation between the source and "
    1916             :            "target SRS. If not specified, the bounding box of the source "
    1917             :            "raster will be used.'/>"
    1918             :            "<Option name='GEOLOC_BACKMAP_OVERSAMPLE_FACTOR' type='float' "
    1919             :            "min='0.1' max='2' description='"
    1920             :            "Oversample factor used to derive the size of the \"backmap\" used "
    1921             :            "for geolocation array transformers.' default='1.3'/>"
    1922             :            "<Option name='GEOLOC_USE_TEMP_DATASETS' type='boolean' "
    1923             :            "description='"
    1924             :            "Whether temporary GeoTIFF datasets should be used to store the "
    1925             :            "backmap. The default is NO, that is to use in-memory arrays, "
    1926             :            "unless the number of pixels of the geolocation array is greater "
    1927             :            "than 16 megapixels.' default='NO'/>"
    1928             :            "<Option name='GEOLOC_ARRAY' alias='SRC_GEOLOC_ARRAY' type='string' "
    1929             :            "description='"
    1930             :            "Name of a GDAL dataset containing a geolocation array and "
    1931             :            "associated metadata. This is an alternative to having geolocation "
    1932             :            "information described in the GEOLOCATION metadata domain of the "
    1933             :            "source dataset. The dataset specified may have a GEOLOCATION "
    1934             :            "metadata domain containing appropriate metadata, however default "
    1935             :            "values are assigned for all omitted items. X_BAND defaults to 1 "
    1936             :            "and Y_BAND to 2, however the dataset must contain exactly 2 bands. "
    1937             :            "PIXEL_OFFSET and LINE_OFFSET default to 0. PIXEL_STEP and "
    1938             :            "LINE_STEP default to the ratio of the width/height of the source "
    1939             :            "dataset divided by the with/height of the geolocation array. "
    1940             :            "SRS defaults to the spatial reference system of the geolocation "
    1941             :            "array dataset, if set, otherwise WGS84 is used. "
    1942             :            "GEOREFERENCING_CONVENTION is selected from the main metadata "
    1943             :            "domain if it is omitted from the GEOLOCATION domain, and if not "
    1944             :            "available TOP_LEFT_CORNER is assigned as a default. "
    1945             :            "If GEOLOC_ARRAY is set SRC_METHOD defaults to GEOLOC_ARRAY.'/>"
    1946             :            "<Option name='DST_GEOLOC_ARRAY' type='string' "
    1947             :            "description='"
    1948             :            "Name of a GDAL dataset that contains at least 2 bands with the X "
    1949             :            "and Y geolocation bands. This is an alternative to having "
    1950             :            "geolocation information described in the GEOLOCATION metadata "
    1951             :            "domain of the destination dataset. See SRC_GEOLOC_ARRAY "
    1952             :            "description for details, assumptions, and defaults. If this "
    1953             :            "option is set, DST_METHOD=GEOLOC_ARRAY will be assumed if not "
    1954             :            "set.'/>"
    1955             :            "<Option name='GEOLOC_NORMALIZE_LONGITUDE_MINUS_180_PLUS_180' "
    1956             :            "type='boolean' "
    1957             :            "description='"
    1958             :            "Force geolocation longitudes into -180,180 when longitude/latitude "
    1959             :            "is the coordinate system of the geolocation arrays' default='NO'>"
    1960             :            "  <Value>YES</Value>"
    1961             :            "  <Value>NO</Value>"
    1962             :            "</Option>"
    1963             :            "<Option name='NUM_THREADS' type='string' "
    1964             :            "description='Number of threads to use'/>"
    1965        1646 :            "</OptionList>";
    1966             : }
    1967             : 
    1968             : /************************************************************************/
    1969             : /*                  GDALCreateGenImgProjTransformer2()                  */
    1970             : /************************************************************************/
    1971             : 
    1972             : /* clang-format off */
    1973             : /**
    1974             :  * Create image to image transformer.
    1975             :  *
    1976             :  * This function creates a transformation object that maps from pixel/line
    1977             :  * coordinates on one image to pixel/line coordinates on another image.  The
    1978             :  * images may potentially be georeferenced in different coordinate systems,
    1979             :  * and may used GCPs to map between their pixel/line coordinates and
    1980             :  * georeferenced coordinates (as opposed to the default assumption that their
    1981             :  * geotransform should be used).
    1982             :  *
    1983             :  * This transformer potentially performs three concatenated transformations.
    1984             :  *
    1985             :  * The first stage is from source image pixel/line coordinates to source
    1986             :  * image georeferenced coordinates, and may be done using the geotransform,
    1987             :  * or if not defined using a polynomial model derived from GCPs.  If GCPs
    1988             :  * are used this stage is accomplished using GDALGCPTransform().
    1989             :  *
    1990             :  * The second stage is to change projections from the source coordinate system
    1991             :  * to the destination coordinate system, assuming they differ.  This is
    1992             :  * accomplished internally using GDALReprojectionTransform().
    1993             :  *
    1994             :  * The third stage is converting from destination image georeferenced
    1995             :  * coordinates to destination image coordinates.  This is done using the
    1996             :  * destination image geotransform, or if not available, using a polynomial
    1997             :  * model derived from GCPs. If GCPs are used this stage is accomplished using
    1998             :  * GDALGCPTransform().  This stage is skipped if hDstDS is NULL when the
    1999             :  * transformation is created.
    2000             :  *
    2001             :  * Supported Options (specified with the -to switch of gdalwarp for example):
    2002             :  * <ul>
    2003             :  * <li> SRC_SRS: WKT SRS, or any string recognized by
    2004             :  * OGRSpatialReference::SetFromUserInput(), to be used as an override for
    2005             :  * hSrcDS.</li>
    2006             :  * <li> DST_SRS: WKT SRS, or any string recognized by
    2007             :  * OGRSpatialReference::SetFromUserInput(),  to be used as an override for
    2008             :  * hDstDS.
    2009             :  * </li>
    2010             :  * <li>PROMOTE_TO_3D=YES/NO: whether to promote SRC_SRS / DST_SRS to 3D.
    2011             :  * Default is NO</li>
    2012             :  * <li> COORDINATE_OPERATION: (GDAL &gt;= 3.0) Coordinate operation, as
    2013             :  * a PROJ or WKT string, used as an override over the normally computed
    2014             :  * pipeline. The pipeline must take into account the axis order of the source
    2015             :  * and target SRS.
    2016             :  * </li>
    2017             :  * <li> ALLOW_BALLPARK=YES/NO: (GDAL &gt;= 3.11) Whether ballpark coordinate
    2018             :  * operations are allowed. Defaults to YES.</li>
    2019             :  * <li> ONLY_BEST=YES/NO/AUTO: (GDAL &gt;= 3.11) By default (at least in the
    2020             :  * PROJ 9.x series), PROJ may use coordinate
    2021             :  * operations that are not the "best" if resources (typically grids) needed
    2022             :  * to use them are missing. It will then fallback to other coordinate operations
    2023             :  * that have a lesser accuracy, for example using Helmert transformations,
    2024             :  * or in the absence of such operations, to ones with potential very rought
    2025             :  * accuracy, using "ballpark" transformations
    2026             :  * (see https://proj.org/glossary.html).
    2027             :  * When calling this method with YES, PROJ will only consider the
    2028             :  * "best" operation, and error out (at Transform() time) if they cannot be
    2029             :  * used.
    2030             :  * This method may be used together with ALLOW_BALLPARK=NO to
    2031             :  * only allow best operations that have a known accuracy.
    2032             :  * Note that this method has no effect on PROJ versions before 9.2.
    2033             :  * The default value for this option can be also set with the
    2034             :  * PROJ_ONLY_BEST_DEFAULT environment variable, or with the "only_best_default"
    2035             :  * setting of proj.ini. Calling SetOnlyBest() overrides such default value.</li>
    2036             :  * <li> COORDINATE_EPOCH: (GDAL &gt;= 3.0) Coordinate epoch,
    2037             :  * expressed as a decimal year. Useful for time-dependent coordinate operations.
    2038             :  * </li>
    2039             :  * <li> SRC_COORDINATE_EPOCH: (GDAL &gt;= 3.4) Coordinate epoch of source CRS,
    2040             :  * expressed as a decimal year. Useful for time-dependent coordinate operations.
    2041             :  * </li>
    2042             :  * <li> DST_COORDINATE_EPOCH: (GDAL &gt;= 3.4) Coordinate epoch of target CRS,
    2043             :  * expressed as a decimal year. Useful for time-dependent coordinate operations.
    2044             :  * </li>
    2045             :  * <li> GCPS_OK: If false, GCPs will not be used, default is TRUE.
    2046             :  * </li>
    2047             :  * <li> REFINE_MINIMUM_GCPS: The minimum amount of GCPs that should be available
    2048             :  * after the refinement.
    2049             :  * </li>
    2050             :  * <li> REFINE_TOLERANCE: The tolerance that specifies when a GCP will be
    2051             :  * eliminated.
    2052             :  * </li>
    2053             :  * <li> MAX_GCP_ORDER: the maximum order to use for GCP derived polynomials if
    2054             :  * possible.  The default is to autoselect based on the number of GCPs.
    2055             :  * A value of -1 triggers use of Thin Plate Spline instead of polynomials.
    2056             :  * </li>
    2057             :  * <li>GCP_ANTIMERIDIAN_UNWRAP=AUTO/YES/NO. (GDAL &gt;= 3.8) Whether to
    2058             :  * "unwrap" longitudes of ground control points that span the antimeridian.
    2059             :  * For datasets with GCPs in longitude/latitude coordinate space spanning the
    2060             :  * antimeridian, longitudes will have a discontinuity on +/- 180 deg, and
    2061             :  * will result in a subset of the GCPs with longitude in the [-180,-170] range
    2062             :  * and another subset in [170, 180]. By default (AUTO), that situation will be
    2063             :  * detected and longitudes in [-180,-170] will be shifted to [180, 190] to get
    2064             :  * a continuous set. This option can be set to YES to force that behavior
    2065             :  * (useful if no SRS information is available), or to NO to disable it.
    2066             :  * </li>
    2067             :  * <li> SRC_METHOD: may have a value which is one of GEOTRANSFORM, GCP_HOMOGRAPHY,
    2068             :  * GCP_POLYNOMIAL, GCP_TPS, GEOLOC_ARRAY, RPC to force only one geolocation
    2069             :  * method to be considered on the source dataset. Will be used for pixel/line
    2070             :  * to georef transformation on the source dataset. NO_GEOTRANSFORM can be
    2071             :  * used to specify the identity geotransform (ungeoreferenced image)
    2072             :  * </li>
    2073             :  * <li> DST_METHOD: may have a value which is one of GEOTRANSFORM,
    2074             :  * GCP_POLYNOMIAL, GCP_HOMOGRAPHY, GCP_TPS, GEOLOC_ARRAY (added in 3.5), RPC to
    2075             :  * force only one
    2076             :  * geolocation method to be considered on the target dataset.  Will be used for
    2077             :  * pixel/line to georef transformation on the destination dataset.
    2078             :  * NO_GEOTRANSFORM can be used to specify the identity geotransform
    2079             :  * (ungeoreferenced image)
    2080             :  * </li>
    2081             :  * <li> RPC_HEIGHT: A fixed height to be used with RPC
    2082             :  * calculations. If RPC_HEIGHT and RPC_DEM are not specified but that the RPC
    2083             :  * metadata domain contains a HEIGHT_DEFAULT item (for example, the DIMAP driver
    2084             :  * may fill it), this value will be used as the RPC_HEIGHT. Otherwise, if none
    2085             :  * of RPC_HEIGHT and RPC_DEM are specified as transformer
    2086             :  * options and if HEIGHT_DEFAULT is no available, a height of 0 will be used.
    2087             :  * </li>
    2088             :  * <li> RPC_DEM: The name of a DEM file to be used with RPC
    2089             :  * calculations. See GDALCreateRPCTransformerV2() for more details.
    2090             :  * </li>
    2091             :  * <li> Other RPC related options. See GDALCreateRPCTransformerV2()
    2092             :  * </li>
    2093             :  * <li>
    2094             :  * INSERT_CENTER_LONG: May be set to FALSE to disable setting up a CENTER_LONG
    2095             :  * value on the coordinate system to rewrap things around the center of the
    2096             :  * image.
    2097             :  * </li>
    2098             :  * <li> SRC_APPROX_ERROR_IN_SRS_UNIT=err_threshold_in_SRS_units. (GDAL
    2099             :  * &gt;= 2.2) Use an approximate transformer for the source transformer. Must be
    2100             :  * defined together with SRC_APPROX_ERROR_IN_PIXEL to be taken into account.
    2101             :  * </li>
    2102             :  * <li> SRC_APPROX_ERROR_IN_PIXEL=err_threshold_in_pixel. (GDAL &gt;= 2.2) Use
    2103             :  * an approximate transformer for the source transformer.. Must be defined
    2104             :  * together with SRC_APPROX_ERROR_IN_SRS_UNIT to be taken into account.
    2105             :  * </li>
    2106             :  * <li>
    2107             :  * DST_APPROX_ERROR_IN_SRS_UNIT=err_threshold_in_SRS_units. (GDAL &gt;= 2.2) Use
    2108             :  * an approximate transformer for the destination transformer. Must be defined
    2109             :  * together with DST_APPROX_ERROR_IN_PIXEL to be taken into account.
    2110             :  * </li>
    2111             :  * <li>
    2112             :  * DST_APPROX_ERROR_IN_PIXEL=err_threshold_in_pixel. (GDAL &gt;= 2.2) Use an
    2113             :  * approximate transformer for the destination transformer. Must be defined
    2114             :  * together with DST_APPROX_ERROR_IN_SRS_UNIT to be taken into account.
    2115             :  * </li>
    2116             :  * <li>
    2117             :  * REPROJECTION_APPROX_ERROR_IN_SRC_SRS_UNIT=err_threshold_in_src_SRS_units.
    2118             :  * (GDAL &gt;= 2.2) Use an approximate transformer for the coordinate
    2119             :  * reprojection. Must be used together with
    2120             :  * REPROJECTION_APPROX_ERROR_IN_DST_SRS_UNIT to be taken into account.
    2121             :  * </li>
    2122             :  * <li>
    2123             :  * REPROJECTION_APPROX_ERROR_IN_DST_SRS_UNIT=err_threshold_in_dst_SRS_units.
    2124             :  * (GDAL &gt;= 2.2) Use an approximate transformer for the coordinate
    2125             :  * reprojection. Must be used together with
    2126             :  * REPROJECTION_APPROX_ERROR_IN_SRC_SRS_UNIT to be taken into account.
    2127             :  * </li>
    2128             :  * <li>
    2129             :  * AREA_OF_INTEREST=west_lon_deg,south_lat_deg,east_lon_deg,north_lat_deg. (GDAL
    2130             :  * &gt;= 3.0) Area of interest, used to compute the best coordinate operation
    2131             :  * between the source and target SRS. If not specified, the bounding box of the
    2132             :  * source raster will be used.
    2133             :  * </li>
    2134             :  * <li> GEOLOC_BACKMAP_OVERSAMPLE_FACTOR=[0.1,2]. (GDAL &gt;= 3.5) Oversample
    2135             :  * factor used to derive the size of the "backmap" used for geolocation array
    2136             :  * transformers. Default value is 1.3.
    2137             :  * </li>
    2138             :  * <li> GEOLOC_USE_TEMP_DATASETS=YES/NO.
    2139             :  * (GDAL &gt;= 3.5) Whether temporary GeoTIFF datasets should be used to store
    2140             :  * the backmap. The default is NO, that is to use in-memory arrays, unless the
    2141             :  * number of pixels of the geolocation array is greater than 16 megapixels.
    2142             :  * </li>
    2143             :  * <li>
    2144             :  * GEOLOC_ARRAY/SRC_GEOLOC_ARRAY=filename. (GDAL &gt;= 3.5.2) Name of a GDAL
    2145             :  * dataset containing a geolocation array and associated metadata. This is an
    2146             :  * alternative to having geolocation information described in the GEOLOCATION
    2147             :  * metadata domain of the source dataset. The dataset specified may have a
    2148             :  * GEOLOCATION metadata domain containing appropriate metadata, however default
    2149             :  * values are assigned for all omitted items. X_BAND defaults to 1 and Y_BAND to
    2150             :  * 2, however the dataset must contain exactly 2 bands. PIXEL_OFFSET and
    2151             :  * LINE_OFFSET default to 0. PIXEL_STEP and LINE_STEP default to the ratio of
    2152             :  * the width/height of the source dataset divided by the with/height of the
    2153             :  * geolocation array. SRS defaults to the geolocation array dataset's spatial
    2154             :  * reference system if set, otherwise WGS84 is used.
    2155             :  * GEOREFERENCING_CONVENTION is selected from the main metadata domain if it
    2156             :  * is omitted from the GEOLOCATION domain, and if not available
    2157             :  * TOP_LEFT_CORNER is assigned as a default.
    2158             :  * If GEOLOC_ARRAY is set SRC_METHOD
    2159             :  * defaults to GEOLOC_ARRAY.
    2160             :  * </li>
    2161             :  * <li>DST_GEOLOC_ARRAY=filename. (GDAL &gt;= 3.5.2) Name of a
    2162             :  * GDAL dataset that contains at least 2 bands with the X and Y geolocation
    2163             :  * bands. This is an alternative to having geolocation information described in
    2164             :  * the GEOLOCATION metadata domain of the destination dataset. See
    2165             :  * SRC_GEOLOC_ARRAY description for details, assumptions, and defaults. If this
    2166             :  * option is set, DST_METHOD=GEOLOC_ARRAY will be assumed if not set.
    2167             :  * </li>
    2168             :  * <li>GEOLOC_NORMALIZE_LONGITUDE_MINUS_180_PLUS_180=YES/NO. (GDAL &gt;= 3.12.0)
    2169             :  * Whether to force geolocation longitudes into -180,180 when longitude/latitude is
    2170             :  * the coordinate system of the geolocation arrays. The default is to enable this mode
    2171             :  * when the values in the geolocation array are in the -180,180, otherwise NO.
    2172             :  * </li>
    2173             :  * </ul>
    2174             :  *
    2175             :  * The use case for the *_APPROX_ERROR_* options is when defining an approximate
    2176             :  * transformer on top of the GenImgProjTransformer globally is not practical.
    2177             :  * Such a use case is when the source dataset has RPC with a RPC DEM. In such
    2178             :  * case we don't want to use the approximate transformer on the RPC
    2179             :  * transformation, as the RPC DEM generally involves non-linearities that the
    2180             :  * approximate transformer will not detect. In such case, we must a
    2181             :  * non-approximated GenImgProjTransformer, but it might be worthwhile to use
    2182             :  * approximate sub- transformers, for example on coordinate reprojection. For
    2183             :  * example if warping from a source dataset with RPC to a destination dataset
    2184             :  * with a UTM projection, since the inverse UTM transformation is rather costly.
    2185             :  * In which case, one can use the REPROJECTION_APPROX_ERROR_IN_SRC_SRS_UNIT and
    2186             :  * REPROJECTION_APPROX_ERROR_IN_DST_SRS_UNIT options.
    2187             :  *
    2188             :  * The list of supported options can also be programmatically obtained with
    2189             :  * GDALGetGenImgProjTranformerOptionList().
    2190             :  *
    2191             :  * @param hSrcDS source dataset, or NULL.
    2192             :  * @param hDstDS destination dataset (or NULL).
    2193             :  * @param papszOptions NULL-terminated list of string options (or NULL).
    2194             :  *
    2195             :  * @return handle suitable for use GDALGenImgProjTransform(), and to be
    2196             :  * deallocated with GDALDestroyGenImgProjTransformer() or NULL on failure.
    2197             :  */
    2198             : /* clang-format on */
    2199             : 
    2200        1645 : void *GDALCreateGenImgProjTransformer2(GDALDatasetH hSrcDS, GDALDatasetH hDstDS,
    2201             :                                        CSLConstList papszOptions)
    2202             : 
    2203             : {
    2204        1645 :     GDALValidateOptions(GDALGetGenImgProjTranformerOptionList(), papszOptions,
    2205             :                         "option", "transformer options");
    2206             : 
    2207        1645 :     double dfWestLongitudeDeg = 0.0;
    2208        1645 :     double dfSouthLatitudeDeg = 0.0;
    2209        1645 :     double dfEastLongitudeDeg = 0.0;
    2210        1645 :     double dfNorthLatitudeDeg = 0.0;
    2211        1645 :     bool bHasAreaOfInterest = false;
    2212        1645 :     if (const char *pszAreaOfInterest =
    2213        1645 :             CSLFetchNameValue(papszOptions, "AREA_OF_INTEREST"))
    2214             :     {
    2215             :         const CPLStringList aosTokens(
    2216           0 :             CSLTokenizeString2(pszAreaOfInterest, ", ", 0));
    2217           0 :         if (aosTokens.size() == 4)
    2218             :         {
    2219           0 :             dfWestLongitudeDeg = CPLAtof(aosTokens[0]);
    2220           0 :             dfSouthLatitudeDeg = CPLAtof(aosTokens[1]);
    2221           0 :             dfEastLongitudeDeg = CPLAtof(aosTokens[2]);
    2222           0 :             dfNorthLatitudeDeg = CPLAtof(aosTokens[3]);
    2223           0 :             bHasAreaOfInterest = true;
    2224             :         }
    2225             :     }
    2226             : 
    2227        1645 :     const char *pszCO = CSLFetchNameValue(papszOptions, "COORDINATE_OPERATION");
    2228             : 
    2229             :     const auto SetAxisMapping =
    2230        3610 :         [papszOptions](OGRSpatialReference &oSRS, const char *pszPrefix)
    2231             :     {
    2232        1204 :         const char *pszMapping = CSLFetchNameValue(
    2233        2408 :             papszOptions, std::string(pszPrefix)
    2234        1204 :                               .append("_DATA_AXIS_TO_SRS_AXIS_MAPPING")
    2235             :                               .c_str());
    2236        1204 :         if (pszMapping)
    2237             :         {
    2238           4 :             CPLStringList aosTokens(CSLTokenizeString2(pszMapping, ",", 0));
    2239           4 :             std::vector<int> anMapping;
    2240           6 :             for (int i = 0; i < aosTokens.size(); ++i)
    2241           4 :                 anMapping.push_back(atoi(aosTokens[i]));
    2242           2 :             oSRS.SetDataAxisToSRSAxisMapping(anMapping);
    2243             :         }
    2244             :         else
    2245             :         {
    2246        1202 :             const char *pszStrategy = CSLFetchNameValueDef(
    2247             :                 papszOptions,
    2248        2404 :                 std::string(pszPrefix).append("_AXIS_MAPPING_STRATEGY").c_str(),
    2249             :                 "TRADITIONAL_GIS_ORDER");
    2250        1202 :             if (EQUAL(pszStrategy, "TRADITIONAL_GIS_ORDER"))
    2251        1201 :                 oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
    2252           1 :             else if (EQUAL(pszStrategy, "AUTHORITY_COMPLIANT"))
    2253           1 :                 oSRS.SetAxisMappingStrategy(OAMS_AUTHORITY_COMPLIANT);
    2254             :             else
    2255             :             {
    2256           0 :                 CPLError(CE_Warning, CPLE_AppDefined,
    2257             :                          "Unrecognized value '%s' for %s", pszStrategy,
    2258           0 :                          std::string(pszPrefix)
    2259           0 :                              .append("_AXIS_MAPPING_STRATEGY")
    2260             :                              .c_str());
    2261           0 :                 return false;
    2262             :             }
    2263             :         }
    2264        1204 :         return true;
    2265        1645 :     };
    2266             : 
    2267             :     /* -------------------------------------------------------------------- */
    2268             :     /*      Initialize the transform info.                                  */
    2269             :     /* -------------------------------------------------------------------- */
    2270             :     GDALGenImgProjTransformInfo *psInfo =
    2271        1645 :         GDALCreateGenImgProjTransformerInternal();
    2272             : 
    2273             :     const auto DealWithForwardOrInverse =
    2274        3276 :         [bHasAreaOfInterest, &dfWestLongitudeDeg, &dfSouthLatitudeDeg,
    2275             :          &dfEastLongitudeDeg, &dfNorthLatitudeDeg, pszCO, papszOptions,
    2276             :          &SetAxisMapping](GDALGenImgProjTransformPart &part, GDALDatasetH hDS,
    2277             :                           const char *pszPrefix, OGRSpatialReference &oSRS,
    2278       31841 :                           bool &bCanUseGeoTransform)
    2279             :     {
    2280             :         const int nOrder =
    2281        3276 :             atoi(CSLFetchNameValueDef(papszOptions, "MAX_GCP_ORDER", "0"));
    2282             : 
    2283             :         const bool bGCPUseOK =
    2284        3276 :             CPLTestBool(CSLFetchNameValueDef(papszOptions, "GCPS_OK", "YES"));
    2285        3276 :         const int nMinimumGcps = atoi(
    2286             :             CSLFetchNameValueDef(papszOptions, "REFINE_MINIMUM_GCPS", "-1"));
    2287             : 
    2288             :         const char *pszRefineTolerance =
    2289        3276 :             CSLFetchNameValue(papszOptions, "REFINE_TOLERANCE");
    2290        3276 :         const bool bRefine = pszRefineTolerance != nullptr;
    2291             :         const double dfTolerance =
    2292        3276 :             pszRefineTolerance ? CPLAtof(pszRefineTolerance) : 0.0;
    2293             : 
    2294             :         const std::string osSRSOptionName =
    2295        9828 :             std::string(pszPrefix).append("_SRS");
    2296             :         const char *pszSRS =
    2297        3276 :             CSLFetchNameValue(papszOptions, osSRSOptionName.c_str());
    2298        3276 :         if (pszSRS)
    2299             :         {
    2300        2404 :             if (pszSRS[0] != '\0' &&
    2301        1200 :                 oSRS.SetFromUserInput(pszSRS) != OGRERR_NONE)
    2302             :             {
    2303           0 :                 CPLError(CE_Failure, CPLE_AppDefined,
    2304             :                          "Failed to import coordinate system `%s'.", pszSRS);
    2305           0 :                 return false;
    2306             :             }
    2307        1204 :             if (!SetAxisMapping(oSRS, osSRSOptionName.c_str()))
    2308           0 :                 return false;
    2309             :         }
    2310             : 
    2311        3276 :         CSLConstList papszMD = nullptr;
    2312             :         GDALRPCInfoV2 sRPCInfo;
    2313             : 
    2314        3276 :         bCanUseGeoTransform = false;
    2315             : 
    2316        3276 :         const char *pszMethod = CSLFetchNameValue(
    2317        6552 :             papszOptions, std::string(pszPrefix).append("_METHOD").c_str());
    2318        3276 :         if (!pszMethod && EQUAL(pszPrefix, "SRC"))
    2319        1608 :             pszMethod = CSLFetchNameValue(papszOptions, "METHOD");
    2320             : 
    2321        3276 :         const char *pszGeolocArray = CSLFetchNameValue(
    2322             :             papszOptions,
    2323        6552 :             std::string(pszPrefix).append("_GEOLOC_ARRAY").c_str());
    2324        3276 :         if (!pszGeolocArray && EQUAL(pszPrefix, "SRC"))
    2325        1644 :             pszGeolocArray = CSLFetchNameValue(papszOptions, "GEOLOC_ARRAY");
    2326        3276 :         if (!pszMethod && pszGeolocArray != nullptr)
    2327           9 :             pszMethod = "GEOLOC_ARRAY";
    2328             : 
    2329             :         /* -------------------------------------------------------------------- */
    2330             :         /*      Get forward and inverse geotransform for the source image.      */
    2331             :         /* -------------------------------------------------------------------- */
    2332        3276 :         if (hDS == nullptr ||
    2333          91 :             (pszMethod != nullptr && EQUAL(pszMethod, "NO_GEOTRANSFORM")))
    2334             :         {
    2335        1322 :             part.adfGeoTransform[0] = 0.0;
    2336        1322 :             part.adfGeoTransform[1] = 1.0;
    2337        1322 :             part.adfGeoTransform[2] = 0.0;
    2338        1322 :             part.adfGeoTransform[3] = 0.0;
    2339        1322 :             part.adfGeoTransform[4] = 0.0;
    2340        1322 :             part.adfGeoTransform[5] = 1.0;
    2341        1322 :             memcpy(part.adfInvGeoTransform, part.adfGeoTransform,
    2342             :                    sizeof(double) * 6);
    2343             :         }
    2344        3841 :         else if ((pszMethod == nullptr || EQUAL(pszMethod, "GEOTRANSFORM")) &&
    2345        1887 :                  GDALGetGeoTransform(hDS, part.adfGeoTransform) == CE_None)
    2346             :         {
    2347        1807 :             if (!GDALInvGeoTransform(part.adfGeoTransform,
    2348        1807 :                                      part.adfInvGeoTransform))
    2349             :             {
    2350           0 :                 CPLError(CE_Failure, CPLE_AppDefined,
    2351             :                          "Cannot invert geotransform");
    2352           0 :                 return false;
    2353             :             }
    2354        1807 :             if (pszSRS == nullptr)
    2355             :             {
    2356        1629 :                 auto hSRS = GDALGetSpatialRef(hDS);
    2357        1629 :                 if (hSRS)
    2358        1314 :                     oSRS = *(OGRSpatialReference::FromHandle(hSRS));
    2359             :             }
    2360        1807 :             if (EQUAL(pszPrefix, "SRC"))
    2361             :             {
    2362        1390 :                 if (!bHasAreaOfInterest && pszCO == nullptr && !oSRS.IsEmpty())
    2363             :                 {
    2364        1177 :                     GDALComputeAreaOfInterest(
    2365        1177 :                         &oSRS, part.adfGeoTransform, GDALGetRasterXSize(hDS),
    2366             :                         GDALGetRasterYSize(hDS), dfWestLongitudeDeg,
    2367             :                         dfSouthLatitudeDeg, dfEastLongitudeDeg,
    2368             :                         dfNorthLatitudeDeg);
    2369             :                 }
    2370        1390 :                 bCanUseGeoTransform = true;
    2371             :             }
    2372             :         }
    2373         147 :         else if (bGCPUseOK &&
    2374          80 :                  ((pszMethod == nullptr && GDALGetGCPCount(hDS) >= 4 &&
    2375         147 :                    GDALGetGCPCount(hDS) < 6) ||
    2376          67 :                   (pszMethod != nullptr &&
    2377         295 :                    EQUAL(pszMethod, "GCP_HOMOGRAPHY"))) &&
    2378          22 :                  GDALGetGCPCount(hDS) > 0)
    2379             :         {
    2380          22 :             if (pszSRS == nullptr)
    2381             :             {
    2382          21 :                 auto hSRS = GDALGetGCPSpatialRef(hDS);
    2383          21 :                 if (hSRS)
    2384          21 :                     oSRS = *(OGRSpatialReference::FromHandle(hSRS));
    2385             :             }
    2386             : 
    2387          22 :             const auto nGCPCount = GDALGetGCPCount(hDS);
    2388          22 :             auto pasGCPList = GDALDuplicateGCPs(nGCPCount, GDALGetGCPs(hDS));
    2389          22 :             GDALGCPAntimeridianUnwrap(nGCPCount, pasGCPList, oSRS,
    2390             :                                       papszOptions);
    2391             : 
    2392          22 :             part.pTransformArg =
    2393          22 :                 GDALCreateHomographyTransformerFromGCPs(nGCPCount, pasGCPList);
    2394             : 
    2395          22 :             GDALDeinitGCPs(nGCPCount, pasGCPList);
    2396          22 :             CPLFree(pasGCPList);
    2397             : 
    2398          22 :             if (part.pTransformArg == nullptr)
    2399             :             {
    2400           0 :                 return false;
    2401             :             }
    2402          22 :             part.pTransformer = GDALHomographyTransform;
    2403             :         }
    2404         125 :         else if (bGCPUseOK &&
    2405          66 :                  (pszMethod == nullptr || EQUAL(pszMethod, "GCP_POLYNOMIAL")) &&
    2406         250 :                  GDALGetGCPCount(hDS) > 0 && nOrder >= 0)
    2407             :         {
    2408           8 :             if (pszSRS == nullptr)
    2409             :             {
    2410           8 :                 auto hSRS = GDALGetGCPSpatialRef(hDS);
    2411           8 :                 if (hSRS)
    2412           8 :                     oSRS = *(OGRSpatialReference::FromHandle(hSRS));
    2413             :             }
    2414             : 
    2415           8 :             const auto nGCPCount = GDALGetGCPCount(hDS);
    2416           8 :             auto pasGCPList = GDALDuplicateGCPs(nGCPCount, GDALGetGCPs(hDS));
    2417           8 :             GDALGCPAntimeridianUnwrap(nGCPCount, pasGCPList, oSRS,
    2418             :                                       papszOptions);
    2419             : 
    2420           8 :             if (bRefine)
    2421             :             {
    2422           0 :                 part.pTransformArg = GDALCreateGCPRefineTransformer(
    2423             :                     nGCPCount, pasGCPList, nOrder, FALSE, dfTolerance,
    2424             :                     nMinimumGcps);
    2425             :             }
    2426             :             else
    2427             :             {
    2428           8 :                 part.pTransformArg = GDALCreateGCPTransformer(
    2429             :                     nGCPCount, pasGCPList, nOrder, FALSE);
    2430             :             }
    2431             : 
    2432           8 :             GDALDeinitGCPs(nGCPCount, pasGCPList);
    2433           8 :             CPLFree(pasGCPList);
    2434             : 
    2435           8 :             if (part.pTransformArg == nullptr)
    2436             :             {
    2437           0 :                 return false;
    2438             :             }
    2439           8 :             part.pTransformer = GDALGCPTransform;
    2440             :         }
    2441             : 
    2442         128 :         else if (bGCPUseOK && GDALGetGCPCount(hDS) > 0 && nOrder <= 0 &&
    2443          11 :                  (pszMethod == nullptr || EQUAL(pszMethod, "GCP_TPS")))
    2444             :         {
    2445          11 :             if (pszSRS == nullptr)
    2446             :             {
    2447          11 :                 auto hSRS = GDALGetGCPSpatialRef(hDS);
    2448          11 :                 if (hSRS)
    2449          10 :                     oSRS = *(OGRSpatialReference::FromHandle(hSRS));
    2450             :             }
    2451             : 
    2452          11 :             const auto nGCPCount = GDALGetGCPCount(hDS);
    2453          11 :             auto pasGCPList = GDALDuplicateGCPs(nGCPCount, GDALGetGCPs(hDS));
    2454          11 :             GDALGCPAntimeridianUnwrap(nGCPCount, pasGCPList, oSRS,
    2455             :                                       papszOptions);
    2456             : 
    2457          11 :             part.pTransformArg = GDALCreateTPSTransformerInt(
    2458             :                 nGCPCount, pasGCPList, FALSE, papszOptions);
    2459             : 
    2460          11 :             GDALDeinitGCPs(nGCPCount, pasGCPList);
    2461          11 :             CPLFree(pasGCPList);
    2462             : 
    2463          11 :             if (part.pTransformArg == nullptr)
    2464             :             {
    2465           2 :                 return false;
    2466             :             }
    2467           9 :             part.pTransformer = GDALTPSTransform;
    2468             :         }
    2469             : 
    2470          52 :         else if ((pszMethod == nullptr || EQUAL(pszMethod, "RPC")) &&
    2471         204 :                  (papszMD = GDALGetMetadata(hDS, "RPC")) != nullptr &&
    2472          46 :                  GDALExtractRPCInfoV2(papszMD, &sRPCInfo))
    2473             :         {
    2474          46 :             CPLStringList aosOptions(papszOptions);
    2475          87 :             if (!CSLFetchNameValue(papszOptions, "RPC_HEIGHT") &&
    2476          41 :                 !CSLFetchNameValue(papszOptions, "RPC_DEM"))
    2477             :             {
    2478           8 :                 if (const char *pszHEIGHT_DEFAULT =
    2479           8 :                         CSLFetchNameValue(papszMD, "HEIGHT_DEFAULT"))
    2480             :                 {
    2481           1 :                     CPLDebug("GDAL",
    2482             :                              "For %s, using RPC_HEIGHT = HEIGHT_DEFAULT = %s",
    2483             :                              pszPrefix, pszHEIGHT_DEFAULT);
    2484           1 :                     aosOptions.SetNameValue("RPC_HEIGHT", pszHEIGHT_DEFAULT);
    2485             :                 }
    2486             :             }
    2487          46 :             part.pTransformArg = GDALCreateRPCTransformerV2(&sRPCInfo, FALSE, 0,
    2488          46 :                                                             aosOptions.List());
    2489          46 :             if (part.pTransformArg == nullptr)
    2490             :             {
    2491           1 :                 return false;
    2492             :             }
    2493          45 :             part.pTransformer = GDALRPCTransform;
    2494          45 :             if (pszSRS == nullptr)
    2495             :             {
    2496          45 :                 oSRS.SetFromUserInput(SRS_WKT_WGS84_LAT_LONG);
    2497          45 :                 oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
    2498             :             }
    2499             :         }
    2500             : 
    2501         119 :         else if ((pszMethod == nullptr || EQUAL(pszMethod, "GEOLOC_ARRAY")) &&
    2502          59 :                  ((papszMD = GDALGetMetadata(hDS, "GEOLOCATION")) != nullptr ||
    2503             :                   pszGeolocArray != nullptr))
    2504             :         {
    2505          52 :             CPLStringList aosGeolocMD;  // keep in this scope
    2506          52 :             if (pszGeolocArray != nullptr)
    2507             :             {
    2508           9 :                 if (papszMD != nullptr)
    2509             :                 {
    2510           0 :                     CPLError(
    2511             :                         CE_Warning, CPLE_AppDefined,
    2512             :                         "Both GEOLOCATION metadata domain on the source "
    2513             :                         "dataset "
    2514             :                         "and [%s_]GEOLOC_ARRAY transformer option are set. "
    2515             :                         "Only using the later.",
    2516             :                         pszPrefix);
    2517             :                 }
    2518           9 :                 aosGeolocMD = GDALCreateGeolocationMetadata(
    2519             :                     hDS, pszGeolocArray,
    2520           9 :                     /* bIsSource= */ EQUAL(pszPrefix, "SRC"));
    2521           9 :                 if (aosGeolocMD.empty())
    2522             :                 {
    2523           3 :                     return false;
    2524             :                 }
    2525           6 :                 papszMD = aosGeolocMD.List();
    2526             :             }
    2527             : 
    2528          49 :             part.pTransformArg = GDALCreateGeoLocTransformerEx(
    2529             :                 hDS, papszMD, FALSE, nullptr, papszOptions);
    2530          49 :             if (part.pTransformArg == nullptr)
    2531             :             {
    2532           2 :                 return false;
    2533             :             }
    2534          47 :             part.pTransformer = GDALGeoLocTransform;
    2535          47 :             if (pszSRS == nullptr)
    2536             :             {
    2537          47 :                 pszSRS = CSLFetchNameValue(papszMD, "SRS");
    2538          47 :                 if (pszSRS)
    2539             :                 {
    2540          45 :                     oSRS.SetFromUserInput(pszSRS);
    2541          45 :                     oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
    2542             :                 }
    2543             :             }
    2544             :         }
    2545             : 
    2546           8 :         else if (pszMethod != nullptr && EQUAL(pszPrefix, "SRC"))
    2547             :         {
    2548           1 :             CPLError(CE_Failure, CPLE_AppDefined,
    2549             :                      "Unable to compute a %s based transformation between "
    2550             :                      "pixel/line and georeferenced coordinates for %s.",
    2551             :                      pszMethod, GDALGetDescription(hDS));
    2552             : 
    2553           1 :             return false;
    2554             :         }
    2555             : 
    2556             :         else
    2557             :         {
    2558           7 :             CPLError(CE_Failure, CPLE_AppDefined,
    2559             :                      "Unable to compute a transformation between pixel/line "
    2560             :                      "and georeferenced coordinates for %s. "
    2561             :                      "There is no affine transformation and no GCPs. "
    2562             :                      "Specify transformation option %s_METHOD=NO_GEOTRANSFORM "
    2563             :                      "to bypass this check.",
    2564             :                      GDALGetDescription(hDS), pszPrefix);
    2565             : 
    2566           7 :             return false;
    2567             :         }
    2568             : 
    2569             :         /* ---------------------------------------------------------------- */
    2570             :         /*      Handle optional source approximation transformer.           */
    2571             :         /* ---------------------------------------------------------------- */
    2572        3260 :         if (part.pTransformer)
    2573             :         {
    2574         131 :             const char *pszApproxErrorFwd = CSLFetchNameValue(
    2575         262 :                 papszOptions, std::string(pszPrefix)
    2576         131 :                                   .append("_APPROX_ERROR_IN_SRS_UNIT")
    2577             :                                   .c_str());
    2578         131 :             const char *pszApproxErrorReverse = CSLFetchNameValue(
    2579         262 :                 papszOptions, std::string(pszPrefix)
    2580         131 :                                   .append("_APPROX_ERROR_IN_PIXEL")
    2581             :                                   .c_str());
    2582         131 :             if (pszApproxErrorFwd && pszApproxErrorReverse)
    2583             :             {
    2584           1 :                 void *pArg = GDALCreateApproxTransformer2(
    2585             :                     part.pTransformer, part.pTransformArg,
    2586             :                     CPLAtof(pszApproxErrorFwd), CPLAtof(pszApproxErrorReverse));
    2587           1 :                 if (pArg == nullptr)
    2588             :                 {
    2589           0 :                     return false;
    2590             :                 }
    2591           1 :                 part.pTransformArg = pArg;
    2592           1 :                 part.pTransformer = GDALApproxTransform;
    2593           1 :                 GDALApproxTransformerOwnsSubtransformer(part.pTransformArg,
    2594             :                                                         TRUE);
    2595             :             }
    2596             :         }
    2597             : 
    2598        3260 :         return true;
    2599        1645 :     };
    2600             : 
    2601             :     /* -------------------------------------------------------------------- */
    2602             :     /*      Get forward and inverse geotransform for the source image.      */
    2603             :     /* -------------------------------------------------------------------- */
    2604        1645 :     bool bCanUseSrcGeoTransform = false;
    2605        3290 :     OGRSpatialReference oSrcSRS;
    2606        1645 :     if (!DealWithForwardOrInverse(psInfo->sSrcParams, hSrcDS, "SRC", oSrcSRS,
    2607             :                                   bCanUseSrcGeoTransform))
    2608             :     {
    2609          14 :         GDALDestroyGenImgProjTransformer(psInfo);
    2610          14 :         return nullptr;
    2611             :     }
    2612             : 
    2613             :     /* -------------------------------------------------------------------- */
    2614             :     /*      Get forward and inverse geotransform for destination image.     */
    2615             :     /*      If we have no destination use a unit transform.                 */
    2616             :     /* -------------------------------------------------------------------- */
    2617        1631 :     bool bIgnored = false;
    2618        3262 :     OGRSpatialReference oDstSRS;
    2619        1631 :     if (!DealWithForwardOrInverse(psInfo->sDstParams, hDstDS, "DST", oDstSRS,
    2620             :                                   bIgnored))
    2621             :     {
    2622           2 :         GDALDestroyGenImgProjTransformer(psInfo);
    2623           2 :         return nullptr;
    2624             :     }
    2625             : 
    2626             :     /* -------------------------------------------------------------------- */
    2627             :     /*      Setup reprojection.                                             */
    2628             :     /* -------------------------------------------------------------------- */
    2629             : 
    2630        1629 :     if (CPLFetchBool(papszOptions, "STRIP_VERT_CS", false))
    2631             :     {
    2632           0 :         if (oSrcSRS.IsCompound())
    2633             :         {
    2634           0 :             oSrcSRS.StripVertical();
    2635             :         }
    2636           0 :         if (oDstSRS.IsCompound())
    2637             :         {
    2638           0 :             oDstSRS.StripVertical();
    2639             :         }
    2640             :     }
    2641             : 
    2642             :     const bool bMayInsertCenterLong =
    2643        2808 :         (bCanUseSrcGeoTransform && !oSrcSRS.IsEmpty() && hSrcDS &&
    2644        1179 :          CPLFetchBool(papszOptions, "INSERT_CENTER_LONG", true));
    2645             :     const char *pszSrcCoordEpoch =
    2646        1629 :         CSLFetchNameValue(papszOptions, "SRC_COORDINATE_EPOCH");
    2647             :     const char *pszDstCoordEpoch =
    2648        1629 :         CSLFetchNameValue(papszOptions, "DST_COORDINATE_EPOCH");
    2649        2973 :     if ((!oSrcSRS.IsEmpty() && !oDstSRS.IsEmpty() &&
    2650        1231 :          (pszSrcCoordEpoch || pszDstCoordEpoch || !oSrcSRS.IsSame(&oDstSRS) ||
    2651        2973 :           (oSrcSRS.IsGeographic() && bMayInsertCenterLong))) ||
    2652             :         pszCO)
    2653             :     {
    2654         880 :         CPLStringList aosOptions;
    2655             : 
    2656         880 :         if (bMayInsertCenterLong)
    2657             :         {
    2658         863 :             InsertCenterLong(hSrcDS, &oSrcSRS, aosOptions);
    2659             :         }
    2660             : 
    2661         880 :         if (CPLFetchBool(papszOptions, "PROMOTE_TO_3D", false))
    2662             :         {
    2663          19 :             oSrcSRS.PromoteTo3D(nullptr);
    2664          19 :             oDstSRS.PromoteTo3D(nullptr);
    2665             :         }
    2666             : 
    2667         880 :         if (!(dfWestLongitudeDeg == 0.0 && dfSouthLatitudeDeg == 0.0 &&
    2668          30 :               dfEastLongitudeDeg == 0.0 && dfNorthLatitudeDeg == 0.0))
    2669             :         {
    2670             :             aosOptions.SetNameValue(
    2671             :                 "AREA_OF_INTEREST",
    2672             :                 CPLSPrintf("%.16g,%.16g,%.16g,%.16g", dfWestLongitudeDeg,
    2673             :                            dfSouthLatitudeDeg, dfEastLongitudeDeg,
    2674         850 :                            dfNorthLatitudeDeg));
    2675             :         }
    2676         880 :         if (pszCO)
    2677             :         {
    2678           7 :             aosOptions.SetNameValue("COORDINATE_OPERATION", pszCO);
    2679             :         }
    2680             : 
    2681             :         const char *pszCoordEpoch =
    2682         880 :             CSLFetchNameValue(papszOptions, "COORDINATE_EPOCH");
    2683         880 :         if (pszCoordEpoch)
    2684             :         {
    2685           1 :             aosOptions.SetNameValue("COORDINATE_EPOCH", pszCoordEpoch);
    2686             :         }
    2687             : 
    2688         880 :         if (pszSrcCoordEpoch)
    2689             :         {
    2690           0 :             aosOptions.SetNameValue("SRC_COORDINATE_EPOCH", pszSrcCoordEpoch);
    2691           0 :             oSrcSRS.SetCoordinateEpoch(CPLAtof(pszSrcCoordEpoch));
    2692             :         }
    2693             : 
    2694         880 :         if (pszDstCoordEpoch)
    2695             :         {
    2696           0 :             aosOptions.SetNameValue("DST_COORDINATE_EPOCH", pszDstCoordEpoch);
    2697           0 :             oDstSRS.SetCoordinateEpoch(CPLAtof(pszDstCoordEpoch));
    2698             :         }
    2699             : 
    2700         880 :         if (const char *pszAllowBallpark =
    2701         880 :                 CSLFetchNameValue(papszOptions, "ALLOW_BALLPARK"))
    2702             :         {
    2703           2 :             aosOptions.SetNameValue("ALLOW_BALLPARK", pszAllowBallpark);
    2704             :         }
    2705             : 
    2706         880 :         if (const char *pszOnlyBest =
    2707         880 :                 CSLFetchNameValue(papszOptions, "ONLY_BEST"))
    2708             :         {
    2709           1 :             aosOptions.SetNameValue("ONLY_BEST", pszOnlyBest);
    2710             :         }
    2711             : 
    2712         884 :         psInfo->pReprojectArg = GDALCreateReprojectionTransformerEx(
    2713         880 :             !oSrcSRS.IsEmpty() ? OGRSpatialReference::ToHandle(&oSrcSRS)
    2714             :                                : nullptr,
    2715         880 :             !oDstSRS.IsEmpty() ? OGRSpatialReference::ToHandle(&oDstSRS)
    2716             :                                : nullptr,
    2717         880 :             aosOptions.List());
    2718             : 
    2719         880 :         if (pszCO)
    2720             :         {
    2721           7 :             psInfo->bHasCustomTransformationPipeline = true;
    2722             :         }
    2723             : 
    2724         880 :         if (psInfo->pReprojectArg == nullptr)
    2725             :         {
    2726           4 :             GDALDestroyGenImgProjTransformer(psInfo);
    2727           4 :             return nullptr;
    2728             :         }
    2729         876 :         psInfo->pReproject = GDALReprojectionTransform;
    2730             : 
    2731             :         /* --------------------------------------------------------------------
    2732             :          */
    2733             :         /*      Handle optional reprojection approximation transformer. */
    2734             :         /* --------------------------------------------------------------------
    2735             :          */
    2736         876 :         const char *psApproxErrorFwd = CSLFetchNameValue(
    2737             :             papszOptions, "REPROJECTION_APPROX_ERROR_IN_DST_SRS_UNIT");
    2738         876 :         const char *psApproxErrorReverse = CSLFetchNameValue(
    2739             :             papszOptions, "REPROJECTION_APPROX_ERROR_IN_SRC_SRS_UNIT");
    2740         876 :         if (psApproxErrorFwd && psApproxErrorReverse)
    2741             :         {
    2742           1 :             void *pArg = GDALCreateApproxTransformer2(
    2743             :                 psInfo->pReproject, psInfo->pReprojectArg,
    2744             :                 CPLAtof(psApproxErrorFwd), CPLAtof(psApproxErrorReverse));
    2745           1 :             if (pArg == nullptr)
    2746             :             {
    2747           0 :                 GDALDestroyGenImgProjTransformer(psInfo);
    2748           0 :                 return nullptr;
    2749             :             }
    2750           1 :             psInfo->pReprojectArg = pArg;
    2751           1 :             psInfo->pReproject = GDALApproxTransform;
    2752           1 :             GDALApproxTransformerOwnsSubtransformer(psInfo->pReprojectArg,
    2753             :                                                     TRUE);
    2754             :         }
    2755             :     }
    2756             : 
    2757        1625 :     return psInfo;
    2758             : }
    2759             : 
    2760             : /************************************************************************/
    2761             : /*                  GDALRefreshGenImgProjTransformer()                  */
    2762             : /************************************************************************/
    2763             : 
    2764        1094 : void GDALRefreshGenImgProjTransformer(void *hTransformArg)
    2765             : {
    2766        1094 :     GDALGenImgProjTransformInfo *psInfo =
    2767             :         static_cast<GDALGenImgProjTransformInfo *>(hTransformArg);
    2768             : 
    2769        1839 :     if (psInfo->pReprojectArg &&
    2770         745 :         psInfo->bCheckWithInvertPROJ != GetCurrentCheckWithInvertPROJ())
    2771             :     {
    2772          68 :         psInfo->bCheckWithInvertPROJ = !psInfo->bCheckWithInvertPROJ;
    2773             : 
    2774             :         CPLXMLNode *psXML =
    2775          68 :             GDALSerializeTransformer(psInfo->pReproject, psInfo->pReprojectArg);
    2776          68 :         GDALDestroyTransformer(psInfo->pReprojectArg);
    2777          68 :         GDALDeserializeTransformer(psXML, &psInfo->pReproject,
    2778             :                                    &psInfo->pReprojectArg);
    2779          68 :         CPLDestroyXMLNode(psXML);
    2780             :     }
    2781        1094 : }
    2782             : 
    2783             : /************************************************************************/
    2784             : /*                  GDALCreateGenImgProjTransformer3()                  */
    2785             : /************************************************************************/
    2786             : 
    2787             : /**
    2788             :  * Create image to image transformer.
    2789             :  *
    2790             :  * This function creates a transformation object that maps from pixel/line
    2791             :  * coordinates on one image to pixel/line coordinates on another image.  The
    2792             :  * images may potentially be georeferenced in different coordinate systems,
    2793             :  * and may used GCPs to map between their pixel/line coordinates and
    2794             :  * georeferenced coordinates (as opposed to the default assumption that their
    2795             :  * geotransform should be used).
    2796             :  *
    2797             :  * This transformer potentially performs three concatenated transformations.
    2798             :  *
    2799             :  * The first stage is from source image pixel/line coordinates to source
    2800             :  * image georeferenced coordinates, and may be done using the geotransform,
    2801             :  * or if not defined using a polynomial model derived from GCPs.  If GCPs
    2802             :  * are used this stage is accomplished using GDALGCPTransform().
    2803             :  *
    2804             :  * The second stage is to change projections from the source coordinate system
    2805             :  * to the destination coordinate system, assuming they differ.  This is
    2806             :  * accomplished internally using GDALReprojectionTransform().
    2807             :  *
    2808             :  * The third stage is converting from destination image georeferenced
    2809             :  * coordinates to destination image coordinates.  This is done using the
    2810             :  * destination image geotransform, or if not available, using a polynomial
    2811             :  * model derived from GCPs. If GCPs are used this stage is accomplished using
    2812             :  * GDALGCPTransform().  This stage is skipped if hDstDS is NULL when the
    2813             :  * transformation is created.
    2814             :  *
    2815             :  * @param pszSrcWKT source WKT (or NULL).
    2816             :  * @param padfSrcGeoTransform source geotransform (or NULL).
    2817             :  * @param pszDstWKT destination WKT (or NULL).
    2818             :  * @param padfDstGeoTransform destination geotransform (or NULL).
    2819             :  *
    2820             :  * @return handle suitable for use GDALGenImgProjTransform(), and to be
    2821             :  * deallocated with GDALDestroyGenImgProjTransformer() or NULL on failure.
    2822             :  */
    2823             : 
    2824           0 : void *GDALCreateGenImgProjTransformer3(const char *pszSrcWKT,
    2825             :                                        const double *padfSrcGeoTransform,
    2826             :                                        const char *pszDstWKT,
    2827             :                                        const double *padfDstGeoTransform)
    2828             : 
    2829             : {
    2830           0 :     OGRSpatialReference oSrcSRS;
    2831           0 :     if (pszSrcWKT)
    2832             :     {
    2833           0 :         oSrcSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
    2834           0 :         if (pszSrcWKT[0] != '\0' &&
    2835           0 :             oSrcSRS.importFromWkt(pszSrcWKT) != OGRERR_NONE)
    2836             :         {
    2837           0 :             CPLError(CE_Failure, CPLE_AppDefined,
    2838             :                      "Failed to import coordinate system `%s'.", pszSrcWKT);
    2839           0 :             return nullptr;
    2840             :         }
    2841             :     }
    2842             : 
    2843           0 :     OGRSpatialReference oDstSRS;
    2844           0 :     if (pszDstWKT)
    2845             :     {
    2846           0 :         oDstSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
    2847           0 :         if (pszDstWKT[0] != '\0' &&
    2848           0 :             oDstSRS.importFromWkt(pszDstWKT) != OGRERR_NONE)
    2849             :         {
    2850           0 :             CPLError(CE_Failure, CPLE_AppDefined,
    2851             :                      "Failed to import coordinate system `%s'.", pszDstWKT);
    2852           0 :             return nullptr;
    2853             :         }
    2854             :     }
    2855           0 :     return GDALCreateGenImgProjTransformer4(
    2856             :         OGRSpatialReference::ToHandle(&oSrcSRS), padfSrcGeoTransform,
    2857           0 :         OGRSpatialReference::ToHandle(&oDstSRS), padfDstGeoTransform, nullptr);
    2858             : }
    2859             : 
    2860             : /************************************************************************/
    2861             : /*                  GDALCreateGenImgProjTransformer4()                  */
    2862             : /************************************************************************/
    2863             : 
    2864             : /**
    2865             :  * Create image to image transformer.
    2866             :  *
    2867             :  * Similar to GDALCreateGenImgProjTransformer3(), except that it takes
    2868             :  * OGRSpatialReferenceH objects and options.
    2869             :  * The options are the ones supported by GDALCreateReprojectionTransformerEx()
    2870             :  *
    2871             :  * @since GDAL 3.0
    2872             :  */
    2873          16 : void *GDALCreateGenImgProjTransformer4(OGRSpatialReferenceH hSrcSRS,
    2874             :                                        const double *padfSrcGeoTransform,
    2875             :                                        OGRSpatialReferenceH hDstSRS,
    2876             :                                        const double *padfDstGeoTransform,
    2877             :                                        const char *const *papszOptions)
    2878             : {
    2879             :     /* -------------------------------------------------------------------- */
    2880             :     /*      Initialize the transform info.                                  */
    2881             :     /* -------------------------------------------------------------------- */
    2882             :     GDALGenImgProjTransformInfo *psInfo =
    2883          16 :         GDALCreateGenImgProjTransformerInternal();
    2884             : 
    2885             :     /* -------------------------------------------------------------------- */
    2886             :     /*      Get forward and inverse geotransform for the source image.      */
    2887             :     /* -------------------------------------------------------------------- */
    2888             : 
    2889             :     const auto SetParams =
    2890          30 :         [](GDALGenImgProjTransformPart &part, const double *padfGT)
    2891             :     {
    2892          30 :         if (padfGT)
    2893             :         {
    2894          30 :             memcpy(part.adfGeoTransform, padfGT, sizeof(part.adfGeoTransform));
    2895          30 :             if (!GDALInvGeoTransform(part.adfGeoTransform,
    2896          30 :                                      part.adfInvGeoTransform))
    2897             :             {
    2898           2 :                 CPLError(CE_Failure, CPLE_AppDefined,
    2899             :                          "Cannot invert geotransform");
    2900           2 :                 return false;
    2901             :             }
    2902             :         }
    2903             :         else
    2904             :         {
    2905           0 :             part.adfGeoTransform[0] = 0.0;
    2906           0 :             part.adfGeoTransform[1] = 1.0;
    2907           0 :             part.adfGeoTransform[2] = 0.0;
    2908           0 :             part.adfGeoTransform[3] = 0.0;
    2909           0 :             part.adfGeoTransform[4] = 0.0;
    2910           0 :             part.adfGeoTransform[5] = 1.0;
    2911           0 :             memcpy(part.adfInvGeoTransform, part.adfGeoTransform,
    2912             :                    sizeof(double) * 6);
    2913             :         }
    2914          28 :         return true;
    2915             :     };
    2916             : 
    2917          16 :     if (!SetParams(psInfo->sSrcParams, padfSrcGeoTransform))
    2918             :     {
    2919           1 :         GDALDestroyGenImgProjTransformer(psInfo);
    2920           1 :         return nullptr;
    2921             :     }
    2922             : 
    2923             :     /* -------------------------------------------------------------------- */
    2924             :     /*      Setup reprojection.                                             */
    2925             :     /* -------------------------------------------------------------------- */
    2926          15 :     OGRSpatialReference *poSrcSRS = OGRSpatialReference::FromHandle(hSrcSRS);
    2927          15 :     OGRSpatialReference *poDstSRS = OGRSpatialReference::FromHandle(hDstSRS);
    2928          30 :     if (!poSrcSRS->IsEmpty() && !poDstSRS->IsEmpty() &&
    2929          15 :         !poSrcSRS->IsSame(poDstSRS))
    2930             :     {
    2931           4 :         psInfo->pReprojectArg =
    2932           4 :             GDALCreateReprojectionTransformerEx(hSrcSRS, hDstSRS, papszOptions);
    2933           4 :         if (psInfo->pReprojectArg == nullptr)
    2934             :         {
    2935           1 :             GDALDestroyGenImgProjTransformer(psInfo);
    2936           1 :             return nullptr;
    2937             :         }
    2938           3 :         psInfo->pReproject = GDALReprojectionTransform;
    2939             :     }
    2940             : 
    2941             :     /* -------------------------------------------------------------------- */
    2942             :     /*      Get forward and inverse geotransform for destination image.     */
    2943             :     /*      If we have no destination matrix use a unit transform.          */
    2944             :     /* -------------------------------------------------------------------- */
    2945          14 :     if (!SetParams(psInfo->sDstParams, padfDstGeoTransform))
    2946             :     {
    2947           1 :         GDALDestroyGenImgProjTransformer(psInfo);
    2948           1 :         return nullptr;
    2949             :     }
    2950             : 
    2951          13 :     return psInfo;
    2952             : }
    2953             : 
    2954             : /************************************************************************/
    2955             : /*            GDALSetGenImgProjTransformerDstGeoTransform()             */
    2956             : /************************************************************************/
    2957             : 
    2958             : /**
    2959             :  * Set GenImgProj output geotransform.
    2960             :  *
    2961             :  * Normally the "destination geotransform", or transformation between
    2962             :  * georeferenced output coordinates and pixel/line coordinates on the
    2963             :  * destination file is extracted from the destination file by
    2964             :  * GDALCreateGenImgProjTransformer() and stored in the GenImgProj private
    2965             :  * info.  However, sometimes it is inconvenient to have an output file
    2966             :  * handle with appropriate geotransform information when creating the
    2967             :  * transformation.  For these cases, this function can be used to apply
    2968             :  * the destination geotransform.
    2969             :  *
    2970             :  * @param hTransformArg the handle to update.
    2971             :  * @param padfGeoTransform the destination geotransform to apply (six doubles).
    2972             :  */
    2973             : 
    2974         875 : void GDALSetGenImgProjTransformerDstGeoTransform(void *hTransformArg,
    2975             :                                                  const double *padfGeoTransform)
    2976             : 
    2977             : {
    2978         875 :     VALIDATE_POINTER0(hTransformArg,
    2979             :                       "GDALSetGenImgProjTransformerDstGeoTransform");
    2980             : 
    2981         875 :     GDALGenImgProjTransformInfo *psInfo =
    2982             :         static_cast<GDALGenImgProjTransformInfo *>(hTransformArg);
    2983             : 
    2984         875 :     memcpy(psInfo->sDstParams.adfGeoTransform, padfGeoTransform,
    2985             :            sizeof(double) * 6);
    2986         875 :     if (!GDALInvGeoTransform(psInfo->sDstParams.adfGeoTransform,
    2987         875 :                              psInfo->sDstParams.adfInvGeoTransform))
    2988             :     {
    2989           0 :         CPLError(CE_Failure, CPLE_AppDefined, "Cannot invert geotransform");
    2990             :     }
    2991             : }
    2992             : 
    2993             : /************************************************************************/
    2994             : /*                  GDALDestroyGenImgProjTransformer()                  */
    2995             : /************************************************************************/
    2996             : 
    2997             : /**
    2998             :  * GenImgProjTransformer deallocator.
    2999             :  *
    3000             :  * This function is used to deallocate the handle created with
    3001             :  * GDALCreateGenImgProjTransformer().
    3002             :  *
    3003             :  * @param hTransformArg the handle to deallocate.
    3004             :  */
    3005             : 
    3006        1900 : void GDALDestroyGenImgProjTransformer(void *hTransformArg)
    3007             : 
    3008             : {
    3009        1900 :     if (hTransformArg == nullptr)
    3010           0 :         return;
    3011             : 
    3012        1900 :     GDALGenImgProjTransformInfo *psInfo =
    3013             :         static_cast<GDALGenImgProjTransformInfo *>(hTransformArg);
    3014             : 
    3015        1900 :     if (psInfo->sSrcParams.pTransformArg != nullptr)
    3016         147 :         GDALDestroyTransformer(psInfo->sSrcParams.pTransformArg);
    3017             : 
    3018        1900 :     if (psInfo->sDstParams.pTransformArg != nullptr)
    3019           7 :         GDALDestroyTransformer(psInfo->sDstParams.pTransformArg);
    3020             : 
    3021        1900 :     if (psInfo->pReprojectArg != nullptr)
    3022         996 :         GDALDestroyTransformer(psInfo->pReprojectArg);
    3023             : 
    3024        1900 :     CPLFree(psInfo);
    3025             : }
    3026             : 
    3027             : /************************************************************************/
    3028             : /*                      GDALGenImgProjTransform()                       */
    3029             : /************************************************************************/
    3030             : 
    3031             : /**
    3032             :  * Perform general image reprojection transformation.
    3033             :  *
    3034             :  * Actually performs the transformation setup in
    3035             :  * GDALCreateGenImgProjTransformer().  This function matches the signature
    3036             :  * required by the GDALTransformerFunc(), and more details on the arguments
    3037             :  * can be found in that topic.
    3038             :  */
    3039             : 
    3040             : #ifdef DEBUG_APPROX_TRANSFORMER
    3041             : int countGDALGenImgProjTransform = 0;
    3042             : #endif
    3043             : 
    3044     2163140 : int GDALGenImgProjTransform(void *pTransformArgIn, int bDstToSrc,
    3045             :                             int nPointCount, double *padfX, double *padfY,
    3046             :                             double *padfZ, int *panSuccess)
    3047             : {
    3048     2163140 :     GDALGenImgProjTransformInfo *psInfo =
    3049             :         static_cast<GDALGenImgProjTransformInfo *>(pTransformArgIn);
    3050             : 
    3051             : #ifdef DEBUG_APPROX_TRANSFORMER
    3052             :     CPLAssert(nPointCount > 0);
    3053             :     countGDALGenImgProjTransform += nPointCount;
    3054             : #endif
    3055             : 
    3056    19851500 :     for (int i = 0; i < nPointCount; i++)
    3057             :     {
    3058    17688400 :         panSuccess[i] = (padfX[i] != HUGE_VAL && padfY[i] != HUGE_VAL);
    3059             :     }
    3060             : 
    3061     2163140 :     int ret = TRUE;
    3062             : 
    3063             :     /* -------------------------------------------------------------------- */
    3064             :     /*      Convert from src (dst) pixel/line to src (dst)                  */
    3065             :     /*      georeferenced coordinates.                                      */
    3066             :     /* -------------------------------------------------------------------- */
    3067             :     {
    3068     2163140 :         const auto params = bDstToSrc ? psInfo->sDstParams : psInfo->sSrcParams;
    3069     2163140 :         const double *padfGeoTransform = params.adfGeoTransform;
    3070     2163140 :         void *pTransformArg = params.pTransformArg;
    3071     2163140 :         GDALTransformerFunc pTransformer = params.pTransformer;
    3072             : 
    3073     2163140 :         if (pTransformArg != nullptr)
    3074             :         {
    3075       41619 :             if (!pTransformer(pTransformArg, FALSE, nPointCount, padfX, padfY,
    3076             :                               padfZ, panSuccess))
    3077        1773 :                 ret = FALSE;
    3078             :         }
    3079             :         else
    3080             :         {
    3081    19730200 :             for (int i = 0; i < nPointCount; i++)
    3082             :             {
    3083    17608700 :                 if (!panSuccess[i])
    3084        2008 :                     continue;
    3085             : 
    3086    17606600 :                 const double dfNewX = padfGeoTransform[0] +
    3087    17606600 :                                       padfX[i] * padfGeoTransform[1] +
    3088    17606600 :                                       padfY[i] * padfGeoTransform[2];
    3089    17606600 :                 const double dfNewY = padfGeoTransform[3] +
    3090    17606600 :                                       padfX[i] * padfGeoTransform[4] +
    3091    17606600 :                                       padfY[i] * padfGeoTransform[5];
    3092             : 
    3093    17606600 :                 padfX[i] = dfNewX;
    3094    17606600 :                 padfY[i] = dfNewY;
    3095             :             }
    3096             :         }
    3097             :     }
    3098             : 
    3099             :     /* -------------------------------------------------------------------- */
    3100             :     /*      Reproject if needed.                                            */
    3101             :     /* -------------------------------------------------------------------- */
    3102     2163140 :     if (psInfo->pReprojectArg)
    3103             :     {
    3104     1812830 :         if (!psInfo->pReproject(psInfo->pReprojectArg, bDstToSrc, nPointCount,
    3105             :                                 padfX, padfY, padfZ, panSuccess))
    3106       20196 :             ret = FALSE;
    3107             :     }
    3108             : 
    3109             :     /* -------------------------------------------------------------------- */
    3110             :     /*      Convert dst (src) georef coordinates back to pixel/line.        */
    3111             :     /* -------------------------------------------------------------------- */
    3112             :     {
    3113     2163140 :         const auto params = bDstToSrc ? psInfo->sSrcParams : psInfo->sDstParams;
    3114     2163140 :         const double *padfInvGeoTransform = params.adfInvGeoTransform;
    3115     2163140 :         void *pTransformArg = params.pTransformArg;
    3116     2163140 :         GDALTransformerFunc pTransformer = params.pTransformer;
    3117             : 
    3118     2163140 :         if (pTransformArg != nullptr)
    3119             :         {
    3120       51490 :             if (!pTransformer(pTransformArg, TRUE, nPointCount, padfX, padfY,
    3121             :                               padfZ, panSuccess))
    3122        1648 :                 ret = FALSE;
    3123             :         }
    3124             :         else
    3125             :         {
    3126    18946000 :             for (int i = 0; i < nPointCount; i++)
    3127             :             {
    3128    16834400 :                 if (!panSuccess[i])
    3129     3555620 :                     continue;
    3130             : 
    3131    13278800 :                 const double dfNewX = padfInvGeoTransform[0] +
    3132    13278800 :                                       padfX[i] * padfInvGeoTransform[1] +
    3133    13278800 :                                       padfY[i] * padfInvGeoTransform[2];
    3134    13278800 :                 const double dfNewY = padfInvGeoTransform[3] +
    3135    13278800 :                                       padfX[i] * padfInvGeoTransform[4] +
    3136    13278800 :                                       padfY[i] * padfInvGeoTransform[5];
    3137             : 
    3138    13278800 :                 padfX[i] = dfNewX;
    3139    13278800 :                 padfY[i] = dfNewY;
    3140             :             }
    3141             :         }
    3142             :     }
    3143             : 
    3144     2163140 :     return ret;
    3145             : }
    3146             : 
    3147             : /************************************************************************/
    3148             : /*              GDALTransformLonLatToDestGenImgProjTransformer()        */
    3149             : /************************************************************************/
    3150             : 
    3151        2752 : int GDALTransformLonLatToDestGenImgProjTransformer(void *hTransformArg,
    3152             :                                                    double *pdfX, double *pdfY)
    3153             : {
    3154        2752 :     GDALGenImgProjTransformInfo *psInfo =
    3155             :         static_cast<GDALGenImgProjTransformInfo *>(hTransformArg);
    3156             : 
    3157        2752 :     if (psInfo->pReprojectArg == nullptr ||
    3158        1564 :         psInfo->pReproject != GDALReprojectionTransform)
    3159        1192 :         return false;
    3160             : 
    3161        1560 :     GDALReprojectionTransformInfo *psReprojInfo =
    3162             :         static_cast<GDALReprojectionTransformInfo *>(psInfo->pReprojectArg);
    3163        3120 :     if (psReprojInfo->poForwardTransform == nullptr ||
    3164        1560 :         psReprojInfo->poForwardTransform->GetSourceCS() == nullptr)
    3165           2 :         return false;
    3166             : 
    3167        1558 :     double z = 0;
    3168        1558 :     int success = true;
    3169        1558 :     auto poSourceCRS = psReprojInfo->poForwardTransform->GetSourceCS();
    3170        2730 :     if (poSourceCRS->IsGeographic() &&
    3171        1172 :         std::fabs(poSourceCRS->GetAngularUnits() -
    3172        1172 :                   CPLAtof(SRS_UA_DEGREE_CONV)) < 1e-9)
    3173             :     {
    3174             :         // Optimization to avoid creating a OGRCoordinateTransformation
    3175        1170 :         OGRAxisOrientation eSourceFirstAxisOrient = OAO_Other;
    3176        1170 :         poSourceCRS->GetAxis(nullptr, 0, &eSourceFirstAxisOrient);
    3177        1170 :         const auto &mapping = poSourceCRS->GetDataAxisToSRSAxisMapping();
    3178        2340 :         if ((mapping[0] == 2 && eSourceFirstAxisOrient == OAO_East) ||
    3179        1170 :             (mapping[0] == 1 && eSourceFirstAxisOrient != OAO_East))
    3180             :         {
    3181           6 :             std::swap(*pdfX, *pdfY);
    3182             :         }
    3183             :     }
    3184             :     else
    3185             :     {
    3186             :         auto poLongLat =
    3187         388 :             std::unique_ptr<OGRSpatialReference>(poSourceCRS->CloneGeogCS());
    3188         388 :         if (poLongLat == nullptr)
    3189           0 :             return false;
    3190         388 :         poLongLat->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
    3191             : 
    3192             :         const bool bCurrentCheckWithInvertProj =
    3193         388 :             GetCurrentCheckWithInvertPROJ();
    3194         388 :         if (!bCurrentCheckWithInvertProj)
    3195         388 :             CPLSetThreadLocalConfigOption("CHECK_WITH_INVERT_PROJ", "YES");
    3196             :         auto poCT = std::unique_ptr<OGRCoordinateTransformation>(
    3197         388 :             OGRCreateCoordinateTransformation(poLongLat.get(), poSourceCRS));
    3198         388 :         if (!bCurrentCheckWithInvertProj)
    3199         388 :             CPLSetThreadLocalConfigOption("CHECK_WITH_INVERT_PROJ", nullptr);
    3200         388 :         if (poCT == nullptr)
    3201           0 :             return false;
    3202             : 
    3203         388 :         poCT->SetEmitErrors(false);
    3204         388 :         if (!poCT->Transform(1, pdfX, pdfY))
    3205           2 :             return false;
    3206             : 
    3207         386 :         if (!psInfo->pReproject(psInfo->pReprojectArg, false, 1, pdfX, pdfY, &z,
    3208         628 :                                 &success) ||
    3209         242 :             !success)
    3210             :         {
    3211         144 :             return false;
    3212             :         }
    3213             :     }
    3214             : 
    3215        1412 :     double *padfGeoTransform = psInfo->sDstParams.adfInvGeoTransform;
    3216        1412 :     void *pTransformArg = psInfo->sDstParams.pTransformArg;
    3217        1412 :     GDALTransformerFunc pTransformer = psInfo->sDstParams.pTransformer;
    3218        1412 :     if (pTransformArg != nullptr)
    3219             :     {
    3220           4 :         if (!pTransformer(pTransformArg, TRUE, 1, pdfX, pdfY, &z, &success) ||
    3221           0 :             !success)
    3222             :         {
    3223           4 :             return false;
    3224             :         }
    3225             :     }
    3226             :     else
    3227             :     {
    3228        1408 :         const double dfNewX = padfGeoTransform[0] +
    3229        1408 :                               pdfX[0] * padfGeoTransform[1] +
    3230        1408 :                               pdfY[0] * padfGeoTransform[2];
    3231        1408 :         const double dfNewY = padfGeoTransform[3] +
    3232        1408 :                               pdfX[0] * padfGeoTransform[4] +
    3233        1408 :                               pdfY[0] * padfGeoTransform[5];
    3234             : 
    3235        1408 :         pdfX[0] = dfNewX;
    3236        1408 :         pdfY[0] = dfNewY;
    3237             :     }
    3238             : 
    3239        1408 :     return true;
    3240             : }
    3241             : 
    3242             : /************************************************************************/
    3243             : /*                 GDALSerializeGenImgProjTransformer()                 */
    3244             : /************************************************************************/
    3245             : 
    3246          78 : static CPLXMLNode *GDALSerializeGenImgProjTransformer(void *pTransformArg)
    3247             : 
    3248             : {
    3249          78 :     GDALGenImgProjTransformInfo *psInfo =
    3250             :         static_cast<GDALGenImgProjTransformInfo *>(pTransformArg);
    3251             : 
    3252             :     CPLXMLNode *psTree =
    3253          78 :         CPLCreateXMLNode(nullptr, CXT_Element, "GenImgProjTransformer");
    3254             : 
    3255             :     const auto SerializePart =
    3256         455 :         [psTree](const char *pszPrefix, const GDALGenImgProjTransformPart &part)
    3257             :     {
    3258         156 :         char szWork[200] = {};
    3259             : 
    3260             :         /* ------------------------------------------------------------- */
    3261             :         /*      Handle transformation.                                   */
    3262             :         /* ------------------------------------------------------------- */
    3263         156 :         if (part.pTransformArg != nullptr)
    3264             :         {
    3265             :             CPLXMLNode *psTransformer =
    3266          13 :                 GDALSerializeTransformer(part.pTransformer, part.pTransformArg);
    3267          13 :             if (psTransformer != nullptr)
    3268             :             {
    3269          13 :                 CPLXMLNode *psTransformerContainer = CPLCreateXMLNode(
    3270             :                     psTree, CXT_Element,
    3271             :                     CPLSPrintf("%s%s", pszPrefix, psTransformer->pszValue));
    3272             : 
    3273          13 :                 CPLAddXMLChild(psTransformerContainer, psTransformer);
    3274             :             }
    3275             :         }
    3276             : 
    3277             :         /* ------------------------------------------------------------- */
    3278             :         /*      Handle geotransforms.                                    */
    3279             :         /* ------------------------------------------------------------- */
    3280             :         else
    3281             :         {
    3282         143 :             CPLsnprintf(szWork, sizeof(szWork),
    3283             :                         "%.17g,%.17g,%.17g,%.17g,%.17g,%.17g",
    3284         143 :                         part.adfGeoTransform[0], part.adfGeoTransform[1],
    3285         143 :                         part.adfGeoTransform[2], part.adfGeoTransform[3],
    3286         143 :                         part.adfGeoTransform[4], part.adfGeoTransform[5]);
    3287         143 :             CPLCreateXMLElementAndValue(
    3288             :                 psTree, CPLSPrintf("%sGeoTransform", pszPrefix), szWork);
    3289             : 
    3290         143 :             CPLsnprintf(szWork, sizeof(szWork),
    3291             :                         "%.17g,%.17g,%.17g,%.17g,%.17g,%.17g",
    3292         143 :                         part.adfInvGeoTransform[0], part.adfInvGeoTransform[1],
    3293         143 :                         part.adfInvGeoTransform[2], part.adfInvGeoTransform[3],
    3294         143 :                         part.adfInvGeoTransform[4], part.adfInvGeoTransform[5]);
    3295         143 :             CPLCreateXMLElementAndValue(
    3296             :                 psTree, CPLSPrintf("%sInvGeoTransform", pszPrefix), szWork);
    3297             :         }
    3298         156 :     };
    3299             : 
    3300          78 :     SerializePart("Src", psInfo->sSrcParams);
    3301          78 :     SerializePart("Dst", psInfo->sDstParams);
    3302             : 
    3303             :     /* -------------------------------------------------------------------- */
    3304             :     /*      Do we have a reprojection transformer?                          */
    3305             :     /* -------------------------------------------------------------------- */
    3306          78 :     if (psInfo->pReprojectArg != nullptr)
    3307             :     {
    3308             : 
    3309             :         CPLXMLNode *psTransformerContainer =
    3310          51 :             CPLCreateXMLNode(psTree, CXT_Element, "ReprojectTransformer");
    3311             : 
    3312             :         CPLXMLNode *psTransformer =
    3313          51 :             GDALSerializeTransformer(psInfo->pReproject, psInfo->pReprojectArg);
    3314          51 :         if (psTransformer != nullptr)
    3315          51 :             CPLAddXMLChild(psTransformerContainer, psTransformer);
    3316             :     }
    3317             : 
    3318          78 :     return psTree;
    3319             : }
    3320             : 
    3321             : /************************************************************************/
    3322             : /*                    GDALDeserializeGeoTransform()                     */
    3323             : /************************************************************************/
    3324             : 
    3325         739 : static void GDALDeserializeGeoTransform(const char *pszGT,
    3326             :                                         double adfGeoTransform[6])
    3327             : {
    3328         739 :     CPLsscanf(pszGT, "%lf,%lf,%lf,%lf,%lf,%lf", adfGeoTransform + 0,
    3329             :               adfGeoTransform + 1, adfGeoTransform + 2, adfGeoTransform + 3,
    3330             :               adfGeoTransform + 4, adfGeoTransform + 5);
    3331         739 : }
    3332             : 
    3333             : /************************************************************************/
    3334             : /*                GDALDeserializeGenImgProjTransformer()                */
    3335             : /************************************************************************/
    3336             : 
    3337         199 : void *GDALDeserializeGenImgProjTransformer(CPLXMLNode *psTree)
    3338             : 
    3339             : {
    3340             :     /* -------------------------------------------------------------------- */
    3341             :     /*      Initialize the transform info.                                  */
    3342             :     /* -------------------------------------------------------------------- */
    3343             :     GDALGenImgProjTransformInfo *psInfo =
    3344         199 :         GDALCreateGenImgProjTransformerInternal();
    3345             : 
    3346             :     const auto DeserializePart =
    3347        1194 :         [psTree](const char *pszPrefix, GDALGenImgProjTransformPart &part)
    3348             :     {
    3349             :         /* ----------------------------------------------------------------- */
    3350             :         /*      Geotransform                                                 */
    3351             :         /* ----------------------------------------------------------------- */
    3352         398 :         if (const auto psGTNode =
    3353         398 :                 CPLGetXMLNode(psTree, CPLSPrintf("%sGeoTransform", pszPrefix)))
    3354             :         {
    3355         382 :             GDALDeserializeGeoTransform(CPLGetXMLValue(psGTNode, "", ""),
    3356         382 :                                         part.adfGeoTransform);
    3357             : 
    3358         382 :             if (const auto psInvGTNode = CPLGetXMLNode(
    3359             :                     psTree, CPLSPrintf("%sInvGeoTransform", pszPrefix)))
    3360             :             {
    3361         357 :                 GDALDeserializeGeoTransform(CPLGetXMLValue(psInvGTNode, "", ""),
    3362         357 :                                             part.adfInvGeoTransform);
    3363             :             }
    3364             :             else
    3365             :             {
    3366          25 :                 if (!GDALInvGeoTransform(part.adfGeoTransform,
    3367          25 :                                          part.adfInvGeoTransform))
    3368             :                 {
    3369           0 :                     CPLError(CE_Failure, CPLE_AppDefined,
    3370             :                              "Cannot invert geotransform");
    3371             :                 }
    3372             :             }
    3373             :         }
    3374             : 
    3375             :         /* ---------------------------------------------------------------- */
    3376             :         /*      Transform                                                   */
    3377             :         /* ---------------------------------------------------------------- */
    3378             :         else
    3379             :         {
    3380          16 :             for (CPLXMLNode *psIter = psTree->psChild; psIter != nullptr;
    3381           0 :                  psIter = psIter->psNext)
    3382             :             {
    3383          16 :                 if (psIter->eType == CXT_Element &&
    3384          16 :                     STARTS_WITH_CI(psIter->pszValue, pszPrefix))
    3385             :                 {
    3386          16 :                     GDALDeserializeTransformer(psIter->psChild,
    3387             :                                                &part.pTransformer,
    3388             :                                                &part.pTransformArg);
    3389          16 :                     break;
    3390             :                 }
    3391             :             }
    3392             :         }
    3393         597 :     };
    3394             : 
    3395         199 :     DeserializePart("Src", psInfo->sSrcParams);
    3396         199 :     DeserializePart("Dst", psInfo->sDstParams);
    3397             : 
    3398             :     /* -------------------------------------------------------------------- */
    3399             :     /*      Reproject transformer                                           */
    3400             :     /* -------------------------------------------------------------------- */
    3401         199 :     CPLXMLNode *psSubtree = CPLGetXMLNode(psTree, "ReprojectTransformer");
    3402         199 :     if (psSubtree != nullptr && psSubtree->psChild != nullptr)
    3403             :     {
    3404          98 :         GDALDeserializeTransformer(psSubtree->psChild, &psInfo->pReproject,
    3405             :                                    &psInfo->pReprojectArg);
    3406             :     }
    3407             : 
    3408         199 :     return psInfo;
    3409             : }
    3410             : 
    3411             : /************************************************************************/
    3412             : /*                 GDALCreateReprojectionTransformer()                  */
    3413             : /************************************************************************/
    3414             : 
    3415             : /**
    3416             :  * Create reprojection transformer.
    3417             :  *
    3418             :  * Creates a callback data structure suitable for use with
    3419             :  * GDALReprojectionTransformation() to represent a transformation from
    3420             :  * one geographic or projected coordinate system to another.  On input
    3421             :  * the coordinate systems are described in OpenGIS WKT format.
    3422             :  *
    3423             :  * Internally the OGRCoordinateTransformation object is used to implement
    3424             :  * the reprojection.
    3425             :  *
    3426             :  * @param pszSrcWKT the coordinate system for the source coordinate system.
    3427             :  * @param pszDstWKT the coordinate system for the destination coordinate
    3428             :  * system.
    3429             :  *
    3430             :  * @return Handle for use with GDALReprojectionTransform(), or NULL if the
    3431             :  * system fails to initialize the reprojection.
    3432             :  **/
    3433             : 
    3434           0 : void *GDALCreateReprojectionTransformer(const char *pszSrcWKT,
    3435             :                                         const char *pszDstWKT)
    3436             : 
    3437             : {
    3438             :     /* -------------------------------------------------------------------- */
    3439             :     /*      Ingest the SRS definitions.                                     */
    3440             :     /* -------------------------------------------------------------------- */
    3441           0 :     OGRSpatialReference oSrcSRS;
    3442           0 :     oSrcSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
    3443           0 :     if (oSrcSRS.importFromWkt(pszSrcWKT) != OGRERR_NONE)
    3444             :     {
    3445           0 :         CPLError(CE_Failure, CPLE_AppDefined,
    3446             :                  "Failed to import coordinate system `%s'.", pszSrcWKT);
    3447           0 :         return nullptr;
    3448             :     }
    3449             : 
    3450           0 :     OGRSpatialReference oDstSRS;
    3451           0 :     oDstSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
    3452           0 :     if (oDstSRS.importFromWkt(pszDstWKT) != OGRERR_NONE)
    3453             :     {
    3454           0 :         CPLError(CE_Failure, CPLE_AppDefined,
    3455             :                  "Failed to import coordinate system `%s'.", pszSrcWKT);
    3456           0 :         return nullptr;
    3457             :     }
    3458             : 
    3459           0 :     return GDALCreateReprojectionTransformerEx(
    3460             :         OGRSpatialReference::ToHandle(&oSrcSRS),
    3461           0 :         OGRSpatialReference::ToHandle(&oDstSRS), nullptr);
    3462             : }
    3463             : 
    3464             : /************************************************************************/
    3465             : /*                 GDALCreateReprojectionTransformerEx()                */
    3466             : /************************************************************************/
    3467             : 
    3468             : /**
    3469             :  * Create reprojection transformer.
    3470             :  *
    3471             :  * Creates a callback data structure suitable for use with
    3472             :  * GDALReprojectionTransformation() to represent a transformation from
    3473             :  * one geographic or projected coordinate system to another.
    3474             :  *
    3475             :  * Internally the OGRCoordinateTransformation object is used to implement
    3476             :  * the reprojection.
    3477             :  *
    3478             :  * @param hSrcSRS the coordinate system for the source coordinate system.
    3479             :  * @param hDstSRS the coordinate system for the destination coordinate
    3480             :  * system.
    3481             :  * @param papszOptions NULL-terminated list of options, or NULL. Currently
    3482             :  * supported options are:
    3483             :  * <ul>
    3484             :  * <li>AREA_OF_INTEREST=west_long,south_lat,east_long,north_lat: Values in
    3485             :  * degrees. longitudes in [-180,180], latitudes in [-90,90].</li>
    3486             :  * <li>COORDINATE_OPERATION=string: PROJ or WKT string representing a
    3487             :  * coordinate operation, overriding the default computed transformation.</li>
    3488             :  * <li>COORDINATE_EPOCH=decimal_year: Coordinate epoch, expressed as a
    3489             :  * decimal year. Useful for time-dependent coordinate operations.</li>
    3490             :  * <li> SRC_COORDINATE_EPOCH: (GDAL &gt;= 3.4) Coordinate epoch of source CRS,
    3491             :  * expressed as a decimal year. Useful for time-dependent coordinate
    3492             :  *operations.</li>
    3493             :  * <li> DST_COORDINATE_EPOCH: (GDAL &gt;= 3.4) Coordinate epoch
    3494             :  *of target CRS, expressed as a decimal year. Useful for time-dependent
    3495             :  *coordinate operations.</li>
    3496             :  * <li> ALLOW_BALLPARK=YES/NO: (GDAL &gt;= 3.11) Whether ballpark coordinate
    3497             :  * operations are allowed. Defaults to YES.</li>
    3498             :  * <li> ONLY_BEST=YES/NO/AUTO: (GDAL &gt;= 3.11) By default (at least in the
    3499             :  * PROJ 9.x series), PROJ may use coordinate
    3500             :  * operations that are not the "best" if resources (typically grids) needed
    3501             :  * to use them are missing. It will then fallback to other coordinate operations
    3502             :  * that have a lesser accuracy, for example using Helmert transformations,
    3503             :  * or in the absence of such operations, to ones with potential very rought
    3504             :  * accuracy, using "ballpark" transformations
    3505             :  * (see https://proj.org/glossary.html).
    3506             :  * When calling this method with YES, PROJ will only consider the
    3507             :  * "best" operation, and error out (at Transform() time) if they cannot be
    3508             :  * used.
    3509             :  * This method may be used together with ALLOW_BALLPARK=NO to
    3510             :  * only allow best operations that have a known accuracy.
    3511             :  * Note that this method has no effect on PROJ versions before 9.2.
    3512             :  * The default value for this option can be also set with the
    3513             :  * PROJ_ONLY_BEST_DEFAULT environment variable, or with the "only_best_default"
    3514             :  * setting of proj.ini. Calling SetOnlyBest() overrides such default value.</li>
    3515             :  * </ul>
    3516             :  *
    3517             :  * @return Handle for use with GDALReprojectionTransform(), or NULL if the
    3518             :  * system fails to initialize the reprojection.
    3519             :  *
    3520             :  * @since GDAL 3.0
    3521             :  **/
    3522             : 
    3523        1069 : void *GDALCreateReprojectionTransformerEx(OGRSpatialReferenceH hSrcSRS,
    3524             :                                           OGRSpatialReferenceH hDstSRS,
    3525             :                                           const char *const *papszOptions)
    3526             : {
    3527        1069 :     OGRSpatialReference *poSrcSRS = OGRSpatialReference::FromHandle(hSrcSRS);
    3528        1069 :     OGRSpatialReference *poDstSRS = OGRSpatialReference::FromHandle(hDstSRS);
    3529             : 
    3530             :     /* -------------------------------------------------------------------- */
    3531             :     /*      Build the forward coordinate transformation.                    */
    3532             :     /* -------------------------------------------------------------------- */
    3533        1069 :     double dfWestLongitudeDeg = 0.0;
    3534        1069 :     double dfSouthLatitudeDeg = 0.0;
    3535        1069 :     double dfEastLongitudeDeg = 0.0;
    3536        1069 :     double dfNorthLatitudeDeg = 0.0;
    3537        1069 :     const char *pszBBOX = CSLFetchNameValue(papszOptions, "AREA_OF_INTEREST");
    3538        1069 :     if (pszBBOX)
    3539             :     {
    3540        1002 :         char **papszTokens = CSLTokenizeString2(pszBBOX, ",", 0);
    3541        1002 :         if (CSLCount(papszTokens) == 4)
    3542             :         {
    3543        1002 :             dfWestLongitudeDeg = CPLAtof(papszTokens[0]);
    3544        1002 :             dfSouthLatitudeDeg = CPLAtof(papszTokens[1]);
    3545        1002 :             dfEastLongitudeDeg = CPLAtof(papszTokens[2]);
    3546        1002 :             dfNorthLatitudeDeg = CPLAtof(papszTokens[3]);
    3547             :         }
    3548        1002 :         CSLDestroy(papszTokens);
    3549             :     }
    3550        1069 :     const char *pszCO = CSLFetchNameValue(papszOptions, "COORDINATE_OPERATION");
    3551             : 
    3552        2138 :     OGRCoordinateTransformationOptions optionsFwd;
    3553        1069 :     if (!(dfWestLongitudeDeg == 0.0 && dfSouthLatitudeDeg == 0.0 &&
    3554             :           dfEastLongitudeDeg == 0.0 && dfNorthLatitudeDeg == 0.0))
    3555             :     {
    3556        1002 :         optionsFwd.SetAreaOfInterest(dfWestLongitudeDeg, dfSouthLatitudeDeg,
    3557             :                                      dfEastLongitudeDeg, dfNorthLatitudeDeg);
    3558             :     }
    3559        1069 :     if (pszCO)
    3560             :     {
    3561           7 :         optionsFwd.SetCoordinateOperation(pszCO, false);
    3562             :     }
    3563             : 
    3564        1069 :     const char *pszCENTER_LONG = CSLFetchNameValue(papszOptions, "CENTER_LONG");
    3565        1069 :     if (pszCENTER_LONG)
    3566             :     {
    3567         703 :         optionsFwd.SetSourceCenterLong(CPLAtof(pszCENTER_LONG));
    3568             :     }
    3569             : 
    3570        1069 :     optionsFwd.SetBallparkAllowed(CPLTestBool(
    3571             :         CSLFetchNameValueDef(papszOptions, "ALLOW_BALLPARK", "YES")));
    3572             : 
    3573             :     const char *pszOnlyBest =
    3574        1069 :         CSLFetchNameValueDef(papszOptions, "ONLY_BEST", "AUTO");
    3575        1069 :     if (!EQUAL(pszOnlyBest, "AUTO"))
    3576             :     {
    3577           1 :         optionsFwd.SetOnlyBest(CPLTestBool(pszOnlyBest));
    3578             :     }
    3579             : 
    3580             :     OGRCoordinateTransformation *poForwardTransform =
    3581        1069 :         OGRCreateCoordinateTransformation(poSrcSRS, poDstSRS, optionsFwd);
    3582             : 
    3583        1069 :     if (poForwardTransform == nullptr)
    3584             :         // OGRCreateCoordinateTransformation() will report errors on its own.
    3585           5 :         return nullptr;
    3586             : 
    3587        1064 :     poForwardTransform->SetEmitErrors(false);
    3588             : 
    3589             :     /* -------------------------------------------------------------------- */
    3590             :     /*      Create a structure to hold the transform info, and also         */
    3591             :     /*      build reverse transform.  We assume that if the forward         */
    3592             :     /*      transform can be created, then so can the reverse one.          */
    3593             :     /* -------------------------------------------------------------------- */
    3594        1064 :     GDALReprojectionTransformInfo *psInfo = new GDALReprojectionTransformInfo();
    3595             : 
    3596        1064 :     psInfo->papszOptions = CSLDuplicate(papszOptions);
    3597        1064 :     psInfo->poForwardTransform = poForwardTransform;
    3598        1064 :     psInfo->dfTime = CPLAtof(CSLFetchNameValueDef(
    3599             :         papszOptions, "COORDINATE_EPOCH",
    3600             :         CSLFetchNameValueDef(
    3601             :             papszOptions, "DST_COORDINATE_EPOCH",
    3602             :             CSLFetchNameValueDef(papszOptions, "SRC_COORDINATE_EPOCH", "0"))));
    3603        1064 :     psInfo->poReverseTransform = poForwardTransform->GetInverse();
    3604             : 
    3605        1064 :     if (psInfo->poReverseTransform)
    3606        1064 :         psInfo->poReverseTransform->SetEmitErrors(false);
    3607             : 
    3608        1064 :     memcpy(psInfo->sTI.abySignature, GDAL_GTI2_SIGNATURE,
    3609             :            strlen(GDAL_GTI2_SIGNATURE));
    3610        1064 :     psInfo->sTI.pszClassName = GDAL_REPROJECTION_TRANSFORMER_CLASS_NAME;
    3611        1064 :     psInfo->sTI.pfnTransform = GDALReprojectionTransform;
    3612        1064 :     psInfo->sTI.pfnCleanup = GDALDestroyReprojectionTransformer;
    3613        1064 :     psInfo->sTI.pfnSerialize = GDALSerializeReprojectionTransformer;
    3614             : 
    3615        1064 :     return psInfo;
    3616             : }
    3617             : 
    3618             : /************************************************************************/
    3619             : /*                 GDALDestroyReprojectionTransformer()                 */
    3620             : /************************************************************************/
    3621             : 
    3622             : /**
    3623             :  * Destroy reprojection transformation.
    3624             :  *
    3625             :  * @param pTransformArg the transformation handle returned by
    3626             :  * GDALCreateReprojectionTransformer().
    3627             :  */
    3628             : 
    3629        1064 : void GDALDestroyReprojectionTransformer(void *pTransformArg)
    3630             : 
    3631             : {
    3632        1064 :     if (pTransformArg == nullptr)
    3633           0 :         return;
    3634             : 
    3635        1064 :     GDALReprojectionTransformInfo *psInfo =
    3636             :         static_cast<GDALReprojectionTransformInfo *>(pTransformArg);
    3637             : 
    3638        1064 :     if (psInfo->poForwardTransform)
    3639        1064 :         OGRCoordinateTransformation::DestroyCT(psInfo->poForwardTransform);
    3640             : 
    3641        1064 :     if (psInfo->poReverseTransform)
    3642        1064 :         OGRCoordinateTransformation::DestroyCT(psInfo->poReverseTransform);
    3643             : 
    3644        1064 :     CSLDestroy(psInfo->papszOptions);
    3645             : 
    3646        1064 :     delete psInfo;
    3647             : }
    3648             : 
    3649             : /************************************************************************/
    3650             : /*                     GDALReprojectionTransform()                      */
    3651             : /************************************************************************/
    3652             : 
    3653             : /**
    3654             :  * Perform reprojection transformation.
    3655             :  *
    3656             :  * Actually performs the reprojection transformation described in
    3657             :  * GDALCreateReprojectionTransformer().  This function matches the
    3658             :  * GDALTransformerFunc() signature.  Details of the arguments are described
    3659             :  * there.
    3660             :  */
    3661             : 
    3662     1813210 : int GDALReprojectionTransform(void *pTransformArg, int bDstToSrc,
    3663             :                               int nPointCount, double *padfX, double *padfY,
    3664             :                               double *padfZ, int *panSuccess)
    3665             : 
    3666             : {
    3667     1813210 :     GDALReprojectionTransformInfo *psInfo =
    3668             :         static_cast<GDALReprojectionTransformInfo *>(pTransformArg);
    3669             :     int bSuccess;
    3670             : 
    3671     1813210 :     std::vector<double> adfTime;
    3672     1813200 :     double *padfT = nullptr;
    3673     1813200 :     if (psInfo->dfTime != 0.0 && nPointCount > 0)
    3674             :     {
    3675           1 :         adfTime.resize(nPointCount, psInfo->dfTime);
    3676           1 :         padfT = &adfTime[0];
    3677             :     }
    3678             : 
    3679     1813200 :     if (bDstToSrc)
    3680             :     {
    3681     1562390 :         if (psInfo->poReverseTransform == nullptr)
    3682             :         {
    3683           0 :             CPLError(
    3684             :                 CE_Failure, CPLE_AppDefined,
    3685             :                 "Inverse coordinate transformation cannot be instantiated");
    3686           0 :             if (panSuccess)
    3687             :             {
    3688           0 :                 for (int i = 0; i < nPointCount; i++)
    3689           0 :                     panSuccess[i] = FALSE;
    3690             :             }
    3691           0 :             bSuccess = false;
    3692             :         }
    3693             :         else
    3694             :         {
    3695     1562400 :             bSuccess = psInfo->poReverseTransform->Transform(
    3696     1562390 :                 nPointCount, padfX, padfY, padfZ, padfT, panSuccess);
    3697             :         }
    3698             :     }
    3699             :     else
    3700      250810 :         bSuccess = psInfo->poForwardTransform->Transform(
    3701      250810 :             nPointCount, padfX, padfY, padfZ, padfT, panSuccess);
    3702             : 
    3703     3626420 :     return bSuccess;
    3704             : }
    3705             : 
    3706             : /************************************************************************/
    3707             : /*                GDALSerializeReprojectionTransformer()                */
    3708             : /************************************************************************/
    3709             : 
    3710         138 : static CPLXMLNode *GDALSerializeReprojectionTransformer(void *pTransformArg)
    3711             : 
    3712             : {
    3713             :     CPLXMLNode *psTree;
    3714         138 :     GDALReprojectionTransformInfo *psInfo =
    3715             :         static_cast<GDALReprojectionTransformInfo *>(pTransformArg);
    3716             : 
    3717         138 :     psTree = CPLCreateXMLNode(nullptr, CXT_Element, "ReprojectionTransformer");
    3718             : 
    3719             :     /* -------------------------------------------------------------------- */
    3720             :     /*      Handle SourceCS.                                                */
    3721             :     /* -------------------------------------------------------------------- */
    3722         276 :     const auto ExportToWkt = [](const OGRSpatialReference *poSRS)
    3723             :     {
    3724             :         // Try first in WKT1 for backward compat
    3725             :         {
    3726         276 :             char *pszWKT = nullptr;
    3727         276 :             const char *const apszOptions[] = {"FORMAT=WKT1", nullptr};
    3728         276 :             CPLErrorHandlerPusher oHandler(CPLQuietErrorHandler);
    3729         276 :             CPLErrorStateBackuper oBackuper;
    3730         276 :             if (poSRS->exportToWkt(&pszWKT, apszOptions) == OGRERR_NONE)
    3731             :             {
    3732         550 :                 std::string osRet(pszWKT);
    3733         275 :                 CPLFree(pszWKT);
    3734         275 :                 return osRet;
    3735             :             }
    3736           1 :             CPLFree(pszWKT);
    3737             :         }
    3738             : 
    3739           1 :         char *pszWKT = nullptr;
    3740           1 :         const char *const apszOptions[] = {"FORMAT=WKT2_2019", nullptr};
    3741           1 :         if (poSRS->exportToWkt(&pszWKT, apszOptions) == OGRERR_NONE)
    3742             :         {
    3743           2 :             std::string osRet(pszWKT);
    3744           1 :             CPLFree(pszWKT);
    3745           1 :             return osRet;
    3746             :         }
    3747           0 :         CPLFree(pszWKT);
    3748           0 :         return std::string();
    3749             :     };
    3750             : 
    3751         138 :     auto poSRS = psInfo->poForwardTransform->GetSourceCS();
    3752         138 :     if (poSRS)
    3753             :     {
    3754         276 :         const auto osWKT = ExportToWkt(poSRS);
    3755         138 :         CPLCreateXMLElementAndValue(psTree, "SourceSRS", osWKT.c_str());
    3756             :     }
    3757             : 
    3758             :     /* -------------------------------------------------------------------- */
    3759             :     /*      Handle DestinationCS.                                           */
    3760             :     /* -------------------------------------------------------------------- */
    3761         138 :     poSRS = psInfo->poForwardTransform->GetTargetCS();
    3762         138 :     if (poSRS)
    3763             :     {
    3764         276 :         const auto osWKT = ExportToWkt(poSRS);
    3765         138 :         CPLCreateXMLElementAndValue(psTree, "TargetSRS", osWKT.c_str());
    3766             :     }
    3767             : 
    3768             :     /* -------------------------------------------------------------------- */
    3769             :     /*      Serialize options.                                              */
    3770             :     /* -------------------------------------------------------------------- */
    3771         138 :     if (psInfo->papszOptions)
    3772             :     {
    3773             :         CPLXMLNode *psOptions =
    3774         125 :             CPLCreateXMLNode(psTree, CXT_Element, "Options");
    3775         308 :         for (auto iter = psInfo->papszOptions; *iter != nullptr; ++iter)
    3776             :         {
    3777         183 :             char *pszKey = nullptr;
    3778         183 :             const char *pszValue = CPLParseNameValue(*iter, &pszKey);
    3779         183 :             if (pszKey && pszValue)
    3780             :             {
    3781             :                 auto elt =
    3782         183 :                     CPLCreateXMLElementAndValue(psOptions, "Option", pszValue);
    3783         183 :                 CPLAddXMLAttributeAndValue(elt, "key", pszKey);
    3784             :             }
    3785         183 :             CPLFree(pszKey);
    3786             :         }
    3787             :     }
    3788             : 
    3789         138 :     return psTree;
    3790             : }
    3791             : 
    3792             : /************************************************************************/
    3793             : /*               GDALDeserializeReprojectionTransformer()               */
    3794             : /************************************************************************/
    3795             : 
    3796         185 : static void *GDALDeserializeReprojectionTransformer(CPLXMLNode *psTree)
    3797             : 
    3798             : {
    3799         185 :     const char *pszSourceSRS = CPLGetXMLValue(psTree, "SourceSRS", nullptr);
    3800         185 :     const char *pszTargetSRS = CPLGetXMLValue(psTree, "TargetSRS", nullptr);
    3801         185 :     char *pszSourceWKT = nullptr, *pszTargetWKT = nullptr;
    3802         185 :     void *pResult = nullptr;
    3803             : 
    3804         370 :     OGRSpatialReference oSrcSRS;
    3805         370 :     OGRSpatialReference oDstSRS;
    3806             : 
    3807         185 :     oSrcSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
    3808         185 :     oDstSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
    3809         185 :     if (pszSourceSRS != nullptr)
    3810             :     {
    3811         185 :         oSrcSRS.SetFromUserInput(pszSourceSRS);
    3812             :     }
    3813             : 
    3814         185 :     if (pszTargetSRS != nullptr)
    3815             :     {
    3816         185 :         oDstSRS.SetFromUserInput(pszTargetSRS);
    3817             :     }
    3818             : 
    3819         185 :     CPLStringList aosList;
    3820         185 :     const CPLXMLNode *psOptions = CPLGetXMLNode(psTree, "Options");
    3821         185 :     if (psOptions)
    3822             :     {
    3823         376 :         for (auto iter = psOptions->psChild; iter; iter = iter->psNext)
    3824             :         {
    3825         227 :             if (iter->eType == CXT_Element &&
    3826         227 :                 strcmp(iter->pszValue, "Option") == 0)
    3827             :             {
    3828         227 :                 const char *pszKey = CPLGetXMLValue(iter, "key", nullptr);
    3829         227 :                 const char *pszValue = CPLGetXMLValue(iter, nullptr, nullptr);
    3830         227 :                 if (pszKey && pszValue)
    3831             :                 {
    3832         227 :                     aosList.SetNameValue(pszKey, pszValue);
    3833             :                 }
    3834             :             }
    3835             :         }
    3836             :     }
    3837             : 
    3838         185 :     pResult = GDALCreateReprojectionTransformerEx(
    3839         185 :         !oSrcSRS.IsEmpty() ? OGRSpatialReference::ToHandle(&oSrcSRS) : nullptr,
    3840         185 :         !oDstSRS.IsEmpty() ? OGRSpatialReference::ToHandle(&oDstSRS) : nullptr,
    3841         185 :         aosList.List());
    3842             : 
    3843         185 :     CPLFree(pszSourceWKT);
    3844         185 :     CPLFree(pszTargetWKT);
    3845             : 
    3846         370 :     return pResult;
    3847             : }
    3848             : 
    3849             : /************************************************************************/
    3850             : /* ==================================================================== */
    3851             : /*      Approximate transformer.                                        */
    3852             : /* ==================================================================== */
    3853             : /************************************************************************/
    3854             : 
    3855             : /************************************************************************/
    3856             : /*                  GDALCreateSimilarApproxTransformer()                */
    3857             : /************************************************************************/
    3858             : 
    3859          25 : static void *GDALCreateSimilarApproxTransformer(void *hTransformArg,
    3860             :                                                 double dfSrcRatioX,
    3861             :                                                 double dfSrcRatioY)
    3862             : {
    3863          25 :     VALIDATE_POINTER1(hTransformArg, "GDALCreateSimilarApproxTransformer",
    3864             :                       nullptr);
    3865             : 
    3866          25 :     GDALApproxTransformInfo *psInfo =
    3867             :         static_cast<GDALApproxTransformInfo *>(hTransformArg);
    3868             : 
    3869          25 :     void *pBaseCBData = GDALCreateSimilarTransformer(psInfo->pBaseCBData,
    3870             :                                                      dfSrcRatioX, dfSrcRatioY);
    3871          25 :     if (pBaseCBData == nullptr)
    3872             :     {
    3873           0 :         return nullptr;
    3874             :     }
    3875             : 
    3876             :     GDALApproxTransformInfo *psClonedInfo =
    3877          25 :         static_cast<GDALApproxTransformInfo *>(GDALCreateApproxTransformer2(
    3878             :             psInfo->pfnBaseTransformer, pBaseCBData, psInfo->dfMaxErrorForward,
    3879             :             psInfo->dfMaxErrorReverse));
    3880          25 :     psClonedInfo->bOwnSubtransformer = TRUE;
    3881             : 
    3882          25 :     return psClonedInfo;
    3883             : }
    3884             : 
    3885             : /************************************************************************/
    3886             : /*                   GDALSerializeApproxTransformer()                   */
    3887             : /************************************************************************/
    3888             : 
    3889          58 : static CPLXMLNode *GDALSerializeApproxTransformer(void *pTransformArg)
    3890             : 
    3891             : {
    3892             :     CPLXMLNode *psTree;
    3893          58 :     GDALApproxTransformInfo *psInfo =
    3894             :         static_cast<GDALApproxTransformInfo *>(pTransformArg);
    3895             : 
    3896          58 :     psTree = CPLCreateXMLNode(nullptr, CXT_Element, "ApproxTransformer");
    3897             : 
    3898             :     /* -------------------------------------------------------------------- */
    3899             :     /*      Attach max error.                                               */
    3900             :     /* -------------------------------------------------------------------- */
    3901          58 :     if (psInfo->dfMaxErrorForward == psInfo->dfMaxErrorReverse)
    3902             :     {
    3903          56 :         CPLCreateXMLElementAndValue(
    3904             :             psTree, "MaxError",
    3905         112 :             CPLString().Printf("%g", psInfo->dfMaxErrorForward));
    3906             :     }
    3907             :     else
    3908             :     {
    3909           2 :         CPLCreateXMLElementAndValue(
    3910             :             psTree, "MaxErrorForward",
    3911           4 :             CPLString().Printf("%g", psInfo->dfMaxErrorForward));
    3912           2 :         CPLCreateXMLElementAndValue(
    3913             :             psTree, "MaxErrorReverse",
    3914           4 :             CPLString().Printf("%g", psInfo->dfMaxErrorReverse));
    3915             :     }
    3916             : 
    3917             :     /* -------------------------------------------------------------------- */
    3918             :     /*      Capture underlying transformer.                                 */
    3919             :     /* -------------------------------------------------------------------- */
    3920             :     CPLXMLNode *psTransformerContainer =
    3921          58 :         CPLCreateXMLNode(psTree, CXT_Element, "BaseTransformer");
    3922             : 
    3923          58 :     CPLXMLNode *psTransformer = GDALSerializeTransformer(
    3924             :         psInfo->pfnBaseTransformer, psInfo->pBaseCBData);
    3925          58 :     if (psTransformer != nullptr)
    3926          58 :         CPLAddXMLChild(psTransformerContainer, psTransformer);
    3927             : 
    3928          58 :     return psTree;
    3929             : }
    3930             : 
    3931             : /************************************************************************/
    3932             : /*                    GDALCreateApproxTransformer()                     */
    3933             : /************************************************************************/
    3934             : 
    3935             : /**
    3936             :  * Create an approximating transformer.
    3937             :  *
    3938             :  * This function creates a context for an approximated transformer.  Basically
    3939             :  * a high precision transformer is supplied as input and internally linear
    3940             :  * approximations are computed to generate results to within a defined
    3941             :  * precision.
    3942             :  *
    3943             :  * The approximation is actually done at the point where GDALApproxTransform()
    3944             :  * calls are made, and depend on the assumption that they are roughly linear.
    3945             :  * The first and last point passed in must be the extreme values and the
    3946             :  * intermediate values should describe a curve between the end points.  The
    3947             :  * approximator transforms and centers using the approximate transformer, and
    3948             :  * then compares the true middle transformed value to a linear approximation
    3949             :  * based on the end points.  If the error is within the supplied threshold then
    3950             :  * the end points are used to linearly approximate all the values otherwise the
    3951             :  * input points are split into two smaller sets, and the function is recursively
    3952             :  * called until a sufficiently small set of points is found that the linear
    3953             :  * approximation is OK, or that all the points are exactly computed.
    3954             :  *
    3955             :  * This function is very suitable for approximating transformation results
    3956             :  * from output pixel/line space to input coordinates for warpers that operate
    3957             :  * on one input scanline at a time.  Care should be taken using it in other
    3958             :  * circumstances as little internal validation is done in order to keep things
    3959             :  * fast.
    3960             :  *
    3961             :  * @param pfnBaseTransformer the high precision transformer which should be
    3962             :  * approximated.
    3963             :  * @param pBaseTransformArg the callback argument for the high precision
    3964             :  * transformer.
    3965             :  * @param dfMaxError the maximum cartesian error in the "output" space that
    3966             :  * is to be accepted in the linear approximation, evaluated as a Manhattan
    3967             :  * distance.
    3968             :  *
    3969             :  * @return callback pointer suitable for use with GDALApproxTransform().  It
    3970             :  * should be deallocated with GDALDestroyApproxTransformer().
    3971             :  */
    3972             : 
    3973        1039 : void *GDALCreateApproxTransformer(GDALTransformerFunc pfnBaseTransformer,
    3974             :                                   void *pBaseTransformArg, double dfMaxError)
    3975             : 
    3976             : {
    3977        1039 :     return GDALCreateApproxTransformer2(pfnBaseTransformer, pBaseTransformArg,
    3978        1039 :                                         dfMaxError, dfMaxError);
    3979             : }
    3980             : 
    3981             : static void *
    3982        1217 : GDALCreateApproxTransformer2(GDALTransformerFunc pfnBaseTransformer,
    3983             :                              void *pBaseTransformArg, double dfMaxErrorForward,
    3984             :                              double dfMaxErrorReverse)
    3985             : 
    3986             : {
    3987        1217 :     GDALApproxTransformInfo *psATInfo = new GDALApproxTransformInfo;
    3988        1217 :     psATInfo->pfnBaseTransformer = pfnBaseTransformer;
    3989        1217 :     psATInfo->pBaseCBData = pBaseTransformArg;
    3990        1217 :     psATInfo->dfMaxErrorForward = dfMaxErrorForward;
    3991        1217 :     psATInfo->dfMaxErrorReverse = dfMaxErrorReverse;
    3992        1217 :     psATInfo->bOwnSubtransformer = FALSE;
    3993             : 
    3994        1217 :     memcpy(psATInfo->sTI.abySignature, GDAL_GTI2_SIGNATURE,
    3995             :            strlen(GDAL_GTI2_SIGNATURE));
    3996        1217 :     psATInfo->sTI.pszClassName = GDAL_APPROX_TRANSFORMER_CLASS_NAME;
    3997        1217 :     psATInfo->sTI.pfnTransform = GDALApproxTransform;
    3998        1217 :     psATInfo->sTI.pfnCleanup = GDALDestroyApproxTransformer;
    3999        1217 :     psATInfo->sTI.pfnSerialize = GDALSerializeApproxTransformer;
    4000        1217 :     psATInfo->sTI.pfnCreateSimilar = GDALCreateSimilarApproxTransformer;
    4001             : 
    4002        1217 :     return psATInfo;
    4003             : }
    4004             : 
    4005             : /************************************************************************/
    4006             : /*              GDALApproxTransformerOwnsSubtransformer()               */
    4007             : /************************************************************************/
    4008             : 
    4009             : /** Set bOwnSubtransformer flag */
    4010        1189 : void GDALApproxTransformerOwnsSubtransformer(void *pCBData, int bOwnFlag)
    4011             : 
    4012             : {
    4013        1189 :     GDALApproxTransformInfo *psATInfo =
    4014             :         static_cast<GDALApproxTransformInfo *>(pCBData);
    4015             : 
    4016        1189 :     psATInfo->bOwnSubtransformer = bOwnFlag;
    4017        1189 : }
    4018             : 
    4019             : /************************************************************************/
    4020             : /*                    GDALDestroyApproxTransformer()                    */
    4021             : /************************************************************************/
    4022             : 
    4023             : /**
    4024             :  * Cleanup approximate transformer.
    4025             :  *
    4026             :  * Deallocates the resources allocated by GDALCreateApproxTransformer().
    4027             :  *
    4028             :  * @param pCBData callback data originally returned by
    4029             :  * GDALCreateApproxTransformer().
    4030             :  */
    4031             : 
    4032        1217 : void GDALDestroyApproxTransformer(void *pCBData)
    4033             : 
    4034             : {
    4035        1217 :     if (pCBData == nullptr)
    4036           0 :         return;
    4037             : 
    4038        1217 :     GDALApproxTransformInfo *psATInfo =
    4039             :         static_cast<GDALApproxTransformInfo *>(pCBData);
    4040             : 
    4041        1217 :     if (psATInfo->bOwnSubtransformer)
    4042        1214 :         GDALDestroyTransformer(psATInfo->pBaseCBData);
    4043             : 
    4044        1217 :     delete psATInfo;
    4045             : }
    4046             : 
    4047             : /************************************************************************/
    4048             : /*                  GDALRefreshApproxTransformer()                      */
    4049             : /************************************************************************/
    4050             : 
    4051          44 : void GDALRefreshApproxTransformer(void *hTransformArg)
    4052             : {
    4053          44 :     GDALApproxTransformInfo *psInfo =
    4054             :         static_cast<GDALApproxTransformInfo *>(hTransformArg);
    4055             : 
    4056          44 :     if (GDALIsTransformer(psInfo->pBaseCBData,
    4057             :                           GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME))
    4058             :     {
    4059          44 :         GDALRefreshGenImgProjTransformer(psInfo->pBaseCBData);
    4060             :     }
    4061          44 : }
    4062             : 
    4063             : /************************************************************************/
    4064             : /*                      GDALApproxTransformInternal()                   */
    4065             : /************************************************************************/
    4066             : 
    4067     1224280 : static int GDALApproxTransformInternal(void *pCBData, int bDstToSrc,
    4068             :                                        int nPoints, double *x, double *y,
    4069             :                                        double *z, int *panSuccess,
    4070             :                                        // SME = Start, Middle, End.
    4071             :                                        const double xSMETransformed[3],
    4072             :                                        const double ySMETransformed[3],
    4073             :                                        const double zSMETransformed[3])
    4074             : {
    4075     1224280 :     GDALApproxTransformInfo *psATInfo =
    4076             :         static_cast<GDALApproxTransformInfo *>(pCBData);
    4077     1224280 :     const int nMiddle = (nPoints - 1) / 2;
    4078             : 
    4079             : #ifdef notdef_sanify_check
    4080             :     {
    4081             :         double x2[3] = {x[0], x[nMiddle], x[nPoints - 1]};
    4082             :         double y2[3] = {y[0], y[nMiddle], y[nPoints - 1]};
    4083             :         double z2[3] = {z[0], z[nMiddle], z[nPoints - 1]};
    4084             :         int anSuccess2[3] = {};
    4085             : 
    4086             :         const int bSuccess = psATInfo->pfnBaseTransformer(
    4087             :             psATInfo->pBaseCBData, bDstToSrc, 3, x2, y2, z2, anSuccess2);
    4088             :         CPLAssert(bSuccess);
    4089             :         CPLAssert(anSuccess2[0]);
    4090             :         CPLAssert(anSuccess2[1]);
    4091             :         CPLAssert(anSuccess2[2]);
    4092             :         CPLAssert(x2[0] == xSMETransformed[0]);
    4093             :         CPLAssert(y2[0] == ySMETransformed[0]);
    4094             :         CPLAssert(z2[0] == zSMETransformed[0]);
    4095             :         CPLAssert(x2[1] == xSMETransformed[1]);
    4096             :         CPLAssert(y2[1] == ySMETransformed[1]);
    4097             :         CPLAssert(z2[1] == zSMETransformed[1]);
    4098             :         CPLAssert(x2[2] == xSMETransformed[2]);
    4099             :         CPLAssert(y2[2] == ySMETransformed[2]);
    4100             :         CPLAssert(z2[2] == zSMETransformed[2]);
    4101             :     }
    4102             : #endif
    4103             : 
    4104             : #ifdef DEBUG_APPROX_TRANSFORMER
    4105             :     fprintf(stderr, "start (%.3f,%.3f) -> (%.3f,%.3f)\n", /*ok*/
    4106             :             x[0], y[0], xSMETransformed[0], ySMETransformed[0]);
    4107             :     fprintf(stderr, "middle (%.3f,%.3f) -> (%.3f,%.3f)\n", /*ok*/
    4108             :             x[nMiddle], y[nMiddle], xSMETransformed[1], ySMETransformed[1]);
    4109             :     fprintf(stderr, "end (%.3f,%.3f) -> (%.3f,%.3f)\n", /*ok*/
    4110             :             x[nPoints - 1], y[nPoints - 1], xSMETransformed[2],
    4111             :             ySMETransformed[2]);
    4112             : #endif
    4113             : 
    4114             :     /* -------------------------------------------------------------------- */
    4115             :     /*      Is the error at the middle acceptable relative to an            */
    4116             :     /*      interpolation of the middle position?                           */
    4117             :     /* -------------------------------------------------------------------- */
    4118     1224280 :     const double dfDeltaX =
    4119     1224280 :         (xSMETransformed[2] - xSMETransformed[0]) / (x[nPoints - 1] - x[0]);
    4120     1224280 :     const double dfDeltaY =
    4121     1224280 :         (ySMETransformed[2] - ySMETransformed[0]) / (x[nPoints - 1] - x[0]);
    4122     1224280 :     const double dfDeltaZ =
    4123     1224280 :         (zSMETransformed[2] - zSMETransformed[0]) / (x[nPoints - 1] - x[0]);
    4124             : 
    4125     1224280 :     const double dfError =
    4126     1224280 :         fabs((xSMETransformed[0] + dfDeltaX * (x[nMiddle] - x[0])) -
    4127     1224280 :              xSMETransformed[1]) +
    4128     1224280 :         fabs((ySMETransformed[0] + dfDeltaY * (x[nMiddle] - x[0])) -
    4129     1224280 :              ySMETransformed[1]);
    4130             : 
    4131     1224280 :     const double dfMaxError =
    4132     1224280 :         (bDstToSrc) ? psATInfo->dfMaxErrorReverse : psATInfo->dfMaxErrorForward;
    4133     1224280 :     if (dfError > dfMaxError)
    4134             :     {
    4135             : #if DEBUG_VERBOSE
    4136             :         CPLDebug("GDAL",
    4137             :                  "ApproxTransformer - "
    4138             :                  "error %g over threshold %g, subdivide %d points.",
    4139             :                  dfError, dfMaxError, nPoints);
    4140             : #endif
    4141             : 
    4142      593363 :         double xMiddle[3] = {x[(nMiddle - 1) / 2], x[nMiddle - 1],
    4143      593363 :                              x[nMiddle + (nPoints - nMiddle - 1) / 2]};
    4144      593363 :         double yMiddle[3] = {y[(nMiddle - 1) / 2], y[nMiddle - 1],
    4145      593363 :                              y[nMiddle + (nPoints - nMiddle - 1) / 2]};
    4146      593363 :         double zMiddle[3] = {z[(nMiddle - 1) / 2], z[nMiddle - 1],
    4147      593363 :                              z[nMiddle + (nPoints - nMiddle - 1) / 2]};
    4148             : 
    4149      593363 :         const bool bUseBaseTransformForHalf1 =
    4150      464934 :             nMiddle <= 5 || y[0] != y[nMiddle - 1] ||
    4151     1523230 :             y[0] != y[(nMiddle - 1) / 2] || x[0] == x[nMiddle - 1] ||
    4152      464934 :             x[0] == x[(nMiddle - 1) / 2];
    4153      593363 :         const bool bUseBaseTransformForHalf2 =
    4154      479069 :             nPoints - nMiddle <= 5 || y[nMiddle] != y[nPoints - 1] ||
    4155      479069 :             y[nMiddle] != y[nMiddle + (nPoints - nMiddle - 1) / 2] ||
    4156     1551500 :             x[nMiddle] == x[nPoints - 1] ||
    4157      479069 :             x[nMiddle] == x[nMiddle + (nPoints - nMiddle - 1) / 2];
    4158             : 
    4159      593363 :         int anSuccess2[3] = {};
    4160      593363 :         int bSuccess = FALSE;
    4161      593363 :         if (!bUseBaseTransformForHalf1 && !bUseBaseTransformForHalf2)
    4162      464934 :             bSuccess = psATInfo->pfnBaseTransformer(
    4163             :                 psATInfo->pBaseCBData, bDstToSrc, 3, xMiddle, yMiddle, zMiddle,
    4164             :                 anSuccess2);
    4165      128429 :         else if (!bUseBaseTransformForHalf1)
    4166             :         {
    4167           0 :             bSuccess = psATInfo->pfnBaseTransformer(
    4168             :                 psATInfo->pBaseCBData, bDstToSrc, 2, xMiddle, yMiddle, zMiddle,
    4169             :                 anSuccess2);
    4170           0 :             anSuccess2[2] = TRUE;
    4171             :         }
    4172      128429 :         else if (!bUseBaseTransformForHalf2)
    4173             :         {
    4174       14135 :             bSuccess = psATInfo->pfnBaseTransformer(
    4175             :                 psATInfo->pBaseCBData, bDstToSrc, 1, xMiddle + 2, yMiddle + 2,
    4176             :                 zMiddle + 2, anSuccess2 + 2);
    4177       14135 :             anSuccess2[0] = TRUE;
    4178       14135 :             anSuccess2[1] = TRUE;
    4179             :         }
    4180             : 
    4181      593363 :         if (!bSuccess || !anSuccess2[0] || !anSuccess2[1] || !anSuccess2[2])
    4182             :         {
    4183      114310 :             bSuccess = psATInfo->pfnBaseTransformer(
    4184             :                 psATInfo->pBaseCBData, bDstToSrc, nMiddle - 1, x + 1, y + 1,
    4185             :                 z + 1, panSuccess + 1);
    4186      228620 :             bSuccess &= psATInfo->pfnBaseTransformer(
    4187      114310 :                 psATInfo->pBaseCBData, bDstToSrc, nPoints - nMiddle - 2,
    4188      114310 :                 x + nMiddle + 1, y + nMiddle + 1, z + nMiddle + 1,
    4189      114310 :                 panSuccess + nMiddle + 1);
    4190             : 
    4191      114310 :             x[0] = xSMETransformed[0];
    4192      114310 :             y[0] = ySMETransformed[0];
    4193      114310 :             z[0] = zSMETransformed[0];
    4194      114310 :             panSuccess[0] = TRUE;
    4195      114310 :             x[nMiddle] = xSMETransformed[1];
    4196      114310 :             y[nMiddle] = ySMETransformed[1];
    4197      114310 :             z[nMiddle] = zSMETransformed[1];
    4198      114310 :             panSuccess[nMiddle] = TRUE;
    4199      114310 :             x[nPoints - 1] = xSMETransformed[2];
    4200      114310 :             y[nPoints - 1] = ySMETransformed[2];
    4201      114310 :             z[nPoints - 1] = zSMETransformed[2];
    4202      114310 :             panSuccess[nPoints - 1] = TRUE;
    4203      114310 :             return bSuccess;
    4204             :         }
    4205             : 
    4206      479053 :         double x2[3] = {};
    4207      479053 :         double y2[3] = {};
    4208      479053 :         double z2[3] = {};
    4209      479053 :         if (!bUseBaseTransformForHalf1)
    4210             :         {
    4211      464918 :             x2[0] = xSMETransformed[0];
    4212      464918 :             y2[0] = ySMETransformed[0];
    4213      464918 :             z2[0] = zSMETransformed[0];
    4214      464918 :             x2[1] = xMiddle[0];
    4215      464918 :             y2[1] = yMiddle[0];
    4216      464918 :             z2[1] = zMiddle[0];
    4217      464918 :             x2[2] = xMiddle[1];
    4218      464918 :             y2[2] = yMiddle[1];
    4219      464918 :             z2[2] = zMiddle[1];
    4220             : 
    4221      464918 :             bSuccess = GDALApproxTransformInternal(
    4222             :                 psATInfo, bDstToSrc, nMiddle, x, y, z, panSuccess, x2, y2, z2);
    4223             :         }
    4224             :         else
    4225             :         {
    4226       14135 :             bSuccess = psATInfo->pfnBaseTransformer(
    4227             :                 psATInfo->pBaseCBData, bDstToSrc, nMiddle - 1, x + 1, y + 1,
    4228             :                 z + 1, panSuccess + 1);
    4229       14135 :             x[0] = xSMETransformed[0];
    4230       14135 :             y[0] = ySMETransformed[0];
    4231       14135 :             z[0] = zSMETransformed[0];
    4232       14135 :             panSuccess[0] = TRUE;
    4233             :         }
    4234             : 
    4235      479053 :         if (!bSuccess)
    4236          24 :             return FALSE;
    4237             : 
    4238      479029 :         if (!bUseBaseTransformForHalf2)
    4239             :         {
    4240      479029 :             x2[0] = xSMETransformed[1];
    4241      479029 :             y2[0] = ySMETransformed[1];
    4242      479029 :             z2[0] = zSMETransformed[1];
    4243      479029 :             x2[1] = xMiddle[2];
    4244      479029 :             y2[1] = yMiddle[2];
    4245      479029 :             z2[1] = zMiddle[2];
    4246      479029 :             x2[2] = xSMETransformed[2];
    4247      479029 :             y2[2] = ySMETransformed[2];
    4248      479029 :             z2[2] = zSMETransformed[2];
    4249             : 
    4250      479029 :             bSuccess = GDALApproxTransformInternal(
    4251      479029 :                 psATInfo, bDstToSrc, nPoints - nMiddle, x + nMiddle,
    4252      479029 :                 y + nMiddle, z + nMiddle, panSuccess + nMiddle, x2, y2, z2);
    4253             :         }
    4254             :         else
    4255             :         {
    4256           0 :             bSuccess = psATInfo->pfnBaseTransformer(
    4257           0 :                 psATInfo->pBaseCBData, bDstToSrc, nPoints - nMiddle - 2,
    4258           0 :                 x + nMiddle + 1, y + nMiddle + 1, z + nMiddle + 1,
    4259           0 :                 panSuccess + nMiddle + 1);
    4260             : 
    4261           0 :             x[nMiddle] = xSMETransformed[1];
    4262           0 :             y[nMiddle] = ySMETransformed[1];
    4263           0 :             z[nMiddle] = zSMETransformed[1];
    4264           0 :             panSuccess[nMiddle] = TRUE;
    4265           0 :             x[nPoints - 1] = xSMETransformed[2];
    4266           0 :             y[nPoints - 1] = ySMETransformed[2];
    4267           0 :             z[nPoints - 1] = zSMETransformed[2];
    4268           0 :             panSuccess[nPoints - 1] = TRUE;
    4269             :         }
    4270             : 
    4271      479029 :         if (!bSuccess)
    4272           2 :             return FALSE;
    4273             : 
    4274      479027 :         return TRUE;
    4275             :     }
    4276             : 
    4277             :     /* -------------------------------------------------------------------- */
    4278             :     /*      Error is OK since this is just used to compute output bounds    */
    4279             :     /*      of newly created file for gdalwarper.  So just use affine       */
    4280             :     /*      approximation of the reverse transform.  Eventually we          */
    4281             :     /*      should implement iterative searching to find a result within    */
    4282             :     /*      our error threshold.                                            */
    4283             :     /*      NOTE: the above comment is not true: gdalwarp uses approximator */
    4284             :     /*      also to compute the source pixel of each target pixel.          */
    4285             :     /* -------------------------------------------------------------------- */
    4286   124694000 :     for (int i = nPoints - 1; i >= 0; i--)
    4287             :     {
    4288             : #ifdef check_error
    4289             :         double xtemp = x[i];
    4290             :         double ytemp = y[i];
    4291             :         double ztemp = z[i];
    4292             :         double x_ori = xtemp;
    4293             :         double y_ori = ytemp;
    4294             :         int btemp = FALSE;
    4295             :         psATInfo->pfnBaseTransformer(psATInfo->pBaseCBData, bDstToSrc, 1,
    4296             :                                      &xtemp, &ytemp, &ztemp, &btemp);
    4297             : #endif
    4298   124063000 :         const double dfDist = (x[i] - x[0]);
    4299   124063000 :         x[i] = xSMETransformed[0] + dfDeltaX * dfDist;
    4300   124063000 :         y[i] = ySMETransformed[0] + dfDeltaY * dfDist;
    4301   124063000 :         z[i] = zSMETransformed[0] + dfDeltaZ * dfDist;
    4302             : #ifdef check_error
    4303             :         const double dfError2 = fabs(x[i] - xtemp) + fabs(y[i] - ytemp);
    4304             :         if (dfError2 > 4 /*10 * dfMaxError*/)
    4305             :         {
    4306             :             /*ok*/ printf("Error = %f on (%f, %f)\n", dfError2, x_ori, y_ori);
    4307             :         }
    4308             : #endif
    4309   124063000 :         panSuccess[i] = TRUE;
    4310             :     }
    4311             : 
    4312      630920 :     return TRUE;
    4313             : }
    4314             : 
    4315             : /************************************************************************/
    4316             : /*                        GDALApproxTransform()                         */
    4317             : /************************************************************************/
    4318             : 
    4319             : /**
    4320             :  * Perform approximate transformation.
    4321             :  *
    4322             :  * Actually performs the approximate transformation described in
    4323             :  * GDALCreateApproxTransformer().  This function matches the
    4324             :  * GDALTransformerFunc() signature.  Details of the arguments are described
    4325             :  * there.
    4326             :  */
    4327             : 
    4328      702683 : int GDALApproxTransform(void *pCBData, int bDstToSrc, int nPoints, double *x,
    4329             :                         double *y, double *z, int *panSuccess)
    4330             : 
    4331             : {
    4332      702683 :     GDALApproxTransformInfo *psATInfo =
    4333             :         static_cast<GDALApproxTransformInfo *>(pCBData);
    4334      702683 :     double x2[3] = {};
    4335      702683 :     double y2[3] = {};
    4336      702683 :     double z2[3] = {};
    4337      702683 :     int anSuccess2[3] = {};
    4338             :     int bSuccess;
    4339             : 
    4340      702683 :     const int nMiddle = (nPoints - 1) / 2;
    4341             : 
    4342             :     /* -------------------------------------------------------------------- */
    4343             :     /*      Bail if our preconditions are not met, or if error is not       */
    4344             :     /*      acceptable.                                                     */
    4345             :     /* -------------------------------------------------------------------- */
    4346      702683 :     int bRet = FALSE;
    4347      702683 :     if (y[0] != y[nPoints - 1] || y[0] != y[nMiddle] ||
    4348      693390 :         x[0] == x[nPoints - 1] || x[0] == x[nMiddle] ||
    4349      286175 :         (psATInfo->dfMaxErrorForward == 0.0 &&
    4350      286175 :          psATInfo->dfMaxErrorReverse == 0.0) ||
    4351             :         nPoints <= 5)
    4352             :     {
    4353      416930 :         bRet = psATInfo->pfnBaseTransformer(psATInfo->pBaseCBData, bDstToSrc,
    4354             :                                             nPoints, x, y, z, panSuccess);
    4355      416928 :         goto end;
    4356             :     }
    4357             : 
    4358             :     /* -------------------------------------------------------------------- */
    4359             :     /*      Transform first, last and middle point.                         */
    4360             :     /* -------------------------------------------------------------------- */
    4361      285753 :     x2[0] = x[0];
    4362      285753 :     y2[0] = y[0];
    4363      285753 :     z2[0] = z[0];
    4364      285753 :     x2[1] = x[nMiddle];
    4365      285753 :     y2[1] = y[nMiddle];
    4366      285753 :     z2[1] = z[nMiddle];
    4367      285753 :     x2[2] = x[nPoints - 1];
    4368      285753 :     y2[2] = y[nPoints - 1];
    4369      285753 :     z2[2] = z[nPoints - 1];
    4370             : 
    4371      285753 :     bSuccess = psATInfo->pfnBaseTransformer(psATInfo->pBaseCBData, bDstToSrc, 3,
    4372             :                                             x2, y2, z2, anSuccess2);
    4373      285754 :     if (!bSuccess || !anSuccess2[0] || !anSuccess2[1] || !anSuccess2[2])
    4374             :     {
    4375        5419 :         bRet = psATInfo->pfnBaseTransformer(psATInfo->pBaseCBData, bDstToSrc,
    4376             :                                             nPoints, x, y, z, panSuccess);
    4377        5416 :         goto end;
    4378             :     }
    4379             : 
    4380      280335 :     bRet = GDALApproxTransformInternal(pCBData, bDstToSrc, nPoints, x, y, z,
    4381             :                                        panSuccess, x2, y2, z2);
    4382             : 
    4383      702683 : end:
    4384             : #ifdef DEBUG_APPROX_TRANSFORMER
    4385             :     for (int i = 0; i < nPoints; i++)
    4386             :         fprintf(stderr, "[%d] (%.10f,%.10f) %d\n", /*ok*/
    4387             :                 i, x[i], y[i], panSuccess[i]);
    4388             : #endif
    4389             : 
    4390      702683 :     return bRet;
    4391             : }
    4392             : 
    4393             : /************************************************************************/
    4394             : /*                  GDALDeserializeApproxTransformer()                  */
    4395             : /************************************************************************/
    4396             : 
    4397         151 : static void *GDALDeserializeApproxTransformer(CPLXMLNode *psTree)
    4398             : 
    4399             : {
    4400         151 :     double dfMaxErrorForward = 0.25;
    4401         151 :     double dfMaxErrorReverse = 0.25;
    4402         151 :     const char *pszMaxError = CPLGetXMLValue(psTree, "MaxError", nullptr);
    4403         151 :     if (pszMaxError != nullptr)
    4404             :     {
    4405         149 :         dfMaxErrorForward = CPLAtof(pszMaxError);
    4406         149 :         dfMaxErrorReverse = dfMaxErrorForward;
    4407             :     }
    4408             :     const char *pszMaxErrorForward =
    4409         151 :         CPLGetXMLValue(psTree, "MaxErrorForward", nullptr);
    4410         151 :     if (pszMaxErrorForward != nullptr)
    4411             :     {
    4412           2 :         dfMaxErrorForward = CPLAtof(pszMaxErrorForward);
    4413             :     }
    4414             :     const char *pszMaxErrorReverse =
    4415         151 :         CPLGetXMLValue(psTree, "MaxErrorReverse", nullptr);
    4416         151 :     if (pszMaxErrorReverse != nullptr)
    4417             :     {
    4418           2 :         dfMaxErrorReverse = CPLAtof(pszMaxErrorReverse);
    4419             :     }
    4420             : 
    4421         151 :     GDALTransformerFunc pfnBaseTransform = nullptr;
    4422         151 :     void *pBaseCBData = nullptr;
    4423             : 
    4424         151 :     CPLXMLNode *psContainer = CPLGetXMLNode(psTree, "BaseTransformer");
    4425             : 
    4426         151 :     if (psContainer != nullptr && psContainer->psChild != nullptr)
    4427             :     {
    4428         151 :         GDALDeserializeTransformer(psContainer->psChild, &pfnBaseTransform,
    4429             :                                    &pBaseCBData);
    4430             :     }
    4431             : 
    4432         151 :     if (pfnBaseTransform == nullptr)
    4433             :     {
    4434           0 :         CPLError(CE_Failure, CPLE_AppDefined,
    4435             :                  "Cannot get base transform for approx transformer.");
    4436           0 :         return nullptr;
    4437             :     }
    4438             : 
    4439         151 :     void *pApproxCBData = GDALCreateApproxTransformer2(
    4440             :         pfnBaseTransform, pBaseCBData, dfMaxErrorForward, dfMaxErrorReverse);
    4441         151 :     GDALApproxTransformerOwnsSubtransformer(pApproxCBData, TRUE);
    4442             : 
    4443         151 :     return pApproxCBData;
    4444             : }
    4445             : 
    4446             : /************************************************************************/
    4447             : /*                 GDALTransformLonLatToDestApproxTransformer()         */
    4448             : /************************************************************************/
    4449             : 
    4450        2398 : int GDALTransformLonLatToDestApproxTransformer(void *hTransformArg,
    4451             :                                                double *pdfX, double *pdfY)
    4452             : {
    4453        2398 :     GDALApproxTransformInfo *psInfo =
    4454             :         static_cast<GDALApproxTransformInfo *>(hTransformArg);
    4455             : 
    4456        2398 :     if (GDALIsTransformer(psInfo->pBaseCBData,
    4457             :                           GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME))
    4458             :     {
    4459        2398 :         return GDALTransformLonLatToDestGenImgProjTransformer(
    4460        2398 :             psInfo->pBaseCBData, pdfX, pdfY);
    4461             :     }
    4462           0 :     return false;
    4463             : }
    4464             : 
    4465             : /************************************************************************/
    4466             : /*                       GDALApplyGeoTransform()                        */
    4467             : /************************************************************************/
    4468             : 
    4469             : /**
    4470             :  * Apply GeoTransform to x/y coordinate.
    4471             :  *
    4472             :  * Applies the following computation, converting a (pixel, line) coordinate
    4473             :  * into a georeferenced (geo_x, geo_y) location.
    4474             :  * \code{.c}
    4475             :  *  *pdfGeoX = padfGeoTransform[0] + dfPixel * padfGeoTransform[1]
    4476             :  *                                 + dfLine  * padfGeoTransform[2];
    4477             :  *  *pdfGeoY = padfGeoTransform[3] + dfPixel * padfGeoTransform[4]
    4478             :  *                                 + dfLine  * padfGeoTransform[5];
    4479             :  * \endcode
    4480             :  *
    4481             :  * @param padfGeoTransform Six coefficient GeoTransform to apply.
    4482             :  * @param dfPixel Input pixel position.
    4483             :  * @param dfLine Input line position.
    4484             :  * @param pdfGeoX output location where geo_x (easting/longitude)
    4485             :  * location is placed.
    4486             :  * @param pdfGeoY output location where geo_y (northing/latitude)
    4487             :  * location is placed.
    4488             :  */
    4489             : 
    4490      306698 : void CPL_STDCALL GDALApplyGeoTransform(const double *padfGeoTransform,
    4491             :                                        double dfPixel, double dfLine,
    4492             :                                        double *pdfGeoX, double *pdfGeoY)
    4493             : {
    4494      306698 :     *pdfGeoX = padfGeoTransform[0] + dfPixel * padfGeoTransform[1] +
    4495      306698 :                dfLine * padfGeoTransform[2];
    4496      306698 :     *pdfGeoY = padfGeoTransform[3] + dfPixel * padfGeoTransform[4] +
    4497      306698 :                dfLine * padfGeoTransform[5];
    4498      306698 : }
    4499             : 
    4500             : /************************************************************************/
    4501             : /*                        GDALInvGeoTransform()                         */
    4502             : /************************************************************************/
    4503             : 
    4504             : /**
    4505             :  * Invert Geotransform.
    4506             :  *
    4507             :  * This function will invert a standard 3x2 set of GeoTransform coefficients.
    4508             :  * This converts the equation from being pixel to geo to being geo to pixel.
    4509             :  *
    4510             :  * @param gt_in Input geotransform (six doubles - unaltered).
    4511             :  * @param gt_out Output geotransform (six doubles - updated).
    4512             :  *
    4513             :  * @return TRUE on success or FALSE if the equation is uninvertable.
    4514             :  */
    4515             : 
    4516        3341 : int CPL_STDCALL GDALInvGeoTransform(const double *gt_in, double *gt_out)
    4517             : 
    4518             : {
    4519             :     // Special case - no rotation - to avoid computing determinate
    4520             :     // and potential precision issues.
    4521        3341 :     if (gt_in[2] == 0.0 && gt_in[4] == 0.0 && gt_in[1] != 0.0 &&
    4522        3284 :         gt_in[5] != 0.0)
    4523             :     {
    4524             :         /*X = gt_in[0] + x * gt_in[1]
    4525             :           Y = gt_in[3] + y * gt_in[5]
    4526             :           -->
    4527             :           x = -gt_in[0] / gt_in[1] + (1 / gt_in[1]) * X
    4528             :           y = -gt_in[3] / gt_in[5] + (1 / gt_in[5]) * Y
    4529             :         */
    4530        3284 :         gt_out[0] = -gt_in[0] / gt_in[1];
    4531        3284 :         gt_out[1] = 1.0 / gt_in[1];
    4532        3284 :         gt_out[2] = 0.0;
    4533        3284 :         gt_out[3] = -gt_in[3] / gt_in[5];
    4534        3284 :         gt_out[4] = 0.0;
    4535        3284 :         gt_out[5] = 1.0 / gt_in[5];
    4536        3284 :         return 1;
    4537             :     }
    4538             : 
    4539             :     // Assume a 3rd row that is [1 0 0].
    4540             : 
    4541             :     // Compute determinate.
    4542             : 
    4543          57 :     const double det = gt_in[1] * gt_in[5] - gt_in[2] * gt_in[4];
    4544         114 :     const double magnitude = std::max(std::max(fabs(gt_in[1]), fabs(gt_in[2])),
    4545          57 :                                       std::max(fabs(gt_in[4]), fabs(gt_in[5])));
    4546             : 
    4547          57 :     if (fabs(det) <= 1e-10 * magnitude * magnitude)
    4548           8 :         return 0;
    4549             : 
    4550          49 :     const double inv_det = 1.0 / det;
    4551             : 
    4552             :     // Compute adjoint, and divide by determinate.
    4553             : 
    4554          49 :     gt_out[1] = gt_in[5] * inv_det;
    4555          49 :     gt_out[4] = -gt_in[4] * inv_det;
    4556             : 
    4557          49 :     gt_out[2] = -gt_in[2] * inv_det;
    4558          49 :     gt_out[5] = gt_in[1] * inv_det;
    4559             : 
    4560          49 :     gt_out[0] = (gt_in[2] * gt_in[3] - gt_in[0] * gt_in[5]) * inv_det;
    4561          49 :     gt_out[3] = (-gt_in[1] * gt_in[3] + gt_in[0] * gt_in[4]) * inv_det;
    4562             : 
    4563          49 :     return 1;
    4564             : }
    4565             : 
    4566             : /************************************************************************/
    4567             : /*                      GDALSerializeTransformer()                      */
    4568             : /************************************************************************/
    4569             : 
    4570         268 : CPLXMLNode *GDALSerializeTransformer(GDALTransformerFunc /* pfnFunc */,
    4571             :                                      void *pTransformArg)
    4572             : {
    4573         268 :     VALIDATE_POINTER1(pTransformArg, "GDALSerializeTransformer", nullptr);
    4574             : 
    4575         268 :     GDALTransformerInfo *psInfo =
    4576             :         static_cast<GDALTransformerInfo *>(pTransformArg);
    4577             : 
    4578         268 :     if (psInfo == nullptr || memcmp(psInfo->abySignature, GDAL_GTI2_SIGNATURE,
    4579             :                                     strlen(GDAL_GTI2_SIGNATURE)) != 0)
    4580             :     {
    4581           0 :         CPLError(CE_Failure, CPLE_AppDefined,
    4582             :                  "Attempt to serialize non-GTI2 transformer.");
    4583           0 :         return nullptr;
    4584             :     }
    4585         268 :     else if (psInfo->pfnSerialize == nullptr)
    4586             :     {
    4587           0 :         CPLError(CE_Failure, CPLE_AppDefined,
    4588             :                  "No serialization function available for this transformer.");
    4589           0 :         return nullptr;
    4590             :     }
    4591             : 
    4592         268 :     return psInfo->pfnSerialize(pTransformArg);
    4593             : }
    4594             : 
    4595             : /************************************************************************/
    4596             : /*                  GDALRegisterTransformDeserializer()                 */
    4597             : /************************************************************************/
    4598             : 
    4599             : static CPLList *psListDeserializer = nullptr;
    4600             : static CPLMutex *hDeserializerMutex = nullptr;
    4601             : 
    4602             : typedef struct
    4603             : {
    4604             :     char *pszTransformName;
    4605             :     GDALTransformerFunc pfnTransformerFunc;
    4606             :     GDALTransformDeserializeFunc pfnDeserializeFunc;
    4607             : } TransformDeserializerInfo;
    4608             : 
    4609           0 : void *GDALRegisterTransformDeserializer(
    4610             :     const char *pszTransformName, GDALTransformerFunc pfnTransformerFunc,
    4611             :     GDALTransformDeserializeFunc pfnDeserializeFunc)
    4612             : {
    4613             :     TransformDeserializerInfo *psInfo =
    4614             :         static_cast<TransformDeserializerInfo *>(
    4615           0 :             CPLMalloc(sizeof(TransformDeserializerInfo)));
    4616           0 :     psInfo->pszTransformName = CPLStrdup(pszTransformName);
    4617           0 :     psInfo->pfnTransformerFunc = pfnTransformerFunc;
    4618           0 :     psInfo->pfnDeserializeFunc = pfnDeserializeFunc;
    4619             : 
    4620           0 :     CPLMutexHolderD(&hDeserializerMutex);
    4621           0 :     psListDeserializer = CPLListInsert(psListDeserializer, psInfo, 0);
    4622             : 
    4623           0 :     return psInfo;
    4624             : }
    4625             : 
    4626             : /************************************************************************/
    4627             : /*                GDALUnregisterTransformDeserializer()                 */
    4628             : /************************************************************************/
    4629             : 
    4630           0 : void GDALUnregisterTransformDeserializer(void *pData)
    4631             : {
    4632           0 :     CPLMutexHolderD(&hDeserializerMutex);
    4633           0 :     CPLList *psList = psListDeserializer;
    4634           0 :     CPLList *psLast = nullptr;
    4635           0 :     while (psList)
    4636             :     {
    4637           0 :         if (psList->pData == pData)
    4638             :         {
    4639           0 :             TransformDeserializerInfo *psInfo =
    4640             :                 static_cast<TransformDeserializerInfo *>(pData);
    4641           0 :             CPLFree(psInfo->pszTransformName);
    4642           0 :             CPLFree(pData);
    4643           0 :             if (psLast)
    4644           0 :                 psLast->psNext = psList->psNext;
    4645             :             else
    4646           0 :                 psListDeserializer = nullptr;
    4647           0 :             CPLFree(psList);
    4648           0 :             break;
    4649             :         }
    4650           0 :         psLast = psList;
    4651           0 :         psList = psList->psNext;
    4652             :     }
    4653           0 : }
    4654             : 
    4655             : /************************************************************************/
    4656             : /*                GDALUnregisterTransformDeserializer()                 */
    4657             : /************************************************************************/
    4658             : 
    4659        1121 : void GDALCleanupTransformDeserializerMutex()
    4660             : {
    4661        1121 :     if (hDeserializerMutex != nullptr)
    4662             :     {
    4663           0 :         CPLDestroyMutex(hDeserializerMutex);
    4664           0 :         hDeserializerMutex = nullptr;
    4665             :     }
    4666        1121 : }
    4667             : 
    4668             : /************************************************************************/
    4669             : /*                     GDALDeserializeTransformer()                     */
    4670             : /************************************************************************/
    4671             : 
    4672         551 : CPLErr GDALDeserializeTransformer(CPLXMLNode *psTree,
    4673             :                                   GDALTransformerFunc *ppfnFunc,
    4674             :                                   void **ppTransformArg)
    4675             : 
    4676             : {
    4677         551 :     *ppfnFunc = nullptr;
    4678         551 :     *ppTransformArg = nullptr;
    4679             : 
    4680         551 :     CPLErrorReset();
    4681             : 
    4682         551 :     if (psTree == nullptr || psTree->eType != CXT_Element)
    4683           0 :         CPLError(CE_Failure, CPLE_AppDefined,
    4684             :                  "Malformed element in GDALDeserializeTransformer");
    4685         551 :     else if (EQUAL(psTree->pszValue, "GenImgProjTransformer"))
    4686             :     {
    4687         199 :         *ppfnFunc = GDALGenImgProjTransform;
    4688         199 :         *ppTransformArg = GDALDeserializeGenImgProjTransformer(psTree);
    4689             :     }
    4690         352 :     else if (EQUAL(psTree->pszValue, "ReprojectionTransformer"))
    4691             :     {
    4692         185 :         *ppfnFunc = GDALReprojectionTransform;
    4693         185 :         *ppTransformArg = GDALDeserializeReprojectionTransformer(psTree);
    4694             :     }
    4695         167 :     else if (EQUAL(psTree->pszValue, "GCPTransformer"))
    4696             :     {
    4697           5 :         *ppfnFunc = GDALGCPTransform;
    4698           5 :         *ppTransformArg = GDALDeserializeGCPTransformer(psTree);
    4699             :     }
    4700         162 :     else if (EQUAL(psTree->pszValue, "TPSTransformer"))
    4701             :     {
    4702           3 :         *ppfnFunc = GDALTPSTransform;
    4703           3 :         *ppTransformArg = GDALDeserializeTPSTransformer(psTree);
    4704             :     }
    4705         159 :     else if (EQUAL(psTree->pszValue, "GeoLocTransformer"))
    4706             :     {
    4707           1 :         *ppfnFunc = GDALGeoLocTransform;
    4708           1 :         *ppTransformArg = GDALDeserializeGeoLocTransformer(psTree);
    4709             :     }
    4710         158 :     else if (EQUAL(psTree->pszValue, "RPCTransformer"))
    4711             :     {
    4712           0 :         *ppfnFunc = GDALRPCTransform;
    4713           0 :         *ppTransformArg = GDALDeserializeRPCTransformer(psTree);
    4714             :     }
    4715         158 :     else if (EQUAL(psTree->pszValue, "ApproxTransformer"))
    4716             :     {
    4717         151 :         *ppfnFunc = GDALApproxTransform;
    4718         151 :         *ppTransformArg = GDALDeserializeApproxTransformer(psTree);
    4719             :     }
    4720           7 :     else if (EQUAL(psTree->pszValue, "HomographyTransformer"))
    4721             :     {
    4722           7 :         *ppfnFunc = GDALHomographyTransform;
    4723           7 :         *ppTransformArg = GDALDeserializeHomographyTransformer(psTree);
    4724             :     }
    4725             :     else
    4726             :     {
    4727           0 :         GDALTransformDeserializeFunc pfnDeserializeFunc = nullptr;
    4728             :         {
    4729           0 :             CPLMutexHolderD(&hDeserializerMutex);
    4730           0 :             CPLList *psList = psListDeserializer;
    4731           0 :             while (psList)
    4732             :             {
    4733           0 :                 TransformDeserializerInfo *psInfo =
    4734             :                     static_cast<TransformDeserializerInfo *>(psList->pData);
    4735           0 :                 if (strcmp(psInfo->pszTransformName, psTree->pszValue) == 0)
    4736             :                 {
    4737           0 :                     *ppfnFunc = psInfo->pfnTransformerFunc;
    4738           0 :                     pfnDeserializeFunc = psInfo->pfnDeserializeFunc;
    4739           0 :                     break;
    4740             :                 }
    4741           0 :                 psList = psList->psNext;
    4742             :             }
    4743             :         }
    4744             : 
    4745           0 :         if (pfnDeserializeFunc != nullptr)
    4746             :         {
    4747           0 :             *ppTransformArg = pfnDeserializeFunc(psTree);
    4748             :         }
    4749             :         else
    4750             :         {
    4751           0 :             CPLError(CE_Failure, CPLE_AppDefined,
    4752             :                      "Unrecognized element '%s' GDALDeserializeTransformer",
    4753             :                      psTree->pszValue);
    4754             :         }
    4755             :     }
    4756             : 
    4757         551 :     return CPLGetLastErrorType();
    4758             : }
    4759             : 
    4760             : /************************************************************************/
    4761             : /*                       GDALDestroyTransformer()                       */
    4762             : /************************************************************************/
    4763             : 
    4764        4151 : void GDALDestroyTransformer(void *pTransformArg)
    4765             : 
    4766             : {
    4767        4151 :     if (pTransformArg == nullptr)
    4768           0 :         return;
    4769             : 
    4770        4151 :     GDALTransformerInfo *psInfo =
    4771             :         static_cast<GDALTransformerInfo *>(pTransformArg);
    4772             : 
    4773        4151 :     if (memcmp(psInfo->abySignature, GDAL_GTI2_SIGNATURE,
    4774             :                strlen(GDAL_GTI2_SIGNATURE)) != 0)
    4775             :     {
    4776           0 :         CPLError(CE_Failure, CPLE_AppDefined,
    4777             :                  "Attempt to destroy non-GTI2 transformer.");
    4778           0 :         return;
    4779             :     }
    4780             : 
    4781        4151 :     psInfo->pfnCleanup(pTransformArg);
    4782             : }
    4783             : 
    4784             : /************************************************************************/
    4785             : /*                         GDALUseTransformer()                         */
    4786             : /************************************************************************/
    4787             : 
    4788        8683 : int GDALUseTransformer(void *pTransformArg, int bDstToSrc, int nPointCount,
    4789             :                        double *x, double *y, double *z, int *panSuccess)
    4790             : {
    4791        8683 :     GDALTransformerInfo *psInfo =
    4792             :         static_cast<GDALTransformerInfo *>(pTransformArg);
    4793             : 
    4794        8683 :     if (psInfo == nullptr || memcmp(psInfo->abySignature, GDAL_GTI2_SIGNATURE,
    4795             :                                     strlen(GDAL_GTI2_SIGNATURE)) != 0)
    4796             :     {
    4797           0 :         CPLError(CE_Failure, CPLE_AppDefined,
    4798             :                  "Attempt to use non-GTI2 transformer.");
    4799           0 :         return FALSE;
    4800             :     }
    4801             : 
    4802        8683 :     return psInfo->pfnTransform(pTransformArg, bDstToSrc, nPointCount, x, y, z,
    4803        8683 :                                 panSuccess);
    4804             : }
    4805             : 
    4806             : /************************************************************************/
    4807             : /*                        GDALCloneTransformer()                        */
    4808             : /************************************************************************/
    4809             : 
    4810          39 : void *GDALCloneTransformer(void *pTransformArg)
    4811             : {
    4812          39 :     GDALTransformerInfo *psInfo =
    4813             :         static_cast<GDALTransformerInfo *>(pTransformArg);
    4814             : 
    4815          39 :     if (psInfo == nullptr || memcmp(psInfo->abySignature, GDAL_GTI2_SIGNATURE,
    4816             :                                     strlen(GDAL_GTI2_SIGNATURE)) != 0)
    4817             :     {
    4818           0 :         CPLError(CE_Failure, CPLE_AppDefined,
    4819             :                  "Attempt to clone non-GTI2 transformer.");
    4820           0 :         return nullptr;
    4821             :     }
    4822             : 
    4823          39 :     if (psInfo->pfnCreateSimilar != nullptr)
    4824             :     {
    4825          20 :         return psInfo->pfnCreateSimilar(psInfo, 1.0, 1.0);
    4826             :     }
    4827             : 
    4828          19 :     if (psInfo->pfnSerialize == nullptr)
    4829             :     {
    4830           0 :         CPLError(CE_Failure, CPLE_AppDefined,
    4831             :                  "No serialization function available for this transformer.");
    4832           0 :         return nullptr;
    4833             :     }
    4834             : 
    4835          19 :     CPLXMLNode *pSerialized = psInfo->pfnSerialize(pTransformArg);
    4836          19 :     if (pSerialized == nullptr)
    4837           0 :         return nullptr;
    4838          19 :     GDALTransformerFunc pfnTransformer = nullptr;
    4839          19 :     void *pClonedTransformArg = nullptr;
    4840          19 :     if (GDALDeserializeTransformer(pSerialized, &pfnTransformer,
    4841          19 :                                    &pClonedTransformArg) != CE_None)
    4842             :     {
    4843           0 :         CPLDestroyXMLNode(pSerialized);
    4844           0 :         CPLFree(pClonedTransformArg);
    4845           0 :         return nullptr;
    4846             :     }
    4847             : 
    4848          19 :     CPLDestroyXMLNode(pSerialized);
    4849          19 :     return pClonedTransformArg;
    4850             : }
    4851             : 
    4852             : /************************************************************************/
    4853             : /*                   GDALCreateSimilarTransformer()                     */
    4854             : /************************************************************************/
    4855             : 
    4856          52 : void *GDALCreateSimilarTransformer(void *pTransformArg, double dfRatioX,
    4857             :                                    double dfRatioY)
    4858             : {
    4859          52 :     GDALTransformerInfo *psInfo =
    4860             :         static_cast<GDALTransformerInfo *>(pTransformArg);
    4861             : 
    4862          52 :     if (psInfo == nullptr || memcmp(psInfo->abySignature, GDAL_GTI2_SIGNATURE,
    4863             :                                     strlen(GDAL_GTI2_SIGNATURE)) != 0)
    4864             :     {
    4865           0 :         CPLError(CE_Failure, CPLE_AppDefined,
    4866             :                  "Attempt to call CreateSimilar on a non-GTI2 transformer.");
    4867           0 :         return nullptr;
    4868             :     }
    4869             : 
    4870          52 :     if (psInfo->pfnCreateSimilar == nullptr)
    4871             :     {
    4872           0 :         CPLError(CE_Failure, CPLE_AppDefined,
    4873             :                  "No CreateSimilar function available for this transformer.");
    4874           0 :         return nullptr;
    4875             :     }
    4876             : 
    4877          52 :     return psInfo->pfnCreateSimilar(psInfo, dfRatioX, dfRatioY);
    4878             : }
    4879             : 
    4880             : /************************************************************************/
    4881             : /*                      GetGenImgProjTransformInfo()                    */
    4882             : /************************************************************************/
    4883             : 
    4884          44 : static GDALTransformerInfo *GetGenImgProjTransformInfo(const char *pszFunc,
    4885             :                                                        void *pTransformArg)
    4886             : {
    4887          44 :     GDALTransformerInfo *psInfo =
    4888             :         static_cast<GDALTransformerInfo *>(pTransformArg);
    4889             : 
    4890          44 :     if (psInfo == nullptr || memcmp(psInfo->abySignature, GDAL_GTI2_SIGNATURE,
    4891             :                                     strlen(GDAL_GTI2_SIGNATURE)) != 0)
    4892             :     {
    4893           0 :         CPLError(CE_Failure, CPLE_AppDefined,
    4894             :                  "Attempt to call %s on "
    4895             :                  "a non-GTI2 transformer.",
    4896             :                  pszFunc);
    4897           0 :         return nullptr;
    4898             :     }
    4899             : 
    4900          44 :     if (EQUAL(psInfo->pszClassName, GDAL_APPROX_TRANSFORMER_CLASS_NAME))
    4901             :     {
    4902          14 :         GDALApproxTransformInfo *psATInfo =
    4903             :             static_cast<GDALApproxTransformInfo *>(pTransformArg);
    4904          14 :         psInfo = static_cast<GDALTransformerInfo *>(psATInfo->pBaseCBData);
    4905             : 
    4906          14 :         if (psInfo == nullptr ||
    4907          14 :             memcmp(psInfo->abySignature, GDAL_GTI2_SIGNATURE,
    4908             :                    strlen(GDAL_GTI2_SIGNATURE)) != 0)
    4909             :         {
    4910           0 :             CPLError(CE_Failure, CPLE_AppDefined,
    4911             :                      "Attempt to call %s on "
    4912             :                      "a non-GTI2 transformer.",
    4913             :                      pszFunc);
    4914           0 :             return nullptr;
    4915             :         }
    4916             :     }
    4917             : 
    4918          44 :     if (EQUAL(psInfo->pszClassName, GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME))
    4919             :     {
    4920          44 :         return psInfo;
    4921             :     }
    4922             : 
    4923           0 :     return nullptr;
    4924             : }
    4925             : 
    4926             : /************************************************************************/
    4927             : /*                 GDALSetTransformerDstGeoTransform()                  */
    4928             : /************************************************************************/
    4929             : 
    4930             : /**
    4931             :  * Set ApproxTransformer or GenImgProj output geotransform.
    4932             :  *
    4933             :  * This is a layer above GDALSetGenImgProjTransformerDstGeoTransform() that
    4934             :  * checks that the passed hTransformArg is compatible.
    4935             :  *
    4936             :  * Normally the "destination geotransform", or transformation between
    4937             :  * georeferenced output coordinates and pixel/line coordinates on the
    4938             :  * destination file is extracted from the destination file by
    4939             :  * GDALCreateGenImgProjTransformer() and stored in the GenImgProj private
    4940             :  * info.  However, sometimes it is inconvenient to have an output file
    4941             :  * handle with appropriate geotransform information when creating the
    4942             :  * transformation.  For these cases, this function can be used to apply
    4943             :  * the destination geotransform.
    4944             :  *
    4945             :  * @param pTransformArg the handle to update.
    4946             :  * @param padfGeoTransform the destination geotransform to apply (six doubles).
    4947             :  */
    4948             : 
    4949          22 : void GDALSetTransformerDstGeoTransform(void *pTransformArg,
    4950             :                                        const double *padfGeoTransform)
    4951             : {
    4952          22 :     VALIDATE_POINTER0(pTransformArg, "GDALSetTransformerDstGeoTransform");
    4953             : 
    4954          22 :     GDALTransformerInfo *psInfo = GetGenImgProjTransformInfo(
    4955             :         "GDALSetTransformerDstGeoTransform", pTransformArg);
    4956          22 :     if (psInfo)
    4957             :     {
    4958          22 :         GDALSetGenImgProjTransformerDstGeoTransform(psInfo, padfGeoTransform);
    4959             :     }
    4960             : }
    4961             : 
    4962             : /************************************************************************/
    4963             : /*                 GDALGetTransformerDstGeoTransform()                  */
    4964             : /************************************************************************/
    4965             : 
    4966             : /**
    4967             :  * Get ApproxTransformer or GenImgProj output geotransform.
    4968             :  *
    4969             :  * @param pTransformArg transformer handle.
    4970             :  * @param padfGeoTransform (output) the destination geotransform to return (six
    4971             :  * doubles).
    4972             :  */
    4973             : 
    4974          22 : void GDALGetTransformerDstGeoTransform(void *pTransformArg,
    4975             :                                        double *padfGeoTransform)
    4976             : {
    4977          22 :     VALIDATE_POINTER0(pTransformArg, "GDALGetTransformerDstGeoTransform");
    4978             : 
    4979          22 :     GDALTransformerInfo *psInfo = GetGenImgProjTransformInfo(
    4980             :         "GDALGetTransformerDstGeoTransform", pTransformArg);
    4981          22 :     if (psInfo)
    4982             :     {
    4983          22 :         GDALGenImgProjTransformInfo *psGenImgProjInfo =
    4984             :             reinterpret_cast<GDALGenImgProjTransformInfo *>(psInfo);
    4985             : 
    4986          22 :         memcpy(padfGeoTransform, psGenImgProjInfo->sDstParams.adfGeoTransform,
    4987             :                sizeof(double) * 6);
    4988             :     }
    4989             : }
    4990             : 
    4991             : /************************************************************************/
    4992             : /*            GDALTransformIsTranslationOnPixelBoundaries()             */
    4993             : /************************************************************************/
    4994             : 
    4995        1529 : bool GDALTransformIsTranslationOnPixelBoundaries(GDALTransformerFunc,
    4996             :                                                  void *pTransformerArg)
    4997             : {
    4998        1529 :     if (GDALIsTransformer(pTransformerArg, GDAL_APPROX_TRANSFORMER_CLASS_NAME))
    4999             :     {
    5000        1199 :         const auto *pApproxInfo =
    5001             :             static_cast<const GDALApproxTransformInfo *>(pTransformerArg);
    5002        1199 :         pTransformerArg = pApproxInfo->pBaseCBData;
    5003             :     }
    5004        1529 :     if (GDALIsTransformer(pTransformerArg, GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME))
    5005             :     {
    5006        1376 :         const auto *pGenImgpProjInfo =
    5007             :             static_cast<GDALGenImgProjTransformInfo *>(pTransformerArg);
    5008         446 :         const auto IsCloseToInteger = [](double dfVal)
    5009         446 :         { return std::fabs(dfVal - std::round(dfVal)) <= 1e-6; };
    5010        2672 :         return pGenImgpProjInfo->sSrcParams.pTransformArg == nullptr &&
    5011        1296 :                pGenImgpProjInfo->sDstParams.pTransformArg == nullptr &&
    5012        1294 :                pGenImgpProjInfo->pReproject == nullptr &&
    5013         526 :                pGenImgpProjInfo->sSrcParams.adfGeoTransform[1] ==
    5014         526 :                    pGenImgpProjInfo->sDstParams.adfGeoTransform[1] &&
    5015         275 :                pGenImgpProjInfo->sSrcParams.adfGeoTransform[5] ==
    5016         275 :                    pGenImgpProjInfo->sDstParams.adfGeoTransform[5] &&
    5017         241 :                pGenImgpProjInfo->sSrcParams.adfGeoTransform[2] ==
    5018         241 :                    pGenImgpProjInfo->sDstParams.adfGeoTransform[2] &&
    5019         241 :                pGenImgpProjInfo->sSrcParams.adfGeoTransform[4] ==
    5020         482 :                    pGenImgpProjInfo->sDstParams.adfGeoTransform[4] &&
    5021             :                // Check that the georeferenced origin of the destination
    5022             :                // geotransform is close to be an integer value when transformed
    5023             :                // to source image coordinates
    5024         241 :                IsCloseToInteger(
    5025         241 :                    pGenImgpProjInfo->sSrcParams.adfInvGeoTransform[0] +
    5026         241 :                    pGenImgpProjInfo->sDstParams.adfGeoTransform[0] *
    5027         241 :                        pGenImgpProjInfo->sSrcParams.adfInvGeoTransform[1] +
    5028         241 :                    pGenImgpProjInfo->sDstParams.adfGeoTransform[3] *
    5029        2913 :                        pGenImgpProjInfo->sSrcParams.adfInvGeoTransform[2]) &&
    5030         205 :                IsCloseToInteger(
    5031         205 :                    pGenImgpProjInfo->sSrcParams.adfInvGeoTransform[3] +
    5032         205 :                    pGenImgpProjInfo->sDstParams.adfGeoTransform[0] *
    5033         205 :                        pGenImgpProjInfo->sSrcParams.adfInvGeoTransform[4] +
    5034         205 :                    pGenImgpProjInfo->sDstParams.adfGeoTransform[3] *
    5035        1581 :                        pGenImgpProjInfo->sSrcParams.adfInvGeoTransform[5]);
    5036             :     }
    5037         153 :     return false;
    5038             : }
    5039             : 
    5040             : /************************************************************************/
    5041             : /*                   GDALTransformIsAffineNoRotation()                  */
    5042             : /************************************************************************/
    5043             : 
    5044          18 : bool GDALTransformIsAffineNoRotation(GDALTransformerFunc, void *pTransformerArg)
    5045             : {
    5046          18 :     if (GDALIsTransformer(pTransformerArg, GDAL_APPROX_TRANSFORMER_CLASS_NAME))
    5047             :     {
    5048          18 :         const auto *pApproxInfo =
    5049             :             static_cast<const GDALApproxTransformInfo *>(pTransformerArg);
    5050          18 :         pTransformerArg = pApproxInfo->pBaseCBData;
    5051             :     }
    5052          18 :     if (GDALIsTransformer(pTransformerArg, GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME))
    5053             :     {
    5054          18 :         const auto *pGenImgpProjInfo =
    5055             :             static_cast<GDALGenImgProjTransformInfo *>(pTransformerArg);
    5056          36 :         return pGenImgpProjInfo->sSrcParams.pTransformArg == nullptr &&
    5057          18 :                pGenImgpProjInfo->sDstParams.pTransformArg == nullptr &&
    5058          18 :                pGenImgpProjInfo->pReproject == nullptr &&
    5059           8 :                pGenImgpProjInfo->sSrcParams.adfGeoTransform[2] == 0 &&
    5060           8 :                pGenImgpProjInfo->sSrcParams.adfGeoTransform[4] == 0 &&
    5061          44 :                pGenImgpProjInfo->sDstParams.adfGeoTransform[2] == 0 &&
    5062          26 :                pGenImgpProjInfo->sDstParams.adfGeoTransform[4] == 0;
    5063             :     }
    5064           0 :     return false;
    5065             : }
    5066             : 
    5067             : /************************************************************************/
    5068             : /*                        GDALTransformHasFastClone()                   */
    5069             : /************************************************************************/
    5070             : 
    5071             : /** Returns whether GDALCloneTransformer() on this transformer is
    5072             :  * "fast"
    5073             :  * Counter-examples are GCPs or TPSs transformers.
    5074             :  */
    5075           2 : bool GDALTransformHasFastClone(void *pTransformerArg)
    5076             : {
    5077           2 :     if (GDALIsTransformer(pTransformerArg, GDAL_APPROX_TRANSFORMER_CLASS_NAME))
    5078             :     {
    5079           1 :         const auto *pApproxInfo =
    5080             :             static_cast<const GDALApproxTransformInfo *>(pTransformerArg);
    5081           1 :         pTransformerArg = pApproxInfo->pBaseCBData;
    5082             :         // Fallback to next lines
    5083             :     }
    5084             : 
    5085           2 :     if (GDALIsTransformer(pTransformerArg, GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME))
    5086             :     {
    5087           2 :         const auto *pGenImgpProjInfo =
    5088             :             static_cast<GDALGenImgProjTransformInfo *>(pTransformerArg);
    5089           2 :         return (pGenImgpProjInfo->sSrcParams.pTransformArg == nullptr ||
    5090           0 :                 GDALTransformHasFastClone(
    5091           4 :                     pGenImgpProjInfo->sSrcParams.pTransformArg)) &&
    5092           2 :                (pGenImgpProjInfo->sDstParams.pTransformArg == nullptr ||
    5093           0 :                 GDALTransformHasFastClone(
    5094           2 :                     pGenImgpProjInfo->sDstParams.pTransformArg));
    5095             :     }
    5096           0 :     else if (GDALIsTransformer(pTransformerArg,
    5097             :                                GDAL_RPC_TRANSFORMER_CLASS_NAME))
    5098             :     {
    5099           0 :         return true;
    5100             :     }
    5101             :     else
    5102             :     {
    5103           0 :         return false;
    5104             :     }
    5105             : }

Generated by: LCOV version 1.14