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

Generated by: LCOV version 1.14