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