Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: Implements Geolocation array based transformer.
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2006, Frank Warmerdam <warmerdam@pobox.com>
9 : * Copyright (c) 2007-2013, Even Rouault <even dot rouault at spatialys.com>
10 : * Copyright (c) 2021, CLS
11 : * Copyright (c) 2022, Planet Labs
12 : *
13 : * SPDX-License-Identifier: MIT
14 : ****************************************************************************/
15 :
16 : #include "cpl_port.h"
17 : #include "gdal_alg.h"
18 : #include "gdal_alg_priv.h"
19 : #include "gdalgeoloc.h"
20 : #include "gdalgeolocquadtree.h"
21 :
22 : #include <climits>
23 : #include <cmath>
24 : #include <cstddef>
25 : #include <cstdlib>
26 : #include <cstring>
27 :
28 : #include <algorithm>
29 : #include <limits>
30 :
31 : #include "cpl_conv.h"
32 : #include "cpl_error.h"
33 : #include "cpl_minixml.h"
34 : #include "cpl_quad_tree.h"
35 : #include "cpl_string.h"
36 : #include "cpl_vsi.h"
37 : #include "gdal.h"
38 : #include "gdal_priv.h"
39 : #include "memdataset.h"
40 :
41 : constexpr float INVALID_BMXY = -10.0f;
42 :
43 : #include "gdalgeoloc_carray_accessor.h"
44 : #include "gdalgeoloc_dataset_accessor.h"
45 :
46 : // #define DEBUG_GEOLOC
47 :
48 : #ifdef DEBUG_GEOLOC
49 : #include "ogrsf_frmts.h"
50 : #endif
51 :
52 : #ifdef DEBUG_GEOLOC
53 : #warning "Remove me before committing"
54 : #endif
55 :
56 : CPL_C_START
57 : CPLXMLNode *GDALSerializeGeoLocTransformer(void *pTransformArg);
58 : void *GDALDeserializeGeoLocTransformer(CPLXMLNode *psTree);
59 : CPL_C_END
60 :
61 : /************************************************************************/
62 : /* ==================================================================== */
63 : /* GDALGeoLocTransformer */
64 : /* ==================================================================== */
65 : /************************************************************************/
66 :
67 : /************************************************************************/
68 : /* UnshiftGeoX() */
69 : /************************************************************************/
70 :
71 : // Renormalize longitudes to [-180,180] range
72 12063200 : static double UnshiftGeoX(const GDALGeoLocTransformInfo *psTransform,
73 : double dfX)
74 : {
75 12063200 : if (!psTransform->bGeographicSRSWithMinus180Plus180LongRange)
76 2747650 : return dfX;
77 9315580 : if (dfX > 180 || dfX < -180)
78 : {
79 35493 : dfX = fmod(dfX + 180.0, 360.0);
80 35493 : if (dfX < 0)
81 372 : dfX += 180.0;
82 : else
83 35121 : dfX -= 180.0;
84 35493 : return dfX;
85 : }
86 9280080 : return dfX;
87 : }
88 :
89 : /************************************************************************/
90 : /* UpdateMinMax() */
91 : /************************************************************************/
92 :
93 253897 : inline void UpdateMinMax(GDALGeoLocTransformInfo *psTransform, double dfGeoLocX,
94 : double dfGeoLocY)
95 : {
96 253897 : if (dfGeoLocX < psTransform->dfMinX)
97 : {
98 1350 : psTransform->dfMinX = dfGeoLocX;
99 1350 : psTransform->dfYAtMinX = dfGeoLocY;
100 : }
101 253897 : if (dfGeoLocX > psTransform->dfMaxX)
102 : {
103 1095 : psTransform->dfMaxX = dfGeoLocX;
104 1095 : psTransform->dfYAtMaxX = dfGeoLocY;
105 : }
106 253897 : if (dfGeoLocY < psTransform->dfMinY)
107 : {
108 1393 : psTransform->dfMinY = dfGeoLocY;
109 1393 : psTransform->dfXAtMinY = dfGeoLocX;
110 : }
111 253897 : if (dfGeoLocY > psTransform->dfMaxY)
112 : {
113 3641 : psTransform->dfMaxY = dfGeoLocY;
114 3641 : psTransform->dfXAtMaxY = dfGeoLocX;
115 : }
116 253897 : }
117 :
118 : /************************************************************************/
119 : /* Clamp() */
120 : /************************************************************************/
121 :
122 3268 : inline double Clamp(double v, double minV, double maxV)
123 : {
124 3268 : return std::min(std::max(v, minV), maxV);
125 : }
126 :
127 : /************************************************************************/
128 : /* START_ITER_PER_BLOCK() */
129 : /************************************************************************/
130 :
131 : #define START_ITER_PER_BLOCK(_rasterXSize, _tileXSize, _rasterYSize, \
132 : _tileYSize, INIT_YBLOCK, _iXStart, _iXEnd, \
133 : _iYStart, _iYEnd) \
134 : { \
135 : const int _nYBlocks = DIV_ROUND_UP(_rasterYSize, _tileYSize); \
136 : const int _nXBlocks = DIV_ROUND_UP(_rasterXSize, _tileXSize); \
137 : for (int _iYBlock = 0; _iYBlock < _nYBlocks; ++_iYBlock) \
138 : { \
139 : const int _iYStart = _iYBlock * _tileYSize; \
140 : const int _iYEnd = _iYBlock == _nYBlocks - 1 \
141 : ? _rasterYSize \
142 : : _iYStart + _tileYSize; \
143 : INIT_YBLOCK; \
144 : for (int _iXBlock = 0; _iXBlock < _nXBlocks; ++_iXBlock) \
145 : { \
146 : const int _iXStart = _iXBlock * _tileXSize; \
147 : const int _iXEnd = _iXBlock == _nXBlocks - 1 \
148 : ? _rasterXSize \
149 : : _iXStart + _tileXSize;
150 :
151 : #define END_ITER_PER_BLOCK \
152 : } \
153 : } \
154 : }
155 :
156 : /************************************************************************/
157 : /* GDALGeoLoc::LoadGeolocFinish() */
158 : /************************************************************************/
159 :
160 : /*! @cond Doxygen_Suppress */
161 :
162 : template <class Accessors>
163 53 : void GDALGeoLoc<Accessors>::LoadGeolocFinish(
164 : GDALGeoLocTransformInfo *psTransform)
165 : {
166 53 : auto pAccessors = static_cast<Accessors *>(psTransform->pAccessors);
167 53 : CSLConstList papszGeolocationInfo = psTransform->papszGeolocationInfo;
168 :
169 : /* -------------------------------------------------------------------- */
170 : /* Scan forward map for lat/long extents. */
171 : /* -------------------------------------------------------------------- */
172 53 : psTransform->dfMinX = std::numeric_limits<double>::max();
173 53 : psTransform->dfMaxX = -std::numeric_limits<double>::max();
174 53 : psTransform->dfMinY = std::numeric_limits<double>::max();
175 53 : psTransform->dfMaxY = -std::numeric_limits<double>::max();
176 :
177 53 : constexpr int TILE_SIZE = GDALGeoLocDatasetAccessors::TILE_SIZE;
178 165 : START_ITER_PER_BLOCK(psTransform->nGeoLocXSize, TILE_SIZE,
179 : psTransform->nGeoLocYSize, TILE_SIZE, (void)0, iXStart,
180 : iXEnd, iYStart, iYEnd)
181 : {
182 2127 : for (int iY = iYStart; iY < iYEnd; ++iY)
183 : {
184 254379 : for (int iX = iXStart; iX < iXEnd; ++iX)
185 : {
186 252309 : double dfX = pAccessors->geolocXAccessor.Get(iX, iY);
187 252309 : if (!psTransform->bHasNoData || dfX != psTransform->dfNoDataX)
188 : {
189 249367 : dfX = UnshiftGeoX(psTransform, dfX);
190 249367 : UpdateMinMax(psTransform, dfX,
191 : pAccessors->geolocYAccessor.Get(iX, iY));
192 : }
193 : }
194 : }
195 : }
196 : END_ITER_PER_BLOCK
197 :
198 : // Check if the SRS is geographic and the geoloc longitudes are in
199 : // [-180,180]
200 53 : const char *pszSRS = CSLFetchNameValue(papszGeolocationInfo, "SRS");
201 53 : if (!psTransform->bGeographicSRSWithMinus180Plus180LongRange && pszSRS &&
202 49 : psTransform->dfMinX >= -180.0 && psTransform->dfMaxX <= 180.0 &&
203 48 : !psTransform->bSwapXY)
204 : {
205 47 : OGRSpatialReference oSRS;
206 47 : psTransform->bGeographicSRSWithMinus180Plus180LongRange =
207 94 : oSRS.importFromWkt(pszSRS) == OGRERR_NONE &&
208 47 : CPL_TO_BOOL(oSRS.IsGeographic());
209 : }
210 :
211 : #ifdef DEBUG_GEOLOC
212 : if (CPLTestBool(CPLGetConfigOption("GEOLOC_DUMP", "NO")))
213 : {
214 : auto poDS = std::unique_ptr<GDALDataset>(
215 : GDALDriver::FromHandle(GDALGetDriverByName("ESRI Shapefile"))
216 : ->Create("/tmp/geoloc_poly.shp", 0, 0, 0, GDT_Unknown,
217 : nullptr));
218 : auto poLayer =
219 : poDS->CreateLayer("geoloc_poly", nullptr, wkbPolygon, nullptr);
220 : auto poLayerDefn = poLayer->GetLayerDefn();
221 : OGRFieldDefn fieldX("x", OFTInteger);
222 : poLayer->CreateField(&fieldX);
223 : OGRFieldDefn fieldY("y", OFTInteger);
224 : poLayer->CreateField(&fieldY);
225 : for (int iY = 0; iY < psTransform->nGeoLocYSize - 1; iY++)
226 : {
227 : for (int iX = 0; iX < psTransform->nGeoLocXSize - 1; iX++)
228 : {
229 : double x0, y0, x1, y1, x2, y2, x3, y3;
230 : if (!PixelLineToXY(psTransform, iX, iY, x0, y0) ||
231 : !PixelLineToXY(psTransform, iX + 1, iY, x2, y2) ||
232 : !PixelLineToXY(psTransform, iX, iY + 1, x1, y1) ||
233 : !PixelLineToXY(psTransform, iX + 1, iY + 1, x3, y3))
234 : {
235 : break;
236 : }
237 : if (psTransform->bGeographicSRSWithMinus180Plus180LongRange &&
238 : std::fabs(x0) > 170 && std::fabs(x1) > 170 &&
239 : std::fabs(x2) > 170 && std::fabs(x3) > 170 &&
240 : (std::fabs(x1 - x0) > 180 || std::fabs(x2 - x0) > 180 ||
241 : std::fabs(x3 - x0) > 180))
242 : {
243 : OGRPolygon *poPoly = new OGRPolygon();
244 : OGRLinearRing *poRing = new OGRLinearRing();
245 : poRing->addPoint(x0 > 0 ? x0 : x0 + 360, y0);
246 : poRing->addPoint(x2 > 0 ? x2 : x2 + 360, y2);
247 : poRing->addPoint(x3 > 0 ? x3 : x3 + 360, y3);
248 : poRing->addPoint(x1 > 0 ? x1 : x1 + 360, y1);
249 : poRing->addPoint(x0 > 0 ? x0 : x0 + 360, y0);
250 : poPoly->addRingDirectly(poRing);
251 : auto poFeature = std::make_unique<OGRFeature>(poLayerDefn);
252 : poFeature->SetField(0, static_cast<int>(iX));
253 : poFeature->SetField(1, static_cast<int>(iY));
254 : poFeature->SetGeometryDirectly(poPoly);
255 : CPL_IGNORE_RET_VAL(poLayer->CreateFeature(poFeature.get()));
256 : if (x0 > 0)
257 : x0 -= 360;
258 : if (x1 > 0)
259 : x1 -= 360;
260 : if (x2 > 0)
261 : x2 -= 360;
262 : if (x3 > 0)
263 : x3 -= 360;
264 : }
265 :
266 : OGRPolygon *poPoly = new OGRPolygon();
267 : OGRLinearRing *poRing = new OGRLinearRing();
268 : poRing->addPoint(x0, y0);
269 : poRing->addPoint(x2, y2);
270 : poRing->addPoint(x3, y3);
271 : poRing->addPoint(x1, y1);
272 : poRing->addPoint(x0, y0);
273 : poPoly->addRingDirectly(poRing);
274 : auto poFeature = std::make_unique<OGRFeature>(poLayerDefn);
275 : poFeature->SetField(0, static_cast<int>(iX));
276 : poFeature->SetField(1, static_cast<int>(iY));
277 : poFeature->SetGeometryDirectly(poPoly);
278 : CPL_IGNORE_RET_VAL(poLayer->CreateFeature(poFeature.get()));
279 : }
280 : }
281 : }
282 : #endif
283 :
284 53 : if (psTransform->bOriginIsTopLeftCorner)
285 : {
286 : // Add "virtual" edge at Y=nGeoLocYSize
287 1983 : for (int iX = 0; iX <= psTransform->nGeoLocXSize; iX++)
288 : {
289 : double dfGeoLocX;
290 : double dfGeoLocY;
291 1954 : if (!PixelLineToXY(psTransform, static_cast<double>(iX),
292 1954 : static_cast<double>(psTransform->nGeoLocYSize),
293 : dfGeoLocX, dfGeoLocY))
294 0 : continue;
295 1954 : if (psTransform->bGeographicSRSWithMinus180Plus180LongRange)
296 1089 : dfGeoLocX = Clamp(dfGeoLocX, -180.0, 180.0);
297 1954 : UpdateMinMax(psTransform, dfGeoLocX, dfGeoLocY);
298 : }
299 :
300 : // Add "virtual" edge at X=nGeoLocXSize
301 1555 : for (int iY = 0; iY <= psTransform->nGeoLocYSize; iY++)
302 : {
303 : double dfGeoLocX;
304 : double dfGeoLocY;
305 1526 : if (!PixelLineToXY(psTransform,
306 1526 : static_cast<double>(psTransform->nGeoLocXSize),
307 : static_cast<double>(iY), dfGeoLocX, dfGeoLocY))
308 0 : continue;
309 1526 : if (psTransform->bGeographicSRSWithMinus180Plus180LongRange)
310 1221 : dfGeoLocX = Clamp(dfGeoLocX, -180.0, 180.0);
311 1526 : UpdateMinMax(psTransform, dfGeoLocX, dfGeoLocY);
312 : }
313 : }
314 : else
315 : {
316 : // Extend by half-pixel on 4 edges for pixel-center convention
317 :
318 401 : for (int iX = 0; iX <= psTransform->nGeoLocXSize; iX++)
319 : {
320 : double dfGeoLocX;
321 : double dfGeoLocY;
322 377 : if (!PixelLineToXY(psTransform, static_cast<double>(iX), -0.5,
323 : dfGeoLocX, dfGeoLocY))
324 114 : continue;
325 263 : if (psTransform->bGeographicSRSWithMinus180Plus180LongRange)
326 242 : dfGeoLocX = Clamp(dfGeoLocX, -180.0, 180.0);
327 263 : UpdateMinMax(psTransform, dfGeoLocX, dfGeoLocY);
328 : }
329 :
330 401 : for (int iX = 0; iX <= psTransform->nGeoLocXSize; iX++)
331 : {
332 : double dfGeoLocX;
333 : double dfGeoLocY;
334 377 : if (!PixelLineToXY(
335 : psTransform, static_cast<double>(iX),
336 377 : static_cast<double>(psTransform->nGeoLocYSize - 1 + 0.5),
337 : dfGeoLocX, dfGeoLocY))
338 118 : continue;
339 259 : if (psTransform->bGeographicSRSWithMinus180Plus180LongRange)
340 238 : dfGeoLocX = Clamp(dfGeoLocX, -180.0, 180.0);
341 259 : UpdateMinMax(psTransform, dfGeoLocX, dfGeoLocY);
342 : }
343 :
344 461 : for (int iY = 0; iY <= psTransform->nGeoLocYSize; iY++)
345 : {
346 : double dfGeoLocX;
347 : double dfGeoLocY;
348 437 : if (!PixelLineToXY(psTransform, -0.5, static_cast<double>(iY),
349 : dfGeoLocX, dfGeoLocY))
350 170 : continue;
351 267 : if (psTransform->bGeographicSRSWithMinus180Plus180LongRange)
352 242 : dfGeoLocX = Clamp(dfGeoLocX, -180.0, 180.0);
353 267 : UpdateMinMax(psTransform, dfGeoLocX, dfGeoLocY);
354 : }
355 :
356 461 : for (int iY = 0; iY <= psTransform->nGeoLocYSize; iY++)
357 : {
358 : double dfGeoLocX;
359 : double dfGeoLocY;
360 437 : if (!PixelLineToXY(psTransform, psTransform->nGeoLocXSize - 1 + 0.5,
361 : static_cast<double>(iY), dfGeoLocX, dfGeoLocY))
362 176 : continue;
363 261 : if (psTransform->bGeographicSRSWithMinus180Plus180LongRange)
364 236 : dfGeoLocX = Clamp(dfGeoLocX, -180.0, 180.0);
365 261 : UpdateMinMax(psTransform, dfGeoLocX, dfGeoLocY);
366 : }
367 : }
368 53 : }
369 :
370 : /************************************************************************/
371 : /* GDALGeoLoc::PixelLineToXY() */
372 : /************************************************************************/
373 :
374 : /** Interpolate a position expressed as (floating point) pixel/line in the
375 : * geolocation array to the corresponding bilinearly-interpolated georeferenced
376 : * position.
377 : *
378 : * The interpolation assumes infinite extension beyond borders of available
379 : * data based on closest grid square.
380 : *
381 : * @param psTransform Transformation info
382 : * @param dfGeoLocPixel Position along the column/pixel axis of the geolocation
383 : * array
384 : * @param dfGeoLocLine Position along the row/line axis of the geolocation
385 : * array
386 : * @param[out] dfX Output X of georeferenced position.
387 : * @param[out] dfY Output Y of georeferenced position.
388 : * @return true if success
389 : */
390 :
391 : template <class Accessors>
392 1460784 : bool GDALGeoLoc<Accessors>::PixelLineToXY(
393 : const GDALGeoLocTransformInfo *psTransform, const double dfGeoLocPixel,
394 : const double dfGeoLocLine, double &dfX, double &dfY)
395 : {
396 1460784 : int iX = static_cast<int>(
397 2921568 : std::min(std::max(0.0, dfGeoLocPixel),
398 1460784 : static_cast<double>(psTransform->nGeoLocXSize - 1)));
399 1460784 : int iY = static_cast<int>(
400 2921568 : std::min(std::max(0.0, dfGeoLocLine),
401 1460784 : static_cast<double>(psTransform->nGeoLocYSize - 1)));
402 :
403 1460784 : auto pAccessors = static_cast<Accessors *>(psTransform->pAccessors);
404 :
405 1852441 : for (int iAttempt = 0; iAttempt < 2; ++iAttempt)
406 : {
407 1852441 : const double dfGLX_0_0 = pAccessors->geolocXAccessor.Get(iX, iY);
408 1852441 : const double dfGLY_0_0 = pAccessors->geolocYAccessor.Get(iX, iY);
409 1852441 : if (psTransform->bHasNoData && dfGLX_0_0 == psTransform->dfNoDataX)
410 : {
411 13386 : return false;
412 : }
413 :
414 : // This assumes infinite extension beyond borders of available
415 : // data based on closest grid square.
416 1839051 : if (iX + 1 < psTransform->nGeoLocXSize &&
417 1615807 : iY + 1 < psTransform->nGeoLocYSize)
418 : {
419 1447394 : const double dfGLX_1_0 =
420 : pAccessors->geolocXAccessor.Get(iX + 1, iY);
421 1447394 : const double dfGLY_1_0 =
422 : pAccessors->geolocYAccessor.Get(iX + 1, iY);
423 1447394 : const double dfGLX_0_1 =
424 : pAccessors->geolocXAccessor.Get(iX, iY + 1);
425 1447394 : const double dfGLY_0_1 =
426 : pAccessors->geolocYAccessor.Get(iX, iY + 1);
427 1447394 : const double dfGLX_1_1 =
428 : pAccessors->geolocXAccessor.Get(iX + 1, iY + 1);
429 1447394 : const double dfGLY_1_1 =
430 : pAccessors->geolocYAccessor.Get(iX + 1, iY + 1);
431 1447394 : if (!psTransform->bHasNoData ||
432 189529 : (dfGLX_1_0 != psTransform->dfNoDataX &&
433 187991 : dfGLX_0_1 != psTransform->dfNoDataX &&
434 186841 : dfGLX_1_1 != psTransform->dfNoDataX))
435 : {
436 : const double dfGLX_1_0_adjusted =
437 1444514 : ShiftGeoX(psTransform, dfGLX_0_0, dfGLX_1_0);
438 : const double dfGLX_0_1_adjusted =
439 1444514 : ShiftGeoX(psTransform, dfGLX_0_0, dfGLX_0_1);
440 : const double dfGLX_1_1_adjusted =
441 1444514 : ShiftGeoX(psTransform, dfGLX_0_0, dfGLX_1_1);
442 1444514 : dfX = (1 - (dfGeoLocLine - iY)) *
443 1444514 : (dfGLX_0_0 + (dfGeoLocPixel - iX) *
444 1444514 : (dfGLX_1_0_adjusted - dfGLX_0_0)) +
445 1444514 : (dfGeoLocLine - iY) *
446 1444514 : (dfGLX_0_1_adjusted +
447 1444514 : (dfGeoLocPixel - iX) *
448 1444514 : (dfGLX_1_1_adjusted - dfGLX_0_1_adjusted));
449 1444514 : dfX = UnshiftGeoX(psTransform, dfX);
450 :
451 1444514 : dfY = (1 - (dfGeoLocLine - iY)) *
452 1444514 : (dfGLY_0_0 +
453 1444514 : (dfGeoLocPixel - iX) * (dfGLY_1_0 - dfGLY_0_0)) +
454 1444514 : (dfGeoLocLine - iY) *
455 1444514 : (dfGLY_0_1 +
456 1444514 : (dfGeoLocPixel - iX) * (dfGLY_1_1 - dfGLY_0_1));
457 1444514 : break;
458 : }
459 : }
460 :
461 394538 : if (iX == psTransform->nGeoLocXSize - 1 && iX >= 1 &&
462 223247 : iY + 1 < psTransform->nGeoLocYSize)
463 : {
464 : // If we are after the right edge, then go one pixel left
465 : // and retry
466 207090 : iX--;
467 207090 : continue;
468 : }
469 187448 : else if (iY == psTransform->nGeoLocYSize - 1 && iY >= 1 &&
470 184568 : iX + 1 < psTransform->nGeoLocXSize)
471 : {
472 : // If we are after the bottom edge, then go one pixel up
473 : // and retry
474 168411 : iY--;
475 168411 : continue;
476 : }
477 19037 : else if (iX == psTransform->nGeoLocXSize - 1 &&
478 16157 : iY == psTransform->nGeoLocYSize - 1 && iX >= 1 && iY >= 1)
479 : {
480 : // If we are after the right and bottom edge, then go one pixel left
481 : // and up and retry
482 16157 : iX--;
483 16157 : iY--;
484 16157 : continue;
485 : }
486 5760 : else if (iX + 1 < psTransform->nGeoLocXSize &&
487 2880 : (!psTransform->bHasNoData ||
488 2880 : pAccessors->geolocXAccessor.Get(iX + 1, iY) !=
489 2880 : psTransform->dfNoDataX))
490 : {
491 1342 : const double dfGLX_1_0 =
492 : pAccessors->geolocXAccessor.Get(iX + 1, iY);
493 1342 : const double dfGLY_1_0 =
494 : pAccessors->geolocYAccessor.Get(iX + 1, iY);
495 1342 : dfX =
496 1342 : dfGLX_0_0 +
497 2684 : (dfGeoLocPixel - iX) *
498 1342 : (ShiftGeoX(psTransform, dfGLX_0_0, dfGLX_1_0) - dfGLX_0_0);
499 1342 : dfX = UnshiftGeoX(psTransform, dfX);
500 1342 : dfY = dfGLY_0_0 + (dfGeoLocPixel - iX) * (dfGLY_1_0 - dfGLY_0_0);
501 : }
502 3076 : else if (iY + 1 < psTransform->nGeoLocYSize &&
503 1538 : (!psTransform->bHasNoData ||
504 1538 : pAccessors->geolocXAccessor.Get(iX, iY + 1) !=
505 1538 : psTransform->dfNoDataX))
506 : {
507 1486 : const double dfGLX_0_1 =
508 : pAccessors->geolocXAccessor.Get(iX, iY + 1);
509 1486 : const double dfGLY_0_1 =
510 : pAccessors->geolocYAccessor.Get(iX, iY + 1);
511 1486 : dfX =
512 1486 : dfGLX_0_0 +
513 2972 : (dfGeoLocLine - iY) *
514 1486 : (ShiftGeoX(psTransform, dfGLX_0_0, dfGLX_0_1) - dfGLX_0_0);
515 1486 : dfX = UnshiftGeoX(psTransform, dfX);
516 1486 : dfY = dfGLY_0_0 + (dfGeoLocLine - iY) * (dfGLY_0_1 - dfGLY_0_0);
517 : }
518 : else
519 : {
520 52 : dfX = UnshiftGeoX(psTransform, dfGLX_0_0);
521 52 : dfY = dfGLY_0_0;
522 : }
523 2880 : break;
524 : }
525 1447394 : return true;
526 : }
527 :
528 : template <class Accessors>
529 11375640 : bool GDALGeoLoc<Accessors>::PixelLineToXY(
530 : const GDALGeoLocTransformInfo *psTransform, const int nGeoLocPixel,
531 : const int nGeoLocLine, double &dfX, double &dfY)
532 : {
533 11375640 : if (nGeoLocPixel >= 0 && nGeoLocPixel < psTransform->nGeoLocXSize &&
534 10563540 : nGeoLocLine >= 0 && nGeoLocLine < psTransform->nGeoLocYSize)
535 : {
536 10417200 : auto pAccessors = static_cast<Accessors *>(psTransform->pAccessors);
537 10417200 : const double dfGLX =
538 : pAccessors->geolocXAccessor.Get(nGeoLocPixel, nGeoLocLine);
539 10417200 : const double dfGLY =
540 : pAccessors->geolocYAccessor.Get(nGeoLocPixel, nGeoLocLine);
541 10417200 : if (psTransform->bHasNoData && dfGLX == psTransform->dfNoDataX)
542 : {
543 50742 : return false;
544 : }
545 10366460 : dfX = UnshiftGeoX(psTransform, dfGLX);
546 10366460 : dfY = dfGLY;
547 10366460 : return true;
548 : }
549 958432 : return PixelLineToXY(psTransform, static_cast<double>(nGeoLocPixel),
550 958432 : static_cast<double>(nGeoLocLine), dfX, dfY);
551 : }
552 :
553 : /************************************************************************/
554 : /* GDALGeoLoc::ExtractSquare() */
555 : /************************************************************************/
556 :
557 : template <class Accessors>
558 3200 : bool GDALGeoLoc<Accessors>::ExtractSquare(
559 : const GDALGeoLocTransformInfo *psTransform, int nX, int nY, double &dfX_0_0,
560 : double &dfY_0_0, double &dfX_1_0, double &dfY_1_0, double &dfX_0_1,
561 : double &dfY_0_1, double &dfX_1_1, double &dfY_1_1)
562 : {
563 6400 : return PixelLineToXY(psTransform, nX, nY, dfX_0_0, dfY_0_0) &&
564 6400 : PixelLineToXY(psTransform, nX + 1, nY, dfX_1_0, dfY_1_0) &&
565 9600 : PixelLineToXY(psTransform, nX, nY + 1, dfX_0_1, dfY_0_1) &&
566 6400 : PixelLineToXY(psTransform, nX + 1, nY + 1, dfX_1_1, dfY_1_1);
567 : }
568 :
569 3200 : bool GDALGeoLocExtractSquare(const GDALGeoLocTransformInfo *psTransform, int nX,
570 : int nY, double &dfX_0_0, double &dfY_0_0,
571 : double &dfX_1_0, double &dfY_1_0, double &dfX_0_1,
572 : double &dfY_0_1, double &dfX_1_1, double &dfY_1_1)
573 : {
574 3200 : if (psTransform->bUseArray)
575 : {
576 3200 : return GDALGeoLoc<GDALGeoLocCArrayAccessors>::ExtractSquare(
577 : psTransform, nX, nY, dfX_0_0, dfY_0_0, dfX_1_0, dfY_1_0, dfX_0_1,
578 3200 : dfY_0_1, dfX_1_1, dfY_1_1);
579 : }
580 : else
581 : {
582 0 : return GDALGeoLoc<GDALGeoLocDatasetAccessors>::ExtractSquare(
583 : psTransform, nX, nY, dfX_0_0, dfY_0_0, dfX_1_0, dfY_1_0, dfX_0_1,
584 0 : dfY_0_1, dfX_1_1, dfY_1_1);
585 : }
586 : }
587 :
588 : /************************************************************************/
589 : /* GDALGeoLocTransform() */
590 : /************************************************************************/
591 :
592 : template <class Accessors>
593 32118 : int GDALGeoLoc<Accessors>::Transform(void *pTransformArg, int bDstToSrc,
594 : int nPointCount, double *padfX,
595 : double *padfY, double * /* padfZ */,
596 : int *panSuccess)
597 : {
598 32118 : int bSuccess = TRUE;
599 32118 : GDALGeoLocTransformInfo *psTransform =
600 : static_cast<GDALGeoLocTransformInfo *>(pTransformArg);
601 :
602 32118 : if (psTransform->bReversed)
603 0 : bDstToSrc = !bDstToSrc;
604 :
605 32118 : const double dfGeorefConventionOffset =
606 32118 : psTransform->bOriginIsTopLeftCorner ? 0 : 0.5;
607 :
608 : /* -------------------------------------------------------------------- */
609 : /* Do original pixel line to target geox/geoy. */
610 : /* -------------------------------------------------------------------- */
611 32118 : if (!bDstToSrc)
612 : {
613 76744 : for (int i = 0; i < nPointCount; i++)
614 : {
615 51795 : if (padfX[i] == HUGE_VAL || padfY[i] == HUGE_VAL)
616 : {
617 3061 : bSuccess = FALSE;
618 3061 : panSuccess[i] = FALSE;
619 3061 : continue;
620 : }
621 :
622 48734 : const double dfGeoLocPixel =
623 48734 : (padfX[i] - psTransform->dfPIXEL_OFFSET) /
624 48734 : psTransform->dfPIXEL_STEP -
625 : dfGeorefConventionOffset;
626 48734 : const double dfGeoLocLine =
627 48734 : (padfY[i] - psTransform->dfLINE_OFFSET) /
628 48734 : psTransform->dfLINE_STEP -
629 : dfGeorefConventionOffset;
630 :
631 48734 : if (!PixelLineToXY(psTransform, dfGeoLocPixel, dfGeoLocLine,
632 48734 : padfX[i], padfY[i]))
633 : {
634 2826 : bSuccess = FALSE;
635 2826 : panSuccess[i] = FALSE;
636 2826 : padfX[i] = HUGE_VAL;
637 2826 : padfY[i] = HUGE_VAL;
638 2826 : continue;
639 : }
640 :
641 45908 : if (psTransform->bSwapXY)
642 : {
643 568 : std::swap(padfX[i], padfY[i]);
644 : }
645 :
646 45908 : panSuccess[i] = TRUE;
647 : }
648 : }
649 :
650 : /* -------------------------------------------------------------------- */
651 : /* geox/geoy to pixel/line using backmap. */
652 : /* -------------------------------------------------------------------- */
653 : else
654 : {
655 7169 : if (psTransform->hQuadTree)
656 : {
657 24 : GDALGeoLocInverseTransformQuadtree(psTransform, nPointCount, padfX,
658 : padfY, panSuccess);
659 24 : return TRUE;
660 : }
661 :
662 7145 : const bool bGeolocMaxAccuracy = CPLTestBool(
663 : CPLGetConfigOption("GDAL_GEOLOC_USE_MAX_ACCURACY", "YES"));
664 :
665 : // Keep those objects in this outer scope, so they are reused, to
666 : // save memory allocations.
667 14290 : OGRPoint oPoint;
668 14290 : OGRLinearRing oRing;
669 7145 : oRing.setNumPoints(5);
670 :
671 7145 : auto pAccessors = static_cast<Accessors *>(psTransform->pAccessors);
672 :
673 243302 : for (int i = 0; i < nPointCount; i++)
674 : {
675 236157 : if (padfX[i] == HUGE_VAL || padfY[i] == HUGE_VAL)
676 : {
677 0 : bSuccess = FALSE;
678 0 : panSuccess[i] = FALSE;
679 0 : continue;
680 : }
681 :
682 236157 : if (psTransform->bSwapXY)
683 : {
684 964 : std::swap(padfX[i], padfY[i]);
685 : }
686 :
687 236157 : const double dfGeoX = padfX[i];
688 236157 : const double dfGeoY = padfY[i];
689 :
690 236157 : const double dfBMX =
691 236157 : ((padfX[i] - psTransform->adfBackMapGeoTransform[0]) /
692 : psTransform->adfBackMapGeoTransform[1]);
693 236157 : const double dfBMY =
694 236157 : ((padfY[i] - psTransform->adfBackMapGeoTransform[3]) /
695 : psTransform->adfBackMapGeoTransform[5]);
696 :
697 236157 : if (!(dfBMX >= 0 && dfBMY >= 0 &&
698 235901 : dfBMX + 1 < psTransform->nBackMapWidth &&
699 234494 : dfBMY + 1 < psTransform->nBackMapHeight))
700 : {
701 2577 : bSuccess = FALSE;
702 2577 : panSuccess[i] = FALSE;
703 2577 : padfX[i] = HUGE_VAL;
704 2577 : padfY[i] = HUGE_VAL;
705 2577 : continue;
706 : }
707 :
708 233580 : const int iBMX = static_cast<int>(dfBMX);
709 233580 : const int iBMY = static_cast<int>(dfBMY);
710 :
711 233580 : const auto fBMX_0_0 = pAccessors->backMapXAccessor.Get(iBMX, iBMY);
712 233580 : const auto fBMY_0_0 = pAccessors->backMapYAccessor.Get(iBMX, iBMY);
713 233580 : if (fBMX_0_0 == INVALID_BMXY)
714 : {
715 92043 : bSuccess = FALSE;
716 92043 : panSuccess[i] = FALSE;
717 92043 : padfX[i] = HUGE_VAL;
718 92043 : padfY[i] = HUGE_VAL;
719 92043 : continue;
720 : }
721 :
722 141537 : const auto fBMX_1_0 =
723 : pAccessors->backMapXAccessor.Get(iBMX + 1, iBMY);
724 141537 : const auto fBMY_1_0 =
725 : pAccessors->backMapYAccessor.Get(iBMX + 1, iBMY);
726 141537 : const auto fBMX_0_1 =
727 : pAccessors->backMapXAccessor.Get(iBMX, iBMY + 1);
728 141537 : const auto fBMY_0_1 =
729 : pAccessors->backMapYAccessor.Get(iBMX, iBMY + 1);
730 141537 : const auto fBMX_1_1 =
731 : pAccessors->backMapXAccessor.Get(iBMX + 1, iBMY + 1);
732 141537 : const auto fBMY_1_1 =
733 : pAccessors->backMapYAccessor.Get(iBMX + 1, iBMY + 1);
734 141537 : if (fBMX_1_0 != INVALID_BMXY && fBMX_0_1 != INVALID_BMXY &&
735 : fBMX_1_1 != INVALID_BMXY)
736 : {
737 138931 : padfX[i] = (1 - (dfBMY - iBMY)) *
738 138931 : (double(fBMX_0_0) +
739 138931 : (dfBMX - iBMX) * double(fBMX_1_0 - fBMX_0_0)) +
740 138931 : (dfBMY - iBMY) *
741 138931 : (double(fBMX_0_1) +
742 138931 : (dfBMX - iBMX) * double(fBMX_1_1 - fBMX_0_1));
743 138931 : padfY[i] = (1 - (dfBMY - iBMY)) *
744 138931 : (double(fBMY_0_0) +
745 138931 : (dfBMX - iBMX) * double(fBMY_1_0 - fBMY_0_0)) +
746 138931 : (dfBMY - iBMY) *
747 138931 : (double(fBMY_0_1) +
748 138931 : (dfBMX - iBMX) * double(fBMY_1_1 - fBMY_0_1));
749 : }
750 2606 : else if (fBMX_1_0 != INVALID_BMXY)
751 : {
752 1897 : padfX[i] = double(fBMX_0_0) +
753 1897 : (dfBMX - iBMX) * double(fBMX_1_0 - fBMX_0_0);
754 1897 : padfY[i] = double(fBMY_0_0) +
755 1897 : (dfBMX - iBMX) * double(fBMY_1_0 - fBMY_0_0);
756 : }
757 709 : else if (fBMX_0_1 != INVALID_BMXY)
758 : {
759 528 : padfX[i] = double(fBMX_0_0) +
760 528 : (dfBMY - iBMY) * double(fBMX_0_1 - fBMX_0_0);
761 528 : padfY[i] = double(fBMY_0_0) +
762 528 : (dfBMY - iBMY) * double(fBMY_0_1 - fBMY_0_0);
763 : }
764 : else
765 : {
766 181 : padfX[i] = double(fBMX_0_0);
767 181 : padfY[i] = double(fBMY_0_0);
768 : }
769 :
770 141537 : const double dfGeoLocPixel =
771 141537 : (padfX[i] - psTransform->dfPIXEL_OFFSET) /
772 141537 : psTransform->dfPIXEL_STEP -
773 : dfGeorefConventionOffset;
774 141537 : const double dfGeoLocLine =
775 141537 : (padfY[i] - psTransform->dfLINE_OFFSET) /
776 141537 : psTransform->dfLINE_STEP -
777 : dfGeorefConventionOffset;
778 : #if 0
779 : CPLDebug("GEOLOC", "%f %f %f %f", padfX[i], padfY[i], dfGeoLocPixel, dfGeoLocLine);
780 : if( !psTransform->bOriginIsTopLeftCorner )
781 : {
782 : if( dfGeoLocPixel + dfGeorefConventionOffset > psTransform->nGeoLocXSize-1 ||
783 : dfGeoLocLine + dfGeorefConventionOffset > psTransform->nGeoLocYSize-1 )
784 : {
785 : panSuccess[i] = FALSE;
786 : padfX[i] = HUGE_VAL;
787 : padfY[i] = HUGE_VAL;
788 : continue;
789 : }
790 : }
791 : #endif
792 141537 : if (!bGeolocMaxAccuracy)
793 : {
794 0 : panSuccess[i] = TRUE;
795 0 : continue;
796 : }
797 :
798 : // Now that we have an approximate solution, identify a matching
799 : // cell in the geolocation array, where we can use inverse bilinear
800 : // interpolation to find the exact solution.
801 :
802 : // NOTE: if the geolocation array is an affine transformation,
803 : // the approximate solution should match the exact one, if the
804 : // backmap has correctly been built.
805 :
806 141537 : oPoint.setX(dfGeoX);
807 141537 : oPoint.setY(dfGeoY);
808 : // The thresholds and radius are rather empirical and have been
809 : // tuned on the product
810 : // S5P_TEST_L2__NO2____20190509T220707_20190509T234837_08137_01_010400_20200220T091343.nc
811 : // that includes the north pole.
812 : // Amended with the test case of
813 : // https://github.com/OSGeo/gdal/issues/5823
814 141537 : const int nSearchRadius =
815 141537 : psTransform->bGeographicSRSWithMinus180Plus180LongRange &&
816 84397 : fabs(dfGeoY) >= 85
817 : ? 5
818 : : 3;
819 141537 : const int nGeoLocPixel =
820 141537 : static_cast<int>(std::floor(dfGeoLocPixel));
821 141537 : const int nGeoLocLine = static_cast<int>(std::floor(dfGeoLocLine));
822 :
823 141537 : bool bDone = false;
824 : // Using the above approximate nGeoLocPixel, nGeoLocLine, try to
825 : // find a forward cell that includes (dfGeoX, dfGeoY), with an
826 : // increasing search radius, up to nSearchRadius.
827 384479 : for (int r = 0; !bDone && r <= nSearchRadius; r++)
828 : {
829 1926383 : for (int iter = 0; !bDone && iter < (r == 0 ? 1 : 8 * r);
830 : ++iter)
831 : {
832 : // For r=1, the below formulas will give the following
833 : // offsets:
834 : // (-1,1), (0,1), (1,1), (1,0), (1,-1), (0,-1), (1,-1)
835 3225340 : const int sx = (r == 0) ? 0
836 1541899 : : (iter < 2 * r) ? -r + iter
837 1130365 : : (iter < 4 * r) ? r
838 741355 : : (iter < 6 * r) ? r - (iter - 4 * r)
839 : : -r;
840 3225340 : const int sy = (r == 0) ? 0
841 1541899 : : (iter < 2 * r) ? r
842 1130365 : : (iter < 4 * r) ? r - (iter - 2 * r)
843 741355 : : (iter < 6 * r) ? -r
844 363473 : : -r + (iter - 6 * r);
845 2089447 : if (nGeoLocPixel >=
846 1683435 : static_cast<int>(psTransform->nGeoLocXSize) - sx ||
847 : nGeoLocLine >=
848 1408282 : static_cast<int>(psTransform->nGeoLocYSize) - sy)
849 : {
850 406013 : continue;
851 : }
852 1277426 : const int iX = nGeoLocPixel + sx;
853 1277426 : const int iY = nGeoLocLine + sy;
854 1277426 : if (iX >= -1 || iY >= -1)
855 : {
856 : double x0, y0, x1, y1, x2, y2, x3, y3;
857 :
858 1277003 : if (!PixelLineToXY(psTransform, iX, iY, x0, y0) ||
859 1233947 : !PixelLineToXY(psTransform, iX + 1, iY, x2, y2) ||
860 3735680 : !PixelLineToXY(psTransform, iX, iY + 1, x1, y1) ||
861 1224727 : !PixelLineToXY(psTransform, iX + 1, iY + 1, x3, y3))
862 : {
863 52754 : continue;
864 : }
865 :
866 1224249 : int nIters = 1;
867 : // For a bounding box crossing the anti-meridian, check
868 : // both around -180 and +180 deg.
869 1224249 : if (psTransform
870 1224249 : ->bGeographicSRSWithMinus180Plus180LongRange &&
871 1035101 : std::fabs(x0) > 170 && std::fabs(x1) > 170 &&
872 53449 : std::fabs(x2) > 170 && std::fabs(x3) > 170 &&
873 38624 : (std::fabs(x1 - x0) > 180 ||
874 38293 : std::fabs(x2 - x0) > 180 ||
875 31416 : std::fabs(x3 - x0) > 180))
876 : {
877 7510 : nIters = 2;
878 7510 : if (x0 > 0)
879 58 : x0 -= 360;
880 7510 : if (x1 > 0)
881 331 : x1 -= 360;
882 7510 : if (x2 > 0)
883 7179 : x2 -= 360;
884 7510 : if (x3 > 0)
885 7510 : x3 -= 360;
886 : }
887 2455875 : for (int iIter = 0; !bDone && iIter < nIters; ++iIter)
888 : {
889 1231621 : if (iIter == 1)
890 : {
891 7372 : x0 += 360;
892 7372 : x1 += 360;
893 7372 : x2 += 360;
894 7372 : x3 += 360;
895 : }
896 1231621 : oRing.setPoint(0, x0, y0);
897 1231621 : oRing.setPoint(1, x2, y2);
898 1231621 : oRing.setPoint(2, x3, y3);
899 1231621 : oRing.setPoint(3, x1, y1);
900 1231621 : oRing.setPoint(4, x0, y0);
901 2344090 : if (oRing.isPointInRing(&oPoint) ||
902 1112472 : oRing.isPointOnRingBoundary(&oPoint))
903 : {
904 120969 : double dfX = static_cast<double>(iX);
905 120969 : double dfY = static_cast<double>(iY);
906 120969 : GDALInverseBilinearInterpolation(
907 : dfGeoX, dfGeoY, x0, y0, x1, y1, x2, y2, x3,
908 : y3, dfX, dfY);
909 :
910 120969 : dfX = (dfX + dfGeorefConventionOffset) *
911 120969 : psTransform->dfPIXEL_STEP +
912 120969 : psTransform->dfPIXEL_OFFSET;
913 120969 : dfY = (dfY + dfGeorefConventionOffset) *
914 120969 : psTransform->dfLINE_STEP +
915 120969 : psTransform->dfLINE_OFFSET;
916 :
917 : #ifdef DEBUG_GEOLOC_REALLY_VERBOSE
918 : CPLDebug("GEOLOC",
919 : "value before adjustment: %f %f, "
920 : "after adjustment: %f %f",
921 : padfX[i], padfY[i], dfX, dfY);
922 : #endif
923 :
924 120969 : padfX[i] = dfX;
925 120969 : padfY[i] = dfY;
926 :
927 120969 : bDone = true;
928 : }
929 : }
930 : }
931 : }
932 : }
933 141537 : if (!bDone)
934 : {
935 20568 : bSuccess = FALSE;
936 20568 : panSuccess[i] = FALSE;
937 20568 : padfX[i] = HUGE_VAL;
938 20568 : padfY[i] = HUGE_VAL;
939 20568 : continue;
940 : }
941 :
942 120969 : panSuccess[i] = TRUE;
943 : }
944 : }
945 :
946 32094 : return bSuccess;
947 : }
948 :
949 : /*! @endcond */
950 :
951 : /************************************************************************/
952 : /* GDALInverseBilinearInterpolation() */
953 : /************************************************************************/
954 :
955 : // (i,j) before the call should correspond to the input coordinates that give
956 : // (x0,y0) as output of the forward interpolation
957 : // After the call it will be updated to the input coordinates that give (x,y)
958 : // This assumes that (x,y) is within the polygon formed by
959 : // (x0, y0), (x2, y2), (x3, y3), (x1, y1), (x0, y0)
960 266866 : void GDALInverseBilinearInterpolation(const double x, const double y,
961 : const double x0, const double y0,
962 : const double x1, const double y1,
963 : const double x2, const double y2,
964 : const double x3, const double y3,
965 : double &i, double &j)
966 : {
967 : // Exact inverse bilinear interpolation method.
968 : // Maths from https://stackoverflow.com/a/812077
969 :
970 266866 : const double A = (x0 - x) * (y0 - y2) - (y0 - y) * (x0 - x2);
971 266866 : const double B = (((x0 - x) * (y1 - y3) - (y0 - y) * (x1 - x3)) +
972 266866 : ((x1 - x) * (y0 - y2) - (y1 - y) * (x0 - x2))) /
973 : 2;
974 266866 : const double C = (x1 - x) * (y1 - y3) - (y1 - y) * (x1 - x3);
975 266866 : const double denom = A - 2 * B + C;
976 : double s;
977 266866 : const double magnitudeOfValues = fabs(A) + fabs(B) + fabs(C);
978 266866 : if (fabs(denom) <= 1e-12 * magnitudeOfValues)
979 : {
980 : // Happens typically when the x_i,y_i points form a rectangle
981 : // Can also happen when they form a triangle.
982 143370 : s = A / (A - C);
983 : }
984 : else
985 : {
986 123496 : const double sqrtTerm = sqrt(B * B - A * C);
987 123496 : const double s1 = ((A - B) + sqrtTerm) / denom;
988 123496 : const double s2 = ((A - B) - sqrtTerm) / denom;
989 123496 : if (s1 < 0 || s1 > 1)
990 102726 : s = s2;
991 : else
992 20770 : s = s1;
993 : }
994 :
995 266866 : const double t_denom_x = (1 - s) * (x0 - x2) + s * (x1 - x3);
996 266866 : if (fabs(t_denom_x) > 1e-12 * magnitudeOfValues)
997 : {
998 265071 : i += ((1 - s) * (x0 - x) + s * (x1 - x)) / t_denom_x;
999 : }
1000 : else
1001 : {
1002 1795 : const double t_denom_y = (1 - s) * (y0 - y2) + s * (y1 - y3);
1003 1795 : if (fabs(t_denom_y) > 1e-12 * magnitudeOfValues)
1004 : {
1005 1795 : i += ((1 - s) * (y0 - y) + s * (y1 - y)) / t_denom_y;
1006 : }
1007 : }
1008 :
1009 266866 : j += s;
1010 266866 : }
1011 :
1012 : /************************************************************************/
1013 : /* GeoLocGenerateBackMap() */
1014 : /************************************************************************/
1015 :
1016 : /*! @cond Doxygen_Suppress */
1017 :
1018 : template <class Accessors>
1019 49 : bool GDALGeoLoc<Accessors>::GenerateBackMap(
1020 : GDALGeoLocTransformInfo *psTransform)
1021 :
1022 : {
1023 49 : CPLDebug("GEOLOC", "Starting backmap generation");
1024 49 : const int nXSize = psTransform->nGeoLocXSize;
1025 49 : const int nYSize = psTransform->nGeoLocYSize;
1026 :
1027 : /* -------------------------------------------------------------------- */
1028 : /* Decide on resolution for backmap. We aim for slightly */
1029 : /* higher resolution than the source but we can't easily */
1030 : /* establish how much dead space there is in the backmap, so it */
1031 : /* is approximate. */
1032 : /* -------------------------------------------------------------------- */
1033 49 : const double dfTargetPixels =
1034 49 : static_cast<double>(nXSize) * nYSize * psTransform->dfOversampleFactor;
1035 : const double dfPixelSizeSquare =
1036 49 : sqrt((psTransform->dfMaxX - psTransform->dfMinX) *
1037 49 : (psTransform->dfMaxY - psTransform->dfMinY) / dfTargetPixels);
1038 49 : if (dfPixelSizeSquare == 0.0)
1039 : {
1040 0 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid pixel size for backmap");
1041 0 : return false;
1042 : }
1043 :
1044 49 : const double dfMinX = psTransform->dfMinX - dfPixelSizeSquare / 2.0;
1045 49 : const double dfMaxX = psTransform->dfMaxX + dfPixelSizeSquare / 2.0;
1046 49 : const double dfMaxY = psTransform->dfMaxY + dfPixelSizeSquare / 2.0;
1047 49 : const double dfMinY = psTransform->dfMinY - dfPixelSizeSquare / 2.0;
1048 49 : const double dfBMXSize = std::ceil((dfMaxX - dfMinX) / dfPixelSizeSquare);
1049 49 : const double dfBMYSize = std::ceil((dfMaxY - dfMinY) / dfPixelSizeSquare);
1050 :
1051 : // +2 : +1 due to afterwards nBMXSize++, and another +1 as security margin
1052 : // for other computations.
1053 49 : if (!(dfBMXSize > 0 && dfBMXSize + 2 < INT_MAX) ||
1054 49 : !(dfBMYSize > 0 && dfBMYSize + 2 < INT_MAX))
1055 : {
1056 0 : CPLError(CE_Failure, CPLE_AppDefined, "Int overflow : %f x %f",
1057 : dfBMXSize, dfBMYSize);
1058 0 : return false;
1059 : }
1060 :
1061 49 : int nBMXSize = static_cast<int>(dfBMXSize);
1062 49 : int nBMYSize = static_cast<int>(dfBMYSize);
1063 :
1064 98 : if (static_cast<size_t>(1 + nBMYSize) >
1065 49 : std::numeric_limits<size_t>::max() / static_cast<size_t>(1 + nBMXSize))
1066 : {
1067 0 : CPLError(CE_Failure, CPLE_AppDefined, "Int overflow : %f x %f",
1068 : dfBMXSize, dfBMYSize);
1069 0 : return false;
1070 : }
1071 :
1072 49 : const double dfPixelXSize = (dfMaxX - dfMinX) / nBMXSize;
1073 49 : const double dfPixelYSize = (dfMaxY - dfMinY) / nBMYSize;
1074 :
1075 : // Extra pixel for right-edge and bottom-edge extensions in TOP_LEFT_CORNER
1076 : // convention.
1077 49 : nBMXSize++;
1078 49 : nBMYSize++;
1079 49 : psTransform->nBackMapWidth = nBMXSize;
1080 49 : psTransform->nBackMapHeight = nBMYSize;
1081 :
1082 49 : psTransform->adfBackMapGeoTransform[0] = dfMinX;
1083 49 : psTransform->adfBackMapGeoTransform[1] = dfPixelXSize;
1084 49 : psTransform->adfBackMapGeoTransform[2] = 0.0;
1085 49 : psTransform->adfBackMapGeoTransform[3] = dfMaxY;
1086 49 : psTransform->adfBackMapGeoTransform[4] = 0.0;
1087 49 : psTransform->adfBackMapGeoTransform[5] = -dfPixelYSize;
1088 :
1089 : /* -------------------------------------------------------------------- */
1090 : /* Allocate backmap. */
1091 : /* -------------------------------------------------------------------- */
1092 49 : auto pAccessors = static_cast<Accessors *>(psTransform->pAccessors);
1093 49 : if (!pAccessors->AllocateBackMap())
1094 0 : return false;
1095 :
1096 49 : const double dfGeorefConventionOffset =
1097 49 : psTransform->bOriginIsTopLeftCorner ? 0 : 0.5;
1098 :
1099 49 : const auto UpdateBackmap =
1100 9857135 : [&](int iBMX, int iBMY, double dfX, double dfY, double tempwt)
1101 : {
1102 870170 : const auto fBMX = pAccessors->backMapXAccessor.Get(iBMX, iBMY);
1103 870170 : const auto fBMY = pAccessors->backMapYAccessor.Get(iBMX, iBMY);
1104 870170 : const float fUpdatedBMX =
1105 : fBMX +
1106 870170 : static_cast<float>(tempwt * ((dfX + dfGeorefConventionOffset) *
1107 870170 : psTransform->dfPIXEL_STEP +
1108 870170 : psTransform->dfPIXEL_OFFSET));
1109 870170 : const float fUpdatedBMY =
1110 : fBMY +
1111 870170 : static_cast<float>(tempwt * ((dfY + dfGeorefConventionOffset) *
1112 870170 : psTransform->dfLINE_STEP +
1113 870170 : psTransform->dfLINE_OFFSET));
1114 870170 : const float fUpdatedWeight =
1115 870170 : pAccessors->backMapWeightAccessor.Get(iBMX, iBMY) +
1116 870170 : static_cast<float>(tempwt);
1117 :
1118 : // Only update the backmap if the updated averaged value results in a
1119 : // geoloc position that isn't too different from the original one.
1120 : // (there's no guarantee that if padfGeoLocX[i] ~= padfGeoLoc[j],
1121 : // padfGeoLoc[alpha * i + (1 - alpha) * j] ~= padfGeoLoc[i] )
1122 870170 : if (fUpdatedWeight > 0)
1123 : {
1124 870167 : const float fX = fUpdatedBMX / fUpdatedWeight;
1125 870167 : const float fY = fUpdatedBMY / fUpdatedWeight;
1126 870167 : const double dfGeoLocPixel =
1127 870167 : (double(fX) - psTransform->dfPIXEL_OFFSET) /
1128 870167 : psTransform->dfPIXEL_STEP -
1129 : dfGeorefConventionOffset;
1130 870167 : const double dfGeoLocLine =
1131 870167 : (double(fY) - psTransform->dfLINE_OFFSET) /
1132 870167 : psTransform->dfLINE_STEP -
1133 : dfGeorefConventionOffset;
1134 870167 : int iXAvg = static_cast<int>(std::max(0.0, dfGeoLocPixel));
1135 870167 : iXAvg = std::min(iXAvg, psTransform->nGeoLocXSize - 1);
1136 870167 : int iYAvg = static_cast<int>(std::max(0.0, dfGeoLocLine));
1137 870167 : iYAvg = std::min(iYAvg, psTransform->nGeoLocYSize - 1);
1138 870167 : const double dfGLX = pAccessors->geolocXAccessor.Get(iXAvg, iYAvg);
1139 870167 : const double dfGLY = pAccessors->geolocYAccessor.Get(iXAvg, iYAvg);
1140 :
1141 870167 : const unsigned iX = static_cast<unsigned>(dfX);
1142 870167 : const unsigned iY = static_cast<unsigned>(dfY);
1143 1740320 : if (!(psTransform->bHasNoData && dfGLX == psTransform->dfNoDataX) &&
1144 870157 : ((iX >= static_cast<unsigned>(nXSize - 1) ||
1145 854459 : iY >= static_cast<unsigned>(nYSize - 1)) ||
1146 843268 : (fabs(dfGLX - pAccessors->geolocXAccessor.Get(iX, iY)) <=
1147 843268 : 2 * dfPixelXSize &&
1148 800820 : fabs(dfGLY - pAccessors->geolocYAccessor.Get(iX, iY)) <=
1149 800820 : 2 * dfPixelYSize)))
1150 : {
1151 827564 : pAccessors->backMapXAccessor.Set(iBMX, iBMY, fUpdatedBMX);
1152 827564 : pAccessors->backMapYAccessor.Set(iBMX, iBMY, fUpdatedBMY);
1153 827564 : pAccessors->backMapWeightAccessor.Set(iBMX, iBMY,
1154 : fUpdatedWeight);
1155 : }
1156 : }
1157 : };
1158 :
1159 : // Keep those objects in this outer scope, so they are reused, to
1160 : // save memory allocations.
1161 98 : OGRPoint oPoint;
1162 98 : OGRLinearRing oRing;
1163 49 : oRing.setNumPoints(5);
1164 :
1165 : /* -------------------------------------------------------------------- */
1166 : /* Run through the whole geoloc array forward projecting and */
1167 : /* pushing into the backmap. */
1168 : /* -------------------------------------------------------------------- */
1169 :
1170 : // Iterate over the (i,j) pixel space of the geolocation array, in a
1171 : // sufficiently dense way that if the geolocation array expressed an affine
1172 : // transformation, we would hit every node of the backmap.
1173 49 : const double dfStep = 1. / psTransform->dfOversampleFactor;
1174 :
1175 49 : constexpr int TILE_SIZE = GDALGeoLocDatasetAccessors::TILE_SIZE;
1176 49 : const int nYBlocks = DIV_ROUND_UP(nYSize, TILE_SIZE);
1177 49 : const int nXBlocks = DIV_ROUND_UP(nXSize, TILE_SIZE);
1178 :
1179 : // First compute for each block the start end ending floating-point
1180 : // pixel/line values
1181 98 : std::vector<std::pair<double, double>> yStartEnd(nYBlocks + 1);
1182 98 : std::vector<std::pair<double, double>> xStartEnd(nXBlocks + 1);
1183 :
1184 : {
1185 49 : int iYBlockLast = -1;
1186 49 : double dfY = -dfStep;
1187 2612 : for (; dfY <= static_cast<double>(nYSize) + 2 * dfStep; dfY += dfStep)
1188 : {
1189 2563 : const int iYBlock = static_cast<int>(dfY / TILE_SIZE);
1190 2563 : if (iYBlock != iYBlockLast)
1191 : {
1192 51 : CPLAssert(iYBlock == iYBlockLast + 1);
1193 51 : if (iYBlockLast >= 0)
1194 2 : yStartEnd[iYBlockLast].second = dfY + dfStep / 10;
1195 51 : yStartEnd[iYBlock].first = dfY;
1196 51 : iYBlockLast = iYBlock;
1197 : }
1198 : }
1199 49 : const int iYBlock = static_cast<int>(dfY / TILE_SIZE);
1200 49 : yStartEnd[iYBlock].second = dfY + dfStep / 10;
1201 : }
1202 :
1203 : {
1204 49 : int iXBlockLast = -1;
1205 49 : double dfX = -dfStep;
1206 3094 : for (; dfX <= static_cast<double>(nXSize) + 2 * dfStep; dfX += dfStep)
1207 : {
1208 3045 : const int iXBlock = static_cast<int>(dfX / TILE_SIZE);
1209 3045 : if (iXBlock != iXBlockLast)
1210 : {
1211 51 : CPLAssert(iXBlock == iXBlockLast + 1);
1212 51 : if (iXBlockLast >= 0)
1213 2 : xStartEnd[iXBlockLast].second = dfX + dfStep / 10;
1214 51 : xStartEnd[iXBlock].first = dfX;
1215 51 : iXBlockLast = iXBlock;
1216 : }
1217 : }
1218 49 : const int iXBlock = static_cast<int>(dfX / TILE_SIZE);
1219 49 : xStartEnd[iXBlock].second = dfX + dfStep / 10;
1220 : }
1221 :
1222 100 : for (int iYBlock = 0; iYBlock < nYBlocks; ++iYBlock)
1223 : {
1224 104 : for (int iXBlock = 0; iXBlock < nXBlocks; ++iXBlock)
1225 : {
1226 : #if 0
1227 : CPLDebug("Process geoloc block (y=%d,x=%d) for y in [%f, %f] and x in [%f, %f]",
1228 : iYBlock, iXBlock,
1229 : yStartEnd[iYBlock].first, yStartEnd[iYBlock].second,
1230 : xStartEnd[iXBlock].first, xStartEnd[iXBlock].second);
1231 : #endif
1232 2883 : for (double dfY = yStartEnd[iYBlock].first;
1233 2883 : dfY < yStartEnd[iYBlock].second; dfY += dfStep)
1234 : {
1235 451339 : for (double dfX = xStartEnd[iXBlock].first;
1236 451339 : dfX < xStartEnd[iXBlock].second; dfX += dfStep)
1237 : {
1238 : // Use forward geolocation array interpolation to compute
1239 : // the georeferenced position corresponding to (dfX, dfY)
1240 : double dfGeoLocX;
1241 : double dfGeoLocY;
1242 448509 : if (!PixelLineToXY(psTransform, dfX, dfY, dfGeoLocX,
1243 : dfGeoLocY))
1244 153033 : continue;
1245 :
1246 : // Compute the floating point coordinates in the pixel space
1247 : // of the backmap
1248 441827 : const double dBMX = static_cast<double>(
1249 441827 : (dfGeoLocX - dfMinX) / dfPixelXSize);
1250 :
1251 441827 : const double dBMY = static_cast<double>(
1252 441827 : (dfMaxY - dfGeoLocY) / dfPixelYSize);
1253 :
1254 : // Get top left index by truncation
1255 441827 : const int iBMX = static_cast<int>(std::floor(dBMX));
1256 441827 : const int iBMY = static_cast<int>(std::floor(dBMY));
1257 :
1258 441827 : if (iBMX >= 0 && iBMX < nBMXSize && iBMY >= 0 &&
1259 : iBMY < nBMYSize)
1260 : {
1261 : // Compute the georeferenced position of the top-left
1262 : // index of the backmap
1263 439907 : double dfGeoX = dfMinX + iBMX * dfPixelXSize;
1264 439907 : const double dfGeoY = dfMaxY - iBMY * dfPixelYSize;
1265 :
1266 439907 : bool bMatchingGeoLocCellFound = false;
1267 :
1268 439907 : const int nOuterIters =
1269 439907 : psTransform->bGeographicSRSWithMinus180Plus180LongRange &&
1270 317622 : fabs(dfGeoX) >= 180
1271 : ? 2
1272 : : 1;
1273 :
1274 879838 : for (int iOuterIter = 0; iOuterIter < nOuterIters;
1275 : ++iOuterIter)
1276 : {
1277 439931 : if (iOuterIter == 1 && dfGeoX >= 180)
1278 0 : dfGeoX -= 360;
1279 439931 : else if (iOuterIter == 1 && dfGeoX <= -180)
1280 24 : dfGeoX += 360;
1281 :
1282 : // Identify a cell (quadrilateral in georeferenced
1283 : // space) in the geolocation array in which dfGeoX,
1284 : // dfGeoY falls into.
1285 439931 : oPoint.setX(dfGeoX);
1286 439931 : oPoint.setY(dfGeoY);
1287 439931 : const int nX = static_cast<int>(std::floor(dfX));
1288 439931 : const int nY = static_cast<int>(std::floor(dfY));
1289 1264706 : for (int sx = -1;
1290 1264706 : !bMatchingGeoLocCellFound && sx <= 0; sx++)
1291 : {
1292 2423953 : for (int sy = -1;
1293 2423953 : !bMatchingGeoLocCellFound && sy <= 0; sy++)
1294 : {
1295 1600470 : const int pixel = nX + sx;
1296 1600470 : const int line = nY + sy;
1297 : double x0, y0, x1, y1, x2, y2, x3, y3;
1298 1600470 : if (!PixelLineToXY(psTransform, pixel, line,
1299 1599902 : x0, y0) ||
1300 1599902 : !PixelLineToXY(psTransform, pixel + 1,
1301 1599682 : line, x2, y2) ||
1302 1599682 : !PixelLineToXY(psTransform, pixel,
1303 3200370 : line + 1, x1, y1) ||
1304 1599252 : !PixelLineToXY(psTransform, pixel + 1,
1305 : line + 1, x3, y3))
1306 : {
1307 1288 : break;
1308 : }
1309 :
1310 1599182 : int nIters = 1;
1311 1599182 : if (psTransform
1312 1599182 : ->bGeographicSRSWithMinus180Plus180LongRange &&
1313 1196821 : std::fabs(x0) > 170 &&
1314 13609 : std::fabs(x1) > 170 &&
1315 13359 : std::fabs(x2) > 170 &&
1316 12098 : std::fabs(x3) > 170 &&
1317 11918 : (std::fabs(x1 - x0) > 180 ||
1318 11861 : std::fabs(x2 - x0) > 180 ||
1319 10499 : std::fabs(x3 - x0) > 180))
1320 : {
1321 1490 : nIters = 2;
1322 1490 : if (x0 > 0)
1323 3 : x0 -= 360;
1324 1490 : if (x1 > 0)
1325 60 : x1 -= 360;
1326 1490 : if (x2 > 0)
1327 1422 : x2 -= 360;
1328 1490 : if (x3 > 0)
1329 1487 : x3 -= 360;
1330 : }
1331 3199850 : for (int iIter = 0; iIter < nIters; ++iIter)
1332 : {
1333 1600672 : if (iIter == 1)
1334 : {
1335 1490 : x0 += 360;
1336 1490 : x1 += 360;
1337 1490 : x2 += 360;
1338 1490 : x3 += 360;
1339 : }
1340 :
1341 1600672 : oRing.setPoint(0, x0, y0);
1342 1600672 : oRing.setPoint(1, x2, y2);
1343 1600672 : oRing.setPoint(2, x3, y3);
1344 1600672 : oRing.setPoint(3, x1, y1);
1345 1600672 : oRing.setPoint(4, x0, y0);
1346 3055670 : if (oRing.isPointInRing(&oPoint) ||
1347 1455002 : oRing.isPointOnRingBoundary(
1348 : &oPoint))
1349 : {
1350 145873 : bMatchingGeoLocCellFound = true;
1351 145873 : double dfBMXValue = pixel;
1352 145873 : double dfBMYValue = line;
1353 145873 : GDALInverseBilinearInterpolation(
1354 : dfGeoX, dfGeoY, x0, y0, x1, y1,
1355 : x2, y2, x3, y3, dfBMXValue,
1356 : dfBMYValue);
1357 :
1358 145873 : dfBMXValue =
1359 145873 : (dfBMXValue +
1360 145873 : dfGeorefConventionOffset) *
1361 145873 : psTransform->dfPIXEL_STEP +
1362 145873 : psTransform->dfPIXEL_OFFSET;
1363 145873 : dfBMYValue =
1364 145873 : (dfBMYValue +
1365 145873 : dfGeorefConventionOffset) *
1366 145873 : psTransform->dfLINE_STEP +
1367 145873 : psTransform->dfLINE_OFFSET;
1368 :
1369 145873 : pAccessors->backMapXAccessor.Set(
1370 : iBMX, iBMY,
1371 : static_cast<float>(dfBMXValue));
1372 145873 : pAccessors->backMapYAccessor.Set(
1373 : iBMX, iBMY,
1374 : static_cast<float>(dfBMYValue));
1375 145873 : pAccessors->backMapWeightAccessor
1376 : .Set(iBMX, iBMY, 1.0f);
1377 : }
1378 : }
1379 : }
1380 : }
1381 : }
1382 439907 : if (bMatchingGeoLocCellFound)
1383 145873 : continue;
1384 : }
1385 :
1386 : // We will end up here in non-nominal cases, with nodata,
1387 : // holes, etc.
1388 :
1389 : // Check if the center is in range
1390 295954 : if (iBMX < -1 || iBMY < -1 || iBMX > nBMXSize ||
1391 : iBMY > nBMYSize)
1392 478 : continue;
1393 :
1394 295476 : const double fracBMX = dBMX - iBMX;
1395 295476 : const double fracBMY = dBMY - iBMY;
1396 :
1397 : // Check logic for top left pixel
1398 295042 : if ((iBMX >= 0) && (iBMY >= 0) && (iBMX < nBMXSize) &&
1399 590518 : (iBMY < nBMYSize) &&
1400 294034 : pAccessors->backMapWeightAccessor.Get(iBMX, iBMY) !=
1401 : 1.0f)
1402 : {
1403 201764 : const double tempwt = (1.0 - fracBMX) * (1.0 - fracBMY);
1404 201764 : UpdateBackmap(iBMX, iBMY, dfX, dfY, tempwt);
1405 : }
1406 :
1407 : // Check logic for top right pixel
1408 295062 : if ((iBMY >= 0) && (iBMX + 1 < nBMXSize) &&
1409 590538 : (iBMY < nBMYSize) &&
1410 294410 : pAccessors->backMapWeightAccessor.Get(iBMX + 1, iBMY) !=
1411 : 1.0f)
1412 : {
1413 202030 : const double tempwt = fracBMX * (1.0 - fracBMY);
1414 202030 : UpdateBackmap(iBMX + 1, iBMY, dfX, dfY, tempwt);
1415 : }
1416 :
1417 : // Check logic for bottom right pixel
1418 590210 : if ((iBMX + 1 < nBMXSize) && (iBMY + 1 < nBMYSize) &&
1419 294734 : pAccessors->backMapWeightAccessor.Get(iBMX + 1,
1420 294734 : iBMY + 1) != 1.0f)
1421 : {
1422 235771 : const double tempwt = fracBMX * fracBMY;
1423 235771 : UpdateBackmap(iBMX + 1, iBMY + 1, dfX, dfY, tempwt);
1424 : }
1425 :
1426 : // Check logic for bottom left pixel
1427 295042 : if ((iBMX >= 0) && (iBMX < nBMXSize) &&
1428 884888 : (iBMY + 1 < nBMYSize) &&
1429 294370 : pAccessors->backMapWeightAccessor.Get(iBMX, iBMY + 1) !=
1430 : 1.0f)
1431 : {
1432 230605 : const double tempwt = (1.0 - fracBMX) * fracBMY;
1433 230605 : UpdateBackmap(iBMX, iBMY + 1, dfX, dfY, tempwt);
1434 : }
1435 : }
1436 : }
1437 : }
1438 : }
1439 :
1440 : // Each pixel in the backmap may have multiple entries.
1441 : // We now go in average it out using the weights
1442 163 : START_ITER_PER_BLOCK(nBMXSize, TILE_SIZE, nBMYSize, TILE_SIZE, (void)0,
1443 : iXStart, iXEnd, iYStart, iYEnd)
1444 : {
1445 2559 : for (int iY = iYStart; iY < iYEnd; ++iY)
1446 : {
1447 346177 : for (int iX = iXStart; iX < iXEnd; ++iX)
1448 : {
1449 : // Check if pixel was only touch during neighbor scan
1450 : // But no real weight was added as source point matched
1451 : // backmap grid node
1452 483432 : const auto weight =
1453 343683 : pAccessors->backMapWeightAccessor.Get(iX, iY);
1454 343683 : if (weight > 0)
1455 : {
1456 247361 : pAccessors->backMapXAccessor.Set(
1457 : iX, iY,
1458 181440 : pAccessors->backMapXAccessor.Get(iX, iY) / weight);
1459 247361 : pAccessors->backMapYAccessor.Set(
1460 : iX, iY,
1461 181440 : pAccessors->backMapYAccessor.Get(iX, iY) / weight);
1462 : }
1463 : else
1464 : {
1465 162243 : pAccessors->backMapXAccessor.Set(iX, iY, INVALID_BMXY);
1466 162243 : pAccessors->backMapYAccessor.Set(iX, iY, INVALID_BMXY);
1467 : }
1468 : }
1469 : }
1470 : }
1471 : END_ITER_PER_BLOCK
1472 :
1473 49 : pAccessors->FreeWghtsBackMap();
1474 :
1475 : // Fill holes in backmap
1476 49 : auto poBackmapDS = pAccessors->GetBackmapDataset();
1477 :
1478 49 : pAccessors->FlushBackmapCaches();
1479 :
1480 : #ifdef DEBUG_GEOLOC
1481 : if (CPLTestBool(CPLGetConfigOption("GEOLOC_DUMP", "NO")))
1482 : {
1483 : poBackmapDS->SetGeoTransform(psTransform->adfBackMapGeoTransform);
1484 : GDALClose(GDALCreateCopy(GDALGetDriverByName("GTiff"),
1485 : "/tmp/geoloc_before_fill.tif", poBackmapDS,
1486 : false, nullptr, nullptr, nullptr));
1487 : }
1488 : #endif
1489 :
1490 49 : constexpr double dfMaxSearchDist = 3.0;
1491 49 : constexpr int nSmoothingIterations = 1;
1492 147 : for (int i = 1; i <= 2; i++)
1493 : {
1494 98 : GDALFillNodata(GDALRasterBand::ToHandle(poBackmapDS->GetRasterBand(i)),
1495 : nullptr, dfMaxSearchDist,
1496 : 0, // unused parameter
1497 : nSmoothingIterations, nullptr, nullptr, nullptr);
1498 : }
1499 :
1500 : #ifdef DEBUG_GEOLOC
1501 : if (CPLTestBool(CPLGetConfigOption("GEOLOC_DUMP", "NO")))
1502 : {
1503 : GDALClose(GDALCreateCopy(GDALGetDriverByName("GTiff"),
1504 : "/tmp/geoloc_after_fill.tif", poBackmapDS,
1505 : false, nullptr, nullptr, nullptr));
1506 : }
1507 : #endif
1508 :
1509 : // A final hole filling logic, proceeding line by line, and filling
1510 : // holes when the backmap values surrounding the hole are close enough.
1511 : struct LastValidStruct
1512 : {
1513 : int iX = -1;
1514 : float bmX = 0;
1515 : };
1516 :
1517 49 : std::vector<LastValidStruct> lastValid(TILE_SIZE);
1518 245 : const auto reinitLine = [&lastValid]()
1519 : {
1520 49 : const size_t nSize = lastValid.size();
1521 49 : lastValid.clear();
1522 49 : lastValid.resize(nSize);
1523 : };
1524 163 : START_ITER_PER_BLOCK(nBMXSize, TILE_SIZE, nBMYSize, TILE_SIZE, reinitLine(),
1525 : iXStart, iXEnd, iYStart, iYEnd)
1526 : {
1527 65 : const int iYCount = iYEnd - iYStart;
1528 2559 : for (int iYIter = 0; iYIter < iYCount; ++iYIter)
1529 : {
1530 2494 : int iLastValidIX = lastValid[iYIter].iX;
1531 2494 : float bmXLastValid = lastValid[iYIter].bmX;
1532 2494 : const int iBMY = iYStart + iYIter;
1533 346177 : for (int iBMX = iXStart; iBMX < iXEnd; ++iBMX)
1534 : {
1535 343683 : const float bmX = pAccessors->backMapXAccessor.Get(iBMX, iBMY);
1536 343683 : if (bmX == INVALID_BMXY)
1537 132704 : continue;
1538 211531 : if (iLastValidIX != -1 && iBMX > iLastValidIX + 1 &&
1539 552 : fabs(bmX - bmXLastValid) <= 2)
1540 : {
1541 354 : const float bmY =
1542 239 : pAccessors->backMapYAccessor.Get(iBMX, iBMY);
1543 354 : const float bmYLastValid =
1544 239 : pAccessors->backMapYAccessor.Get(iLastValidIX, iBMY);
1545 239 : if (fabs(bmY - bmYLastValid) <= 2)
1546 : {
1547 266 : for (int iBMXInner = iLastValidIX + 1; iBMXInner < iBMX;
1548 : ++iBMXInner)
1549 : {
1550 145 : const float alpha =
1551 145 : static_cast<float>(iBMXInner - iLastValidIX) /
1552 145 : (iBMX - iLastValidIX);
1553 145 : pAccessors->backMapXAccessor.Set(
1554 : iBMXInner, iBMY,
1555 145 : (1.0f - alpha) * bmXLastValid + alpha * bmX);
1556 145 : pAccessors->backMapYAccessor.Set(
1557 : iBMXInner, iBMY,
1558 145 : (1.0f - alpha) * bmYLastValid + alpha * bmY);
1559 : }
1560 : }
1561 : }
1562 210979 : iLastValidIX = iBMX;
1563 210979 : bmXLastValid = bmX;
1564 : }
1565 2494 : lastValid[iYIter].iX = iLastValidIX;
1566 2494 : lastValid[iYIter].bmX = bmXLastValid;
1567 : }
1568 : }
1569 : END_ITER_PER_BLOCK
1570 :
1571 : #ifdef DEBUG_GEOLOC
1572 : if (CPLTestBool(CPLGetConfigOption("GEOLOC_DUMP", "NO")))
1573 : {
1574 : pAccessors->FlushBackmapCaches();
1575 :
1576 : GDALClose(GDALCreateCopy(GDALGetDriverByName("GTiff"),
1577 : "/tmp/geoloc_after_line_fill.tif", poBackmapDS,
1578 : false, nullptr, nullptr, nullptr));
1579 : }
1580 : #endif
1581 :
1582 49 : pAccessors->ReleaseBackmapDataset(poBackmapDS);
1583 49 : CPLDebug("GEOLOC", "Ending backmap generation");
1584 :
1585 49 : return true;
1586 : }
1587 :
1588 : /*! @endcond */
1589 :
1590 : /************************************************************************/
1591 : /* GDALGeoLocRescale() */
1592 : /************************************************************************/
1593 :
1594 4 : static void GDALGeoLocRescale(char **&papszMD, const char *pszItem,
1595 : double dfRatio, double dfDefaultVal)
1596 : {
1597 : const double dfVal =
1598 4 : dfRatio * CPLAtofM(CSLFetchNameValueDef(
1599 4 : papszMD, pszItem, CPLSPrintf("%.17g", dfDefaultVal)));
1600 :
1601 4 : papszMD = CSLSetNameValue(papszMD, pszItem, CPLSPrintf("%.17g", dfVal));
1602 4 : }
1603 :
1604 : /************************************************************************/
1605 : /* GDALCreateSimilarGeoLocTransformer() */
1606 : /************************************************************************/
1607 :
1608 1 : static void *GDALCreateSimilarGeoLocTransformer(void *hTransformArg,
1609 : double dfRatioX,
1610 : double dfRatioY)
1611 : {
1612 1 : VALIDATE_POINTER1(hTransformArg, "GDALCreateSimilarGeoLocTransformer",
1613 : nullptr);
1614 :
1615 1 : GDALGeoLocTransformInfo *psInfo =
1616 : static_cast<GDALGeoLocTransformInfo *>(hTransformArg);
1617 :
1618 1 : char **papszGeolocationInfo = CSLDuplicate(psInfo->papszGeolocationInfo);
1619 :
1620 1 : if (dfRatioX != 1.0 || dfRatioY != 1.0)
1621 : {
1622 1 : GDALGeoLocRescale(papszGeolocationInfo, "PIXEL_OFFSET", dfRatioX, 0.0);
1623 1 : GDALGeoLocRescale(papszGeolocationInfo, "LINE_OFFSET", dfRatioY, 0.0);
1624 1 : GDALGeoLocRescale(papszGeolocationInfo, "PIXEL_STEP", 1.0 / dfRatioX,
1625 : 1.0);
1626 1 : GDALGeoLocRescale(papszGeolocationInfo, "LINE_STEP", 1.0 / dfRatioY,
1627 : 1.0);
1628 : }
1629 :
1630 : auto psInfoNew =
1631 2 : static_cast<GDALGeoLocTransformInfo *>(GDALCreateGeoLocTransformer(
1632 1 : nullptr, papszGeolocationInfo, psInfo->bReversed));
1633 1 : psInfoNew->dfOversampleFactor = psInfo->dfOversampleFactor;
1634 :
1635 1 : CSLDestroy(papszGeolocationInfo);
1636 :
1637 1 : return psInfoNew;
1638 : }
1639 :
1640 : /************************************************************************/
1641 : /* GDALCreateGeolocationMetadata() */
1642 : /************************************************************************/
1643 :
1644 : /** Synthesize the content of a GEOLOCATION metadata domain from a
1645 : * geolocation dataset.
1646 : *
1647 : * This is used when doing gdalwarp -to GEOLOC_ARRAY=some.tif
1648 : */
1649 13 : CPLStringList GDALCreateGeolocationMetadata(GDALDatasetH hBaseDS,
1650 : const char *pszGeolocationDataset,
1651 : bool bIsSource)
1652 : {
1653 26 : CPLStringList aosMD;
1654 :
1655 : // Open geolocation dataset
1656 : auto poGeolocDS = std::unique_ptr<GDALDataset>(
1657 26 : GDALDataset::Open(pszGeolocationDataset, GDAL_OF_RASTER));
1658 13 : if (poGeolocDS == nullptr)
1659 : {
1660 2 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid dataset: %s",
1661 : pszGeolocationDataset);
1662 2 : return CPLStringList();
1663 : }
1664 11 : const int nGeoLocXSize = poGeolocDS->GetRasterXSize();
1665 11 : const int nGeoLocYSize = poGeolocDS->GetRasterYSize();
1666 11 : if (nGeoLocXSize == 0 || nGeoLocYSize == 0)
1667 : {
1668 0 : CPLError(CE_Failure, CPLE_AppDefined,
1669 : "Invalid dataset dimension for %s: %dx%d",
1670 : pszGeolocationDataset, nGeoLocXSize, nGeoLocYSize);
1671 0 : return CPLStringList();
1672 : }
1673 :
1674 : // Import the GEOLOCATION metadata from the geolocation dataset, if existing
1675 : auto papszGeolocMDFromGeolocDS =
1676 11 : poGeolocDS->GetMetadata(GDAL_MDD_GEOLOCATION);
1677 11 : if (papszGeolocMDFromGeolocDS)
1678 3 : aosMD = CSLDuplicate(papszGeolocMDFromGeolocDS);
1679 :
1680 11 : aosMD.SetNameValue("X_DATASET", pszGeolocationDataset);
1681 11 : aosMD.SetNameValue("Y_DATASET", pszGeolocationDataset);
1682 :
1683 : // Set X_BAND, Y_BAND to 1, 2 if they are not specified in the initial
1684 : // GEOLOC metadata domain.and the geolocation dataset has 2 bands.
1685 19 : if (aosMD.FetchNameValue("X_BAND") == nullptr &&
1686 8 : aosMD.FetchNameValue("Y_BAND") == nullptr)
1687 : {
1688 8 : if (poGeolocDS->GetRasterCount() != 2)
1689 : {
1690 1 : CPLError(CE_Failure, CPLE_AppDefined,
1691 : "Expected 2 bands for %s. Got %d", pszGeolocationDataset,
1692 : poGeolocDS->GetRasterCount());
1693 1 : return CPLStringList();
1694 : }
1695 7 : aosMD.SetNameValue("X_BAND", "1");
1696 7 : aosMD.SetNameValue("Y_BAND", "2");
1697 : }
1698 :
1699 : // Set the geoloc SRS from the geolocation dataset SRS if there's no
1700 : // explicit one in the initial GEOLOC metadata domain.
1701 10 : if (aosMD.FetchNameValue("SRS") == nullptr)
1702 : {
1703 10 : auto poSRS = poGeolocDS->GetSpatialRef();
1704 10 : if (poSRS)
1705 : {
1706 3 : char *pszWKT = nullptr;
1707 3 : poSRS->exportToWkt(&pszWKT);
1708 3 : aosMD.SetNameValue("SRS", pszWKT);
1709 3 : CPLFree(pszWKT);
1710 : }
1711 : }
1712 10 : if (aosMD.FetchNameValue("SRS") == nullptr)
1713 : {
1714 7 : aosMD.SetNameValue("SRS", SRS_WKT_WGS84_LAT_LONG);
1715 : }
1716 :
1717 : // Set default values for PIXEL/LINE_OFFSET/STEP if not present.
1718 10 : if (aosMD.FetchNameValue("PIXEL_OFFSET") == nullptr)
1719 7 : aosMD.SetNameValue("PIXEL_OFFSET", "0");
1720 :
1721 10 : if (aosMD.FetchNameValue("LINE_OFFSET") == nullptr)
1722 7 : aosMD.SetNameValue("LINE_OFFSET", "0");
1723 :
1724 10 : if (aosMD.FetchNameValue("PIXEL_STEP") == nullptr)
1725 : {
1726 : aosMD.SetNameValue(
1727 7 : "PIXEL_STEP", CPLSPrintf("%.17g", static_cast<double>(
1728 7 : GDALGetRasterXSize(hBaseDS)) /
1729 7 : nGeoLocXSize));
1730 : }
1731 :
1732 10 : if (aosMD.FetchNameValue("LINE_STEP") == nullptr)
1733 : {
1734 : aosMD.SetNameValue(
1735 7 : "LINE_STEP", CPLSPrintf("%.17g", static_cast<double>(
1736 7 : GDALGetRasterYSize(hBaseDS)) /
1737 7 : nGeoLocYSize));
1738 : }
1739 :
1740 10 : if (aosMD.FetchNameValue("GEOREFERENCING_CONVENTION") == nullptr)
1741 : {
1742 : const char *pszConvention =
1743 10 : poGeolocDS->GetMetadataItem("GEOREFERENCING_CONVENTION");
1744 10 : if (pszConvention)
1745 2 : aosMD.SetNameValue("GEOREFERENCING_CONVENTION", pszConvention);
1746 : }
1747 :
1748 20 : std::string osDebugMsg;
1749 10 : osDebugMsg = "Synthesized GEOLOCATION metadata for ";
1750 10 : osDebugMsg += bIsSource ? "source" : "target";
1751 10 : osDebugMsg += ":\n";
1752 102 : for (int i = 0; i < aosMD.size(); ++i)
1753 : {
1754 92 : osDebugMsg += " ";
1755 92 : osDebugMsg += aosMD[i];
1756 92 : osDebugMsg += '\n';
1757 : }
1758 10 : CPLDebug("GEOLOC", "%s", osDebugMsg.c_str());
1759 :
1760 10 : return aosMD;
1761 : }
1762 :
1763 : /************************************************************************/
1764 : /* GDALCreateGeoLocTransformer() */
1765 : /************************************************************************/
1766 :
1767 55 : void *GDALCreateGeoLocTransformerEx(GDALDatasetH hBaseDS,
1768 : CSLConstList papszGeolocationInfo,
1769 : int bReversed, const char *pszSourceDataset,
1770 : CSLConstList papszTransformOptions)
1771 :
1772 : {
1773 :
1774 55 : if (CSLFetchNameValue(papszGeolocationInfo, "PIXEL_OFFSET") == nullptr ||
1775 53 : CSLFetchNameValue(papszGeolocationInfo, "LINE_OFFSET") == nullptr ||
1776 53 : CSLFetchNameValue(papszGeolocationInfo, "PIXEL_STEP") == nullptr ||
1777 53 : CSLFetchNameValue(papszGeolocationInfo, "LINE_STEP") == nullptr ||
1778 161 : CSLFetchNameValue(papszGeolocationInfo, "X_BAND") == nullptr ||
1779 53 : CSLFetchNameValue(papszGeolocationInfo, "Y_BAND") == nullptr)
1780 : {
1781 2 : CPLError(CE_Failure, CPLE_AppDefined,
1782 : "Missing some geolocation fields in "
1783 : "GDALCreateGeoLocTransformer()");
1784 2 : return nullptr;
1785 : }
1786 :
1787 : /* -------------------------------------------------------------------- */
1788 : /* Initialize core info. */
1789 : /* -------------------------------------------------------------------- */
1790 : GDALGeoLocTransformInfo *psTransform =
1791 : static_cast<GDALGeoLocTransformInfo *>(
1792 53 : CPLCalloc(sizeof(GDALGeoLocTransformInfo), 1));
1793 :
1794 53 : psTransform->bReversed = CPL_TO_BOOL(bReversed);
1795 53 : psTransform->dfOversampleFactor = std::max(
1796 106 : 0.1,
1797 106 : std::min(2.0,
1798 53 : CPLAtof(CSLFetchNameValueDef(
1799 : papszTransformOptions, "GEOLOC_BACKMAP_OVERSAMPLE_FACTOR",
1800 : CPLGetConfigOption("GDAL_GEOLOC_BACKMAP_OVERSAMPLE_FACTOR",
1801 53 : "1.3")))));
1802 :
1803 53 : memcpy(psTransform->sTI.abySignature, GDAL_GTI2_SIGNATURE,
1804 : strlen(GDAL_GTI2_SIGNATURE));
1805 53 : psTransform->sTI.pszClassName = "GDALGeoLocTransformer";
1806 53 : psTransform->sTI.pfnTransform = GDALGeoLocTransform;
1807 53 : psTransform->sTI.pfnCleanup = GDALDestroyGeoLocTransformer;
1808 53 : psTransform->sTI.pfnSerialize = GDALSerializeGeoLocTransformer;
1809 53 : psTransform->sTI.pfnCreateSimilar = GDALCreateSimilarGeoLocTransformer;
1810 :
1811 53 : psTransform->papszGeolocationInfo = CSLDuplicate(papszGeolocationInfo);
1812 :
1813 53 : psTransform->bGeographicSRSWithMinus180Plus180LongRange =
1814 53 : CPLTestBool(CSLFetchNameValueDef(
1815 : papszTransformOptions,
1816 : "GEOLOC_NORMALIZE_LONGITUDE_MINUS_180_PLUS_180", "NO"));
1817 :
1818 : /* -------------------------------------------------------------------- */
1819 : /* Pull geolocation info from the options/metadata. */
1820 : /* -------------------------------------------------------------------- */
1821 53 : psTransform->dfPIXEL_OFFSET =
1822 53 : CPLAtof(CSLFetchNameValue(papszGeolocationInfo, "PIXEL_OFFSET"));
1823 53 : psTransform->dfLINE_OFFSET =
1824 53 : CPLAtof(CSLFetchNameValue(papszGeolocationInfo, "LINE_OFFSET"));
1825 53 : psTransform->dfPIXEL_STEP =
1826 53 : CPLAtof(CSLFetchNameValue(papszGeolocationInfo, "PIXEL_STEP"));
1827 53 : psTransform->dfLINE_STEP =
1828 53 : CPLAtof(CSLFetchNameValue(papszGeolocationInfo, "LINE_STEP"));
1829 :
1830 53 : psTransform->bOriginIsTopLeftCorner = EQUAL(
1831 : CSLFetchNameValueDef(papszGeolocationInfo, "GEOREFERENCING_CONVENTION",
1832 : "TOP_LEFT_CORNER"),
1833 : "TOP_LEFT_CORNER");
1834 :
1835 : /* -------------------------------------------------------------------- */
1836 : /* Establish access to geolocation dataset(s). */
1837 : /* -------------------------------------------------------------------- */
1838 : const char *pszDSName =
1839 53 : CSLFetchNameValue(papszGeolocationInfo, "X_DATASET");
1840 53 : if (pszDSName != nullptr)
1841 : {
1842 106 : CPLConfigOptionSetter oSetter("CPL_ALLOW_VSISTDIN", "NO", true);
1843 53 : if (CPLTestBool(CSLFetchNameValueDef(
1844 54 : papszGeolocationInfo, "X_DATASET_RELATIVE_TO_SOURCE", "NO")) &&
1845 1 : (hBaseDS != nullptr || pszSourceDataset))
1846 : {
1847 8 : const CPLString osFilename = CPLProjectRelativeFilenameSafe(
1848 7 : CPLGetDirnameSafe(pszSourceDataset
1849 : ? pszSourceDataset
1850 3 : : GDALGetDescription(hBaseDS))
1851 : .c_str(),
1852 4 : pszDSName);
1853 4 : psTransform->hDS_X =
1854 4 : GDALOpenShared(osFilename.c_str(), GA_ReadOnly);
1855 : }
1856 : else
1857 : {
1858 49 : psTransform->hDS_X = GDALOpenShared(pszDSName, GA_ReadOnly);
1859 : }
1860 : }
1861 : else
1862 : {
1863 0 : psTransform->hDS_X = hBaseDS;
1864 0 : if (hBaseDS)
1865 : {
1866 0 : GDALReferenceDataset(psTransform->hDS_X);
1867 0 : psTransform->papszGeolocationInfo =
1868 0 : CSLSetNameValue(psTransform->papszGeolocationInfo, "X_DATASET",
1869 : GDALGetDescription(hBaseDS));
1870 : }
1871 : }
1872 :
1873 53 : pszDSName = CSLFetchNameValue(papszGeolocationInfo, "Y_DATASET");
1874 53 : if (pszDSName != nullptr)
1875 : {
1876 106 : CPLConfigOptionSetter oSetter("CPL_ALLOW_VSISTDIN", "NO", true);
1877 53 : if (CPLTestBool(CSLFetchNameValueDef(
1878 54 : papszGeolocationInfo, "Y_DATASET_RELATIVE_TO_SOURCE", "NO")) &&
1879 1 : (hBaseDS != nullptr || pszSourceDataset))
1880 : {
1881 8 : const CPLString osFilename = CPLProjectRelativeFilenameSafe(
1882 7 : CPLGetDirnameSafe(pszSourceDataset
1883 : ? pszSourceDataset
1884 3 : : GDALGetDescription(hBaseDS))
1885 : .c_str(),
1886 4 : pszDSName);
1887 4 : psTransform->hDS_Y =
1888 4 : GDALOpenShared(osFilename.c_str(), GA_ReadOnly);
1889 : }
1890 : else
1891 : {
1892 49 : psTransform->hDS_Y = GDALOpenShared(pszDSName, GA_ReadOnly);
1893 : }
1894 : }
1895 : else
1896 : {
1897 0 : psTransform->hDS_Y = hBaseDS;
1898 0 : if (hBaseDS)
1899 : {
1900 0 : GDALReferenceDataset(psTransform->hDS_Y);
1901 0 : psTransform->papszGeolocationInfo =
1902 0 : CSLSetNameValue(psTransform->papszGeolocationInfo, "Y_DATASET",
1903 : GDALGetDescription(hBaseDS));
1904 : }
1905 : }
1906 :
1907 53 : if (psTransform->hDS_X == nullptr || psTransform->hDS_Y == nullptr)
1908 : {
1909 0 : GDALDestroyGeoLocTransformer(psTransform);
1910 0 : return nullptr;
1911 : }
1912 :
1913 : /* -------------------------------------------------------------------- */
1914 : /* Get the band handles. */
1915 : /* -------------------------------------------------------------------- */
1916 : const int nXBand =
1917 53 : std::max(1, atoi(CSLFetchNameValue(papszGeolocationInfo, "X_BAND")));
1918 53 : psTransform->hBand_X = GDALGetRasterBand(psTransform->hDS_X, nXBand);
1919 :
1920 53 : psTransform->dfNoDataX = GDALGetRasterNoDataValue(
1921 : psTransform->hBand_X, &(psTransform->bHasNoData));
1922 :
1923 : const int nYBand =
1924 53 : std::max(1, atoi(CSLFetchNameValue(papszGeolocationInfo, "Y_BAND")));
1925 53 : psTransform->hBand_Y = GDALGetRasterBand(psTransform->hDS_Y, nYBand);
1926 :
1927 53 : if (psTransform->hBand_X == nullptr || psTransform->hBand_Y == nullptr)
1928 : {
1929 0 : GDALDestroyGeoLocTransformer(psTransform);
1930 0 : return nullptr;
1931 : }
1932 :
1933 53 : psTransform->bSwapXY = CPLTestBool(
1934 : CSLFetchNameValueDef(papszGeolocationInfo, "SWAP_XY", "NO"));
1935 :
1936 : /* -------------------------------------------------------------------- */
1937 : /* Check that X and Y bands have the same dimensions */
1938 : /* -------------------------------------------------------------------- */
1939 53 : const int nXSize_XBand = GDALGetRasterXSize(psTransform->hDS_X);
1940 53 : const int nYSize_XBand = GDALGetRasterYSize(psTransform->hDS_X);
1941 53 : const int nXSize_YBand = GDALGetRasterXSize(psTransform->hDS_Y);
1942 53 : const int nYSize_YBand = GDALGetRasterYSize(psTransform->hDS_Y);
1943 53 : if (nYSize_XBand == 1 || nYSize_YBand == 1)
1944 : {
1945 15 : if (nYSize_XBand != 1 || nYSize_YBand != 1)
1946 : {
1947 0 : CPLError(CE_Failure, CPLE_AppDefined,
1948 : "X_BAND and Y_BAND should have both nYSize == 1");
1949 0 : GDALDestroyGeoLocTransformer(psTransform);
1950 0 : return nullptr;
1951 : }
1952 : }
1953 38 : else if (nXSize_XBand != nXSize_YBand || nYSize_XBand != nYSize_YBand)
1954 : {
1955 0 : CPLError(CE_Failure, CPLE_AppDefined,
1956 : "X_BAND and Y_BAND do not have the same dimensions");
1957 0 : GDALDestroyGeoLocTransformer(psTransform);
1958 0 : return nullptr;
1959 : }
1960 :
1961 53 : if (nXSize_XBand <= 0 || nYSize_XBand <= 0 || nXSize_YBand <= 0 ||
1962 : nYSize_YBand <= 0)
1963 : {
1964 0 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid X_BAND / Y_BAND size");
1965 0 : GDALDestroyGeoLocTransformer(psTransform);
1966 0 : return nullptr;
1967 : }
1968 :
1969 : // Is it a regular grid ? That is:
1970 : // The XBAND contains the x coordinates for all lines.
1971 : // The YBAND contains the y coordinates for all columns.
1972 53 : const bool bIsRegularGrid = (nYSize_XBand == 1 && nYSize_YBand == 1);
1973 :
1974 53 : const int nXSize = nXSize_XBand;
1975 53 : const int nYSize = bIsRegularGrid ? nXSize_YBand : nYSize_XBand;
1976 :
1977 106 : if (static_cast<size_t>(nXSize) >
1978 53 : std::numeric_limits<size_t>::max() / static_cast<size_t>(nYSize))
1979 : {
1980 0 : CPLError(CE_Failure, CPLE_AppDefined, "Int overflow : %d x %d", nXSize,
1981 : nYSize);
1982 0 : GDALDestroyGeoLocTransformer(psTransform);
1983 0 : return nullptr;
1984 : }
1985 :
1986 53 : psTransform->nGeoLocXSize = nXSize;
1987 53 : psTransform->nGeoLocYSize = nYSize;
1988 :
1989 53 : if (hBaseDS && psTransform->dfPIXEL_OFFSET == 0 &&
1990 49 : psTransform->dfLINE_OFFSET == 0 && psTransform->dfPIXEL_STEP == 1 &&
1991 41 : psTransform->dfLINE_STEP == 1)
1992 : {
1993 81 : if (GDALGetRasterXSize(hBaseDS) > nXSize ||
1994 40 : GDALGetRasterYSize(hBaseDS) > nYSize)
1995 : {
1996 2 : CPLError(CE_Warning, CPLE_AppDefined,
1997 : "Geolocation array is %d x %d large, "
1998 : "whereas dataset is %d x %d large. Result might be "
1999 : "incorrect due to lack of values in geolocation array.",
2000 : nXSize, nYSize, GDALGetRasterXSize(hBaseDS),
2001 : GDALGetRasterYSize(hBaseDS));
2002 : }
2003 : }
2004 :
2005 : /* -------------------------------------------------------------------- */
2006 : /* Load the geolocation array. */
2007 : /* -------------------------------------------------------------------- */
2008 :
2009 : // The quadtree method is experimental. It simplifies the code
2010 : // significantly, but unfortunately burns more RAM and is slower.
2011 : const bool bUseQuadtree =
2012 53 : EQUAL(CPLGetConfigOption("GDAL_GEOLOC_INVERSE_METHOD", "BACKMAP"),
2013 : "QUADTREE");
2014 :
2015 : // Decide if we should C-arrays for geoloc and backmap, or on-disk
2016 : // temporary datasets.
2017 53 : const char *pszUseTempDatasets = CSLFetchNameValueDef(
2018 : papszTransformOptions, "GEOLOC_USE_TEMP_DATASETS",
2019 : CPLGetConfigOption("GDAL_GEOLOC_USE_TEMP_DATASETS", nullptr));
2020 53 : if (pszUseTempDatasets)
2021 4 : psTransform->bUseArray = !CPLTestBool(pszUseTempDatasets);
2022 : else
2023 : {
2024 49 : constexpr int MEGAPIXEL_LIMIT = 24;
2025 49 : psTransform->bUseArray =
2026 49 : nXSize < MEGAPIXEL_LIMIT * 1000 * 1000 / nYSize;
2027 49 : if (!psTransform->bUseArray)
2028 : {
2029 0 : CPLDebug("GEOLOC",
2030 : "Using temporary GTiff backing to store backmap, because "
2031 : "geoloc arrays require %d megapixels, exceeding the %d "
2032 : "megapixels limit. You can set the "
2033 : "GDAL_GEOLOC_USE_TEMP_DATASETS configuration option to "
2034 : "NO to force RAM storage of backmap",
2035 0 : static_cast<int>(static_cast<int64_t>(nXSize) * nYSize /
2036 : (1000 * 1000)),
2037 : MEGAPIXEL_LIMIT);
2038 : }
2039 : }
2040 :
2041 53 : if (psTransform->bUseArray)
2042 : {
2043 51 : auto pAccessors = new GDALGeoLocCArrayAccessors(psTransform);
2044 51 : psTransform->pAccessors = pAccessors;
2045 51 : if (!pAccessors->Load(bIsRegularGrid, bUseQuadtree))
2046 : {
2047 0 : GDALDestroyGeoLocTransformer(psTransform);
2048 0 : return nullptr;
2049 : }
2050 : }
2051 : else
2052 : {
2053 2 : auto pAccessors = new GDALGeoLocDatasetAccessors(psTransform);
2054 2 : psTransform->pAccessors = pAccessors;
2055 2 : if (!pAccessors->Load(bIsRegularGrid, bUseQuadtree))
2056 : {
2057 0 : GDALDestroyGeoLocTransformer(psTransform);
2058 0 : return nullptr;
2059 : }
2060 : }
2061 :
2062 53 : return psTransform;
2063 : }
2064 :
2065 : /** Create GeoLocation transformer */
2066 1 : void *GDALCreateGeoLocTransformer(GDALDatasetH hBaseDS,
2067 : char **papszGeolocationInfo, int bReversed)
2068 :
2069 : {
2070 1 : return GDALCreateGeoLocTransformerEx(hBaseDS, papszGeolocationInfo,
2071 1 : bReversed, nullptr, nullptr);
2072 : }
2073 :
2074 : /************************************************************************/
2075 : /* GDALDestroyGeoLocTransformer() */
2076 : /************************************************************************/
2077 :
2078 : /** Destroy GeoLocation transformer */
2079 53 : void GDALDestroyGeoLocTransformer(void *pTransformAlg)
2080 :
2081 : {
2082 53 : if (pTransformAlg == nullptr)
2083 0 : return;
2084 :
2085 53 : GDALGeoLocTransformInfo *psTransform =
2086 : static_cast<GDALGeoLocTransformInfo *>(pTransformAlg);
2087 :
2088 53 : CSLDestroy(psTransform->papszGeolocationInfo);
2089 :
2090 53 : if (psTransform->bUseArray)
2091 : delete static_cast<GDALGeoLocCArrayAccessors *>(
2092 51 : psTransform->pAccessors);
2093 : else
2094 : delete static_cast<GDALGeoLocDatasetAccessors *>(
2095 2 : psTransform->pAccessors);
2096 :
2097 106 : if (psTransform->hDS_X != nullptr &&
2098 53 : GDALDereferenceDataset(psTransform->hDS_X) == 0)
2099 31 : GDALClose(psTransform->hDS_X);
2100 :
2101 106 : if (psTransform->hDS_Y != nullptr &&
2102 53 : GDALDereferenceDataset(psTransform->hDS_Y) == 0)
2103 45 : GDALClose(psTransform->hDS_Y);
2104 :
2105 53 : if (psTransform->hQuadTree != nullptr)
2106 4 : CPLQuadTreeDestroy(psTransform->hQuadTree);
2107 :
2108 53 : CPLFree(pTransformAlg);
2109 : }
2110 :
2111 : /** Use GeoLocation transformer */
2112 32118 : int GDALGeoLocTransform(void *pTransformArg, int bDstToSrc, int nPointCount,
2113 : double *padfX, double *padfY, double *padfZ,
2114 : int *panSuccess)
2115 : {
2116 32118 : GDALGeoLocTransformInfo *psTransform =
2117 : static_cast<GDALGeoLocTransformInfo *>(pTransformArg);
2118 32118 : if (psTransform->bUseArray)
2119 : {
2120 29995 : return GDALGeoLoc<GDALGeoLocCArrayAccessors>::Transform(
2121 : pTransformArg, bDstToSrc, nPointCount, padfX, padfY, padfZ,
2122 29995 : panSuccess);
2123 : }
2124 : else
2125 : {
2126 2123 : return GDALGeoLoc<GDALGeoLocDatasetAccessors>::Transform(
2127 : pTransformArg, bDstToSrc, nPointCount, padfX, padfY, padfZ,
2128 2123 : panSuccess);
2129 : }
2130 : }
2131 :
2132 : /************************************************************************/
2133 : /* GDALSerializeGeoLocTransformer() */
2134 : /************************************************************************/
2135 :
2136 1 : CPLXMLNode *GDALSerializeGeoLocTransformer(void *pTransformArg)
2137 :
2138 : {
2139 1 : VALIDATE_POINTER1(pTransformArg, "GDALSerializeGeoLocTransformer", nullptr);
2140 :
2141 1 : GDALGeoLocTransformInfo *psInfo =
2142 : static_cast<GDALGeoLocTransformInfo *>(pTransformArg);
2143 :
2144 : CPLXMLNode *psTree =
2145 1 : CPLCreateXMLNode(nullptr, CXT_Element, "GeoLocTransformer");
2146 :
2147 : /* -------------------------------------------------------------------- */
2148 : /* Serialize bReversed. */
2149 : /* -------------------------------------------------------------------- */
2150 1 : CPLCreateXMLElementAndValue(
2151 : psTree, "Reversed",
2152 2 : CPLString().Printf("%d", static_cast<int>(psInfo->bReversed)));
2153 :
2154 : /* -------------------------------------------------------------------- */
2155 : /* geoloc metadata. */
2156 : /* -------------------------------------------------------------------- */
2157 1 : char **papszMD = psInfo->papszGeolocationInfo;
2158 1 : CPLXMLNode *psMD = CPLCreateXMLNode(psTree, CXT_Element, "Metadata");
2159 :
2160 10 : for (int i = 0; papszMD != nullptr && papszMD[i] != nullptr; i++)
2161 : {
2162 9 : char *pszKey = nullptr;
2163 9 : const char *pszRawValue = CPLParseNameValue(papszMD[i], &pszKey);
2164 :
2165 9 : CPLXMLNode *psMDI = CPLCreateXMLNode(psMD, CXT_Element, "MDI");
2166 9 : CPLSetXMLValue(psMDI, "#key", pszKey);
2167 9 : CPLCreateXMLNode(psMDI, CXT_Text, pszRawValue);
2168 :
2169 9 : CPLFree(pszKey);
2170 : }
2171 :
2172 1 : return psTree;
2173 : }
2174 :
2175 : /************************************************************************/
2176 : /* GDALDeserializeGeoLocTransformer() */
2177 : /************************************************************************/
2178 :
2179 1 : void *GDALDeserializeGeoLocTransformer(CPLXMLNode *psTree)
2180 :
2181 : {
2182 : /* -------------------------------------------------------------------- */
2183 : /* Collect metadata. */
2184 : /* -------------------------------------------------------------------- */
2185 1 : CPLXMLNode *psMetadata = CPLGetXMLNode(psTree, "Metadata");
2186 :
2187 1 : if (psMetadata == nullptr || psMetadata->eType != CXT_Element ||
2188 1 : !EQUAL(psMetadata->pszValue, "Metadata"))
2189 0 : return nullptr;
2190 :
2191 1 : char **papszMD = nullptr;
2192 :
2193 12 : for (CPLXMLNode *psMDI = psMetadata->psChild; psMDI != nullptr;
2194 11 : psMDI = psMDI->psNext)
2195 : {
2196 11 : if (!EQUAL(psMDI->pszValue, "MDI") || psMDI->eType != CXT_Element ||
2197 11 : psMDI->psChild == nullptr || psMDI->psChild->psNext == nullptr ||
2198 11 : psMDI->psChild->eType != CXT_Attribute ||
2199 11 : psMDI->psChild->psChild == nullptr)
2200 0 : continue;
2201 :
2202 11 : papszMD = CSLSetNameValue(papszMD, psMDI->psChild->psChild->pszValue,
2203 11 : psMDI->psChild->psNext->pszValue);
2204 : }
2205 :
2206 : /* -------------------------------------------------------------------- */
2207 : /* Get other flags. */
2208 : /* -------------------------------------------------------------------- */
2209 1 : const int bReversed = atoi(CPLGetXMLValue(psTree, "Reversed", "0"));
2210 :
2211 : /* -------------------------------------------------------------------- */
2212 : /* Generate transformation. */
2213 : /* -------------------------------------------------------------------- */
2214 :
2215 : const char *pszSourceDataset =
2216 1 : CPLGetXMLValue(psTree, "SourceDataset", nullptr);
2217 :
2218 1 : void *pResult = GDALCreateGeoLocTransformerEx(nullptr, papszMD, bReversed,
2219 : pszSourceDataset, nullptr);
2220 :
2221 : /* -------------------------------------------------------------------- */
2222 : /* Cleanup GCP copy. */
2223 : /* -------------------------------------------------------------------- */
2224 1 : CSLDestroy(papszMD);
2225 :
2226 1 : return pResult;
2227 : }
|
