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 10418 : bool GDALIsTransformer(void *hTransformerArg, const char *pszClassName)
67 : {
68 10418 : if (!hTransformerArg)
69 0 : return false;
70 : // All transformers should have a GDALTransformerInfo member as their first members
71 10418 : GDALTransformerInfo *psInfo =
72 : static_cast<GDALTransformerInfo *>(hTransformerArg);
73 10418 : return memcmp(psInfo->abySignature, GDAL_GTI2_SIGNATURE,
74 19942 : strlen(GDAL_GTI2_SIGNATURE)) == 0 &&
75 19942 : 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 40 : CPLErr CPL_STDCALL GDALSuggestedWarpOutput(GDALDatasetH hSrcDS,
173 : GDALTransformerFunc pfnTransformer,
174 : void *pTransformArg,
175 : double *padfGeoTransformOut,
176 : int *pnPixels, int *pnLines)
177 :
178 : {
179 40 : VALIDATE_POINTER1(hSrcDS, "GDALSuggestedWarpOutput", CE_Failure);
180 :
181 40 : double adfExtent[4] = {};
182 :
183 40 : return GDALSuggestedWarpOutput2(hSrcDS, pfnTransformer, pTransformArg,
184 : padfGeoTransformOut, pnPixels, pnLines,
185 40 : adfExtent, 0);
186 : }
187 :
188 488 : 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 488 : double adfX[21] = {};
194 488 : double adfY[21] = {};
195 :
196 488 : const double dfMaxXOut = padfExtent[2];
197 488 : const double dfMaxYOut = padfExtent[3];
198 :
199 : // Take 20 steps.
200 488 : int nSamplePoints = 0;
201 10736 : for (double dfRatio = 0.0; dfRatio <= 1.01; dfRatio += 0.05)
202 : {
203 : // Ensure we end exactly at the end.
204 10248 : if (dfRatio > 0.99)
205 488 : dfRatio = 1.0;
206 :
207 : // Along right.
208 10248 : adfX[nSamplePoints] = dfMaxXOut;
209 10248 : adfY[nSamplePoints] = dfMaxYOut - dfPixelSizeY * dfRatio * nLines;
210 10248 : nSamplePoints++;
211 : }
212 488 : double adfZ[21] = {};
213 :
214 488 : int abSuccess[21] = {};
215 :
216 488 : bool bErr = false;
217 488 : if (!pfnTransformer(pTransformArg, TRUE, nSamplePoints, adfX, adfY, adfZ,
218 : abSuccess))
219 : {
220 0 : bErr = true;
221 : }
222 :
223 488 : if (!bErr && !pfnTransformer(pTransformArg, FALSE, nSamplePoints, adfX,
224 : adfY, adfZ, abSuccess))
225 : {
226 0 : bErr = true;
227 : }
228 :
229 488 : nSamplePoints = 0;
230 488 : int nBadCount = 0;
231 10736 : for (double dfRatio = 0.0; !bErr && dfRatio <= 1.01; dfRatio += 0.05)
232 : {
233 10248 : const double expected_x = dfMaxXOut;
234 10248 : const double expected_y = dfMaxYOut - dfPixelSizeY * dfRatio * nLines;
235 10248 : if (fabs(adfX[nSamplePoints] - expected_x) > dfPixelSizeX ||
236 6018 : fabs(adfY[nSamplePoints] - expected_y) > dfPixelSizeY)
237 4230 : nBadCount++;
238 10248 : nSamplePoints++;
239 : }
240 :
241 488 : return nBadCount == nSamplePoints;
242 : }
243 :
244 387 : 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 387 : double adfX[21] = {};
250 387 : double adfY[21] = {};
251 :
252 387 : const double dfMinXOut = padfExtent[0];
253 387 : const double dfMinYOut = padfExtent[1];
254 :
255 : // Take 20 steps.
256 387 : int nSamplePoints = 0;
257 8514 : for (double dfRatio = 0.0; dfRatio <= 1.01; dfRatio += 0.05)
258 : {
259 : // Ensure we end exactly at the end.
260 8127 : if (dfRatio > 0.99)
261 387 : dfRatio = 1.0;
262 :
263 : // Along right.
264 8127 : adfX[nSamplePoints] = dfMinXOut + dfPixelSizeX * dfRatio * nPixels;
265 8127 : adfY[nSamplePoints] = dfMinYOut;
266 8127 : nSamplePoints++;
267 : }
268 387 : double adfZ[21] = {};
269 :
270 387 : int abSuccess[21] = {};
271 :
272 387 : bool bErr = false;
273 387 : if (!pfnTransformer(pTransformArg, TRUE, nSamplePoints, adfX, adfY, adfZ,
274 : abSuccess))
275 : {
276 0 : bErr = true;
277 : }
278 :
279 387 : if (!bErr && !pfnTransformer(pTransformArg, FALSE, nSamplePoints, adfX,
280 : adfY, adfZ, abSuccess))
281 : {
282 0 : bErr = true;
283 : }
284 :
285 387 : nSamplePoints = 0;
286 387 : int nBadCount = 0;
287 8514 : for (double dfRatio = 0.0; !bErr && dfRatio <= 1.01; dfRatio += 0.05)
288 : {
289 8127 : const double expected_x = dfMinXOut + dfPixelSizeX * dfRatio * nPixels;
290 8127 : const double expected_y = dfMinYOut;
291 8127 : if (fabs(adfX[nSamplePoints] - expected_x) > dfPixelSizeX ||
292 5952 : fabs(adfY[nSamplePoints] - expected_y) > dfPixelSizeY)
293 2175 : nBadCount++;
294 8127 : nSamplePoints++;
295 : }
296 :
297 387 : 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 896 : 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 896 : VALIDATE_POINTER1(hSrcDS, "GDALSuggestedWarpOutput2", CE_Failure);
355 :
356 : const bool bIsGDALGenImgProjTransform{
357 1792 : pTransformArg &&
358 896 : GDALIsTransformer(pTransformArg, GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME)};
359 :
360 : /* -------------------------------------------------------------------- */
361 : /* Setup sample points all around the edge of the input raster. */
362 : /* -------------------------------------------------------------------- */
363 896 : if (bIsGDALGenImgProjTransform)
364 : {
365 : // In case CHECK_WITH_INVERT_PROJ has been modified.
366 896 : 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 896 : const int nInXSize = GDALGetRasterXSize(hSrcDS);
376 896 : const int nInYSize = GDALGetRasterYSize(hSrcDS);
377 :
378 : /* ------------------------------------------------------------- */
379 : /* Special case for warping on the same (or null) CRS. */
380 : /* ------------------------------------------------------------- */
381 896 : if ((!nOptions || (nOptions & GDAL_SWO_FORCE_SQUARE_PIXEL) == 0) &&
382 895 : pTransformArg && bIsGDALGenImgProjTransform)
383 : {
384 895 : const GDALGenImgProjTransformInfo *psInfo =
385 : static_cast<const GDALGenImgProjTransformInfo *>(pTransformArg);
386 :
387 895 : if (!psInfo->pSrcTransformer &&
388 822 : !psInfo->bHasCustomTransformationPipeline &&
389 818 : !psInfo->pDstTransformer && psInfo->adfSrcGeoTransform[2] == 0 &&
390 818 : psInfo->adfSrcGeoTransform[4] == 0 &&
391 818 : psInfo->adfDstGeoTransform[0] == 0 &&
392 806 : psInfo->adfDstGeoTransform[1] == 1 &&
393 806 : psInfo->adfDstGeoTransform[2] == 0 &&
394 806 : psInfo->adfDstGeoTransform[3] == 0 &&
395 806 : psInfo->adfDstGeoTransform[4] == 0 &&
396 806 : psInfo->adfDstGeoTransform[5] == 1)
397 : {
398 806 : const OGRSpatialReference *poSourceCRS = nullptr;
399 806 : const OGRSpatialReference *poTargetCRS = nullptr;
400 :
401 806 : if (psInfo->pReprojectArg)
402 : {
403 576 : const GDALReprojectionTransformInfo *psRTI =
404 : static_cast<const GDALReprojectionTransformInfo *>(
405 : psInfo->pReprojectArg);
406 576 : poSourceCRS = psRTI->poForwardTransform->GetSourceCS();
407 576 : poTargetCRS = psRTI->poForwardTransform->GetTargetCS();
408 : }
409 :
410 1382 : if ((!poSourceCRS && !poTargetCRS) ||
411 576 : (poSourceCRS && poTargetCRS &&
412 576 : poSourceCRS->IsSame(poTargetCRS)))
413 : {
414 :
415 575 : const bool bNorthUp{psInfo->adfSrcGeoTransform[5] < 0.0};
416 :
417 575 : memcpy(padfGeoTransformOut, psInfo->adfSrcGeoTransform,
418 : sizeof(double) * 6);
419 :
420 575 : if (!bNorthUp)
421 : {
422 38 : padfGeoTransformOut[3] = padfGeoTransformOut[3] +
423 38 : nInYSize * padfGeoTransformOut[5];
424 38 : padfGeoTransformOut[5] = -padfGeoTransformOut[5];
425 : }
426 :
427 575 : *pnPixels = nInXSize;
428 575 : *pnLines = nInYSize;
429 :
430 : // Calculate extent from hSrcDS
431 575 : if (padfExtent)
432 : {
433 575 : padfExtent[0] = psInfo->adfSrcGeoTransform[0];
434 575 : padfExtent[1] = psInfo->adfSrcGeoTransform[3] +
435 575 : nInYSize * psInfo->adfSrcGeoTransform[5];
436 575 : padfExtent[2] = psInfo->adfSrcGeoTransform[0] +
437 575 : nInXSize * psInfo->adfSrcGeoTransform[1];
438 575 : padfExtent[3] = psInfo->adfSrcGeoTransform[3];
439 575 : if (!bNorthUp)
440 : {
441 38 : std::swap(padfExtent[1], padfExtent[3]);
442 : }
443 : }
444 575 : return CE_None;
445 : }
446 : }
447 : }
448 :
449 321 : const int N_PIXELSTEP = 50;
450 : int nSteps = static_cast<int>(
451 321 : static_cast<double>(std::min(nInYSize, nInXSize)) / N_PIXELSTEP + 0.5);
452 321 : if (nSteps < 20)
453 296 : nSteps = 20;
454 25 : else if (nSteps > 100)
455 10 : 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 321 : retry:
463 321 : int nStepsPlusOne = nSteps + 1;
464 321 : int nSampleMax = nStepsPlusOne * nStepsPlusOne;
465 :
466 321 : double dfStep = 1.0 / nSteps;
467 321 : double *padfY = nullptr;
468 321 : double *padfZ = nullptr;
469 321 : double *padfYRevert = nullptr;
470 321 : double *padfZRevert = nullptr;
471 :
472 : int *pabSuccess = static_cast<int *>(
473 321 : VSI_MALLOC3_VERBOSE(sizeof(int), nStepsPlusOne, nStepsPlusOne));
474 : double *padfX = static_cast<double *>(
475 321 : VSI_MALLOC3_VERBOSE(sizeof(double) * 3, nStepsPlusOne, nStepsPlusOne));
476 : double *padfXRevert = static_cast<double *>(
477 321 : VSI_MALLOC3_VERBOSE(sizeof(double) * 3, nStepsPlusOne, nStepsPlusOne));
478 321 : 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 321 : padfY = padfX + nSampleMax;
492 321 : padfZ = padfX + nSampleMax * 2;
493 321 : padfYRevert = padfXRevert + nSampleMax;
494 321 : padfZRevert = padfXRevert + nSampleMax * 2;
495 :
496 : // Take N_STEPS steps.
497 7965 : for (int iStep = 0; iStep <= nSteps; iStep++)
498 : {
499 7644 : double dfRatio = (iStep == nSteps) ? 1.0 : iStep * dfStep;
500 7644 : int iStep2 = iStep;
501 :
502 : // Along top.
503 7644 : padfX[iStep2] = dfRatio * nInXSize;
504 7644 : padfY[iStep2] = 0.0;
505 7644 : padfZ[iStep2] = 0.0;
506 :
507 : // Along bottom.
508 7644 : iStep2 += nStepsPlusOne;
509 7644 : padfX[iStep2] = dfRatio * nInXSize;
510 7644 : padfY[iStep2] = nInYSize;
511 7644 : padfZ[iStep2] = 0.0;
512 :
513 : // Along left.
514 7644 : iStep2 += nStepsPlusOne;
515 7644 : padfX[iStep2] = 0.0;
516 7644 : padfY[iStep2] = dfRatio * nInYSize;
517 7644 : padfZ[iStep2] = 0.0;
518 :
519 : // Along right.
520 7644 : iStep2 += nStepsPlusOne;
521 7644 : padfX[iStep2] = nInXSize;
522 7644 : padfY[iStep2] = dfRatio * nInYSize;
523 7644 : padfZ[iStep2] = 0.0;
524 : }
525 :
526 321 : int nSamplePoints = 4 * nStepsPlusOne;
527 :
528 321 : memset(pabSuccess, 1, sizeof(int) * nSampleMax);
529 :
530 : /* -------------------------------------------------------------------- */
531 : /* Transform them to the output coordinate system. */
532 : /* -------------------------------------------------------------------- */
533 321 : 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 321 : constexpr int SIGN_FINAL_UNINIT = -2;
546 321 : constexpr int SIGN_FINAL_INVALID = 0;
547 321 : int iSignDiscontinuity = SIGN_FINAL_UNINIT;
548 321 : int nFailedCount = 0;
549 321 : const int iSignArray[2] = {-1, 1};
550 30897 : for (int i = 0; i < nSamplePoints; i++)
551 : {
552 30576 : 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 27612 : if (iSignDiscontinuity == 1 || iSignDiscontinuity == -1)
560 : {
561 9777 : if (!((iSignDiscontinuity * padfX[i] > 0 &&
562 9728 : 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 17835 : else if (iSignDiscontinuity == SIGN_FINAL_UNINIT)
569 : {
570 741 : for (const auto &iSign : iSignArray)
571 : {
572 554 : if ((iSign * padfX[i] > 0 && iSign * padfX[i] <= 180.0) ||
573 421 : (fabs(padfX[i] - iSign * -180.0) < 1e-8))
574 : {
575 133 : iSignDiscontinuity = iSign;
576 133 : break;
577 : }
578 : }
579 320 : if (iSignDiscontinuity == SIGN_FINAL_UNINIT)
580 : {
581 187 : iSignDiscontinuity = SIGN_FINAL_INVALID;
582 : }
583 : }
584 : }
585 : else
586 : {
587 2964 : nFailedCount++;
588 : }
589 : }
590 :
591 321 : if (iSignDiscontinuity == 1 || iSignDiscontinuity == -1)
592 : {
593 10338 : for (int i = 0; i < nSamplePoints; i++)
594 : {
595 10232 : if (pabSuccess[i])
596 : {
597 9615 : 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 321 : if (nFailedCount)
621 : {
622 49 : CPLDebug("WARP", "At least one point failed after direct transform");
623 : }
624 : else
625 : {
626 272 : memcpy(padfXRevert, padfX, nSamplePoints * sizeof(double));
627 272 : memcpy(padfYRevert, padfY, nSamplePoints * sizeof(double));
628 272 : memcpy(padfZRevert, padfZ, nSamplePoints * sizeof(double));
629 272 : if (pfnTransformer(pTransformArg, TRUE, nSamplePoints, padfXRevert,
630 272 : padfYRevert, padfZRevert, pabSuccess))
631 : {
632 24863 : for (int i = 0; nFailedCount == 0 && i < nSamplePoints; i++)
633 : {
634 24607 : if (!pabSuccess[i])
635 : {
636 16 : nFailedCount++;
637 16 : break;
638 : }
639 :
640 24591 : double dfRatio = (i % nStepsPlusOne) * dfStep;
641 24591 : if (dfRatio > 0.99)
642 1022 : dfRatio = 1.0;
643 :
644 24591 : double dfExpectedX = 0.0;
645 24591 : double dfExpectedY = 0.0;
646 24591 : if (i < nStepsPlusOne)
647 : {
648 6279 : dfExpectedX = dfRatio * nInXSize;
649 : }
650 18312 : else if (i < 2 * nStepsPlusOne)
651 : {
652 6232 : dfExpectedX = dfRatio * nInXSize;
653 6232 : dfExpectedY = nInYSize;
654 : }
655 12080 : else if (i < 3 * nStepsPlusOne)
656 : {
657 6102 : dfExpectedY = dfRatio * nInYSize;
658 : }
659 : else
660 : {
661 5978 : dfExpectedX = nInXSize;
662 5978 : dfExpectedY = dfRatio * nInYSize;
663 : }
664 :
665 24591 : if (fabs(padfXRevert[i] - dfExpectedX) >
666 24591 : nInXSize / static_cast<double>(nSteps) ||
667 24582 : fabs(padfYRevert[i] - dfExpectedY) >
668 24582 : nInYSize / static_cast<double>(nSteps))
669 9 : nFailedCount++;
670 : }
671 272 : 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 321 : 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 321 : double dfMinXOut = 0.0;
729 321 : double dfMinYOut = 0.0;
730 321 : double dfMaxXOut = 0.0;
731 321 : double dfMaxYOut = 0.0;
732 321 : bool bGotInitialPoint = false;
733 :
734 321 : nFailedCount = 0;
735 73459 : for (int i = 0; i < nSamplePoints; i++)
736 : {
737 73138 : int x_i = 0;
738 73138 : int y_i = 0;
739 :
740 73138 : if (nSamplePoints == nSampleMax)
741 : {
742 49394 : x_i = i % nStepsPlusOne;
743 49394 : y_i = i / nStepsPlusOne;
744 : }
745 : else
746 : {
747 23744 : if (i < 2 * nStepsPlusOne)
748 : {
749 11872 : x_i = i % nStepsPlusOne;
750 11872 : y_i = (i < nStepsPlusOne) ? 0 : nSteps;
751 : }
752 : }
753 :
754 73138 : if (x_i > 0 && (pabSuccess[i - 1] || pabSuccess[i]))
755 : {
756 47905 : double x_out_before = padfX[i - 1];
757 47905 : double x_out_after = padfX[i];
758 47905 : int nIter = 0;
759 47905 : double x_in_before =
760 47905 : static_cast<double>(x_i - 1) * nInXSize / nSteps;
761 47905 : double x_in_after = static_cast<double>(x_i) * nInXSize / nSteps;
762 47905 : int invalid_before = !(pabSuccess[i - 1]);
763 47905 : 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 63581 : while ((invalid_before || invalid_after ||
770 133758 : 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 73138 : if (!pabSuccess[i])
830 : {
831 12036 : nFailedCount++;
832 12036 : continue;
833 : }
834 :
835 61102 : if (bGotInitialPoint)
836 : {
837 60800 : dfMinXOut = std::min(dfMinXOut, padfX[i]);
838 60800 : dfMinYOut = std::min(dfMinYOut, padfY[i]);
839 60800 : dfMaxXOut = std::max(dfMaxXOut, padfX[i]);
840 60800 : dfMaxYOut = std::max(dfMaxYOut, padfY[i]);
841 : }
842 : else
843 : {
844 302 : bGotInitialPoint = true;
845 302 : dfMinXOut = padfX[i];
846 302 : dfMaxXOut = padfX[i];
847 302 : dfMinYOut = padfY[i];
848 302 : dfMaxYOut = padfY[i];
849 : }
850 : }
851 :
852 321 : 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 317 : if (nFailedCount)
867 45 : CPLDebug("GDAL",
868 : "GDALSuggestedWarpOutput(): %d out of %d points failed to "
869 : "transform.",
870 : nFailedCount, nSamplePoints);
871 :
872 317 : bool bIsGeographicCoordsDeg = false;
873 317 : if (bIsGDALGenImgProjTransform)
874 : {
875 317 : const GDALGenImgProjTransformInfo *pGIPTI =
876 : static_cast<const GDALGenImgProjTransformInfo *>(pTransformArg);
877 317 : 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 291 : else if (pGIPTI->pSrcTransformer == nullptr &&
945 248 : pGIPTI->pDstTransformer == nullptr &&
946 248 : pGIPTI->pReproject == GDALReprojectionTransform &&
947 237 : pGIPTI->adfDstGeoTransform[0] == 0 &&
948 235 : pGIPTI->adfDstGeoTransform[1] == 1 &&
949 235 : pGIPTI->adfDstGeoTransform[2] == 0 &&
950 235 : pGIPTI->adfDstGeoTransform[3] == 0 &&
951 235 : pGIPTI->adfDstGeoTransform[4] == 0 &&
952 235 : pGIPTI->adfDstGeoTransform[5] == 1)
953 : {
954 : /* ------------------------------------------------------------- */
955 : /* Special case for warping using source geotransform and */
956 : /* reprojection to deal with the poles. */
957 : /* ------------------------------------------------------------- */
958 235 : const GDALReprojectionTransformInfo *psRTI =
959 : static_cast<const GDALReprojectionTransformInfo *>(
960 : pGIPTI->pReprojectArg);
961 : const OGRSpatialReference *poSourceCRS =
962 235 : psRTI->poForwardTransform->GetSourceCS();
963 : const OGRSpatialReference *poTargetCRS =
964 235 : psRTI->poForwardTransform->GetTargetCS();
965 469 : if (poTargetCRS != nullptr &&
966 234 : psRTI->poReverseTransform != nullptr &&
967 234 : poTargetCRS->IsGeographic() &&
968 90 : fabs(poTargetCRS->GetAngularUnits() -
969 559 : CPLAtof(SRS_UA_DEGREE_CONV)) < 1e-9 &&
970 90 : (!poSourceCRS || !poSourceCRS->IsGeographic()))
971 : {
972 83 : bIsGeographicCoordsDeg = true;
973 :
974 83 : std::unique_ptr<CPLConfigOptionSetter> poSetter;
975 83 : 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 249 : for (const auto &sign : iSignArray)
988 : {
989 166 : double X = 0.0;
990 166 : const double Yinit = 90.0 * sign;
991 166 : double Y = Yinit;
992 166 : if (psRTI->poReverseTransform->Transform(1, &X, &Y))
993 : {
994 112 : const auto invGT = pGIPTI->adfSrcInvGeoTransform;
995 112 : const double x = invGT[0] + X * invGT[1] + Y * invGT[2];
996 112 : const double y = invGT[3] + X * invGT[4] + Y * invGT[5];
997 112 : constexpr double EPSILON = 1e-5;
998 112 : 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 83 : 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 235 : if (poSourceCRS != nullptr && poTargetCRS != nullptr &&
1051 234 : pGIPTI->adfSrcGeoTransform[1] != 0 &&
1052 234 : pGIPTI->adfSrcGeoTransform[2] == 0 &&
1053 234 : pGIPTI->adfSrcGeoTransform[4] == 0 &&
1054 234 : pGIPTI->adfSrcGeoTransform[5] != 0)
1055 : {
1056 234 : const double dfULX = pGIPTI->adfSrcGeoTransform[0];
1057 234 : const double dfULY = pGIPTI->adfSrcGeoTransform[3];
1058 234 : const double dfLRX =
1059 234 : dfULX + pGIPTI->adfSrcGeoTransform[1] * nInXSize;
1060 234 : const double dfLRY =
1061 234 : dfULY + pGIPTI->adfSrcGeoTransform[5] * nInYSize;
1062 234 : const double dfMinSrcX = std::min(dfULX, dfLRX);
1063 234 : const double dfMinSrcY = std::min(dfULY, dfLRY);
1064 234 : const double dfMaxSrcX = std::max(dfULX, dfLRX);
1065 234 : const double dfMaxSrcY = std::max(dfULY, dfLRY);
1066 234 : double dfTmpMinXOut = std::numeric_limits<double>::max();
1067 234 : double dfTmpMinYOut = std::numeric_limits<double>::max();
1068 234 : double dfTmpMaxXOut = std::numeric_limits<double>::min();
1069 234 : double dfTmpMaxYOut = std::numeric_limits<double>::min();
1070 468 : if (psRTI->poForwardTransform->TransformBounds(
1071 : dfMinSrcX, dfMinSrcY, dfMaxSrcX, dfMaxSrcY,
1072 : &dfTmpMinXOut, &dfTmpMinYOut, &dfTmpMaxXOut,
1073 : &dfTmpMaxYOut,
1074 234 : 2)) // minimum number of points as we already have a
1075 : // logic above to sample
1076 : {
1077 223 : dfMinXOut = std::min(dfMinXOut, dfTmpMinXOut);
1078 223 : dfMinYOut = std::min(dfMinYOut, dfTmpMinYOut);
1079 223 : dfMaxXOut = std::max(dfMaxXOut, dfTmpMaxXOut);
1080 223 : 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 317 : double dfDiagonalDist = 0.0;
1094 317 : double dfDeltaX = 0.0;
1095 317 : double dfDeltaY = 0.0;
1096 :
1097 317 : if (pabSuccess[0] && pabSuccess[nSamplePoints - 1])
1098 : {
1099 275 : dfDeltaX = padfX[nSamplePoints - 1] - padfX[0];
1100 275 : 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 317 : if (dfDeltaX == 0.0 || dfDeltaY == 0.0)
1105 : {
1106 47 : dfDeltaX = dfMaxXOut - dfMinXOut;
1107 47 : dfDeltaY = dfMaxYOut - dfMinYOut;
1108 : }
1109 :
1110 317 : dfDiagonalDist = sqrt(dfDeltaX * dfDeltaX + dfDeltaY * dfDeltaY);
1111 :
1112 : /* -------------------------------------------------------------------- */
1113 : /* Compute a pixel size from this. */
1114 : /* -------------------------------------------------------------------- */
1115 : const double dfPixelSize =
1116 317 : dfDiagonalDist / sqrt(static_cast<double>(nInXSize) * nInXSize +
1117 317 : static_cast<double>(nInYSize) * nInYSize);
1118 :
1119 317 : const double dfPixels = (dfMaxXOut - dfMinXOut) / dfPixelSize;
1120 317 : const double dfLines = (dfMaxYOut - dfMinYOut) / dfPixelSize;
1121 :
1122 317 : const int knIntMaxMinusOne = std::numeric_limits<int>::max() - 1;
1123 317 : 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 317 : 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 309 : *pnPixels = static_cast<int>(dfPixels + 0.5);
1145 309 : *pnLines = static_cast<int>(dfLines + 0.5);
1146 : }
1147 :
1148 317 : double dfPixelSizeX = dfPixelSize;
1149 317 : double dfPixelSizeY = dfPixelSize;
1150 :
1151 317 : 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 496 : for (const auto &dfRatio : adfRatioArray)
1159 : {
1160 488 : const double dfTryPixelSizeX =
1161 488 : dfPixelSizeX - dfPixelSizeX * dfRatio / *pnPixels;
1162 488 : double adfExtent[4] = {dfMinXOut, dfMaxYOut - (*pnLines) * dfPixelSizeY,
1163 488 : dfMinXOut + (*pnPixels) * dfTryPixelSizeX,
1164 488 : dfMaxYOut};
1165 488 : if (!GDALSuggestedWarpOutput2_MustAdjustForRightBorder(
1166 : pfnTransformer, pTransformArg, adfExtent, *pnPixels, *pnLines,
1167 : dfTryPixelSizeX, dfPixelSizeY))
1168 : {
1169 309 : dfPixelSizeX = dfTryPixelSizeX;
1170 309 : 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 397 : for (const auto &dfRatio : adfRatioArray)
1180 : {
1181 387 : const double dfTryPixelSizeY =
1182 387 : dfPixelSizeY - dfPixelSizeY * dfRatio / *pnLines;
1183 : double adfExtent[4] = {
1184 387 : dfMinXOut, dfMaxYOut - (*pnLines) * dfTryPixelSizeY,
1185 387 : dfMinXOut + (*pnPixels) * dfPixelSizeX, dfMaxYOut};
1186 387 : if (!GDALSuggestedWarpOutput2_MustAdjustForBottomBorder(
1187 : pfnTransformer, pTransformArg, adfExtent, *pnPixels, *pnLines,
1188 : dfPixelSizeX, dfTryPixelSizeY))
1189 : {
1190 307 : dfPixelSizeY = dfTryPixelSizeY;
1191 307 : break;
1192 : }
1193 : }
1194 :
1195 : /* -------------------------------------------------------------------- */
1196 : /* Recompute some bounds so that all return values are consistent */
1197 : /* -------------------------------------------------------------------- */
1198 317 : double dfMaxXOutNew = dfMinXOut + (*pnPixels) * dfPixelSizeX;
1199 317 : if (bIsGeographicCoordsDeg &&
1200 107 : ((dfMaxXOut <= 180 && dfMaxXOutNew > 180) || dfMaxXOut == 180))
1201 : {
1202 3 : dfMaxXOut = 180;
1203 3 : dfPixelSizeX = (dfMaxXOut - dfMinXOut) / *pnPixels;
1204 : }
1205 : else
1206 : {
1207 314 : dfMaxXOut = dfMaxXOutNew;
1208 : }
1209 :
1210 317 : double dfMinYOutNew = dfMaxYOut - (*pnLines) * dfPixelSizeY;
1211 317 : if (bIsGeographicCoordsDeg && dfMinYOut >= -90 && dfMinYOutNew < -90)
1212 : {
1213 0 : dfMinYOut = -90;
1214 0 : dfPixelSizeY = (dfMaxYOut - dfMinYOut) / *pnLines;
1215 : }
1216 : else
1217 : {
1218 317 : dfMinYOut = dfMinYOutNew;
1219 : }
1220 :
1221 : /* -------------------------------------------------------------------- */
1222 : /* Return raw extents. */
1223 : /* -------------------------------------------------------------------- */
1224 317 : padfExtent[0] = dfMinXOut;
1225 317 : padfExtent[1] = dfMinYOut;
1226 317 : padfExtent[2] = dfMaxXOut;
1227 317 : padfExtent[3] = dfMaxYOut;
1228 :
1229 : /* -------------------------------------------------------------------- */
1230 : /* Set the output geotransform. */
1231 : /* -------------------------------------------------------------------- */
1232 317 : padfGeoTransformOut[0] = dfMinXOut;
1233 317 : padfGeoTransformOut[1] = dfPixelSizeX;
1234 317 : padfGeoTransformOut[2] = 0.0;
1235 317 : padfGeoTransformOut[3] = dfMaxYOut;
1236 317 : padfGeoTransformOut[4] = 0.0;
1237 317 : padfGeoTransformOut[5] = -dfPixelSizeY;
1238 :
1239 317 : CPLFree(padfX);
1240 317 : CPLFree(padfXRevert);
1241 317 : CPLFree(pabSuccess);
1242 :
1243 317 : return CE_None;
1244 : }
1245 :
1246 : /************************************************************************/
1247 : /* GetCurrentCheckWithInvertPROJ() */
1248 : /************************************************************************/
1249 :
1250 2632 : static bool GetCurrentCheckWithInvertPROJ()
1251 : {
1252 2632 : 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 1637 : static GDALGenImgProjTransformInfo *GDALCreateGenImgProjTransformerInternal()
1264 : {
1265 : /* -------------------------------------------------------------------- */
1266 : /* Initialize the transform info. */
1267 : /* -------------------------------------------------------------------- */
1268 : GDALGenImgProjTransformInfo *psInfo =
1269 : static_cast<GDALGenImgProjTransformInfo *>(
1270 1637 : CPLCalloc(sizeof(GDALGenImgProjTransformInfo), 1));
1271 :
1272 1637 : memcpy(psInfo->sTI.abySignature, GDAL_GTI2_SIGNATURE,
1273 : strlen(GDAL_GTI2_SIGNATURE));
1274 1637 : psInfo->sTI.pszClassName = GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME;
1275 1637 : psInfo->sTI.pfnTransform = GDALGenImgProjTransform;
1276 1637 : psInfo->sTI.pfnCleanup = GDALDestroyGenImgProjTransformer;
1277 1637 : psInfo->sTI.pfnSerialize = GDALSerializeGenImgProjTransformer;
1278 1637 : psInfo->sTI.pfnCreateSimilar = GDALCreateSimilarGenImgProjTransformer;
1279 :
1280 1637 : psInfo->bCheckWithInvertPROJ = GetCurrentCheckWithInvertPROJ();
1281 1637 : psInfo->bHasCustomTransformationPipeline = false;
1282 :
1283 1637 : 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 715 : static void InsertCenterLong(GDALDatasetH hDS, OGRSpatialReference *poSRS,
1428 : CPLStringList &aosOptions)
1429 :
1430 : {
1431 1252 : if (!poSRS->IsGeographic() || std::fabs(poSRS->GetAngularUnits() -
1432 537 : CPLAtof(SRS_UA_DEGREE_CONV)) > 1e-9)
1433 : {
1434 179 : 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 1006 : 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 1006 : bool ret = false;
1501 :
1502 1006 : if (!poSRS)
1503 0 : return false;
1504 :
1505 1006 : OGRSpatialReference oSrcSRSHoriz(*poSRS);
1506 1006 : if (oSrcSRSHoriz.IsCompound())
1507 : {
1508 17 : oSrcSRSHoriz.StripVertical();
1509 : }
1510 :
1511 1006 : OGRSpatialReference *poGeog = oSrcSRSHoriz.CloneGeogCS();
1512 1006 : if (poGeog)
1513 : {
1514 1006 : poGeog->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1515 1006 : poGeog->SetAngularUnits(SRS_UA_DEGREE, CPLAtof(SRS_UA_DEGREE_CONV));
1516 :
1517 1006 : auto poCT = OGRCreateCoordinateTransformation(&oSrcSRSHoriz, poGeog);
1518 1006 : if (poCT)
1519 : {
1520 1006 : poCT->SetEmitErrors(false);
1521 :
1522 : double x[4], y[4];
1523 1006 : x[0] = adfGT[0];
1524 1006 : y[0] = adfGT[3];
1525 1006 : x[1] = adfGT[0] + nXSize * adfGT[1];
1526 1006 : y[1] = adfGT[3];
1527 1006 : x[2] = adfGT[0];
1528 1006 : y[2] = adfGT[3] + nYSize * adfGT[5];
1529 1006 : x[3] = x[1];
1530 1006 : y[3] = y[2];
1531 1006 : int validity[4] = {false, false, false, false};
1532 1006 : poCT->Transform(4, x, y, nullptr, validity);
1533 1006 : dfWestLongitudeDeg = std::numeric_limits<double>::max();
1534 1006 : dfSouthLatitudeDeg = std::numeric_limits<double>::max();
1535 1006 : dfEastLongitudeDeg = -std::numeric_limits<double>::max();
1536 1006 : dfNorthLatitudeDeg = -std::numeric_limits<double>::max();
1537 5030 : for (int i = 0; i < 4; i++)
1538 : {
1539 4024 : if (validity[i])
1540 : {
1541 4012 : ret = true;
1542 4012 : dfWestLongitudeDeg = std::min(dfWestLongitudeDeg, x[i]);
1543 4012 : dfSouthLatitudeDeg = std::min(dfSouthLatitudeDeg, y[i]);
1544 4012 : dfEastLongitudeDeg = std::max(dfEastLongitudeDeg, x[i]);
1545 4012 : dfNorthLatitudeDeg = std::max(dfNorthLatitudeDeg, y[i]);
1546 : }
1547 : }
1548 1006 : if (validity[0] && validity[1] && x[0] > x[1])
1549 : {
1550 3 : dfWestLongitudeDeg = x[0];
1551 3 : dfEastLongitudeDeg = x[1];
1552 : }
1553 1006 : if (ret && std::fabs(dfWestLongitudeDeg) <= 180 &&
1554 1002 : std::fabs(dfEastLongitudeDeg) <= 180 &&
1555 998 : std::fabs(dfSouthLatitudeDeg) <= 90 &&
1556 994 : std::fabs(dfNorthLatitudeDeg) <= 90)
1557 : {
1558 994 : 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 1006 : OGRCoordinateTransformation::DestroyCT(poCT);
1571 : }
1572 :
1573 1006 : delete poGeog;
1574 : }
1575 :
1576 1006 : 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 1391 : void *GDALCreateGenImgProjTransformer2(GDALDatasetH hSrcDS, GDALDatasetH hDstDS,
1908 : char **papszOptions)
1909 :
1910 : {
1911 1391 : CSLConstList papszMD = nullptr;
1912 : GDALRPCInfoV2 sRPCInfo;
1913 1391 : const char *pszMethod = CSLFetchNameValue(papszOptions, "SRC_METHOD");
1914 1391 : if (pszMethod == nullptr)
1915 1380 : pszMethod = CSLFetchNameValue(papszOptions, "METHOD");
1916 1391 : const char *pszSrcSRS = CSLFetchNameValue(papszOptions, "SRC_SRS");
1917 1391 : const char *pszDstSRS = CSLFetchNameValue(papszOptions, "DST_SRS");
1918 :
1919 1391 : const char *pszValue = CSLFetchNameValue(papszOptions, "MAX_GCP_ORDER");
1920 1391 : const int nOrder = pszValue ? atoi(pszValue) : 0;
1921 :
1922 1391 : pszValue = CSLFetchNameValue(papszOptions, "GCPS_OK");
1923 1391 : const bool bGCPUseOK = pszValue ? CPLTestBool(pszValue) : true;
1924 :
1925 1391 : pszValue = CSLFetchNameValue(papszOptions, "REFINE_MINIMUM_GCPS");
1926 1391 : const int nMinimumGcps = pszValue ? atoi(pszValue) : -1;
1927 :
1928 1391 : pszValue = CSLFetchNameValue(papszOptions, "REFINE_TOLERANCE");
1929 1391 : const bool bRefine = pszValue != nullptr;
1930 1391 : const double dfTolerance = pszValue ? CPLAtof(pszValue) : 0.0;
1931 :
1932 1391 : double dfWestLongitudeDeg = 0.0;
1933 1391 : double dfSouthLatitudeDeg = 0.0;
1934 1391 : double dfEastLongitudeDeg = 0.0;
1935 1391 : double dfNorthLatitudeDeg = 0.0;
1936 1391 : bool bHasAreaOfInterest = false;
1937 1391 : pszValue = CSLFetchNameValue(papszOptions, "AREA_OF_INTEREST");
1938 1391 : 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 1391 : const char *pszCO = CSLFetchNameValue(papszOptions, "COORDINATE_OPERATION");
1953 :
1954 : const auto SetAxisMapping =
1955 3160 : [papszOptions](OGRSpatialReference &oSRS, const char *pszPrefix)
1956 : {
1957 1054 : const char *pszMapping = CSLFetchNameValue(
1958 2108 : papszOptions, std::string(pszPrefix)
1959 1054 : .append("_DATA_AXIS_TO_SRS_AXIS_MAPPING")
1960 : .c_str());
1961 1054 : if (pszMapping)
1962 : {
1963 4 : CPLStringList aosTokens(CSLTokenizeString2(pszMapping, ",", 0));
1964 4 : std::vector<int> anMapping;
1965 6 : for (int i = 0; i < aosTokens.size(); ++i)
1966 4 : anMapping.push_back(atoi(aosTokens[i]));
1967 2 : oSRS.SetDataAxisToSRSAxisMapping(anMapping);
1968 : }
1969 : else
1970 : {
1971 1052 : const char *pszStrategy = CSLFetchNameValueDef(
1972 : papszOptions,
1973 2104 : std::string(pszPrefix).append("_AXIS_MAPPING_STRATEGY").c_str(),
1974 : "TRADITIONAL_GIS_ORDER");
1975 1052 : if (EQUAL(pszStrategy, "TRADITIONAL_GIS_ORDER"))
1976 1051 : oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1977 1 : else if (EQUAL(pszStrategy, "AUTHORITY_COMPLIANT"))
1978 1 : oSRS.SetAxisMappingStrategy(OAMS_AUTHORITY_COMPLIANT);
1979 : else
1980 : {
1981 0 : CPLError(CE_Warning, CPLE_AppDefined,
1982 : "Unrecognized value '%s' for %s", pszStrategy,
1983 0 : std::string(pszPrefix)
1984 0 : .append("_AXIS_MAPPING_STRATEGY")
1985 : .c_str());
1986 0 : return false;
1987 : }
1988 : }
1989 1054 : return true;
1990 1391 : };
1991 :
1992 2782 : OGRSpatialReference oSrcSRS;
1993 1391 : if (pszSrcSRS)
1994 : {
1995 234 : if (pszSrcSRS[0] != '\0' &&
1996 116 : oSrcSRS.SetFromUserInput(pszSrcSRS) != OGRERR_NONE)
1997 : {
1998 0 : CPLError(CE_Failure, CPLE_AppDefined,
1999 : "Failed to import coordinate system `%s'.", pszSrcSRS);
2000 0 : return nullptr;
2001 : }
2002 118 : if (!SetAxisMapping(oSrcSRS, "SRC_SRS"))
2003 0 : return nullptr;
2004 : }
2005 :
2006 2782 : OGRSpatialReference oDstSRS;
2007 1391 : if (pszDstSRS)
2008 : {
2009 1870 : if (pszDstSRS[0] != '\0' &&
2010 934 : oDstSRS.SetFromUserInput(pszDstSRS) != OGRERR_NONE)
2011 : {
2012 0 : CPLError(CE_Failure, CPLE_AppDefined,
2013 : "Failed to import coordinate system `%s'.", pszDstSRS);
2014 0 : return nullptr;
2015 : }
2016 936 : if (!SetAxisMapping(oDstSRS, "DST_SRS"))
2017 0 : return nullptr;
2018 : }
2019 :
2020 : const char *pszSrcGeolocArray =
2021 1391 : CSLFetchNameValueDef(papszOptions, "SRC_GEOLOC_ARRAY",
2022 : CSLFetchNameValue(papszOptions, "GEOLOC_ARRAY"));
2023 1391 : if (pszMethod == nullptr && pszSrcGeolocArray != nullptr)
2024 7 : pszMethod = "GEOLOC_ARRAY";
2025 :
2026 : /* -------------------------------------------------------------------- */
2027 : /* Initialize the transform info. */
2028 : /* -------------------------------------------------------------------- */
2029 : GDALGenImgProjTransformInfo *psInfo =
2030 1391 : GDALCreateGenImgProjTransformerInternal();
2031 :
2032 1391 : bool bCanUseSrcGeoTransform = false;
2033 :
2034 : /* -------------------------------------------------------------------- */
2035 : /* Get forward and inverse geotransform for the source image. */
2036 : /* -------------------------------------------------------------------- */
2037 1391 : if (hSrcDS == nullptr ||
2038 73 : (pszMethod != nullptr && EQUAL(pszMethod, "NO_GEOTRANSFORM")))
2039 : {
2040 87 : psInfo->adfSrcGeoTransform[0] = 0.0;
2041 87 : psInfo->adfSrcGeoTransform[1] = 1.0;
2042 87 : psInfo->adfSrcGeoTransform[2] = 0.0;
2043 87 : psInfo->adfSrcGeoTransform[3] = 0.0;
2044 87 : psInfo->adfSrcGeoTransform[4] = 0.0;
2045 87 : psInfo->adfSrcGeoTransform[5] = 1.0;
2046 87 : memcpy(psInfo->adfSrcInvGeoTransform, psInfo->adfSrcGeoTransform,
2047 : sizeof(double) * 6);
2048 : }
2049 2546 : else if ((pszMethod == nullptr || EQUAL(pszMethod, "GEOTRANSFORM")) &&
2050 1242 : GDALGetGeoTransform(hSrcDS, psInfo->adfSrcGeoTransform) == CE_None)
2051 : {
2052 1170 : if (!GDALInvGeoTransform(psInfo->adfSrcGeoTransform,
2053 1170 : psInfo->adfSrcInvGeoTransform))
2054 : {
2055 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot invert geotransform");
2056 0 : GDALDestroyGenImgProjTransformer(psInfo);
2057 0 : return nullptr;
2058 : }
2059 1170 : if (pszSrcSRS == nullptr)
2060 : {
2061 1097 : auto hSRS = GDALGetSpatialRef(hSrcDS);
2062 1097 : if (hSRS)
2063 922 : oSrcSRS = *(OGRSpatialReference::FromHandle(hSRS));
2064 : }
2065 1170 : if (!bHasAreaOfInterest && pszCO == nullptr && !oSrcSRS.IsEmpty())
2066 : {
2067 990 : GDALComputeAreaOfInterest(&oSrcSRS, psInfo->adfSrcGeoTransform,
2068 : GDALGetRasterXSize(hSrcDS),
2069 : GDALGetRasterYSize(hSrcDS),
2070 : dfWestLongitudeDeg, dfSouthLatitudeDeg,
2071 : dfEastLongitudeDeg, dfNorthLatitudeDeg);
2072 : }
2073 1170 : bCanUseSrcGeoTransform = true;
2074 : }
2075 134 : else if (bGCPUseOK &&
2076 62 : (pszMethod == nullptr || EQUAL(pszMethod, "GCP_POLYNOMIAL")) &&
2077 268 : GDALGetGCPCount(hSrcDS) > 0 && nOrder >= 0)
2078 : {
2079 29 : if (pszSrcSRS == nullptr)
2080 : {
2081 28 : auto hSRS = GDALGetGCPSpatialRef(hSrcDS);
2082 28 : if (hSRS)
2083 28 : oSrcSRS = *(OGRSpatialReference::FromHandle(hSRS));
2084 : }
2085 :
2086 29 : const auto nGCPCount = GDALGetGCPCount(hSrcDS);
2087 29 : auto pasGCPList = GDALDuplicateGCPs(nGCPCount, GDALGetGCPs(hSrcDS));
2088 29 : GDALGCPAntimeridianUnwrap(nGCPCount, pasGCPList, oSrcSRS, papszOptions);
2089 :
2090 29 : if (bRefine)
2091 : {
2092 0 : psInfo->pSrcTransformArg = GDALCreateGCPRefineTransformer(
2093 : nGCPCount, pasGCPList, nOrder, FALSE, dfTolerance,
2094 : nMinimumGcps);
2095 : }
2096 : else
2097 : {
2098 29 : psInfo->pSrcTransformArg =
2099 29 : GDALCreateGCPTransformer(nGCPCount, pasGCPList, nOrder, FALSE);
2100 : }
2101 :
2102 29 : GDALDeinitGCPs(nGCPCount, pasGCPList);
2103 29 : CPLFree(pasGCPList);
2104 :
2105 29 : if (psInfo->pSrcTransformArg == nullptr)
2106 : {
2107 0 : GDALDestroyGenImgProjTransformer(psInfo);
2108 0 : return nullptr;
2109 : }
2110 29 : psInfo->pSrcTransformer = GDALGCPTransform;
2111 : }
2112 :
2113 116 : else if (bGCPUseOK && GDALGetGCPCount(hSrcDS) > 0 && nOrder <= 0 &&
2114 11 : (pszMethod == nullptr || EQUAL(pszMethod, "GCP_TPS")))
2115 : {
2116 11 : if (pszSrcSRS == nullptr)
2117 : {
2118 11 : auto hSRS = GDALGetGCPSpatialRef(hSrcDS);
2119 11 : if (hSRS)
2120 10 : oSrcSRS = *(OGRSpatialReference::FromHandle(hSRS));
2121 : }
2122 :
2123 11 : const auto nGCPCount = GDALGetGCPCount(hSrcDS);
2124 11 : auto pasGCPList = GDALDuplicateGCPs(nGCPCount, GDALGetGCPs(hSrcDS));
2125 11 : GDALGCPAntimeridianUnwrap(nGCPCount, pasGCPList, oSrcSRS, papszOptions);
2126 :
2127 11 : psInfo->pSrcTransformArg = GDALCreateTPSTransformerInt(
2128 : nGCPCount, pasGCPList, FALSE, papszOptions);
2129 :
2130 11 : GDALDeinitGCPs(nGCPCount, pasGCPList);
2131 11 : CPLFree(pasGCPList);
2132 :
2133 11 : if (psInfo->pSrcTransformArg == nullptr)
2134 : {
2135 2 : GDALDestroyGenImgProjTransformer(psInfo);
2136 2 : return nullptr;
2137 : }
2138 9 : psInfo->pSrcTransformer = GDALTPSTransform;
2139 : }
2140 :
2141 48 : else if ((pszMethod == nullptr || EQUAL(pszMethod, "RPC")) &&
2142 185 : (papszMD = GDALGetMetadata(hSrcDS, "RPC")) != nullptr &&
2143 43 : GDALExtractRPCInfoV2(papszMD, &sRPCInfo))
2144 : {
2145 43 : psInfo->pSrcTransformArg =
2146 43 : GDALCreateRPCTransformerV2(&sRPCInfo, FALSE, 0, papszOptions);
2147 43 : if (psInfo->pSrcTransformArg == nullptr)
2148 : {
2149 1 : GDALDestroyGenImgProjTransformer(psInfo);
2150 1 : return nullptr;
2151 : }
2152 42 : psInfo->pSrcTransformer = GDALRPCTransform;
2153 42 : if (pszSrcSRS == nullptr)
2154 : {
2155 42 : oSrcSRS.SetFromUserInput(SRS_WKT_WGS84_LAT_LONG);
2156 42 : oSrcSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
2157 : }
2158 : }
2159 :
2160 102 : else if ((pszMethod == nullptr || EQUAL(pszMethod, "GEOLOC_ARRAY")) &&
2161 51 : ((papszMD = GDALGetMetadata(hSrcDS, "GEOLOCATION")) != nullptr ||
2162 : pszSrcGeolocArray != nullptr))
2163 : {
2164 44 : CPLStringList aosGeolocMD; // keep in this scope
2165 44 : if (pszSrcGeolocArray != nullptr)
2166 : {
2167 7 : if (papszMD != nullptr)
2168 : {
2169 0 : CPLError(
2170 : CE_Warning, CPLE_AppDefined,
2171 : "Both GEOLOCATION metadata domain on the source dataset "
2172 : "and [SRC_]GEOLOC_ARRAY transformer option are set. "
2173 : "Only using the later.");
2174 : }
2175 : aosGeolocMD =
2176 14 : GDALCreateGeolocationMetadata(hSrcDS, pszSrcGeolocArray,
2177 7 : /* bIsSource= */ true);
2178 7 : if (aosGeolocMD.empty())
2179 : {
2180 2 : GDALDestroyGenImgProjTransformer(psInfo);
2181 2 : return nullptr;
2182 : }
2183 5 : papszMD = aosGeolocMD.List();
2184 : }
2185 :
2186 42 : psInfo->pSrcTransformArg = GDALCreateGeoLocTransformerEx(
2187 : hSrcDS, papszMD, FALSE, nullptr, papszOptions);
2188 42 : if (psInfo->pSrcTransformArg == nullptr)
2189 : {
2190 1 : GDALDestroyGenImgProjTransformer(psInfo);
2191 1 : return nullptr;
2192 : }
2193 41 : psInfo->pSrcTransformer = GDALGeoLocTransform;
2194 41 : if (pszSrcSRS == nullptr)
2195 : {
2196 41 : pszSrcSRS = CSLFetchNameValue(papszMD, "SRS");
2197 41 : if (pszSrcSRS)
2198 : {
2199 39 : oSrcSRS.SetFromUserInput(pszSrcSRS);
2200 39 : oSrcSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
2201 : }
2202 : }
2203 : }
2204 :
2205 7 : else if (pszMethod != nullptr)
2206 : {
2207 0 : CPLError(CE_Failure, CPLE_AppDefined,
2208 : "Unable to compute a %s based transformation between "
2209 : "pixel/line and georeferenced coordinates for %s.",
2210 : pszMethod, GDALGetDescription(hSrcDS));
2211 :
2212 0 : GDALDestroyGenImgProjTransformer(psInfo);
2213 0 : return nullptr;
2214 : }
2215 :
2216 : else
2217 : {
2218 7 : CPLError(CE_Failure, CPLE_AppDefined,
2219 : "Unable to compute a transformation between pixel/line "
2220 : "and georeferenced coordinates for %s. "
2221 : "There is no affine transformation and no GCPs. "
2222 : "Specify transformation option SRC_METHOD=NO_GEOTRANSFORM "
2223 : "to bypass this check.",
2224 : GDALGetDescription(hSrcDS));
2225 :
2226 7 : GDALDestroyGenImgProjTransformer(psInfo);
2227 7 : return nullptr;
2228 : }
2229 :
2230 : /* -------------------------------------------------------------------- */
2231 : /* Handle optional source approximation transformer. */
2232 : /* -------------------------------------------------------------------- */
2233 1378 : if (psInfo->pSrcTransformer)
2234 : {
2235 : const char *pszSrcApproxErrorFwd =
2236 121 : CSLFetchNameValue(papszOptions, "SRC_APPROX_ERROR_IN_SRS_UNIT");
2237 : const char *pszSrcApproxErrorReverse =
2238 121 : CSLFetchNameValue(papszOptions, "SRC_APPROX_ERROR_IN_PIXEL");
2239 121 : if (pszSrcApproxErrorFwd && pszSrcApproxErrorReverse)
2240 : {
2241 1 : void *pArg = GDALCreateApproxTransformer2(
2242 : psInfo->pSrcTransformer, psInfo->pSrcTransformArg,
2243 : CPLAtof(pszSrcApproxErrorFwd),
2244 : CPLAtof(pszSrcApproxErrorReverse));
2245 1 : if (pArg == nullptr)
2246 : {
2247 0 : GDALDestroyGenImgProjTransformer(psInfo);
2248 0 : return nullptr;
2249 : }
2250 1 : psInfo->pSrcTransformArg = pArg;
2251 1 : psInfo->pSrcTransformer = GDALApproxTransform;
2252 1 : GDALApproxTransformerOwnsSubtransformer(psInfo->pSrcTransformArg,
2253 : TRUE);
2254 : }
2255 : }
2256 :
2257 : /* -------------------------------------------------------------------- */
2258 : /* Get forward and inverse geotransform for destination image. */
2259 : /* If we have no destination use a unit transform. */
2260 : /* -------------------------------------------------------------------- */
2261 1378 : const char *pszDstMethod = CSLFetchNameValue(papszOptions, "DST_METHOD");
2262 : const char *pszDstGeolocArray =
2263 1378 : CSLFetchNameValue(papszOptions, "DST_GEOLOC_ARRAY");
2264 1378 : if (pszDstMethod == nullptr && pszDstGeolocArray != nullptr)
2265 2 : pszDstMethod = "GEOLOC_ARRAY";
2266 :
2267 1378 : if (!hDstDS ||
2268 5 : (pszDstMethod != nullptr && EQUAL(pszDstMethod, "NO_GEOTRANSFORM")))
2269 : {
2270 1031 : psInfo->adfDstGeoTransform[0] = 0.0;
2271 1031 : psInfo->adfDstGeoTransform[1] = 1.0;
2272 1031 : psInfo->adfDstGeoTransform[2] = 0.0;
2273 1031 : psInfo->adfDstGeoTransform[3] = 0.0;
2274 1031 : psInfo->adfDstGeoTransform[4] = 0.0;
2275 1031 : psInfo->adfDstGeoTransform[5] = 1.0;
2276 1031 : memcpy(psInfo->adfDstInvGeoTransform, psInfo->adfDstGeoTransform,
2277 : sizeof(double) * 6);
2278 : }
2279 692 : else if ((pszDstMethod == nullptr || EQUAL(pszDstMethod, "GEOTRANSFORM")) &&
2280 345 : GDALGetGeoTransform(hDstDS, psInfo->adfDstGeoTransform) == CE_None)
2281 : {
2282 338 : if (pszDstSRS == nullptr)
2283 : {
2284 256 : auto hSRS = GDALGetSpatialRef(hDstDS);
2285 256 : if (hSRS)
2286 175 : oDstSRS = *(OGRSpatialReference::FromHandle(hSRS));
2287 : }
2288 338 : if (!GDALInvGeoTransform(psInfo->adfDstGeoTransform,
2289 338 : psInfo->adfDstInvGeoTransform))
2290 : {
2291 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot invert geotransform");
2292 0 : GDALDestroyGenImgProjTransformer(psInfo);
2293 0 : return nullptr;
2294 : }
2295 : }
2296 9 : else if (bGCPUseOK &&
2297 2 : (pszDstMethod == nullptr ||
2298 2 : EQUAL(pszDstMethod, "GCP_POLYNOMIAL")) &&
2299 18 : GDALGetGCPCount(hDstDS) > 0 && nOrder >= 0)
2300 : {
2301 0 : if (pszDstSRS == nullptr)
2302 : {
2303 0 : auto hSRS = GDALGetGCPSpatialRef(hDstDS);
2304 0 : if (hSRS)
2305 0 : oDstSRS = *(OGRSpatialReference::FromHandle(hSRS));
2306 : }
2307 :
2308 0 : const auto nGCPCount = GDALGetGCPCount(hDstDS);
2309 0 : auto pasGCPList = GDALDuplicateGCPs(nGCPCount, GDALGetGCPs(hDstDS));
2310 0 : GDALGCPAntimeridianUnwrap(nGCPCount, pasGCPList, oDstSRS, papszOptions);
2311 :
2312 0 : if (bRefine)
2313 : {
2314 0 : psInfo->pDstTransformArg = GDALCreateGCPRefineTransformer(
2315 : nGCPCount, pasGCPList, nOrder, FALSE, dfTolerance,
2316 : nMinimumGcps);
2317 : }
2318 : else
2319 : {
2320 0 : psInfo->pDstTransformArg =
2321 0 : GDALCreateGCPTransformer(nGCPCount, pasGCPList, nOrder, FALSE);
2322 : }
2323 :
2324 0 : GDALDeinitGCPs(nGCPCount, pasGCPList);
2325 0 : CPLFree(pasGCPList);
2326 :
2327 0 : if (psInfo->pDstTransformArg == nullptr)
2328 : {
2329 0 : GDALDestroyGenImgProjTransformer(psInfo);
2330 0 : return nullptr;
2331 : }
2332 0 : psInfo->pDstTransformer = GDALGCPTransform;
2333 : }
2334 9 : else if (bGCPUseOK && GDALGetGCPCount(hDstDS) > 0 && nOrder <= 0 &&
2335 0 : (pszDstMethod == nullptr || EQUAL(pszDstMethod, "GCP_TPS")))
2336 : {
2337 0 : if (pszDstSRS == nullptr)
2338 : {
2339 0 : auto hSRS = GDALGetGCPSpatialRef(hDstDS);
2340 0 : if (hSRS)
2341 0 : oDstSRS = *(OGRSpatialReference::FromHandle(hSRS));
2342 : }
2343 :
2344 0 : const auto nGCPCount = GDALGetGCPCount(hDstDS);
2345 0 : auto pasGCPList = GDALDuplicateGCPs(nGCPCount, GDALGetGCPs(hDstDS));
2346 0 : GDALGCPAntimeridianUnwrap(nGCPCount, pasGCPList, oDstSRS, papszOptions);
2347 :
2348 0 : psInfo->pDstTransformArg = GDALCreateTPSTransformerInt(
2349 : nGCPCount, pasGCPList, FALSE, papszOptions);
2350 :
2351 0 : GDALDeinitGCPs(nGCPCount, pasGCPList);
2352 0 : CPLFree(pasGCPList);
2353 :
2354 0 : if (psInfo->pDstTransformArg == nullptr)
2355 : {
2356 0 : GDALDestroyGenImgProjTransformer(psInfo);
2357 0 : return nullptr;
2358 : }
2359 0 : psInfo->pDstTransformer = GDALTPSTransform;
2360 : }
2361 2 : else if ((pszDstMethod == nullptr || EQUAL(pszDstMethod, "RPC")) &&
2362 13 : (papszMD = GDALGetMetadata(hDstDS, "RPC")) != nullptr &&
2363 2 : GDALExtractRPCInfoV2(papszMD, &sRPCInfo))
2364 : {
2365 2 : psInfo->pDstTransformArg =
2366 2 : GDALCreateRPCTransformerV2(&sRPCInfo, FALSE, 0, papszOptions);
2367 2 : if (psInfo->pDstTransformArg == nullptr)
2368 : {
2369 0 : GDALDestroyGenImgProjTransformer(psInfo);
2370 0 : return nullptr;
2371 : }
2372 2 : psInfo->pDstTransformer = GDALRPCTransform;
2373 2 : if (pszDstSRS == nullptr)
2374 : {
2375 2 : oDstSRS.SetFromUserInput(SRS_WKT_WGS84_LAT_LONG);
2376 2 : oDstSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
2377 : }
2378 : }
2379 14 : else if ((pszDstMethod == nullptr || EQUAL(pszDstMethod, "GEOLOC_ARRAY")) &&
2380 7 : ((papszMD = GDALGetMetadata(hDstDS, "GEOLOCATION")) != nullptr ||
2381 : pszDstGeolocArray != nullptr))
2382 : {
2383 7 : CPLStringList aosGeolocMD; // keep in this scope
2384 7 : if (pszDstGeolocArray != nullptr)
2385 : {
2386 2 : if (papszMD != nullptr)
2387 : {
2388 0 : CPLError(
2389 : CE_Warning, CPLE_AppDefined,
2390 : "Both GEOLOCATION metadata domain on the target dataset "
2391 : "and DST_GEOLOC_ARRAY transformer option are set. "
2392 : "Only using the later.");
2393 : }
2394 : aosGeolocMD =
2395 4 : GDALCreateGeolocationMetadata(hDstDS, pszDstGeolocArray,
2396 2 : /* bIsSource= */ false);
2397 2 : if (aosGeolocMD.empty())
2398 : {
2399 1 : GDALDestroyGenImgProjTransformer(psInfo);
2400 1 : return nullptr;
2401 : }
2402 1 : papszMD = aosGeolocMD.List();
2403 : }
2404 :
2405 6 : psInfo->pDstTransformArg = GDALCreateGeoLocTransformerEx(
2406 : hDstDS, papszMD, FALSE, nullptr, papszOptions);
2407 6 : if (psInfo->pDstTransformArg == nullptr)
2408 : {
2409 1 : GDALDestroyGenImgProjTransformer(psInfo);
2410 1 : return nullptr;
2411 : }
2412 5 : psInfo->pDstTransformer = GDALGeoLocTransform;
2413 5 : if (pszDstSRS == nullptr)
2414 : {
2415 5 : pszDstSRS = CSLFetchNameValue(papszMD, "SRS");
2416 5 : if (pszDstSRS)
2417 : {
2418 5 : oDstSRS.SetFromUserInput(pszDstSRS);
2419 5 : oDstSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
2420 : }
2421 : }
2422 : }
2423 : else
2424 : {
2425 0 : CPLError(CE_Failure, CPLE_AppDefined,
2426 : "Unable to compute a transformation between pixel/line "
2427 : "and georeferenced coordinates for %s. "
2428 : "There is no affine transformation and no GCPs. "
2429 : "Specify transformation option DST_METHOD=NO_GEOTRANSFORM "
2430 : "to bypass this check.",
2431 : GDALGetDescription(hDstDS));
2432 :
2433 0 : GDALDestroyGenImgProjTransformer(psInfo);
2434 0 : return nullptr;
2435 : }
2436 :
2437 : /* -------------------------------------------------------------------- */
2438 : /* Handle optional destination approximation transformer. */
2439 : /* -------------------------------------------------------------------- */
2440 1376 : if (psInfo->pDstTransformer)
2441 : {
2442 : const char *pszDstApproxErrorFwd =
2443 7 : CSLFetchNameValue(papszOptions, "DST_APPROX_ERROR_IN_PIXEL");
2444 : const char *pszDstApproxErrorReverse =
2445 7 : CSLFetchNameValue(papszOptions, "DST_APPROX_ERROR_IN_SRS_UNIT");
2446 7 : if (pszDstApproxErrorFwd && pszDstApproxErrorReverse)
2447 : {
2448 0 : void *pArg = GDALCreateApproxTransformer2(
2449 : psInfo->pDstTransformer, psInfo->pDstTransformArg,
2450 : CPLAtof(pszDstApproxErrorFwd),
2451 : CPLAtof(pszDstApproxErrorReverse));
2452 0 : if (pArg == nullptr)
2453 : {
2454 0 : GDALDestroyGenImgProjTransformer(psInfo);
2455 0 : return nullptr;
2456 : }
2457 0 : psInfo->pDstTransformArg = pArg;
2458 0 : psInfo->pDstTransformer = GDALApproxTransform;
2459 0 : GDALApproxTransformerOwnsSubtransformer(psInfo->pDstTransformArg,
2460 : TRUE);
2461 : }
2462 : }
2463 :
2464 : /* -------------------------------------------------------------------- */
2465 : /* Setup reprojection. */
2466 : /* -------------------------------------------------------------------- */
2467 :
2468 1376 : if (CPLFetchBool(papszOptions, "STRIP_VERT_CS", false))
2469 : {
2470 0 : if (oSrcSRS.IsCompound())
2471 : {
2472 0 : oSrcSRS.StripVertical();
2473 : }
2474 0 : if (oDstSRS.IsCompound())
2475 : {
2476 0 : oDstSRS.StripVertical();
2477 : }
2478 : }
2479 :
2480 : const bool bMayInsertCenterLong =
2481 2368 : (bCanUseSrcGeoTransform && !oSrcSRS.IsEmpty() && hSrcDS &&
2482 992 : CPLFetchBool(papszOptions, "INSERT_CENTER_LONG", true));
2483 : const char *pszSrcCoordEpoch =
2484 1376 : CSLFetchNameValue(papszOptions, "SRC_COORDINATE_EPOCH");
2485 : const char *pszDstCoordEpoch =
2486 1376 : CSLFetchNameValue(papszOptions, "DST_COORDINATE_EPOCH");
2487 2529 : if ((!oSrcSRS.IsEmpty() && !oDstSRS.IsEmpty() &&
2488 1052 : (pszSrcCoordEpoch || pszDstCoordEpoch || !oSrcSRS.IsSame(&oDstSRS) ||
2489 2529 : (oSrcSRS.IsGeographic() && bMayInsertCenterLong))) ||
2490 : pszCO)
2491 : {
2492 732 : CPLStringList aosOptions;
2493 :
2494 732 : if (bMayInsertCenterLong)
2495 : {
2496 715 : InsertCenterLong(hSrcDS, &oSrcSRS, aosOptions);
2497 : }
2498 :
2499 732 : if (CPLFetchBool(papszOptions, "PROMOTE_TO_3D", false))
2500 : {
2501 19 : oSrcSRS.PromoteTo3D(nullptr);
2502 19 : oDstSRS.PromoteTo3D(nullptr);
2503 : }
2504 :
2505 732 : if (!(dfWestLongitudeDeg == 0.0 && dfSouthLatitudeDeg == 0.0 &&
2506 29 : dfEastLongitudeDeg == 0.0 && dfNorthLatitudeDeg == 0.0))
2507 : {
2508 : aosOptions.SetNameValue(
2509 : "AREA_OF_INTEREST",
2510 : CPLSPrintf("%.16g,%.16g,%.16g,%.16g", dfWestLongitudeDeg,
2511 : dfSouthLatitudeDeg, dfEastLongitudeDeg,
2512 703 : dfNorthLatitudeDeg));
2513 : }
2514 732 : if (pszCO)
2515 : {
2516 7 : aosOptions.SetNameValue("COORDINATE_OPERATION", pszCO);
2517 : }
2518 :
2519 : const char *pszCoordEpoch =
2520 732 : CSLFetchNameValue(papszOptions, "COORDINATE_EPOCH");
2521 732 : if (pszCoordEpoch)
2522 : {
2523 1 : aosOptions.SetNameValue("COORDINATE_EPOCH", pszCoordEpoch);
2524 : }
2525 :
2526 732 : if (pszSrcCoordEpoch)
2527 : {
2528 0 : aosOptions.SetNameValue("SRC_COORDINATE_EPOCH", pszSrcCoordEpoch);
2529 0 : oSrcSRS.SetCoordinateEpoch(CPLAtof(pszSrcCoordEpoch));
2530 : }
2531 :
2532 732 : if (pszDstCoordEpoch)
2533 : {
2534 0 : aosOptions.SetNameValue("DST_COORDINATE_EPOCH", pszDstCoordEpoch);
2535 0 : oDstSRS.SetCoordinateEpoch(CPLAtof(pszDstCoordEpoch));
2536 : }
2537 :
2538 736 : psInfo->pReprojectArg = GDALCreateReprojectionTransformerEx(
2539 732 : !oSrcSRS.IsEmpty() ? OGRSpatialReference::ToHandle(&oSrcSRS)
2540 : : nullptr,
2541 732 : !oDstSRS.IsEmpty() ? OGRSpatialReference::ToHandle(&oDstSRS)
2542 : : nullptr,
2543 732 : aosOptions.List());
2544 :
2545 732 : if (pszCO)
2546 : {
2547 7 : psInfo->bHasCustomTransformationPipeline = true;
2548 : }
2549 :
2550 732 : if (psInfo->pReprojectArg == nullptr)
2551 : {
2552 1 : GDALDestroyGenImgProjTransformer(psInfo);
2553 1 : return nullptr;
2554 : }
2555 731 : psInfo->pReproject = GDALReprojectionTransform;
2556 :
2557 : /* --------------------------------------------------------------------
2558 : */
2559 : /* Handle optional reprojection approximation transformer. */
2560 : /* --------------------------------------------------------------------
2561 : */
2562 731 : const char *psApproxErrorFwd = CSLFetchNameValue(
2563 : papszOptions, "REPROJECTION_APPROX_ERROR_IN_DST_SRS_UNIT");
2564 731 : const char *psApproxErrorReverse = CSLFetchNameValue(
2565 : papszOptions, "REPROJECTION_APPROX_ERROR_IN_SRC_SRS_UNIT");
2566 731 : if (psApproxErrorFwd && psApproxErrorReverse)
2567 : {
2568 1 : void *pArg = GDALCreateApproxTransformer2(
2569 : psInfo->pReproject, psInfo->pReprojectArg,
2570 : CPLAtof(psApproxErrorFwd), CPLAtof(psApproxErrorReverse));
2571 1 : if (pArg == nullptr)
2572 : {
2573 0 : GDALDestroyGenImgProjTransformer(psInfo);
2574 0 : return nullptr;
2575 : }
2576 1 : psInfo->pReprojectArg = pArg;
2577 1 : psInfo->pReproject = GDALApproxTransform;
2578 1 : GDALApproxTransformerOwnsSubtransformer(psInfo->pReprojectArg,
2579 : TRUE);
2580 : }
2581 : }
2582 :
2583 1375 : return psInfo;
2584 : }
2585 :
2586 : /************************************************************************/
2587 : /* GDALRefreshGenImgProjTransformer() */
2588 : /************************************************************************/
2589 :
2590 948 : void GDALRefreshGenImgProjTransformer(void *hTransformArg)
2591 : {
2592 948 : GDALGenImgProjTransformInfo *psInfo =
2593 : static_cast<GDALGenImgProjTransformInfo *>(hTransformArg);
2594 :
2595 1589 : if (psInfo->pReprojectArg &&
2596 641 : psInfo->bCheckWithInvertPROJ != GetCurrentCheckWithInvertPROJ())
2597 : {
2598 66 : psInfo->bCheckWithInvertPROJ = !psInfo->bCheckWithInvertPROJ;
2599 :
2600 : CPLXMLNode *psXML =
2601 66 : GDALSerializeTransformer(psInfo->pReproject, psInfo->pReprojectArg);
2602 66 : GDALDestroyTransformer(psInfo->pReprojectArg);
2603 66 : GDALDeserializeTransformer(psXML, &psInfo->pReproject,
2604 : &psInfo->pReprojectArg);
2605 66 : CPLDestroyXMLNode(psXML);
2606 : }
2607 948 : }
2608 :
2609 : /************************************************************************/
2610 : /* GDALCreateGenImgProjTransformer3() */
2611 : /************************************************************************/
2612 :
2613 : /**
2614 : * Create image to image transformer.
2615 : *
2616 : * This function creates a transformation object that maps from pixel/line
2617 : * coordinates on one image to pixel/line coordinates on another image. The
2618 : * images may potentially be georeferenced in different coordinate systems,
2619 : * and may used GCPs to map between their pixel/line coordinates and
2620 : * georeferenced coordinates (as opposed to the default assumption that their
2621 : * geotransform should be used).
2622 : *
2623 : * This transformer potentially performs three concatenated transformations.
2624 : *
2625 : * The first stage is from source image pixel/line coordinates to source
2626 : * image georeferenced coordinates, and may be done using the geotransform,
2627 : * or if not defined using a polynomial model derived from GCPs. If GCPs
2628 : * are used this stage is accomplished using GDALGCPTransform().
2629 : *
2630 : * The second stage is to change projections from the source coordinate system
2631 : * to the destination coordinate system, assuming they differ. This is
2632 : * accomplished internally using GDALReprojectionTransform().
2633 : *
2634 : * The third stage is converting from destination image georeferenced
2635 : * coordinates to destination image coordinates. This is done using the
2636 : * destination image geotransform, or if not available, using a polynomial
2637 : * model derived from GCPs. If GCPs are used this stage is accomplished using
2638 : * GDALGCPTransform(). This stage is skipped if hDstDS is NULL when the
2639 : * transformation is created.
2640 : *
2641 : * @param pszSrcWKT source WKT (or NULL).
2642 : * @param padfSrcGeoTransform source geotransform (or NULL).
2643 : * @param pszDstWKT destination WKT (or NULL).
2644 : * @param padfDstGeoTransform destination geotransform (or NULL).
2645 : *
2646 : * @return handle suitable for use GDALGenImgProjTransform(), and to be
2647 : * deallocated with GDALDestroyGenImgProjTransformer() or NULL on failure.
2648 : */
2649 :
2650 0 : void *GDALCreateGenImgProjTransformer3(const char *pszSrcWKT,
2651 : const double *padfSrcGeoTransform,
2652 : const char *pszDstWKT,
2653 : const double *padfDstGeoTransform)
2654 :
2655 : {
2656 0 : OGRSpatialReference oSrcSRS;
2657 0 : if (pszSrcWKT)
2658 : {
2659 0 : oSrcSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
2660 0 : if (pszSrcWKT[0] != '\0' &&
2661 0 : oSrcSRS.importFromWkt(pszSrcWKT) != OGRERR_NONE)
2662 : {
2663 0 : CPLError(CE_Failure, CPLE_AppDefined,
2664 : "Failed to import coordinate system `%s'.", pszSrcWKT);
2665 0 : return nullptr;
2666 : }
2667 : }
2668 :
2669 0 : OGRSpatialReference oDstSRS;
2670 0 : if (pszDstWKT)
2671 : {
2672 0 : oDstSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
2673 0 : if (pszDstWKT[0] != '\0' &&
2674 0 : oDstSRS.importFromWkt(pszDstWKT) != OGRERR_NONE)
2675 : {
2676 0 : CPLError(CE_Failure, CPLE_AppDefined,
2677 : "Failed to import coordinate system `%s'.", pszDstWKT);
2678 0 : return nullptr;
2679 : }
2680 : }
2681 0 : return GDALCreateGenImgProjTransformer4(
2682 : OGRSpatialReference::ToHandle(&oSrcSRS), padfSrcGeoTransform,
2683 0 : OGRSpatialReference::ToHandle(&oDstSRS), padfDstGeoTransform, nullptr);
2684 : }
2685 :
2686 : /************************************************************************/
2687 : /* GDALCreateGenImgProjTransformer4() */
2688 : /************************************************************************/
2689 :
2690 : /**
2691 : * Create image to image transformer.
2692 : *
2693 : * Similar to GDALCreateGenImgProjTransformer3(), except that it takes
2694 : * OGRSpatialReferenceH objects and options.
2695 : * The options are the ones supported by GDALCreateReprojectionTransformerEx()
2696 : *
2697 : * @since GDAL 3.0
2698 : */
2699 16 : void *GDALCreateGenImgProjTransformer4(OGRSpatialReferenceH hSrcSRS,
2700 : const double *padfSrcGeoTransform,
2701 : OGRSpatialReferenceH hDstSRS,
2702 : const double *padfDstGeoTransform,
2703 : const char *const *papszOptions)
2704 : {
2705 : /* -------------------------------------------------------------------- */
2706 : /* Initialize the transform info. */
2707 : /* -------------------------------------------------------------------- */
2708 : GDALGenImgProjTransformInfo *psInfo =
2709 16 : GDALCreateGenImgProjTransformerInternal();
2710 :
2711 : /* -------------------------------------------------------------------- */
2712 : /* Get forward and inverse geotransform for the source image. */
2713 : /* -------------------------------------------------------------------- */
2714 16 : if (padfSrcGeoTransform)
2715 : {
2716 16 : memcpy(psInfo->adfSrcGeoTransform, padfSrcGeoTransform,
2717 : sizeof(psInfo->adfSrcGeoTransform));
2718 16 : if (!GDALInvGeoTransform(psInfo->adfSrcGeoTransform,
2719 16 : psInfo->adfSrcInvGeoTransform))
2720 : {
2721 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot invert geotransform");
2722 1 : GDALDestroyGenImgProjTransformer(psInfo);
2723 1 : return nullptr;
2724 : }
2725 : }
2726 : else
2727 : {
2728 0 : psInfo->adfSrcGeoTransform[0] = 0.0;
2729 0 : psInfo->adfSrcGeoTransform[1] = 1.0;
2730 0 : psInfo->adfSrcGeoTransform[2] = 0.0;
2731 0 : psInfo->adfSrcGeoTransform[3] = 0.0;
2732 0 : psInfo->adfSrcGeoTransform[4] = 0.0;
2733 0 : psInfo->adfSrcGeoTransform[5] = 1.0;
2734 0 : memcpy(psInfo->adfSrcInvGeoTransform, psInfo->adfSrcGeoTransform,
2735 : sizeof(double) * 6);
2736 : }
2737 :
2738 : /* -------------------------------------------------------------------- */
2739 : /* Setup reprojection. */
2740 : /* -------------------------------------------------------------------- */
2741 15 : OGRSpatialReference *poSrcSRS = OGRSpatialReference::FromHandle(hSrcSRS);
2742 15 : OGRSpatialReference *poDstSRS = OGRSpatialReference::FromHandle(hDstSRS);
2743 30 : if (!poSrcSRS->IsEmpty() && !poDstSRS->IsEmpty() &&
2744 15 : !poSrcSRS->IsSame(poDstSRS))
2745 : {
2746 4 : psInfo->pReprojectArg =
2747 4 : GDALCreateReprojectionTransformerEx(hSrcSRS, hDstSRS, papszOptions);
2748 4 : if (psInfo->pReprojectArg == nullptr)
2749 : {
2750 1 : GDALDestroyGenImgProjTransformer(psInfo);
2751 1 : return nullptr;
2752 : }
2753 3 : psInfo->pReproject = GDALReprojectionTransform;
2754 : }
2755 :
2756 : /* -------------------------------------------------------------------- */
2757 : /* Get forward and inverse geotransform for destination image. */
2758 : /* If we have no destination matrix use a unit transform. */
2759 : /* -------------------------------------------------------------------- */
2760 14 : if (padfDstGeoTransform)
2761 : {
2762 14 : memcpy(psInfo->adfDstGeoTransform, padfDstGeoTransform,
2763 : sizeof(psInfo->adfDstGeoTransform));
2764 14 : if (!GDALInvGeoTransform(psInfo->adfDstGeoTransform,
2765 14 : psInfo->adfDstInvGeoTransform))
2766 : {
2767 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot invert geotransform");
2768 1 : GDALDestroyGenImgProjTransformer(psInfo);
2769 1 : return nullptr;
2770 : }
2771 : }
2772 : else
2773 : {
2774 0 : psInfo->adfDstGeoTransform[0] = 0.0;
2775 0 : psInfo->adfDstGeoTransform[1] = 1.0;
2776 0 : psInfo->adfDstGeoTransform[2] = 0.0;
2777 0 : psInfo->adfDstGeoTransform[3] = 0.0;
2778 0 : psInfo->adfDstGeoTransform[4] = 0.0;
2779 0 : psInfo->adfDstGeoTransform[5] = 1.0;
2780 0 : memcpy(psInfo->adfDstInvGeoTransform, psInfo->adfDstGeoTransform,
2781 : sizeof(double) * 6);
2782 : }
2783 :
2784 13 : return psInfo;
2785 : }
2786 :
2787 : /************************************************************************/
2788 : /* GDALSetGenImgProjTransformerDstGeoTransform() */
2789 : /************************************************************************/
2790 :
2791 : /**
2792 : * Set GenImgProj output geotransform.
2793 : *
2794 : * Normally the "destination geotransform", or transformation between
2795 : * georeferenced output coordinates and pixel/line coordinates on the
2796 : * destination file is extracted from the destination file by
2797 : * GDALCreateGenImgProjTransformer() and stored in the GenImgProj private
2798 : * info. However, sometimes it is inconvenient to have an output file
2799 : * handle with appropriate geotransform information when creating the
2800 : * transformation. For these cases, this function can be used to apply
2801 : * the destination geotransform.
2802 : *
2803 : * @param hTransformArg the handle to update.
2804 : * @param padfGeoTransform the destination geotransform to apply (six doubles).
2805 : */
2806 :
2807 808 : void GDALSetGenImgProjTransformerDstGeoTransform(void *hTransformArg,
2808 : const double *padfGeoTransform)
2809 :
2810 : {
2811 808 : VALIDATE_POINTER0(hTransformArg,
2812 : "GDALSetGenImgProjTransformerDstGeoTransform");
2813 :
2814 808 : GDALGenImgProjTransformInfo *psInfo =
2815 : static_cast<GDALGenImgProjTransformInfo *>(hTransformArg);
2816 :
2817 808 : memcpy(psInfo->adfDstGeoTransform, padfGeoTransform, sizeof(double) * 6);
2818 808 : if (!GDALInvGeoTransform(psInfo->adfDstGeoTransform,
2819 808 : psInfo->adfDstInvGeoTransform))
2820 : {
2821 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot invert geotransform");
2822 : }
2823 : }
2824 :
2825 : /************************************************************************/
2826 : /* GDALDestroyGenImgProjTransformer() */
2827 : /************************************************************************/
2828 :
2829 : /**
2830 : * GenImgProjTransformer deallocator.
2831 : *
2832 : * This function is used to deallocate the handle created with
2833 : * GDALCreateGenImgProjTransformer().
2834 : *
2835 : * @param hTransformArg the handle to deallocate.
2836 : */
2837 :
2838 1637 : void GDALDestroyGenImgProjTransformer(void *hTransformArg)
2839 :
2840 : {
2841 1637 : if (hTransformArg == nullptr)
2842 0 : return;
2843 :
2844 1637 : GDALGenImgProjTransformInfo *psInfo =
2845 : static_cast<GDALGenImgProjTransformInfo *>(hTransformArg);
2846 :
2847 1637 : if (psInfo->pSrcTransformArg != nullptr)
2848 144 : GDALDestroyTransformer(psInfo->pSrcTransformArg);
2849 :
2850 1637 : if (psInfo->pDstTransformArg != nullptr)
2851 7 : GDALDestroyTransformer(psInfo->pDstTransformArg);
2852 :
2853 1637 : if (psInfo->pReprojectArg != nullptr)
2854 842 : GDALDestroyTransformer(psInfo->pReprojectArg);
2855 :
2856 1637 : CPLFree(psInfo);
2857 : }
2858 :
2859 : /************************************************************************/
2860 : /* GDALGenImgProjTransform() */
2861 : /************************************************************************/
2862 :
2863 : /**
2864 : * Perform general image reprojection transformation.
2865 : *
2866 : * Actually performs the transformation setup in
2867 : * GDALCreateGenImgProjTransformer(). This function matches the signature
2868 : * required by the GDALTransformerFunc(), and more details on the arguments
2869 : * can be found in that topic.
2870 : */
2871 :
2872 : #ifdef DEBUG_APPROX_TRANSFORMER
2873 : int countGDALGenImgProjTransform = 0;
2874 : #endif
2875 :
2876 1922290 : int GDALGenImgProjTransform(void *pTransformArgIn, int bDstToSrc,
2877 : int nPointCount, double *padfX, double *padfY,
2878 : double *padfZ, int *panSuccess)
2879 : {
2880 1922290 : GDALGenImgProjTransformInfo *psInfo =
2881 : static_cast<GDALGenImgProjTransformInfo *>(pTransformArgIn);
2882 :
2883 : #ifdef DEBUG_APPROX_TRANSFORMER
2884 : CPLAssert(nPointCount > 0);
2885 : countGDALGenImgProjTransform += nPointCount;
2886 : #endif
2887 :
2888 21553700 : for (int i = 0; i < nPointCount; i++)
2889 : {
2890 19631400 : panSuccess[i] = (padfX[i] != HUGE_VAL && padfY[i] != HUGE_VAL);
2891 : }
2892 :
2893 : /* -------------------------------------------------------------------- */
2894 : /* Convert from src (dst) pixel/line to src (dst) */
2895 : /* georeferenced coordinates. */
2896 : /* -------------------------------------------------------------------- */
2897 1922290 : double *padfGeoTransform = nullptr;
2898 1922290 : void *pTransformArg = nullptr;
2899 1922290 : GDALTransformerFunc pTransformer = nullptr;
2900 1922290 : if (bDstToSrc)
2901 : {
2902 1550370 : padfGeoTransform = psInfo->adfDstGeoTransform;
2903 1550370 : pTransformArg = psInfo->pDstTransformArg;
2904 1550370 : pTransformer = psInfo->pDstTransformer;
2905 : }
2906 : else
2907 : {
2908 371922 : padfGeoTransform = psInfo->adfSrcGeoTransform;
2909 371922 : pTransformArg = psInfo->pSrcTransformArg;
2910 371922 : pTransformer = psInfo->pSrcTransformer;
2911 : }
2912 :
2913 1922290 : if (pTransformArg != nullptr)
2914 : {
2915 41607 : if (!pTransformer(pTransformArg, FALSE, nPointCount, padfX, padfY,
2916 : padfZ, panSuccess))
2917 0 : return FALSE;
2918 : }
2919 : else
2920 : {
2921 21433400 : for (int i = 0; i < nPointCount; i++)
2922 : {
2923 19552800 : if (!panSuccess[i])
2924 1758 : continue;
2925 :
2926 19551000 : const double dfNewX = padfGeoTransform[0] +
2927 19551000 : padfX[i] * padfGeoTransform[1] +
2928 19551000 : padfY[i] * padfGeoTransform[2];
2929 19551000 : const double dfNewY = padfGeoTransform[3] +
2930 19551000 : padfX[i] * padfGeoTransform[4] +
2931 19551000 : padfY[i] * padfGeoTransform[5];
2932 :
2933 19551000 : padfX[i] = dfNewX;
2934 19551000 : padfY[i] = dfNewY;
2935 : }
2936 : }
2937 :
2938 : /* -------------------------------------------------------------------- */
2939 : /* Reproject if needed. */
2940 : /* -------------------------------------------------------------------- */
2941 1922290 : if (psInfo->pReprojectArg)
2942 : {
2943 1609160 : if (!psInfo->pReproject(psInfo->pReprojectArg, bDstToSrc, nPointCount,
2944 : padfX, padfY, padfZ, panSuccess))
2945 0 : return FALSE;
2946 : }
2947 :
2948 : /* -------------------------------------------------------------------- */
2949 : /* Convert dst (src) georef coordinates back to pixel/line. */
2950 : /* -------------------------------------------------------------------- */
2951 1922280 : if (bDstToSrc)
2952 : {
2953 1550350 : padfGeoTransform = psInfo->adfSrcInvGeoTransform;
2954 1550350 : pTransformArg = psInfo->pSrcTransformArg;
2955 1550350 : pTransformer = psInfo->pSrcTransformer;
2956 : }
2957 : else
2958 : {
2959 371925 : padfGeoTransform = psInfo->adfDstInvGeoTransform;
2960 371925 : pTransformArg = psInfo->pDstTransformArg;
2961 371925 : pTransformer = psInfo->pDstTransformer;
2962 : }
2963 :
2964 1922280 : if (pTransformArg != nullptr)
2965 : {
2966 51665 : if (!pTransformer(pTransformArg, TRUE, nPointCount, padfX, padfY, padfZ,
2967 : panSuccess))
2968 0 : return FALSE;
2969 : }
2970 : else
2971 : {
2972 20648500 : for (int i = 0; i < nPointCount; i++)
2973 : {
2974 18777900 : if (!panSuccess[i])
2975 3523210 : continue;
2976 :
2977 15254700 : const double dfNewX = padfGeoTransform[0] +
2978 15254700 : padfX[i] * padfGeoTransform[1] +
2979 15254700 : padfY[i] * padfGeoTransform[2];
2980 15254700 : const double dfNewY = padfGeoTransform[3] +
2981 15254700 : padfX[i] * padfGeoTransform[4] +
2982 15254700 : padfY[i] * padfGeoTransform[5];
2983 :
2984 15254700 : padfX[i] = dfNewX;
2985 15254700 : padfY[i] = dfNewY;
2986 : }
2987 : }
2988 :
2989 1922280 : return TRUE;
2990 : }
2991 :
2992 : /************************************************************************/
2993 : /* GDALTransformLonLatToDestGenImgProjTransformer() */
2994 : /************************************************************************/
2995 :
2996 2496 : int GDALTransformLonLatToDestGenImgProjTransformer(void *hTransformArg,
2997 : double *pdfX, double *pdfY)
2998 : {
2999 2496 : GDALGenImgProjTransformInfo *psInfo =
3000 : static_cast<GDALGenImgProjTransformInfo *>(hTransformArg);
3001 :
3002 2496 : if (psInfo->pReprojectArg == nullptr ||
3003 1418 : psInfo->pReproject != GDALReprojectionTransform)
3004 1082 : return false;
3005 :
3006 1414 : GDALReprojectionTransformInfo *psReprojInfo =
3007 : static_cast<GDALReprojectionTransformInfo *>(psInfo->pReprojectArg);
3008 2828 : if (psReprojInfo->poForwardTransform == nullptr ||
3009 1414 : psReprojInfo->poForwardTransform->GetSourceCS() == nullptr)
3010 2 : return false;
3011 :
3012 1412 : double z = 0;
3013 1412 : int success = true;
3014 1412 : auto poSourceCRS = psReprojInfo->poForwardTransform->GetSourceCS();
3015 2504 : if (poSourceCRS->IsGeographic() &&
3016 1092 : std::fabs(poSourceCRS->GetAngularUnits() -
3017 1092 : CPLAtof(SRS_UA_DEGREE_CONV)) < 1e-9)
3018 : {
3019 : // Optimization to avoid creating a OGRCoordinateTransformation
3020 1090 : OGRAxisOrientation eSourceFirstAxisOrient = OAO_Other;
3021 1090 : poSourceCRS->GetAxis(nullptr, 0, &eSourceFirstAxisOrient);
3022 1090 : const auto &mapping = poSourceCRS->GetDataAxisToSRSAxisMapping();
3023 2180 : if ((mapping[0] == 2 && eSourceFirstAxisOrient == OAO_East) ||
3024 1090 : (mapping[0] == 1 && eSourceFirstAxisOrient != OAO_East))
3025 : {
3026 6 : std::swap(*pdfX, *pdfY);
3027 : }
3028 : }
3029 : else
3030 : {
3031 : auto poLongLat =
3032 322 : std::unique_ptr<OGRSpatialReference>(poSourceCRS->CloneGeogCS());
3033 322 : if (poLongLat == nullptr)
3034 0 : return false;
3035 322 : poLongLat->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3036 :
3037 : const bool bCurrentCheckWithInvertProj =
3038 322 : GetCurrentCheckWithInvertPROJ();
3039 322 : if (!bCurrentCheckWithInvertProj)
3040 322 : CPLSetThreadLocalConfigOption("CHECK_WITH_INVERT_PROJ", "YES");
3041 : auto poCT = std::unique_ptr<OGRCoordinateTransformation>(
3042 322 : OGRCreateCoordinateTransformation(poLongLat.get(), poSourceCRS));
3043 322 : if (!bCurrentCheckWithInvertProj)
3044 322 : CPLSetThreadLocalConfigOption("CHECK_WITH_INVERT_PROJ", nullptr);
3045 322 : if (poCT == nullptr)
3046 0 : return false;
3047 :
3048 322 : poCT->SetEmitErrors(false);
3049 322 : if (!poCT->Transform(1, pdfX, pdfY))
3050 2 : return false;
3051 :
3052 320 : if (!psInfo->pReproject(psInfo->pReprojectArg, false, 1, pdfX, pdfY, &z,
3053 640 : &success) ||
3054 320 : !success)
3055 : {
3056 138 : return false;
3057 : }
3058 : }
3059 :
3060 1272 : double *padfGeoTransform = psInfo->adfDstInvGeoTransform;
3061 1272 : void *pTransformArg = psInfo->pDstTransformArg;
3062 1272 : GDALTransformerFunc pTransformer = psInfo->pDstTransformer;
3063 1272 : if (pTransformArg != nullptr)
3064 : {
3065 8 : if (!pTransformer(pTransformArg, TRUE, 1, pdfX, pdfY, &z, &success) ||
3066 4 : !success)
3067 : {
3068 4 : return false;
3069 : }
3070 : }
3071 : else
3072 : {
3073 1268 : const double dfNewX = padfGeoTransform[0] +
3074 1268 : pdfX[0] * padfGeoTransform[1] +
3075 1268 : pdfY[0] * padfGeoTransform[2];
3076 1268 : const double dfNewY = padfGeoTransform[3] +
3077 1268 : pdfX[0] * padfGeoTransform[4] +
3078 1268 : pdfY[0] * padfGeoTransform[5];
3079 :
3080 1268 : pdfX[0] = dfNewX;
3081 1268 : pdfY[0] = dfNewY;
3082 : }
3083 :
3084 1268 : return true;
3085 : }
3086 :
3087 : /************************************************************************/
3088 : /* GDALSerializeGenImgProjTransformer() */
3089 : /************************************************************************/
3090 :
3091 77 : static CPLXMLNode *GDALSerializeGenImgProjTransformer(void *pTransformArg)
3092 :
3093 : {
3094 77 : GDALGenImgProjTransformInfo *psInfo =
3095 : static_cast<GDALGenImgProjTransformInfo *>(pTransformArg);
3096 :
3097 : CPLXMLNode *psTree =
3098 77 : CPLCreateXMLNode(nullptr, CXT_Element, "GenImgProjTransformer");
3099 :
3100 77 : char szWork[200] = {};
3101 :
3102 : /* -------------------------------------------------------------------- */
3103 : /* Handle source transformation. */
3104 : /* -------------------------------------------------------------------- */
3105 77 : if (psInfo->pSrcTransformArg != nullptr)
3106 : {
3107 13 : CPLXMLNode *psTransformer = GDALSerializeTransformer(
3108 : psInfo->pSrcTransformer, psInfo->pSrcTransformArg);
3109 13 : if (psTransformer != nullptr)
3110 : {
3111 : CPLXMLNode *psTransformerContainer =
3112 13 : CPLCreateXMLNode(psTree, CXT_Element,
3113 : CPLSPrintf("Src%s", psTransformer->pszValue));
3114 :
3115 13 : CPLAddXMLChild(psTransformerContainer, psTransformer);
3116 : }
3117 : }
3118 :
3119 : /* -------------------------------------------------------------------- */
3120 : /* Handle source geotransforms. */
3121 : /* -------------------------------------------------------------------- */
3122 : else
3123 : {
3124 64 : CPLsnprintf(
3125 : szWork, sizeof(szWork), "%.17g,%.17g,%.17g,%.17g,%.17g,%.17g",
3126 : psInfo->adfSrcGeoTransform[0], psInfo->adfSrcGeoTransform[1],
3127 : psInfo->adfSrcGeoTransform[2], psInfo->adfSrcGeoTransform[3],
3128 : psInfo->adfSrcGeoTransform[4], psInfo->adfSrcGeoTransform[5]);
3129 64 : CPLCreateXMLElementAndValue(psTree, "SrcGeoTransform", szWork);
3130 :
3131 64 : CPLsnprintf(
3132 : szWork, sizeof(szWork), "%.17g,%.17g,%.17g,%.17g,%.17g,%.17g",
3133 : psInfo->adfSrcInvGeoTransform[0], psInfo->adfSrcInvGeoTransform[1],
3134 : psInfo->adfSrcInvGeoTransform[2], psInfo->adfSrcInvGeoTransform[3],
3135 : psInfo->adfSrcInvGeoTransform[4], psInfo->adfSrcInvGeoTransform[5]);
3136 64 : CPLCreateXMLElementAndValue(psTree, "SrcInvGeoTransform", szWork);
3137 : }
3138 :
3139 : /* -------------------------------------------------------------------- */
3140 : /* Handle dest transformation. */
3141 : /* -------------------------------------------------------------------- */
3142 77 : if (psInfo->pDstTransformArg != nullptr)
3143 : {
3144 0 : CPLXMLNode *psTransformer = GDALSerializeTransformer(
3145 : psInfo->pDstTransformer, psInfo->pDstTransformArg);
3146 0 : if (psTransformer != nullptr)
3147 : {
3148 : CPLXMLNode *psTransformerContainer =
3149 0 : CPLCreateXMLNode(psTree, CXT_Element,
3150 : CPLSPrintf("Dst%s", psTransformer->pszValue));
3151 :
3152 0 : CPLAddXMLChild(psTransformerContainer, psTransformer);
3153 : }
3154 : }
3155 :
3156 : /* -------------------------------------------------------------------- */
3157 : /* Handle destination geotransforms. */
3158 : /* -------------------------------------------------------------------- */
3159 : else
3160 : {
3161 77 : CPLsnprintf(
3162 : szWork, sizeof(szWork), "%.17g,%.17g,%.17g,%.17g,%.17g,%.17g",
3163 : psInfo->adfDstGeoTransform[0], psInfo->adfDstGeoTransform[1],
3164 : psInfo->adfDstGeoTransform[2], psInfo->adfDstGeoTransform[3],
3165 : psInfo->adfDstGeoTransform[4], psInfo->adfDstGeoTransform[5]);
3166 77 : CPLCreateXMLElementAndValue(psTree, "DstGeoTransform", szWork);
3167 :
3168 77 : CPLsnprintf(
3169 : szWork, sizeof(szWork), "%.17g,%.17g,%.17g,%.17g,%.17g,%.17g",
3170 : psInfo->adfDstInvGeoTransform[0], psInfo->adfDstInvGeoTransform[1],
3171 : psInfo->adfDstInvGeoTransform[2], psInfo->adfDstInvGeoTransform[3],
3172 : psInfo->adfDstInvGeoTransform[4], psInfo->adfDstInvGeoTransform[5]);
3173 77 : CPLCreateXMLElementAndValue(psTree, "DstInvGeoTransform", szWork);
3174 : }
3175 :
3176 : /* -------------------------------------------------------------------- */
3177 : /* Do we have a reprojection transformer? */
3178 : /* -------------------------------------------------------------------- */
3179 77 : if (psInfo->pReprojectArg != nullptr)
3180 : {
3181 :
3182 : CPLXMLNode *psTransformerContainer =
3183 50 : CPLCreateXMLNode(psTree, CXT_Element, "ReprojectTransformer");
3184 :
3185 : CPLXMLNode *psTransformer =
3186 50 : GDALSerializeTransformer(psInfo->pReproject, psInfo->pReprojectArg);
3187 50 : if (psTransformer != nullptr)
3188 50 : CPLAddXMLChild(psTransformerContainer, psTransformer);
3189 : }
3190 :
3191 77 : return psTree;
3192 : }
3193 :
3194 : /************************************************************************/
3195 : /* GDALDeserializeGeoTransform() */
3196 : /************************************************************************/
3197 :
3198 735 : static void GDALDeserializeGeoTransform(const char *pszGT,
3199 : double adfGeoTransform[6])
3200 : {
3201 735 : CPLsscanf(pszGT, "%lf,%lf,%lf,%lf,%lf,%lf", adfGeoTransform + 0,
3202 : adfGeoTransform + 1, adfGeoTransform + 2, adfGeoTransform + 3,
3203 : adfGeoTransform + 4, adfGeoTransform + 5);
3204 735 : }
3205 :
3206 : /************************************************************************/
3207 : /* GDALDeserializeGenImgProjTransformer() */
3208 : /************************************************************************/
3209 :
3210 198 : void *GDALDeserializeGenImgProjTransformer(CPLXMLNode *psTree)
3211 :
3212 : {
3213 : /* -------------------------------------------------------------------- */
3214 : /* Initialize the transform info. */
3215 : /* -------------------------------------------------------------------- */
3216 : GDALGenImgProjTransformInfo *psInfo =
3217 198 : GDALCreateGenImgProjTransformerInternal();
3218 :
3219 : /* -------------------------------------------------------------------- */
3220 : /* SrcGeotransform */
3221 : /* -------------------------------------------------------------------- */
3222 198 : if (CPLGetXMLNode(psTree, "SrcGeoTransform") != nullptr)
3223 : {
3224 182 : GDALDeserializeGeoTransform(
3225 : CPLGetXMLValue(psTree, "SrcGeoTransform", ""),
3226 182 : psInfo->adfSrcGeoTransform);
3227 :
3228 182 : if (CPLGetXMLNode(psTree, "SrcInvGeoTransform") != nullptr)
3229 : {
3230 182 : GDALDeserializeGeoTransform(
3231 : CPLGetXMLValue(psTree, "SrcInvGeoTransform", ""),
3232 182 : psInfo->adfSrcInvGeoTransform);
3233 : }
3234 : else
3235 : {
3236 0 : if (!GDALInvGeoTransform(psInfo->adfSrcGeoTransform,
3237 0 : psInfo->adfSrcInvGeoTransform))
3238 : {
3239 0 : CPLError(CE_Failure, CPLE_AppDefined,
3240 : "Cannot invert geotransform");
3241 : }
3242 : }
3243 : }
3244 :
3245 : /* -------------------------------------------------------------------- */
3246 : /* Src Transform */
3247 : /* -------------------------------------------------------------------- */
3248 : else
3249 : {
3250 16 : for (CPLXMLNode *psIter = psTree->psChild; psIter != nullptr;
3251 0 : psIter = psIter->psNext)
3252 : {
3253 16 : if (psIter->eType == CXT_Element &&
3254 16 : STARTS_WITH_CI(psIter->pszValue, "Src"))
3255 : {
3256 16 : GDALDeserializeTransformer(psIter->psChild,
3257 : &psInfo->pSrcTransformer,
3258 : &psInfo->pSrcTransformArg);
3259 16 : break;
3260 : }
3261 : }
3262 : }
3263 :
3264 : /* -------------------------------------------------------------------- */
3265 : /* DstGeotransform */
3266 : /* -------------------------------------------------------------------- */
3267 198 : if (CPLGetXMLNode(psTree, "DstGeoTransform") != nullptr)
3268 : {
3269 198 : GDALDeserializeGeoTransform(
3270 : CPLGetXMLValue(psTree, "DstGeoTransform", ""),
3271 198 : psInfo->adfDstGeoTransform);
3272 :
3273 198 : if (CPLGetXMLNode(psTree, "DstInvGeoTransform") != nullptr)
3274 : {
3275 173 : GDALDeserializeGeoTransform(
3276 : CPLGetXMLValue(psTree, "DstInvGeoTransform", ""),
3277 173 : psInfo->adfDstInvGeoTransform);
3278 : }
3279 : else
3280 : {
3281 25 : if (!GDALInvGeoTransform(psInfo->adfDstGeoTransform,
3282 25 : psInfo->adfDstInvGeoTransform))
3283 : {
3284 0 : CPLError(CE_Failure, CPLE_AppDefined,
3285 : "Cannot invert geotransform");
3286 : }
3287 : }
3288 : }
3289 :
3290 : /* -------------------------------------------------------------------- */
3291 : /* Dst Transform */
3292 : /* -------------------------------------------------------------------- */
3293 : else
3294 : {
3295 0 : for (CPLXMLNode *psIter = psTree->psChild; psIter != nullptr;
3296 0 : psIter = psIter->psNext)
3297 : {
3298 0 : if (psIter->eType == CXT_Element &&
3299 0 : STARTS_WITH_CI(psIter->pszValue, "Dst"))
3300 : {
3301 0 : GDALDeserializeTransformer(psIter->psChild,
3302 : &psInfo->pDstTransformer,
3303 : &psInfo->pDstTransformArg);
3304 0 : break;
3305 : }
3306 : }
3307 : }
3308 :
3309 : /* -------------------------------------------------------------------- */
3310 : /* Reproject transformer */
3311 : /* -------------------------------------------------------------------- */
3312 198 : CPLXMLNode *psSubtree = CPLGetXMLNode(psTree, "ReprojectTransformer");
3313 198 : if (psSubtree != nullptr && psSubtree->psChild != nullptr)
3314 : {
3315 97 : GDALDeserializeTransformer(psSubtree->psChild, &psInfo->pReproject,
3316 : &psInfo->pReprojectArg);
3317 : }
3318 :
3319 198 : return psInfo;
3320 : }
3321 :
3322 : /************************************************************************/
3323 : /* GDALCreateReprojectionTransformer() */
3324 : /************************************************************************/
3325 :
3326 : /**
3327 : * Create reprojection transformer.
3328 : *
3329 : * Creates a callback data structure suitable for use with
3330 : * GDALReprojectionTransformation() to represent a transformation from
3331 : * one geographic or projected coordinate system to another. On input
3332 : * the coordinate systems are described in OpenGIS WKT format.
3333 : *
3334 : * Internally the OGRCoordinateTransformation object is used to implement
3335 : * the reprojection.
3336 : *
3337 : * @param pszSrcWKT the coordinate system for the source coordinate system.
3338 : * @param pszDstWKT the coordinate system for the destination coordinate
3339 : * system.
3340 : *
3341 : * @return Handle for use with GDALReprojectionTransform(), or NULL if the
3342 : * system fails to initialize the reprojection.
3343 : **/
3344 :
3345 0 : void *GDALCreateReprojectionTransformer(const char *pszSrcWKT,
3346 : const char *pszDstWKT)
3347 :
3348 : {
3349 : /* -------------------------------------------------------------------- */
3350 : /* Ingest the SRS definitions. */
3351 : /* -------------------------------------------------------------------- */
3352 0 : OGRSpatialReference oSrcSRS;
3353 0 : oSrcSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3354 0 : if (oSrcSRS.importFromWkt(pszSrcWKT) != OGRERR_NONE)
3355 : {
3356 0 : CPLError(CE_Failure, CPLE_AppDefined,
3357 : "Failed to import coordinate system `%s'.", pszSrcWKT);
3358 0 : return nullptr;
3359 : }
3360 :
3361 0 : OGRSpatialReference oDstSRS;
3362 0 : oDstSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3363 0 : if (oDstSRS.importFromWkt(pszDstWKT) != OGRERR_NONE)
3364 : {
3365 0 : CPLError(CE_Failure, CPLE_AppDefined,
3366 : "Failed to import coordinate system `%s'.", pszSrcWKT);
3367 0 : return nullptr;
3368 : }
3369 :
3370 0 : return GDALCreateReprojectionTransformerEx(
3371 : OGRSpatialReference::ToHandle(&oSrcSRS),
3372 0 : OGRSpatialReference::ToHandle(&oDstSRS), nullptr);
3373 : }
3374 :
3375 : /************************************************************************/
3376 : /* GDALCreateReprojectionTransformerEx() */
3377 : /************************************************************************/
3378 :
3379 : /**
3380 : * Create reprojection transformer.
3381 : *
3382 : * Creates a callback data structure suitable for use with
3383 : * GDALReprojectionTransformation() to represent a transformation from
3384 : * one geographic or projected coordinate system to another.
3385 : *
3386 : * Internally the OGRCoordinateTransformation object is used to implement
3387 : * the reprojection.
3388 : *
3389 : * @param hSrcSRS the coordinate system for the source coordinate system.
3390 : * @param hDstSRS the coordinate system for the destination coordinate
3391 : * system.
3392 : * @param papszOptions NULL-terminated list of options, or NULL. Currently
3393 : * supported options are:
3394 : * <ul>
3395 : * <li>AREA_OF_INTEREST=west_long,south_lat,east_long,north_lat: Values in
3396 : * degrees. longitudes in [-180,180], latitudes in [-90,90].</li>
3397 : * <li>COORDINATE_OPERATION=string: PROJ or WKT string representing a
3398 : * coordinate operation, overriding the default computed transformation.</li>
3399 : * <li>COORDINATE_EPOCH=decimal_year: Coordinate epoch, expressed as a
3400 : * decimal year. Useful for time-dependent coordinate operations.</li>
3401 : * <li> SRC_COORDINATE_EPOCH: (GDAL >= 3.4) Coordinate epoch of source CRS,
3402 : * expressed as a decimal year. Useful for time-dependent coordinate
3403 : *operations.</li> <li> DST_COORDINATE_EPOCH: (GDAL >= 3.4) Coordinate epoch
3404 : *of target CRS, expressed as a decimal year. Useful for time-dependent
3405 : *coordinate operations.</li>
3406 : * </ul>
3407 : *
3408 : * @return Handle for use with GDALReprojectionTransform(), or NULL if the
3409 : * system fails to initialize the reprojection.
3410 : *
3411 : * @since GDAL 3.0
3412 : **/
3413 :
3414 910 : void *GDALCreateReprojectionTransformerEx(OGRSpatialReferenceH hSrcSRS,
3415 : OGRSpatialReferenceH hDstSRS,
3416 : const char *const *papszOptions)
3417 : {
3418 910 : OGRSpatialReference *poSrcSRS = OGRSpatialReference::FromHandle(hSrcSRS);
3419 910 : OGRSpatialReference *poDstSRS = OGRSpatialReference::FromHandle(hDstSRS);
3420 :
3421 : /* -------------------------------------------------------------------- */
3422 : /* Build the forward coordinate transformation. */
3423 : /* -------------------------------------------------------------------- */
3424 910 : double dfWestLongitudeDeg = 0.0;
3425 910 : double dfSouthLatitudeDeg = 0.0;
3426 910 : double dfEastLongitudeDeg = 0.0;
3427 910 : double dfNorthLatitudeDeg = 0.0;
3428 910 : const char *pszBBOX = CSLFetchNameValue(papszOptions, "AREA_OF_INTEREST");
3429 910 : if (pszBBOX)
3430 : {
3431 844 : char **papszTokens = CSLTokenizeString2(pszBBOX, ",", 0);
3432 844 : if (CSLCount(papszTokens) == 4)
3433 : {
3434 844 : dfWestLongitudeDeg = CPLAtof(papszTokens[0]);
3435 844 : dfSouthLatitudeDeg = CPLAtof(papszTokens[1]);
3436 844 : dfEastLongitudeDeg = CPLAtof(papszTokens[2]);
3437 844 : dfNorthLatitudeDeg = CPLAtof(papszTokens[3]);
3438 : }
3439 844 : CSLDestroy(papszTokens);
3440 : }
3441 910 : const char *pszCO = CSLFetchNameValue(papszOptions, "COORDINATE_OPERATION");
3442 :
3443 1820 : OGRCoordinateTransformationOptions optionsFwd;
3444 910 : if (!(dfWestLongitudeDeg == 0.0 && dfSouthLatitudeDeg == 0.0 &&
3445 : dfEastLongitudeDeg == 0.0 && dfNorthLatitudeDeg == 0.0))
3446 : {
3447 844 : optionsFwd.SetAreaOfInterest(dfWestLongitudeDeg, dfSouthLatitudeDeg,
3448 : dfEastLongitudeDeg, dfNorthLatitudeDeg);
3449 : }
3450 910 : if (pszCO)
3451 : {
3452 7 : optionsFwd.SetCoordinateOperation(pszCO, false);
3453 : }
3454 :
3455 910 : const char *pszCENTER_LONG = CSLFetchNameValue(papszOptions, "CENTER_LONG");
3456 910 : if (pszCENTER_LONG)
3457 : {
3458 605 : optionsFwd.SetSourceCenterLong(CPLAtof(pszCENTER_LONG));
3459 : }
3460 :
3461 : OGRCoordinateTransformation *poForwardTransform =
3462 910 : OGRCreateCoordinateTransformation(poSrcSRS, poDstSRS, optionsFwd);
3463 :
3464 910 : if (poForwardTransform == nullptr)
3465 : // OGRCreateCoordinateTransformation() will report errors on its own.
3466 2 : return nullptr;
3467 :
3468 908 : poForwardTransform->SetEmitErrors(false);
3469 :
3470 : /* -------------------------------------------------------------------- */
3471 : /* Create a structure to hold the transform info, and also */
3472 : /* build reverse transform. We assume that if the forward */
3473 : /* transform can be created, then so can the reverse one. */
3474 : /* -------------------------------------------------------------------- */
3475 908 : GDALReprojectionTransformInfo *psInfo = new GDALReprojectionTransformInfo();
3476 :
3477 908 : psInfo->papszOptions = CSLDuplicate(papszOptions);
3478 908 : psInfo->poForwardTransform = poForwardTransform;
3479 908 : psInfo->dfTime = CPLAtof(CSLFetchNameValueDef(
3480 : papszOptions, "COORDINATE_EPOCH",
3481 : CSLFetchNameValueDef(
3482 : papszOptions, "DST_COORDINATE_EPOCH",
3483 : CSLFetchNameValueDef(papszOptions, "SRC_COORDINATE_EPOCH", "0"))));
3484 908 : psInfo->poReverseTransform = poForwardTransform->GetInverse();
3485 :
3486 908 : if (psInfo->poReverseTransform)
3487 908 : psInfo->poReverseTransform->SetEmitErrors(false);
3488 :
3489 908 : memcpy(psInfo->sTI.abySignature, GDAL_GTI2_SIGNATURE,
3490 : strlen(GDAL_GTI2_SIGNATURE));
3491 908 : psInfo->sTI.pszClassName = "GDALReprojectionTransformer";
3492 908 : psInfo->sTI.pfnTransform = GDALReprojectionTransform;
3493 908 : psInfo->sTI.pfnCleanup = GDALDestroyReprojectionTransformer;
3494 908 : psInfo->sTI.pfnSerialize = GDALSerializeReprojectionTransformer;
3495 :
3496 908 : return psInfo;
3497 : }
3498 :
3499 : /************************************************************************/
3500 : /* GDALDestroyReprojectionTransformer() */
3501 : /************************************************************************/
3502 :
3503 : /**
3504 : * Destroy reprojection transformation.
3505 : *
3506 : * @param pTransformArg the transformation handle returned by
3507 : * GDALCreateReprojectionTransformer().
3508 : */
3509 :
3510 908 : void GDALDestroyReprojectionTransformer(void *pTransformArg)
3511 :
3512 : {
3513 908 : if (pTransformArg == nullptr)
3514 0 : return;
3515 :
3516 908 : GDALReprojectionTransformInfo *psInfo =
3517 : static_cast<GDALReprojectionTransformInfo *>(pTransformArg);
3518 :
3519 908 : if (psInfo->poForwardTransform)
3520 908 : OGRCoordinateTransformation::DestroyCT(psInfo->poForwardTransform);
3521 :
3522 908 : if (psInfo->poReverseTransform)
3523 908 : OGRCoordinateTransformation::DestroyCT(psInfo->poReverseTransform);
3524 :
3525 908 : CSLDestroy(psInfo->papszOptions);
3526 :
3527 908 : delete psInfo;
3528 : }
3529 :
3530 : /************************************************************************/
3531 : /* GDALReprojectionTransform() */
3532 : /************************************************************************/
3533 :
3534 : /**
3535 : * Perform reprojection transformation.
3536 : *
3537 : * Actually performs the reprojection transformation described in
3538 : * GDALCreateReprojectionTransformer(). This function matches the
3539 : * GDALTransformerFunc() signature. Details of the arguments are described
3540 : * there.
3541 : */
3542 :
3543 1609480 : int GDALReprojectionTransform(void *pTransformArg, int bDstToSrc,
3544 : int nPointCount, double *padfX, double *padfY,
3545 : double *padfZ, int *panSuccess)
3546 :
3547 : {
3548 1609480 : GDALReprojectionTransformInfo *psInfo =
3549 : static_cast<GDALReprojectionTransformInfo *>(pTransformArg);
3550 : int bSuccess;
3551 :
3552 1609480 : std::vector<double> adfTime;
3553 1609470 : double *padfT = nullptr;
3554 1609470 : if (psInfo->dfTime != 0.0 && nPointCount > 0)
3555 : {
3556 1 : adfTime.resize(nPointCount, psInfo->dfTime);
3557 1 : padfT = &adfTime[0];
3558 : }
3559 :
3560 1609470 : if (bDstToSrc)
3561 : {
3562 1371370 : if (psInfo->poReverseTransform == nullptr)
3563 : {
3564 0 : CPLError(
3565 : CE_Failure, CPLE_AppDefined,
3566 : "Inverse coordinate transformation cannot be instantiated");
3567 0 : if (panSuccess)
3568 : {
3569 0 : for (int i = 0; i < nPointCount; i++)
3570 0 : panSuccess[i] = FALSE;
3571 : }
3572 0 : bSuccess = false;
3573 : }
3574 : else
3575 : {
3576 1371370 : bSuccess = psInfo->poReverseTransform->Transform(
3577 1371370 : nPointCount, padfX, padfY, padfZ, padfT, panSuccess);
3578 : }
3579 : }
3580 : else
3581 238099 : bSuccess = psInfo->poForwardTransform->Transform(
3582 238099 : nPointCount, padfX, padfY, padfZ, padfT, panSuccess);
3583 :
3584 3218920 : return bSuccess;
3585 : }
3586 :
3587 : /************************************************************************/
3588 : /* GDALSerializeReprojectionTransformer() */
3589 : /************************************************************************/
3590 :
3591 127 : static CPLXMLNode *GDALSerializeReprojectionTransformer(void *pTransformArg)
3592 :
3593 : {
3594 : CPLXMLNode *psTree;
3595 127 : GDALReprojectionTransformInfo *psInfo =
3596 : static_cast<GDALReprojectionTransformInfo *>(pTransformArg);
3597 :
3598 127 : psTree = CPLCreateXMLNode(nullptr, CXT_Element, "ReprojectionTransformer");
3599 :
3600 : /* -------------------------------------------------------------------- */
3601 : /* Handle SourceCS. */
3602 : /* -------------------------------------------------------------------- */
3603 254 : const auto ExportToWkt = [](const OGRSpatialReference *poSRS)
3604 : {
3605 : // Try first in WKT1 for backward compat
3606 : {
3607 254 : char *pszWKT = nullptr;
3608 254 : const char *const apszOptions[] = {"FORMAT=WKT1", nullptr};
3609 254 : CPLErrorHandlerPusher oHandler(CPLQuietErrorHandler);
3610 254 : CPLErrorStateBackuper oBackuper;
3611 254 : if (poSRS->exportToWkt(&pszWKT, apszOptions) == OGRERR_NONE)
3612 : {
3613 506 : std::string osRet(pszWKT);
3614 253 : CPLFree(pszWKT);
3615 253 : return osRet;
3616 : }
3617 1 : CPLFree(pszWKT);
3618 : }
3619 :
3620 1 : char *pszWKT = nullptr;
3621 1 : const char *const apszOptions[] = {"FORMAT=WKT2_2019", nullptr};
3622 1 : if (poSRS->exportToWkt(&pszWKT, apszOptions) == OGRERR_NONE)
3623 : {
3624 2 : std::string osRet(pszWKT);
3625 1 : CPLFree(pszWKT);
3626 1 : return osRet;
3627 : }
3628 0 : CPLFree(pszWKT);
3629 0 : return std::string();
3630 : };
3631 :
3632 127 : auto poSRS = psInfo->poForwardTransform->GetSourceCS();
3633 127 : if (poSRS)
3634 : {
3635 254 : const auto osWKT = ExportToWkt(poSRS);
3636 127 : CPLCreateXMLElementAndValue(psTree, "SourceSRS", osWKT.c_str());
3637 : }
3638 :
3639 : /* -------------------------------------------------------------------- */
3640 : /* Handle DestinationCS. */
3641 : /* -------------------------------------------------------------------- */
3642 127 : poSRS = psInfo->poForwardTransform->GetTargetCS();
3643 127 : if (poSRS)
3644 : {
3645 254 : const auto osWKT = ExportToWkt(poSRS);
3646 127 : CPLCreateXMLElementAndValue(psTree, "TargetSRS", osWKT.c_str());
3647 : }
3648 :
3649 : /* -------------------------------------------------------------------- */
3650 : /* Serialize options. */
3651 : /* -------------------------------------------------------------------- */
3652 127 : if (psInfo->papszOptions)
3653 : {
3654 : CPLXMLNode *psOptions =
3655 114 : CPLCreateXMLNode(psTree, CXT_Element, "Options");
3656 276 : for (auto iter = psInfo->papszOptions; *iter != nullptr; ++iter)
3657 : {
3658 162 : char *pszKey = nullptr;
3659 162 : const char *pszValue = CPLParseNameValue(*iter, &pszKey);
3660 162 : if (pszKey && pszValue)
3661 : {
3662 : auto elt =
3663 162 : CPLCreateXMLElementAndValue(psOptions, "Option", pszValue);
3664 162 : CPLAddXMLAttributeAndValue(elt, "key", pszKey);
3665 : }
3666 162 : CPLFree(pszKey);
3667 : }
3668 : }
3669 :
3670 127 : return psTree;
3671 : }
3672 :
3673 : /************************************************************************/
3674 : /* GDALDeserializeReprojectionTransformer() */
3675 : /************************************************************************/
3676 :
3677 174 : static void *GDALDeserializeReprojectionTransformer(CPLXMLNode *psTree)
3678 :
3679 : {
3680 174 : const char *pszSourceSRS = CPLGetXMLValue(psTree, "SourceSRS", nullptr);
3681 174 : const char *pszTargetSRS = CPLGetXMLValue(psTree, "TargetSRS", nullptr);
3682 174 : char *pszSourceWKT = nullptr, *pszTargetWKT = nullptr;
3683 174 : void *pResult = nullptr;
3684 :
3685 348 : OGRSpatialReference oSrcSRS;
3686 348 : OGRSpatialReference oDstSRS;
3687 :
3688 174 : oSrcSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3689 174 : oDstSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3690 174 : if (pszSourceSRS != nullptr)
3691 : {
3692 174 : oSrcSRS.SetFromUserInput(pszSourceSRS);
3693 : }
3694 :
3695 174 : if (pszTargetSRS != nullptr)
3696 : {
3697 174 : oDstSRS.SetFromUserInput(pszTargetSRS);
3698 : }
3699 :
3700 174 : CPLStringList aosList;
3701 174 : const CPLXMLNode *psOptions = CPLGetXMLNode(psTree, "Options");
3702 174 : if (psOptions)
3703 : {
3704 344 : for (auto iter = psOptions->psChild; iter; iter = iter->psNext)
3705 : {
3706 206 : if (iter->eType == CXT_Element &&
3707 206 : strcmp(iter->pszValue, "Option") == 0)
3708 : {
3709 206 : const char *pszKey = CPLGetXMLValue(iter, "key", nullptr);
3710 206 : const char *pszValue = CPLGetXMLValue(iter, nullptr, nullptr);
3711 206 : if (pszKey && pszValue)
3712 : {
3713 206 : aosList.SetNameValue(pszKey, pszValue);
3714 : }
3715 : }
3716 : }
3717 : }
3718 :
3719 174 : pResult = GDALCreateReprojectionTransformerEx(
3720 174 : !oSrcSRS.IsEmpty() ? OGRSpatialReference::ToHandle(&oSrcSRS) : nullptr,
3721 174 : !oDstSRS.IsEmpty() ? OGRSpatialReference::ToHandle(&oDstSRS) : nullptr,
3722 174 : aosList.List());
3723 :
3724 174 : CPLFree(pszSourceWKT);
3725 174 : CPLFree(pszTargetWKT);
3726 :
3727 348 : return pResult;
3728 : }
3729 :
3730 : /************************************************************************/
3731 : /* ==================================================================== */
3732 : /* Approximate transformer. */
3733 : /* ==================================================================== */
3734 : /************************************************************************/
3735 :
3736 : typedef struct
3737 : {
3738 : GDALTransformerInfo sTI;
3739 :
3740 : GDALTransformerFunc pfnBaseTransformer;
3741 : void *pBaseCBData;
3742 : double dfMaxErrorForward;
3743 : double dfMaxErrorReverse;
3744 :
3745 : int bOwnSubtransformer;
3746 : } ApproxTransformInfo;
3747 :
3748 : /************************************************************************/
3749 : /* GDALCreateSimilarApproxTransformer() */
3750 : /************************************************************************/
3751 :
3752 17 : static void *GDALCreateSimilarApproxTransformer(void *hTransformArg,
3753 : double dfSrcRatioX,
3754 : double dfSrcRatioY)
3755 : {
3756 17 : VALIDATE_POINTER1(hTransformArg, "GDALCreateSimilarApproxTransformer",
3757 : nullptr);
3758 :
3759 17 : ApproxTransformInfo *psInfo =
3760 : static_cast<ApproxTransformInfo *>(hTransformArg);
3761 :
3762 : ApproxTransformInfo *psClonedInfo = static_cast<ApproxTransformInfo *>(
3763 17 : CPLMalloc(sizeof(ApproxTransformInfo)));
3764 :
3765 17 : memcpy(psClonedInfo, psInfo, sizeof(ApproxTransformInfo));
3766 17 : if (psClonedInfo->pBaseCBData)
3767 : {
3768 17 : psClonedInfo->pBaseCBData = GDALCreateSimilarTransformer(
3769 : psInfo->pBaseCBData, dfSrcRatioX, dfSrcRatioY);
3770 17 : if (psClonedInfo->pBaseCBData == nullptr)
3771 : {
3772 0 : CPLFree(psClonedInfo);
3773 0 : return nullptr;
3774 : }
3775 : }
3776 17 : psClonedInfo->bOwnSubtransformer = TRUE;
3777 :
3778 17 : return psClonedInfo;
3779 : }
3780 :
3781 : /************************************************************************/
3782 : /* GDALSerializeApproxTransformer() */
3783 : /************************************************************************/
3784 :
3785 57 : static CPLXMLNode *GDALSerializeApproxTransformer(void *pTransformArg)
3786 :
3787 : {
3788 : CPLXMLNode *psTree;
3789 57 : ApproxTransformInfo *psInfo =
3790 : static_cast<ApproxTransformInfo *>(pTransformArg);
3791 :
3792 57 : psTree = CPLCreateXMLNode(nullptr, CXT_Element, "ApproxTransformer");
3793 :
3794 : /* -------------------------------------------------------------------- */
3795 : /* Attach max error. */
3796 : /* -------------------------------------------------------------------- */
3797 57 : if (psInfo->dfMaxErrorForward == psInfo->dfMaxErrorReverse)
3798 : {
3799 55 : CPLCreateXMLElementAndValue(
3800 : psTree, "MaxError",
3801 110 : CPLString().Printf("%g", psInfo->dfMaxErrorForward));
3802 : }
3803 : else
3804 : {
3805 2 : CPLCreateXMLElementAndValue(
3806 : psTree, "MaxErrorForward",
3807 4 : CPLString().Printf("%g", psInfo->dfMaxErrorForward));
3808 2 : CPLCreateXMLElementAndValue(
3809 : psTree, "MaxErrorReverse",
3810 4 : CPLString().Printf("%g", psInfo->dfMaxErrorReverse));
3811 : }
3812 :
3813 : /* -------------------------------------------------------------------- */
3814 : /* Capture underlying transformer. */
3815 : /* -------------------------------------------------------------------- */
3816 : CPLXMLNode *psTransformerContainer =
3817 57 : CPLCreateXMLNode(psTree, CXT_Element, "BaseTransformer");
3818 :
3819 57 : CPLXMLNode *psTransformer = GDALSerializeTransformer(
3820 : psInfo->pfnBaseTransformer, psInfo->pBaseCBData);
3821 57 : if (psTransformer != nullptr)
3822 57 : CPLAddXMLChild(psTransformerContainer, psTransformer);
3823 :
3824 57 : return psTree;
3825 : }
3826 :
3827 : /************************************************************************/
3828 : /* GDALCreateApproxTransformer() */
3829 : /************************************************************************/
3830 :
3831 : /**
3832 : * Create an approximating transformer.
3833 : *
3834 : * This function creates a context for an approximated transformer. Basically
3835 : * a high precision transformer is supplied as input and internally linear
3836 : * approximations are computed to generate results to within a defined
3837 : * precision.
3838 : *
3839 : * The approximation is actually done at the point where GDALApproxTransform()
3840 : * calls are made, and depend on the assumption that they are roughly linear.
3841 : * The first and last point passed in must be the extreme values and the
3842 : * intermediate values should describe a curve between the end points. The
3843 : * approximator transforms and centers using the approximate transformer, and
3844 : * then compares the true middle transformed value to a linear approximation
3845 : * based on the end points. If the error is within the supplied threshold then
3846 : * the end points are used to linearly approximate all the values otherwise the
3847 : * input points are split into two smaller sets, and the function is recursively
3848 : * called until a sufficiently small set of points is found that the linear
3849 : * approximation is OK, or that all the points are exactly computed.
3850 : *
3851 : * This function is very suitable for approximating transformation results
3852 : * from output pixel/line space to input coordinates for warpers that operate
3853 : * on one input scanline at a time. Care should be taken using it in other
3854 : * circumstances as little internal validation is done in order to keep things
3855 : * fast.
3856 : *
3857 : * @param pfnBaseTransformer the high precision transformer which should be
3858 : * approximated.
3859 : * @param pBaseTransformArg the callback argument for the high precision
3860 : * transformer.
3861 : * @param dfMaxError the maximum cartesian error in the "output" space that
3862 : * is to be accepted in the linear approximation, evaluated as a Manhattan
3863 : * distance.
3864 : *
3865 : * @return callback pointer suitable for use with GDALApproxTransform(). It
3866 : * should be deallocated with GDALDestroyApproxTransformer().
3867 : */
3868 :
3869 894 : void *GDALCreateApproxTransformer(GDALTransformerFunc pfnBaseTransformer,
3870 : void *pBaseTransformArg, double dfMaxError)
3871 :
3872 : {
3873 894 : return GDALCreateApproxTransformer2(pfnBaseTransformer, pBaseTransformArg,
3874 894 : dfMaxError, dfMaxError);
3875 : }
3876 :
3877 : static void *
3878 1046 : GDALCreateApproxTransformer2(GDALTransformerFunc pfnBaseTransformer,
3879 : void *pBaseTransformArg, double dfMaxErrorForward,
3880 : double dfMaxErrorReverse)
3881 :
3882 : {
3883 : ApproxTransformInfo *psATInfo = static_cast<ApproxTransformInfo *>(
3884 1046 : CPLMalloc(sizeof(ApproxTransformInfo)));
3885 1046 : psATInfo->pfnBaseTransformer = pfnBaseTransformer;
3886 1046 : psATInfo->pBaseCBData = pBaseTransformArg;
3887 1046 : psATInfo->dfMaxErrorForward = dfMaxErrorForward;
3888 1046 : psATInfo->dfMaxErrorReverse = dfMaxErrorReverse;
3889 1046 : psATInfo->bOwnSubtransformer = FALSE;
3890 :
3891 1046 : memcpy(psATInfo->sTI.abySignature, GDAL_GTI2_SIGNATURE,
3892 : strlen(GDAL_GTI2_SIGNATURE));
3893 1046 : psATInfo->sTI.pszClassName = GDAL_APPROX_TRANSFORMER_CLASS_NAME;
3894 1046 : psATInfo->sTI.pfnTransform = GDALApproxTransform;
3895 1046 : psATInfo->sTI.pfnCleanup = GDALDestroyApproxTransformer;
3896 1046 : psATInfo->sTI.pfnSerialize = GDALSerializeApproxTransformer;
3897 1046 : psATInfo->sTI.pfnCreateSimilar = GDALCreateSimilarApproxTransformer;
3898 :
3899 1046 : return psATInfo;
3900 : }
3901 :
3902 : /************************************************************************/
3903 : /* GDALApproxTransformerOwnsSubtransformer() */
3904 : /************************************************************************/
3905 :
3906 : /** Set bOwnSubtransformer flag */
3907 1043 : void GDALApproxTransformerOwnsSubtransformer(void *pCBData, int bOwnFlag)
3908 :
3909 : {
3910 1043 : ApproxTransformInfo *psATInfo = static_cast<ApproxTransformInfo *>(pCBData);
3911 :
3912 1043 : psATInfo->bOwnSubtransformer = bOwnFlag;
3913 1043 : }
3914 :
3915 : /************************************************************************/
3916 : /* GDALDestroyApproxTransformer() */
3917 : /************************************************************************/
3918 :
3919 : /**
3920 : * Cleanup approximate transformer.
3921 : *
3922 : * Deallocates the resources allocated by GDALCreateApproxTransformer().
3923 : *
3924 : * @param pCBData callback data originally returned by
3925 : * GDALCreateApproxTransformer().
3926 : */
3927 :
3928 1063 : void GDALDestroyApproxTransformer(void *pCBData)
3929 :
3930 : {
3931 1063 : if (pCBData == nullptr)
3932 0 : return;
3933 :
3934 1063 : ApproxTransformInfo *psATInfo = static_cast<ApproxTransformInfo *>(pCBData);
3935 :
3936 1063 : if (psATInfo->bOwnSubtransformer)
3937 1060 : GDALDestroyTransformer(psATInfo->pBaseCBData);
3938 :
3939 1063 : CPLFree(pCBData);
3940 : }
3941 :
3942 : /************************************************************************/
3943 : /* GDALRefreshApproxTransformer() */
3944 : /************************************************************************/
3945 :
3946 44 : void GDALRefreshApproxTransformer(void *hTransformArg)
3947 : {
3948 44 : ApproxTransformInfo *psInfo =
3949 : static_cast<ApproxTransformInfo *>(hTransformArg);
3950 :
3951 44 : if (GDALIsTransformer(psInfo->pBaseCBData,
3952 : GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME))
3953 : {
3954 44 : GDALRefreshGenImgProjTransformer(psInfo->pBaseCBData);
3955 : }
3956 44 : }
3957 :
3958 : /************************************************************************/
3959 : /* GDALApproxTransformInternal() */
3960 : /************************************************************************/
3961 :
3962 1111940 : static int GDALApproxTransformInternal(void *pCBData, int bDstToSrc,
3963 : int nPoints, double *x, double *y,
3964 : double *z, int *panSuccess,
3965 : // SME = Start, Middle, End.
3966 : const double xSMETransformed[3],
3967 : const double ySMETransformed[3],
3968 : const double zSMETransformed[3])
3969 : {
3970 1111940 : ApproxTransformInfo *psATInfo = static_cast<ApproxTransformInfo *>(pCBData);
3971 1111940 : const int nMiddle = (nPoints - 1) / 2;
3972 :
3973 : #ifdef notdef_sanify_check
3974 : {
3975 : double x2[3] = {x[0], x[nMiddle], x[nPoints - 1]};
3976 : double y2[3] = {y[0], y[nMiddle], y[nPoints - 1]};
3977 : double z2[3] = {z[0], z[nMiddle], z[nPoints - 1]};
3978 : int anSuccess2[3] = {};
3979 :
3980 : const int bSuccess = psATInfo->pfnBaseTransformer(
3981 : psATInfo->pBaseCBData, bDstToSrc, 3, x2, y2, z2, anSuccess2);
3982 : CPLAssert(bSuccess);
3983 : CPLAssert(anSuccess2[0]);
3984 : CPLAssert(anSuccess2[1]);
3985 : CPLAssert(anSuccess2[2]);
3986 : CPLAssert(x2[0] == xSMETransformed[0]);
3987 : CPLAssert(y2[0] == ySMETransformed[0]);
3988 : CPLAssert(z2[0] == zSMETransformed[0]);
3989 : CPLAssert(x2[1] == xSMETransformed[1]);
3990 : CPLAssert(y2[1] == ySMETransformed[1]);
3991 : CPLAssert(z2[1] == zSMETransformed[1]);
3992 : CPLAssert(x2[2] == xSMETransformed[2]);
3993 : CPLAssert(y2[2] == ySMETransformed[2]);
3994 : CPLAssert(z2[2] == zSMETransformed[2]);
3995 : }
3996 : #endif
3997 :
3998 : #ifdef DEBUG_APPROX_TRANSFORMER
3999 : fprintf(stderr, "start (%.3f,%.3f) -> (%.3f,%.3f)\n", /*ok*/
4000 : x[0], y[0], xSMETransformed[0], ySMETransformed[0]);
4001 : fprintf(stderr, "middle (%.3f,%.3f) -> (%.3f,%.3f)\n", /*ok*/
4002 : x[nMiddle], y[nMiddle], xSMETransformed[1], ySMETransformed[1]);
4003 : fprintf(stderr, "end (%.3f,%.3f) -> (%.3f,%.3f)\n", /*ok*/
4004 : x[nPoints - 1], y[nPoints - 1], xSMETransformed[2],
4005 : ySMETransformed[2]);
4006 : #endif
4007 :
4008 : /* -------------------------------------------------------------------- */
4009 : /* Is the error at the middle acceptable relative to an */
4010 : /* interpolation of the middle position? */
4011 : /* -------------------------------------------------------------------- */
4012 1111940 : const double dfDeltaX =
4013 1111940 : (xSMETransformed[2] - xSMETransformed[0]) / (x[nPoints - 1] - x[0]);
4014 1111940 : const double dfDeltaY =
4015 1111940 : (ySMETransformed[2] - ySMETransformed[0]) / (x[nPoints - 1] - x[0]);
4016 1111940 : const double dfDeltaZ =
4017 1111940 : (zSMETransformed[2] - zSMETransformed[0]) / (x[nPoints - 1] - x[0]);
4018 :
4019 1111940 : const double dfError =
4020 1111940 : fabs((xSMETransformed[0] + dfDeltaX * (x[nMiddle] - x[0])) -
4021 1111940 : xSMETransformed[1]) +
4022 1111940 : fabs((ySMETransformed[0] + dfDeltaY * (x[nMiddle] - x[0])) -
4023 1111940 : ySMETransformed[1]);
4024 :
4025 1111940 : const double dfMaxError =
4026 1111940 : (bDstToSrc) ? psATInfo->dfMaxErrorReverse : psATInfo->dfMaxErrorForward;
4027 1111940 : if (dfError > dfMaxError)
4028 : {
4029 : #if DEBUG_VERBOSE
4030 : CPLDebug("GDAL",
4031 : "ApproxTransformer - "
4032 : "error %g over threshold %g, subdivide %d points.",
4033 : dfError, dfMaxError, nPoints);
4034 : #endif
4035 :
4036 590163 : double xMiddle[3] = {x[(nMiddle - 1) / 2], x[nMiddle - 1],
4037 590163 : x[nMiddle + (nPoints - nMiddle - 1) / 2]};
4038 590163 : double yMiddle[3] = {y[(nMiddle - 1) / 2], y[nMiddle - 1],
4039 590163 : y[nMiddle + (nPoints - nMiddle - 1) / 2]};
4040 590163 : double zMiddle[3] = {z[(nMiddle - 1) / 2], z[nMiddle - 1],
4041 590163 : z[nMiddle + (nPoints - nMiddle - 1) / 2]};
4042 :
4043 590163 : const bool bUseBaseTransformForHalf1 =
4044 462280 : nMiddle <= 5 || y[0] != y[nMiddle - 1] ||
4045 1514720 : y[0] != y[(nMiddle - 1) / 2] || x[0] == x[nMiddle - 1] ||
4046 462280 : x[0] == x[(nMiddle - 1) / 2];
4047 590163 : const bool bUseBaseTransformForHalf2 =
4048 476268 : nPoints - nMiddle <= 5 || y[nMiddle] != y[nPoints - 1] ||
4049 476268 : y[nMiddle] != y[nMiddle + (nPoints - nMiddle - 1) / 2] ||
4050 1542700 : x[nMiddle] == x[nPoints - 1] ||
4051 476268 : x[nMiddle] == x[nMiddle + (nPoints - nMiddle - 1) / 2];
4052 :
4053 590163 : int anSuccess2[3] = {};
4054 590163 : int bSuccess = FALSE;
4055 590163 : if (!bUseBaseTransformForHalf1 && !bUseBaseTransformForHalf2)
4056 462280 : bSuccess = psATInfo->pfnBaseTransformer(
4057 : psATInfo->pBaseCBData, bDstToSrc, 3, xMiddle, yMiddle, zMiddle,
4058 : anSuccess2);
4059 127883 : else if (!bUseBaseTransformForHalf1)
4060 : {
4061 0 : bSuccess = psATInfo->pfnBaseTransformer(
4062 : psATInfo->pBaseCBData, bDstToSrc, 2, xMiddle, yMiddle, zMiddle,
4063 : anSuccess2);
4064 0 : anSuccess2[2] = TRUE;
4065 : }
4066 127883 : else if (!bUseBaseTransformForHalf2)
4067 : {
4068 13988 : bSuccess = psATInfo->pfnBaseTransformer(
4069 : psATInfo->pBaseCBData, bDstToSrc, 1, xMiddle + 2, yMiddle + 2,
4070 : zMiddle + 2, anSuccess2 + 2);
4071 13988 : anSuccess2[0] = TRUE;
4072 13988 : anSuccess2[1] = TRUE;
4073 : }
4074 :
4075 590163 : if (!bSuccess || !anSuccess2[0] || !anSuccess2[1] || !anSuccess2[2])
4076 : {
4077 113911 : bSuccess = psATInfo->pfnBaseTransformer(
4078 : psATInfo->pBaseCBData, bDstToSrc, nMiddle - 1, x + 1, y + 1,
4079 : z + 1, panSuccess + 1);
4080 227822 : bSuccess &= psATInfo->pfnBaseTransformer(
4081 113911 : psATInfo->pBaseCBData, bDstToSrc, nPoints - nMiddle - 2,
4082 113911 : x + nMiddle + 1, y + nMiddle + 1, z + nMiddle + 1,
4083 113911 : panSuccess + nMiddle + 1);
4084 :
4085 113911 : x[0] = xSMETransformed[0];
4086 113911 : y[0] = ySMETransformed[0];
4087 113911 : z[0] = zSMETransformed[0];
4088 113911 : panSuccess[0] = TRUE;
4089 113911 : x[nMiddle] = xSMETransformed[1];
4090 113911 : y[nMiddle] = ySMETransformed[1];
4091 113911 : z[nMiddle] = zSMETransformed[1];
4092 113911 : panSuccess[nMiddle] = TRUE;
4093 113911 : x[nPoints - 1] = xSMETransformed[2];
4094 113911 : y[nPoints - 1] = ySMETransformed[2];
4095 113911 : z[nPoints - 1] = zSMETransformed[2];
4096 113911 : panSuccess[nPoints - 1] = TRUE;
4097 113911 : return bSuccess;
4098 : }
4099 :
4100 476252 : double x2[3] = {};
4101 476252 : double y2[3] = {};
4102 476252 : double z2[3] = {};
4103 476252 : if (!bUseBaseTransformForHalf1)
4104 : {
4105 462264 : x2[0] = xSMETransformed[0];
4106 462264 : y2[0] = ySMETransformed[0];
4107 462264 : z2[0] = zSMETransformed[0];
4108 462264 : x2[1] = xMiddle[0];
4109 462264 : y2[1] = yMiddle[0];
4110 462264 : z2[1] = zMiddle[0];
4111 462264 : x2[2] = xMiddle[1];
4112 462264 : y2[2] = yMiddle[1];
4113 462264 : z2[2] = zMiddle[1];
4114 :
4115 462264 : bSuccess = GDALApproxTransformInternal(
4116 : psATInfo, bDstToSrc, nMiddle, x, y, z, panSuccess, x2, y2, z2);
4117 : }
4118 : else
4119 : {
4120 13988 : bSuccess = psATInfo->pfnBaseTransformer(
4121 : psATInfo->pBaseCBData, bDstToSrc, nMiddle - 1, x + 1, y + 1,
4122 : z + 1, panSuccess + 1);
4123 13988 : x[0] = xSMETransformed[0];
4124 13988 : y[0] = ySMETransformed[0];
4125 13988 : z[0] = zSMETransformed[0];
4126 13988 : panSuccess[0] = TRUE;
4127 : }
4128 :
4129 476252 : if (!bSuccess)
4130 0 : return FALSE;
4131 :
4132 476252 : if (!bUseBaseTransformForHalf2)
4133 : {
4134 476252 : x2[0] = xSMETransformed[1];
4135 476252 : y2[0] = ySMETransformed[1];
4136 476252 : z2[0] = zSMETransformed[1];
4137 476252 : x2[1] = xMiddle[2];
4138 476252 : y2[1] = yMiddle[2];
4139 476252 : z2[1] = zMiddle[2];
4140 476252 : x2[2] = xSMETransformed[2];
4141 476252 : y2[2] = ySMETransformed[2];
4142 476252 : z2[2] = zSMETransformed[2];
4143 :
4144 476252 : bSuccess = GDALApproxTransformInternal(
4145 476252 : psATInfo, bDstToSrc, nPoints - nMiddle, x + nMiddle,
4146 476252 : y + nMiddle, z + nMiddle, panSuccess + nMiddle, x2, y2, z2);
4147 : }
4148 : else
4149 : {
4150 0 : bSuccess = psATInfo->pfnBaseTransformer(
4151 0 : psATInfo->pBaseCBData, bDstToSrc, nPoints - nMiddle - 2,
4152 0 : x + nMiddle + 1, y + nMiddle + 1, z + nMiddle + 1,
4153 0 : panSuccess + nMiddle + 1);
4154 :
4155 0 : x[nMiddle] = xSMETransformed[1];
4156 0 : y[nMiddle] = ySMETransformed[1];
4157 0 : z[nMiddle] = zSMETransformed[1];
4158 0 : panSuccess[nMiddle] = TRUE;
4159 0 : x[nPoints - 1] = xSMETransformed[2];
4160 0 : y[nPoints - 1] = ySMETransformed[2];
4161 0 : z[nPoints - 1] = zSMETransformed[2];
4162 0 : panSuccess[nPoints - 1] = TRUE;
4163 : }
4164 :
4165 476252 : if (!bSuccess)
4166 0 : return FALSE;
4167 :
4168 476252 : return TRUE;
4169 : }
4170 :
4171 : /* -------------------------------------------------------------------- */
4172 : /* Error is OK since this is just used to compute output bounds */
4173 : /* of newly created file for gdalwarper. So just use affine */
4174 : /* approximation of the reverse transform. Eventually we */
4175 : /* should implement iterative searching to find a result within */
4176 : /* our error threshold. */
4177 : /* NOTE: the above comment is not true: gdalwarp uses approximator */
4178 : /* also to compute the source pixel of each target pixel. */
4179 : /* -------------------------------------------------------------------- */
4180 97365100 : for (int i = nPoints - 1; i >= 0; i--)
4181 : {
4182 : #ifdef check_error
4183 : double xtemp = x[i];
4184 : double ytemp = y[i];
4185 : double ztemp = z[i];
4186 : double x_ori = xtemp;
4187 : double y_ori = ytemp;
4188 : int btemp = FALSE;
4189 : psATInfo->pfnBaseTransformer(psATInfo->pBaseCBData, bDstToSrc, 1,
4190 : &xtemp, &ytemp, &ztemp, &btemp);
4191 : #endif
4192 96843300 : const double dfDist = (x[i] - x[0]);
4193 96843300 : x[i] = xSMETransformed[0] + dfDeltaX * dfDist;
4194 96843300 : y[i] = ySMETransformed[0] + dfDeltaY * dfDist;
4195 96843300 : z[i] = zSMETransformed[0] + dfDeltaZ * dfDist;
4196 : #ifdef check_error
4197 : const double dfError2 = fabs(x[i] - xtemp) + fabs(y[i] - ytemp);
4198 : if (dfError2 > 4 /*10 * dfMaxError*/)
4199 : {
4200 : /*ok*/ printf("Error = %f on (%f, %f)\n", dfError2, x_ori, y_ori);
4201 : }
4202 : #endif
4203 96843300 : panSuccess[i] = TRUE;
4204 : }
4205 :
4206 521776 : return TRUE;
4207 : }
4208 :
4209 : /************************************************************************/
4210 : /* GDALApproxTransform() */
4211 : /************************************************************************/
4212 :
4213 : /**
4214 : * Perform approximate transformation.
4215 : *
4216 : * Actually performs the approximate transformation described in
4217 : * GDALCreateApproxTransformer(). This function matches the
4218 : * GDALTransformerFunc() signature. Details of the arguments are described
4219 : * there.
4220 : */
4221 :
4222 515328 : int GDALApproxTransform(void *pCBData, int bDstToSrc, int nPoints, double *x,
4223 : double *y, double *z, int *panSuccess)
4224 :
4225 : {
4226 515328 : ApproxTransformInfo *psATInfo = static_cast<ApproxTransformInfo *>(pCBData);
4227 515328 : double x2[3] = {};
4228 515328 : double y2[3] = {};
4229 515328 : double z2[3] = {};
4230 515328 : int anSuccess2[3] = {};
4231 : int bSuccess;
4232 :
4233 515328 : const int nMiddle = (nPoints - 1) / 2;
4234 :
4235 : /* -------------------------------------------------------------------- */
4236 : /* Bail if our preconditions are not met, or if error is not */
4237 : /* acceptable. */
4238 : /* -------------------------------------------------------------------- */
4239 515328 : int bRet = FALSE;
4240 515328 : if (y[0] != y[nPoints - 1] || y[0] != y[nMiddle] ||
4241 507144 : x[0] == x[nPoints - 1] || x[0] == x[nMiddle] ||
4242 178881 : (psATInfo->dfMaxErrorForward == 0.0 &&
4243 178881 : psATInfo->dfMaxErrorReverse == 0.0) ||
4244 : nPoints <= 5)
4245 : {
4246 336781 : bRet = psATInfo->pfnBaseTransformer(psATInfo->pBaseCBData, bDstToSrc,
4247 : nPoints, x, y, z, panSuccess);
4248 336770 : goto end;
4249 : }
4250 :
4251 : /* -------------------------------------------------------------------- */
4252 : /* Transform first, last and middle point. */
4253 : /* -------------------------------------------------------------------- */
4254 178547 : x2[0] = x[0];
4255 178547 : y2[0] = y[0];
4256 178547 : z2[0] = z[0];
4257 178547 : x2[1] = x[nMiddle];
4258 178547 : y2[1] = y[nMiddle];
4259 178547 : z2[1] = z[nMiddle];
4260 178547 : x2[2] = x[nPoints - 1];
4261 178547 : y2[2] = y[nPoints - 1];
4262 178547 : z2[2] = z[nPoints - 1];
4263 :
4264 178547 : bSuccess = psATInfo->pfnBaseTransformer(psATInfo->pBaseCBData, bDstToSrc, 3,
4265 : x2, y2, z2, anSuccess2);
4266 178547 : if (!bSuccess || !anSuccess2[0] || !anSuccess2[1] || !anSuccess2[2])
4267 : {
4268 5130 : bRet = psATInfo->pfnBaseTransformer(psATInfo->pBaseCBData, bDstToSrc,
4269 : nPoints, x, y, z, panSuccess);
4270 5123 : goto end;
4271 : }
4272 :
4273 173417 : bRet = GDALApproxTransformInternal(pCBData, bDstToSrc, nPoints, x, y, z,
4274 : panSuccess, x2, y2, z2);
4275 :
4276 515318 : end:
4277 : #ifdef DEBUG_APPROX_TRANSFORMER
4278 : for (int i = 0; i < nPoints; i++)
4279 : fprintf(stderr, "[%d] (%.10f,%.10f) %d\n", /*ok*/
4280 : i, x[i], y[i], panSuccess[i]);
4281 : #endif
4282 :
4283 515318 : return bRet;
4284 : }
4285 :
4286 : /************************************************************************/
4287 : /* GDALDeserializeApproxTransformer() */
4288 : /************************************************************************/
4289 :
4290 150 : static void *GDALDeserializeApproxTransformer(CPLXMLNode *psTree)
4291 :
4292 : {
4293 150 : double dfMaxErrorForward = 0.25;
4294 150 : double dfMaxErrorReverse = 0.25;
4295 150 : const char *pszMaxError = CPLGetXMLValue(psTree, "MaxError", nullptr);
4296 150 : if (pszMaxError != nullptr)
4297 : {
4298 148 : dfMaxErrorForward = CPLAtof(pszMaxError);
4299 148 : dfMaxErrorReverse = dfMaxErrorForward;
4300 : }
4301 : const char *pszMaxErrorForward =
4302 150 : CPLGetXMLValue(psTree, "MaxErrorForward", nullptr);
4303 150 : if (pszMaxErrorForward != nullptr)
4304 : {
4305 2 : dfMaxErrorForward = CPLAtof(pszMaxErrorForward);
4306 : }
4307 : const char *pszMaxErrorReverse =
4308 150 : CPLGetXMLValue(psTree, "MaxErrorReverse", nullptr);
4309 150 : if (pszMaxErrorReverse != nullptr)
4310 : {
4311 2 : dfMaxErrorReverse = CPLAtof(pszMaxErrorReverse);
4312 : }
4313 :
4314 150 : GDALTransformerFunc pfnBaseTransform = nullptr;
4315 150 : void *pBaseCBData = nullptr;
4316 :
4317 150 : CPLXMLNode *psContainer = CPLGetXMLNode(psTree, "BaseTransformer");
4318 :
4319 150 : if (psContainer != nullptr && psContainer->psChild != nullptr)
4320 : {
4321 150 : GDALDeserializeTransformer(psContainer->psChild, &pfnBaseTransform,
4322 : &pBaseCBData);
4323 : }
4324 :
4325 150 : if (pfnBaseTransform == nullptr)
4326 : {
4327 0 : CPLError(CE_Failure, CPLE_AppDefined,
4328 : "Cannot get base transform for approx transformer.");
4329 0 : return nullptr;
4330 : }
4331 :
4332 150 : void *pApproxCBData = GDALCreateApproxTransformer2(
4333 : pfnBaseTransform, pBaseCBData, dfMaxErrorForward, dfMaxErrorReverse);
4334 150 : GDALApproxTransformerOwnsSubtransformer(pApproxCBData, TRUE);
4335 :
4336 150 : return pApproxCBData;
4337 : }
4338 :
4339 : /************************************************************************/
4340 : /* GDALTransformLonLatToDestApproxTransformer() */
4341 : /************************************************************************/
4342 :
4343 2090 : int GDALTransformLonLatToDestApproxTransformer(void *hTransformArg,
4344 : double *pdfX, double *pdfY)
4345 : {
4346 2090 : ApproxTransformInfo *psInfo =
4347 : static_cast<ApproxTransformInfo *>(hTransformArg);
4348 :
4349 2090 : if (GDALIsTransformer(psInfo->pBaseCBData,
4350 : GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME))
4351 : {
4352 2090 : return GDALTransformLonLatToDestGenImgProjTransformer(
4353 2090 : psInfo->pBaseCBData, pdfX, pdfY);
4354 : }
4355 0 : return false;
4356 : }
4357 :
4358 : /************************************************************************/
4359 : /* GDALApplyGeoTransform() */
4360 : /************************************************************************/
4361 :
4362 : /**
4363 : * Apply GeoTransform to x/y coordinate.
4364 : *
4365 : * Applies the following computation, converting a (pixel, line) coordinate
4366 : * into a georeferenced (geo_x, geo_y) location.
4367 : * \code{.c}
4368 : * *pdfGeoX = padfGeoTransform[0] + dfPixel * padfGeoTransform[1]
4369 : * + dfLine * padfGeoTransform[2];
4370 : * *pdfGeoY = padfGeoTransform[3] + dfPixel * padfGeoTransform[4]
4371 : * + dfLine * padfGeoTransform[5];
4372 : * \endcode
4373 : *
4374 : * @param padfGeoTransform Six coefficient GeoTransform to apply.
4375 : * @param dfPixel Input pixel position.
4376 : * @param dfLine Input line position.
4377 : * @param pdfGeoX output location where geo_x (easting/longitude)
4378 : * location is placed.
4379 : * @param pdfGeoY output location where geo_y (northing/latitude)
4380 : * location is placed.
4381 : */
4382 :
4383 304732 : void CPL_STDCALL GDALApplyGeoTransform(const double *padfGeoTransform,
4384 : double dfPixel, double dfLine,
4385 : double *pdfGeoX, double *pdfGeoY)
4386 : {
4387 304732 : *pdfGeoX = padfGeoTransform[0] + dfPixel * padfGeoTransform[1] +
4388 304732 : dfLine * padfGeoTransform[2];
4389 304732 : *pdfGeoY = padfGeoTransform[3] + dfPixel * padfGeoTransform[4] +
4390 304732 : dfLine * padfGeoTransform[5];
4391 304732 : }
4392 :
4393 : /************************************************************************/
4394 : /* GDALInvGeoTransform() */
4395 : /************************************************************************/
4396 :
4397 : /**
4398 : * Invert Geotransform.
4399 : *
4400 : * This function will invert a standard 3x2 set of GeoTransform coefficients.
4401 : * This converts the equation from being pixel to geo to being geo to pixel.
4402 : *
4403 : * @param gt_in Input geotransform (six doubles - unaltered).
4404 : * @param gt_out Output geotransform (six doubles - updated).
4405 : *
4406 : * @return TRUE on success or FALSE if the equation is uninvertable.
4407 : */
4408 :
4409 2942 : int CPL_STDCALL GDALInvGeoTransform(const double *gt_in, double *gt_out)
4410 :
4411 : {
4412 : // Special case - no rotation - to avoid computing determinate
4413 : // and potential precision issues.
4414 2942 : if (gt_in[2] == 0.0 && gt_in[4] == 0.0 && gt_in[1] != 0.0 &&
4415 2888 : gt_in[5] != 0.0)
4416 : {
4417 : /*X = gt_in[0] + x * gt_in[1]
4418 : Y = gt_in[3] + y * gt_in[5]
4419 : -->
4420 : x = -gt_in[0] / gt_in[1] + (1 / gt_in[1]) * X
4421 : y = -gt_in[3] / gt_in[5] + (1 / gt_in[5]) * Y
4422 : */
4423 2888 : gt_out[0] = -gt_in[0] / gt_in[1];
4424 2888 : gt_out[1] = 1.0 / gt_in[1];
4425 2888 : gt_out[2] = 0.0;
4426 2888 : gt_out[3] = -gt_in[3] / gt_in[5];
4427 2888 : gt_out[4] = 0.0;
4428 2888 : gt_out[5] = 1.0 / gt_in[5];
4429 2888 : return 1;
4430 : }
4431 :
4432 : // Assume a 3rd row that is [1 0 0].
4433 :
4434 : // Compute determinate.
4435 :
4436 54 : const double det = gt_in[1] * gt_in[5] - gt_in[2] * gt_in[4];
4437 108 : const double magnitude = std::max(std::max(fabs(gt_in[1]), fabs(gt_in[2])),
4438 54 : std::max(fabs(gt_in[4]), fabs(gt_in[5])));
4439 :
4440 54 : if (fabs(det) <= 1e-10 * magnitude * magnitude)
4441 5 : return 0;
4442 :
4443 49 : const double inv_det = 1.0 / det;
4444 :
4445 : // Compute adjoint, and divide by determinate.
4446 :
4447 49 : gt_out[1] = gt_in[5] * inv_det;
4448 49 : gt_out[4] = -gt_in[4] * inv_det;
4449 :
4450 49 : gt_out[2] = -gt_in[2] * inv_det;
4451 49 : gt_out[5] = gt_in[1] * inv_det;
4452 :
4453 49 : gt_out[0] = (gt_in[2] * gt_in[3] - gt_in[0] * gt_in[5]) * inv_det;
4454 49 : gt_out[3] = (-gt_in[1] * gt_in[3] + gt_in[0] * gt_in[4]) * inv_det;
4455 :
4456 49 : return 1;
4457 : }
4458 :
4459 : /************************************************************************/
4460 : /* GDALSerializeTransformer() */
4461 : /************************************************************************/
4462 :
4463 263 : CPLXMLNode *GDALSerializeTransformer(GDALTransformerFunc /* pfnFunc */,
4464 : void *pTransformArg)
4465 : {
4466 263 : VALIDATE_POINTER1(pTransformArg, "GDALSerializeTransformer", nullptr);
4467 :
4468 263 : GDALTransformerInfo *psInfo =
4469 : static_cast<GDALTransformerInfo *>(pTransformArg);
4470 :
4471 263 : if (psInfo == nullptr || memcmp(psInfo->abySignature, GDAL_GTI2_SIGNATURE,
4472 : strlen(GDAL_GTI2_SIGNATURE)) != 0)
4473 : {
4474 0 : CPLError(CE_Failure, CPLE_AppDefined,
4475 : "Attempt to serialize non-GTI2 transformer.");
4476 0 : return nullptr;
4477 : }
4478 263 : else if (psInfo->pfnSerialize == nullptr)
4479 : {
4480 0 : CPLError(CE_Failure, CPLE_AppDefined,
4481 : "No serialization function available for this transformer.");
4482 0 : return nullptr;
4483 : }
4484 :
4485 263 : return psInfo->pfnSerialize(pTransformArg);
4486 : }
4487 :
4488 : /************************************************************************/
4489 : /* GDALRegisterTransformDeserializer() */
4490 : /************************************************************************/
4491 :
4492 : static CPLList *psListDeserializer = nullptr;
4493 : static CPLMutex *hDeserializerMutex = nullptr;
4494 :
4495 : typedef struct
4496 : {
4497 : char *pszTransformName;
4498 : GDALTransformerFunc pfnTransformerFunc;
4499 : GDALTransformDeserializeFunc pfnDeserializeFunc;
4500 : } TransformDeserializerInfo;
4501 :
4502 0 : void *GDALRegisterTransformDeserializer(
4503 : const char *pszTransformName, GDALTransformerFunc pfnTransformerFunc,
4504 : GDALTransformDeserializeFunc pfnDeserializeFunc)
4505 : {
4506 : TransformDeserializerInfo *psInfo =
4507 : static_cast<TransformDeserializerInfo *>(
4508 0 : CPLMalloc(sizeof(TransformDeserializerInfo)));
4509 0 : psInfo->pszTransformName = CPLStrdup(pszTransformName);
4510 0 : psInfo->pfnTransformerFunc = pfnTransformerFunc;
4511 0 : psInfo->pfnDeserializeFunc = pfnDeserializeFunc;
4512 :
4513 0 : CPLMutexHolderD(&hDeserializerMutex);
4514 0 : psListDeserializer = CPLListInsert(psListDeserializer, psInfo, 0);
4515 :
4516 0 : return psInfo;
4517 : }
4518 :
4519 : /************************************************************************/
4520 : /* GDALUnregisterTransformDeserializer() */
4521 : /************************************************************************/
4522 :
4523 0 : void GDALUnregisterTransformDeserializer(void *pData)
4524 : {
4525 0 : CPLMutexHolderD(&hDeserializerMutex);
4526 0 : CPLList *psList = psListDeserializer;
4527 0 : CPLList *psLast = nullptr;
4528 0 : while (psList)
4529 : {
4530 0 : if (psList->pData == pData)
4531 : {
4532 0 : TransformDeserializerInfo *psInfo =
4533 : static_cast<TransformDeserializerInfo *>(pData);
4534 0 : CPLFree(psInfo->pszTransformName);
4535 0 : CPLFree(pData);
4536 0 : if (psLast)
4537 0 : psLast->psNext = psList->psNext;
4538 : else
4539 0 : psListDeserializer = nullptr;
4540 0 : CPLFree(psList);
4541 0 : break;
4542 : }
4543 0 : psLast = psList;
4544 0 : psList = psList->psNext;
4545 : }
4546 0 : }
4547 :
4548 : /************************************************************************/
4549 : /* GDALUnregisterTransformDeserializer() */
4550 : /************************************************************************/
4551 :
4552 941 : void GDALCleanupTransformDeserializerMutex()
4553 : {
4554 941 : if (hDeserializerMutex != nullptr)
4555 : {
4556 0 : CPLDestroyMutex(hDeserializerMutex);
4557 0 : hDeserializerMutex = nullptr;
4558 : }
4559 941 : }
4560 :
4561 : /************************************************************************/
4562 : /* GDALDeserializeTransformer() */
4563 : /************************************************************************/
4564 :
4565 538 : CPLErr GDALDeserializeTransformer(CPLXMLNode *psTree,
4566 : GDALTransformerFunc *ppfnFunc,
4567 : void **ppTransformArg)
4568 :
4569 : {
4570 538 : *ppfnFunc = nullptr;
4571 538 : *ppTransformArg = nullptr;
4572 :
4573 538 : CPLErrorReset();
4574 :
4575 538 : if (psTree == nullptr || psTree->eType != CXT_Element)
4576 0 : CPLError(CE_Failure, CPLE_AppDefined,
4577 : "Malformed element in GDALDeserializeTransformer");
4578 538 : else if (EQUAL(psTree->pszValue, "GenImgProjTransformer"))
4579 : {
4580 198 : *ppfnFunc = GDALGenImgProjTransform;
4581 198 : *ppTransformArg = GDALDeserializeGenImgProjTransformer(psTree);
4582 : }
4583 340 : else if (EQUAL(psTree->pszValue, "ReprojectionTransformer"))
4584 : {
4585 174 : *ppfnFunc = GDALReprojectionTransform;
4586 174 : *ppTransformArg = GDALDeserializeReprojectionTransformer(psTree);
4587 : }
4588 166 : else if (EQUAL(psTree->pszValue, "GCPTransformer"))
4589 : {
4590 12 : *ppfnFunc = GDALGCPTransform;
4591 12 : *ppTransformArg = GDALDeserializeGCPTransformer(psTree);
4592 : }
4593 154 : else if (EQUAL(psTree->pszValue, "TPSTransformer"))
4594 : {
4595 3 : *ppfnFunc = GDALTPSTransform;
4596 3 : *ppTransformArg = GDALDeserializeTPSTransformer(psTree);
4597 : }
4598 151 : else if (EQUAL(psTree->pszValue, "GeoLocTransformer"))
4599 : {
4600 1 : *ppfnFunc = GDALGeoLocTransform;
4601 1 : *ppTransformArg = GDALDeserializeGeoLocTransformer(psTree);
4602 : }
4603 150 : else if (EQUAL(psTree->pszValue, "RPCTransformer"))
4604 : {
4605 0 : *ppfnFunc = GDALRPCTransform;
4606 0 : *ppTransformArg = GDALDeserializeRPCTransformer(psTree);
4607 : }
4608 150 : else if (EQUAL(psTree->pszValue, "ApproxTransformer"))
4609 : {
4610 150 : *ppfnFunc = GDALApproxTransform;
4611 150 : *ppTransformArg = GDALDeserializeApproxTransformer(psTree);
4612 : }
4613 : else
4614 : {
4615 0 : GDALTransformDeserializeFunc pfnDeserializeFunc = nullptr;
4616 : {
4617 0 : CPLMutexHolderD(&hDeserializerMutex);
4618 0 : CPLList *psList = psListDeserializer;
4619 0 : while (psList)
4620 : {
4621 0 : TransformDeserializerInfo *psInfo =
4622 : static_cast<TransformDeserializerInfo *>(psList->pData);
4623 0 : if (strcmp(psInfo->pszTransformName, psTree->pszValue) == 0)
4624 : {
4625 0 : *ppfnFunc = psInfo->pfnTransformerFunc;
4626 0 : pfnDeserializeFunc = psInfo->pfnDeserializeFunc;
4627 0 : break;
4628 : }
4629 0 : psList = psList->psNext;
4630 : }
4631 : }
4632 :
4633 0 : if (pfnDeserializeFunc != nullptr)
4634 : {
4635 0 : *ppTransformArg = pfnDeserializeFunc(psTree);
4636 : }
4637 : else
4638 : {
4639 0 : CPLError(CE_Failure, CPLE_AppDefined,
4640 : "Unrecognized element '%s' GDALDeserializeTransformer",
4641 : psTree->pszValue);
4642 : }
4643 : }
4644 :
4645 538 : return CPLGetLastErrorType();
4646 : }
4647 :
4648 : /************************************************************************/
4649 : /* GDALDestroyTransformer() */
4650 : /************************************************************************/
4651 :
4652 3493 : void GDALDestroyTransformer(void *pTransformArg)
4653 :
4654 : {
4655 3493 : if (pTransformArg == nullptr)
4656 10 : return;
4657 :
4658 3483 : GDALTransformerInfo *psInfo =
4659 : static_cast<GDALTransformerInfo *>(pTransformArg);
4660 :
4661 3483 : if (memcmp(psInfo->abySignature, GDAL_GTI2_SIGNATURE,
4662 : strlen(GDAL_GTI2_SIGNATURE)) != 0)
4663 : {
4664 0 : CPLError(CE_Failure, CPLE_AppDefined,
4665 : "Attempt to destroy non-GTI2 transformer.");
4666 0 : return;
4667 : }
4668 :
4669 3483 : psInfo->pfnCleanup(pTransformArg);
4670 : }
4671 :
4672 : /************************************************************************/
4673 : /* GDALUseTransformer() */
4674 : /************************************************************************/
4675 :
4676 8680 : int GDALUseTransformer(void *pTransformArg, int bDstToSrc, int nPointCount,
4677 : double *x, double *y, double *z, int *panSuccess)
4678 : {
4679 8680 : GDALTransformerInfo *psInfo =
4680 : static_cast<GDALTransformerInfo *>(pTransformArg);
4681 :
4682 8680 : if (psInfo == nullptr || memcmp(psInfo->abySignature, GDAL_GTI2_SIGNATURE,
4683 : strlen(GDAL_GTI2_SIGNATURE)) != 0)
4684 : {
4685 0 : CPLError(CE_Failure, CPLE_AppDefined,
4686 : "Attempt to use non-GTI2 transformer.");
4687 0 : return FALSE;
4688 : }
4689 :
4690 8680 : return psInfo->pfnTransform(pTransformArg, bDstToSrc, nPointCount, x, y, z,
4691 8680 : panSuccess);
4692 : }
4693 :
4694 : /************************************************************************/
4695 : /* GDALCloneTransformer() */
4696 : /************************************************************************/
4697 :
4698 23 : void *GDALCloneTransformer(void *pTransformArg)
4699 : {
4700 23 : GDALTransformerInfo *psInfo =
4701 : static_cast<GDALTransformerInfo *>(pTransformArg);
4702 :
4703 23 : if (psInfo == nullptr || memcmp(psInfo->abySignature, GDAL_GTI2_SIGNATURE,
4704 : strlen(GDAL_GTI2_SIGNATURE)) != 0)
4705 : {
4706 0 : CPLError(CE_Failure, CPLE_AppDefined,
4707 : "Attempt to clone non-GTI2 transformer.");
4708 0 : return nullptr;
4709 : }
4710 :
4711 23 : if (psInfo->pfnCreateSimilar != nullptr)
4712 : {
4713 12 : return psInfo->pfnCreateSimilar(psInfo, 1.0, 1.0);
4714 : }
4715 :
4716 11 : if (psInfo->pfnSerialize == nullptr)
4717 : {
4718 0 : CPLError(CE_Failure, CPLE_AppDefined,
4719 : "No serialization function available for this transformer.");
4720 0 : return nullptr;
4721 : }
4722 :
4723 11 : CPLXMLNode *pSerialized = psInfo->pfnSerialize(pTransformArg);
4724 11 : if (pSerialized == nullptr)
4725 0 : return nullptr;
4726 11 : GDALTransformerFunc pfnTransformer = nullptr;
4727 11 : void *pClonedTransformArg = nullptr;
4728 11 : if (GDALDeserializeTransformer(pSerialized, &pfnTransformer,
4729 11 : &pClonedTransformArg) != CE_None)
4730 : {
4731 0 : CPLDestroyXMLNode(pSerialized);
4732 0 : CPLFree(pClonedTransformArg);
4733 0 : return nullptr;
4734 : }
4735 :
4736 11 : CPLDestroyXMLNode(pSerialized);
4737 11 : return pClonedTransformArg;
4738 : }
4739 :
4740 : /************************************************************************/
4741 : /* GDALCreateSimilarTransformer() */
4742 : /************************************************************************/
4743 :
4744 44 : void *GDALCreateSimilarTransformer(void *pTransformArg, double dfRatioX,
4745 : double dfRatioY)
4746 : {
4747 44 : GDALTransformerInfo *psInfo =
4748 : static_cast<GDALTransformerInfo *>(pTransformArg);
4749 :
4750 44 : if (psInfo == nullptr || memcmp(psInfo->abySignature, GDAL_GTI2_SIGNATURE,
4751 : strlen(GDAL_GTI2_SIGNATURE)) != 0)
4752 : {
4753 0 : CPLError(CE_Failure, CPLE_AppDefined,
4754 : "Attempt to call CreateSimilar on a non-GTI2 transformer.");
4755 0 : return nullptr;
4756 : }
4757 :
4758 44 : if (psInfo->pfnCreateSimilar == nullptr)
4759 : {
4760 0 : CPLError(CE_Failure, CPLE_AppDefined,
4761 : "No CreateSimilar function available for this transformer.");
4762 0 : return nullptr;
4763 : }
4764 :
4765 44 : return psInfo->pfnCreateSimilar(psInfo, dfRatioX, dfRatioY);
4766 : }
4767 :
4768 : /************************************************************************/
4769 : /* GetGenImgProjTransformInfo() */
4770 : /************************************************************************/
4771 :
4772 44 : static GDALTransformerInfo *GetGenImgProjTransformInfo(const char *pszFunc,
4773 : void *pTransformArg)
4774 : {
4775 44 : GDALTransformerInfo *psInfo =
4776 : static_cast<GDALTransformerInfo *>(pTransformArg);
4777 :
4778 44 : if (psInfo == nullptr || memcmp(psInfo->abySignature, GDAL_GTI2_SIGNATURE,
4779 : strlen(GDAL_GTI2_SIGNATURE)) != 0)
4780 : {
4781 0 : CPLError(CE_Failure, CPLE_AppDefined,
4782 : "Attempt to call %s on "
4783 : "a non-GTI2 transformer.",
4784 : pszFunc);
4785 0 : return nullptr;
4786 : }
4787 :
4788 44 : if (EQUAL(psInfo->pszClassName, GDAL_APPROX_TRANSFORMER_CLASS_NAME))
4789 : {
4790 14 : ApproxTransformInfo *psATInfo =
4791 : static_cast<ApproxTransformInfo *>(pTransformArg);
4792 14 : psInfo = static_cast<GDALTransformerInfo *>(psATInfo->pBaseCBData);
4793 :
4794 14 : if (psInfo == nullptr ||
4795 14 : memcmp(psInfo->abySignature, GDAL_GTI2_SIGNATURE,
4796 : strlen(GDAL_GTI2_SIGNATURE)) != 0)
4797 : {
4798 0 : CPLError(CE_Failure, CPLE_AppDefined,
4799 : "Attempt to call %s on "
4800 : "a non-GTI2 transformer.",
4801 : pszFunc);
4802 0 : return nullptr;
4803 : }
4804 : }
4805 :
4806 44 : if (EQUAL(psInfo->pszClassName, GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME))
4807 : {
4808 44 : return psInfo;
4809 : }
4810 :
4811 0 : return nullptr;
4812 : }
4813 :
4814 : /************************************************************************/
4815 : /* GDALSetTransformerDstGeoTransform() */
4816 : /************************************************************************/
4817 :
4818 : /**
4819 : * Set ApproxTransformer or GenImgProj output geotransform.
4820 : *
4821 : * This is a layer above GDALSetGenImgProjTransformerDstGeoTransform() that
4822 : * checks that the passed hTransformArg is compatible.
4823 : *
4824 : * Normally the "destination geotransform", or transformation between
4825 : * georeferenced output coordinates and pixel/line coordinates on the
4826 : * destination file is extracted from the destination file by
4827 : * GDALCreateGenImgProjTransformer() and stored in the GenImgProj private
4828 : * info. However, sometimes it is inconvenient to have an output file
4829 : * handle with appropriate geotransform information when creating the
4830 : * transformation. For these cases, this function can be used to apply
4831 : * the destination geotransform.
4832 : *
4833 : * @param pTransformArg the handle to update.
4834 : * @param padfGeoTransform the destination geotransform to apply (six doubles).
4835 : */
4836 :
4837 22 : void GDALSetTransformerDstGeoTransform(void *pTransformArg,
4838 : const double *padfGeoTransform)
4839 : {
4840 22 : VALIDATE_POINTER0(pTransformArg, "GDALSetTransformerDstGeoTransform");
4841 :
4842 22 : GDALTransformerInfo *psInfo = GetGenImgProjTransformInfo(
4843 : "GDALSetTransformerDstGeoTransform", pTransformArg);
4844 22 : if (psInfo)
4845 : {
4846 22 : GDALSetGenImgProjTransformerDstGeoTransform(psInfo, padfGeoTransform);
4847 : }
4848 : }
4849 :
4850 : /************************************************************************/
4851 : /* GDALGetTransformerDstGeoTransform() */
4852 : /************************************************************************/
4853 :
4854 : /**
4855 : * Get ApproxTransformer or GenImgProj output geotransform.
4856 : *
4857 : * @param pTransformArg transformer handle.
4858 : * @param padfGeoTransform (output) the destination geotransform to return (six
4859 : * doubles).
4860 : */
4861 :
4862 22 : void GDALGetTransformerDstGeoTransform(void *pTransformArg,
4863 : double *padfGeoTransform)
4864 : {
4865 22 : VALIDATE_POINTER0(pTransformArg, "GDALGetTransformerDstGeoTransform");
4866 :
4867 22 : GDALTransformerInfo *psInfo = GetGenImgProjTransformInfo(
4868 : "GDALGetTransformerDstGeoTransform", pTransformArg);
4869 22 : if (psInfo)
4870 : {
4871 22 : GDALGenImgProjTransformInfo *psGenImgProjInfo =
4872 : reinterpret_cast<GDALGenImgProjTransformInfo *>(psInfo);
4873 :
4874 22 : memcpy(padfGeoTransform, psGenImgProjInfo->adfDstGeoTransform,
4875 : sizeof(double) * 6);
4876 : }
4877 : }
4878 :
4879 : /************************************************************************/
4880 : /* GDALTransformIsTranslationOnPixelBoundaries() */
4881 : /************************************************************************/
4882 :
4883 1401 : bool GDALTransformIsTranslationOnPixelBoundaries(GDALTransformerFunc,
4884 : void *pTransformerArg)
4885 : {
4886 1401 : if (GDALIsTransformer(pTransformerArg, GDAL_APPROX_TRANSFORMER_CLASS_NAME))
4887 : {
4888 1045 : const auto *pApproxInfo =
4889 : static_cast<const ApproxTransformInfo *>(pTransformerArg);
4890 1045 : pTransformerArg = pApproxInfo->pBaseCBData;
4891 : }
4892 1401 : if (GDALIsTransformer(pTransformerArg, GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME))
4893 : {
4894 1248 : const auto *pGenImgpProjInfo =
4895 : static_cast<GDALGenImgProjTransformInfo *>(pTransformerArg);
4896 394 : const auto IsCloseToInteger = [](double dfVal)
4897 394 : { return std::fabs(dfVal - std::round(dfVal)) <= 1e-6; };
4898 2417 : return pGenImgpProjInfo->pSrcTransformArg == nullptr &&
4899 1169 : pGenImgpProjInfo->pDstTransformArg == nullptr &&
4900 1167 : pGenImgpProjInfo->pReproject == nullptr &&
4901 472 : pGenImgpProjInfo->adfSrcGeoTransform[1] ==
4902 472 : pGenImgpProjInfo->adfDstGeoTransform[1] &&
4903 229 : pGenImgpProjInfo->adfSrcGeoTransform[5] ==
4904 229 : pGenImgpProjInfo->adfDstGeoTransform[5] &&
4905 215 : pGenImgpProjInfo->adfSrcGeoTransform[2] ==
4906 215 : pGenImgpProjInfo->adfDstGeoTransform[2] &&
4907 215 : pGenImgpProjInfo->adfSrcGeoTransform[4] ==
4908 430 : pGenImgpProjInfo->adfDstGeoTransform[4] &&
4909 : // Check that the georeferenced origin of the destination
4910 : // geotransform is close to be an integer value when transformed
4911 : // to source image coordinates
4912 215 : IsCloseToInteger(
4913 215 : pGenImgpProjInfo->adfSrcInvGeoTransform[0] +
4914 215 : pGenImgpProjInfo->adfDstGeoTransform[0] *
4915 215 : pGenImgpProjInfo->adfSrcInvGeoTransform[1] +
4916 215 : pGenImgpProjInfo->adfDstGeoTransform[3] *
4917 2632 : pGenImgpProjInfo->adfSrcInvGeoTransform[2]) &&
4918 179 : IsCloseToInteger(pGenImgpProjInfo->adfSrcInvGeoTransform[3] +
4919 179 : pGenImgpProjInfo->adfDstGeoTransform[0] *
4920 179 : pGenImgpProjInfo->adfSrcInvGeoTransform[4] +
4921 179 : pGenImgpProjInfo->adfDstGeoTransform[3] *
4922 1427 : pGenImgpProjInfo->adfSrcInvGeoTransform[5]);
4923 : }
4924 153 : return false;
4925 : }
4926 :
4927 : /************************************************************************/
4928 : /* GDALTransformIsAffineNoRotation() */
4929 : /************************************************************************/
4930 :
4931 18 : bool GDALTransformIsAffineNoRotation(GDALTransformerFunc, void *pTransformerArg)
4932 : {
4933 18 : if (GDALIsTransformer(pTransformerArg, GDAL_APPROX_TRANSFORMER_CLASS_NAME))
4934 : {
4935 18 : const auto *pApproxInfo =
4936 : static_cast<const ApproxTransformInfo *>(pTransformerArg);
4937 18 : pTransformerArg = pApproxInfo->pBaseCBData;
4938 : }
4939 18 : if (GDALIsTransformer(pTransformerArg, GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME))
4940 : {
4941 18 : const auto *pGenImgpProjInfo =
4942 : static_cast<GDALGenImgProjTransformInfo *>(pTransformerArg);
4943 36 : return pGenImgpProjInfo->pSrcTransformArg == nullptr &&
4944 18 : pGenImgpProjInfo->pDstTransformArg == nullptr &&
4945 18 : pGenImgpProjInfo->pReproject == nullptr &&
4946 8 : pGenImgpProjInfo->adfSrcGeoTransform[2] == 0 &&
4947 8 : pGenImgpProjInfo->adfSrcGeoTransform[4] == 0 &&
4948 44 : pGenImgpProjInfo->adfDstGeoTransform[2] == 0 &&
4949 26 : pGenImgpProjInfo->adfDstGeoTransform[4] == 0;
4950 : }
4951 0 : return false;
4952 : }
4953 :
4954 : /************************************************************************/
4955 : /* GDALTransformHasFastClone() */
4956 : /************************************************************************/
4957 :
4958 : /** Returns whether GDALCloneTransformer() on this transformer is
4959 : * "fast"
4960 : * Counter-examples are GCPs or TPSs transformers.
4961 : */
4962 2 : bool GDALTransformHasFastClone(void *pTransformerArg)
4963 : {
4964 2 : if (GDALIsTransformer(pTransformerArg, GDAL_APPROX_TRANSFORMER_CLASS_NAME))
4965 : {
4966 1 : const auto *pApproxInfo =
4967 : static_cast<const ApproxTransformInfo *>(pTransformerArg);
4968 1 : pTransformerArg = pApproxInfo->pBaseCBData;
4969 : // Fallback to next lines
4970 : }
4971 :
4972 2 : if (GDALIsTransformer(pTransformerArg, GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME))
4973 : {
4974 2 : const auto *pGenImgpProjInfo =
4975 : static_cast<GDALGenImgProjTransformInfo *>(pTransformerArg);
4976 2 : return (pGenImgpProjInfo->pSrcTransformArg == nullptr ||
4977 0 : GDALTransformHasFastClone(
4978 4 : pGenImgpProjInfo->pSrcTransformArg)) &&
4979 2 : (pGenImgpProjInfo->pDstTransformArg == nullptr ||
4980 2 : GDALTransformHasFastClone(pGenImgpProjInfo->pDstTransformArg));
4981 : }
4982 0 : else if (GDALIsTransformer(pTransformerArg,
4983 : GDAL_RPC_TRANSFORMER_CLASS_NAME))
4984 : {
4985 0 : return true;
4986 : }
4987 : else
4988 : {
4989 0 : return false;
4990 : }
4991 : }
|
