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