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