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