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 reused, 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] = (1 - (dfBMY - iBMY)) *
738 120417 : (double(fBMX_0_0) +
739 120417 : (dfBMX - iBMX) * double(fBMX_1_0 - fBMX_0_0)) +
740 120417 : (dfBMY - iBMY) *
741 120417 : (double(fBMX_0_1) +
742 120417 : (dfBMX - iBMX) * double(fBMX_1_1 - fBMX_0_1));
743 120417 : padfY[i] = (1 - (dfBMY - iBMY)) *
744 120417 : (double(fBMY_0_0) +
745 120417 : (dfBMX - iBMX) * double(fBMY_1_0 - fBMY_0_0)) +
746 120417 : (dfBMY - iBMY) *
747 120417 : (double(fBMY_0_1) +
748 120417 : (dfBMX - iBMX) * double(fBMY_1_1 - fBMY_0_1));
749 : }
750 2231 : else if (fBMX_1_0 != INVALID_BMXY)
751 : {
752 1686 : padfX[i] = double(fBMX_0_0) +
753 1686 : (dfBMX - iBMX) * double(fBMX_1_0 - fBMX_0_0);
754 1686 : padfY[i] = double(fBMY_0_0) +
755 1686 : (dfBMX - iBMX) * double(fBMY_1_0 - fBMY_0_0);
756 : }
757 545 : else if (fBMX_0_1 != INVALID_BMXY)
758 : {
759 410 : padfX[i] = double(fBMX_0_0) +
760 410 : (dfBMY - iBMY) * double(fBMX_0_1 - fBMX_0_0);
761 410 : padfY[i] = double(fBMY_0_0) +
762 410 : (dfBMY - iBMY) * double(fBMY_0_1 - fBMY_0_0);
763 : }
764 : else
765 : {
766 135 : padfX[i] = double(fBMX_0_0);
767 135 : padfY[i] = double(fBMY_0_0);
768 : }
769 :
770 122648 : const double dfGeoLocPixel =
771 122648 : (padfX[i] - psTransform->dfPIXEL_OFFSET) /
772 122648 : psTransform->dfPIXEL_STEP -
773 : dfGeorefConventionOffset;
774 122648 : const double dfGeoLocLine =
775 122648 : (padfY[i] - psTransform->dfLINE_OFFSET) /
776 122648 : 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 122648 : 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 122648 : oPoint.setX(dfGeoX);
807 122648 : 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 122648 : const int nSearchRadius =
815 122648 : psTransform->bGeographicSRSWithMinus180Plus180LongRange &&
816 72891 : fabs(dfGeoY) >= 85
817 : ? 5
818 : : 3;
819 122648 : const int nGeoLocPixel =
820 122648 : static_cast<int>(std::floor(dfGeoLocPixel));
821 122648 : const int nGeoLocLine = static_cast<int>(std::floor(dfGeoLocLine));
822 :
823 122648 : 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 332038 : for (int r = 0; !bDone && r <= nSearchRadius; r++)
828 : {
829 1693997 : 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 2846570 : const int sx = (r == 0) ? 0
836 1361964 : : (iter < 2 * r) ? -r + iter
837 1000168 : : (iter < 4 * r) ? r
838 656311 : : (iter < 6 * r) ? r - (iter - 4 * r)
839 : : -r;
840 2846570 : const int sy = (r == 0) ? 0
841 1361964 : : (iter < 2 * r) ? r
842 1000168 : : (iter < 4 * r) ? r - (iter - 2 * r)
843 656311 : : (iter < 6 * r) ? -r
844 321886 : : -r + (iter - 6 * r);
845 1809036 : if (nGeoLocPixel >=
846 1484612 : static_cast<int>(psTransform->nGeoLocXSize) - sx ||
847 : nGeoLocLine >=
848 1249792 : static_cast<int>(psTransform->nGeoLocYSize) - sy)
849 : {
850 324429 : continue;
851 : }
852 1160183 : const int iX = nGeoLocPixel + sx;
853 1160183 : const int iY = nGeoLocLine + sy;
854 1160183 : if (iX >= -1 || iY >= -1)
855 : {
856 : double x0, y0, x1, y1, x2, y2, x3, y3;
857 :
858 1159874 : if (!PixelLineToXY(psTransform, iX, iY, x0, y0) ||
859 1120228 : !PixelLineToXY(psTransform, iX + 1, iY, x2, y2) ||
860 3391710 : !PixelLineToXY(psTransform, iX, iY + 1, x1, y1) ||
861 1111616 : !PixelLineToXY(psTransform, iX + 1, iY + 1, x3, y3))
862 : {
863 48686 : continue;
864 : }
865 :
866 1111188 : int nIters = 1;
867 : // For a bounding box crossing the anti-meridian, check
868 : // both around -180 and +180 deg.
869 1111188 : if (psTransform
870 1111188 : ->bGeographicSRSWithMinus180Plus180LongRange &&
871 965179 : std::fabs(x0) > 170 && std::fabs(x1) > 170 &&
872 51767 : std::fabs(x2) > 170 && std::fabs(x3) > 170 &&
873 37344 : (std::fabs(x1 - x0) > 180 ||
874 37032 : std::fabs(x2 - x0) > 180 ||
875 30492 : std::fabs(x3 - x0) > 180))
876 : {
877 7144 : nIters = 2;
878 7144 : if (x0 > 0)
879 40 : x0 -= 360;
880 7144 : if (x1 > 0)
881 312 : x1 -= 360;
882 7144 : if (x2 > 0)
883 6832 : x2 -= 360;
884 7144 : if (x3 > 0)
885 7144 : x3 -= 360;
886 : }
887 2229399 : for (int iIter = 0; !bDone && iIter < nIters; ++iIter)
888 : {
889 1118214 : if (iIter == 1)
890 : {
891 7026 : x0 += 360;
892 7026 : x1 += 360;
893 7026 : x2 += 360;
894 7026 : x3 += 360;
895 : }
896 1118214 : oRing.setPoint(0, x0, y0);
897 1118214 : oRing.setPoint(1, x2, y2);
898 1118214 : oRing.setPoint(2, x3, y3);
899 1118214 : oRing.setPoint(3, x1, y1);
900 1118214 : oRing.setPoint(4, x0, y0);
901 2132786 : if (oRing.isPointInRing(&oPoint) ||
902 1014572 : oRing.isPointOnRingBoundary(&oPoint))
903 : {
904 105074 : double dfX = static_cast<double>(iX);
905 105074 : double dfY = static_cast<double>(iY);
906 105074 : GDALInverseBilinearInterpolation(
907 : dfGeoX, dfGeoY, x0, y0, x1, y1, x2, y2, x3,
908 : y3, dfX, dfY);
909 :
910 105074 : dfX = (dfX + dfGeorefConventionOffset) *
911 105074 : psTransform->dfPIXEL_STEP +
912 105074 : psTransform->dfPIXEL_OFFSET;
913 105074 : dfY = (dfY + dfGeorefConventionOffset) *
914 105074 : psTransform->dfLINE_STEP +
915 105074 : 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 105074 : padfX[i] = dfX;
925 105074 : padfY[i] = dfY;
926 :
927 105074 : bDone = true;
928 : }
929 : }
930 : }
931 : }
932 : }
933 122648 : if (!bDone)
934 : {
935 17574 : bSuccess = FALSE;
936 17574 : panSuccess[i] = FALSE;
937 17574 : padfX[i] = HUGE_VAL;
938 17574 : padfY[i] = HUGE_VAL;
939 17574 : continue;
940 : }
941 :
942 105074 : panSuccess[i] = TRUE;
943 : }
944 : }
945 :
946 31156 : 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 243336 : 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 243336 : const double A = (x0 - x) * (y0 - y2) - (y0 - y) * (x0 - x2);
971 243336 : const double B = (((x0 - x) * (y1 - y3) - (y0 - y) * (x1 - x3)) +
972 243336 : ((x1 - x) * (y0 - y2) - (y1 - y) * (x0 - x2))) /
973 : 2;
974 243336 : const double C = (x1 - x) * (y1 - y3) - (y1 - y) * (x1 - x3);
975 243336 : const double denom = A - 2 * B + C;
976 : double s;
977 243336 : const double magnitudeOfValues = fabs(A) + fabs(B) + fabs(C);
978 243336 : 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 130885 : s = A / (A - C);
983 : }
984 : else
985 : {
986 112451 : const double sqrtTerm = sqrt(B * B - A * C);
987 112451 : const double s1 = ((A - B) + sqrtTerm) / denom;
988 112451 : const double s2 = ((A - B) - sqrtTerm) / denom;
989 112451 : if (s1 < 0 || s1 > 1)
990 95924 : s = s2;
991 : else
992 16527 : s = s1;
993 : }
994 :
995 243336 : const double t_denom_x = (1 - s) * (x0 - x2) + s * (x1 - x3);
996 243336 : if (fabs(t_denom_x) > 1e-12 * magnitudeOfValues)
997 : {
998 241826 : i += ((1 - s) * (x0 - x) + s * (x1 - x)) / t_denom_x;
999 : }
1000 : else
1001 : {
1002 1510 : const double t_denom_y = (1 - s) * (y0 - y2) + s * (y1 - y3);
1003 1510 : if (fabs(t_denom_y) > 1e-12 * magnitudeOfValues)
1004 : {
1005 1510 : i += ((1 - s) * (y0 - y) + s * (y1 - y)) / t_denom_y;
1006 : }
1007 : }
1008 :
1009 243336 : j += s;
1010 243336 : }
1011 :
1012 : /************************************************************************/
1013 : /* GeoLocGenerateBackMap() */
1014 : /************************************************************************/
1015 :
1016 : /*! @cond Doxygen_Suppress */
1017 :
1018 : template <class Accessors>
1019 45 : bool GDALGeoLoc<Accessors>::GenerateBackMap(
1020 : GDALGeoLocTransformInfo *psTransform)
1021 :
1022 : {
1023 45 : CPLDebug("GEOLOC", "Starting backmap generation");
1024 45 : const int nXSize = psTransform->nGeoLocXSize;
1025 45 : 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 45 : const double dfTargetPixels =
1034 45 : static_cast<double>(nXSize) * nYSize * psTransform->dfOversampleFactor;
1035 : const double dfPixelSizeSquare =
1036 45 : sqrt((psTransform->dfMaxX - psTransform->dfMinX) *
1037 45 : (psTransform->dfMaxY - psTransform->dfMinY) / dfTargetPixels);
1038 45 : if (dfPixelSizeSquare == 0.0)
1039 : {
1040 0 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid pixel size for backmap");
1041 0 : return false;
1042 : }
1043 :
1044 45 : const double dfMinX = psTransform->dfMinX - dfPixelSizeSquare / 2.0;
1045 45 : const double dfMaxX = psTransform->dfMaxX + dfPixelSizeSquare / 2.0;
1046 45 : const double dfMaxY = psTransform->dfMaxY + dfPixelSizeSquare / 2.0;
1047 45 : const double dfMinY = psTransform->dfMinY - dfPixelSizeSquare / 2.0;
1048 45 : const double dfBMXSize = std::ceil((dfMaxX - dfMinX) / dfPixelSizeSquare);
1049 45 : 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 45 : if (!(dfBMXSize > 0 && dfBMXSize + 2 < INT_MAX) ||
1054 45 : !(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 45 : int nBMXSize = static_cast<int>(dfBMXSize);
1062 45 : int nBMYSize = static_cast<int>(dfBMYSize);
1063 :
1064 90 : if (static_cast<size_t>(1 + nBMYSize) >
1065 45 : 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 45 : const double dfPixelXSize = (dfMaxX - dfMinX) / nBMXSize;
1073 45 : const double dfPixelYSize = (dfMaxY - dfMinY) / nBMYSize;
1074 :
1075 : // Extra pixel for right-edge and bottom-edge extensions in TOP_LEFT_CORNER
1076 : // convention.
1077 45 : nBMXSize++;
1078 45 : nBMYSize++;
1079 45 : psTransform->nBackMapWidth = nBMXSize;
1080 45 : psTransform->nBackMapHeight = nBMYSize;
1081 :
1082 45 : psTransform->adfBackMapGeoTransform[0] = dfMinX;
1083 45 : psTransform->adfBackMapGeoTransform[1] = dfPixelXSize;
1084 45 : psTransform->adfBackMapGeoTransform[2] = 0.0;
1085 45 : psTransform->adfBackMapGeoTransform[3] = dfMaxY;
1086 45 : psTransform->adfBackMapGeoTransform[4] = 0.0;
1087 45 : psTransform->adfBackMapGeoTransform[5] = -dfPixelYSize;
1088 :
1089 : /* -------------------------------------------------------------------- */
1090 : /* Allocate backmap. */
1091 : /* -------------------------------------------------------------------- */
1092 45 : auto pAccessors = static_cast<Accessors *>(psTransform->pAccessors);
1093 45 : if (!pAccessors->AllocateBackMap())
1094 0 : return false;
1095 :
1096 45 : const double dfGeorefConventionOffset =
1097 45 : psTransform->bOriginIsTopLeftCorner ? 0 : 0.5;
1098 :
1099 45 : const auto UpdateBackmap =
1100 8962995 : [&](int iBMX, int iBMY, double dfX, double dfY, double tempwt)
1101 : {
1102 822655 : const auto fBMX = pAccessors->backMapXAccessor.Get(iBMX, iBMY);
1103 822655 : const auto fBMY = pAccessors->backMapYAccessor.Get(iBMX, iBMY);
1104 822655 : const float fUpdatedBMX =
1105 : fBMX +
1106 822655 : static_cast<float>(tempwt * ((dfX + dfGeorefConventionOffset) *
1107 822655 : psTransform->dfPIXEL_STEP +
1108 822655 : psTransform->dfPIXEL_OFFSET));
1109 822655 : const float fUpdatedBMY =
1110 : fBMY +
1111 822655 : static_cast<float>(tempwt * ((dfY + dfGeorefConventionOffset) *
1112 822655 : psTransform->dfLINE_STEP +
1113 822655 : psTransform->dfLINE_OFFSET));
1114 822655 : const float fUpdatedWeight =
1115 822655 : pAccessors->backMapWeightAccessor.Get(iBMX, iBMY) +
1116 822655 : 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 822655 : if (fUpdatedWeight > 0)
1123 : {
1124 822652 : const float fX = fUpdatedBMX / fUpdatedWeight;
1125 822652 : const float fY = fUpdatedBMY / fUpdatedWeight;
1126 822652 : const double dfGeoLocPixel =
1127 822652 : (double(fX) - psTransform->dfPIXEL_OFFSET) /
1128 822652 : psTransform->dfPIXEL_STEP -
1129 : dfGeorefConventionOffset;
1130 822652 : const double dfGeoLocLine =
1131 822652 : (double(fY) - psTransform->dfLINE_OFFSET) /
1132 822652 : psTransform->dfLINE_STEP -
1133 : dfGeorefConventionOffset;
1134 822652 : int iXAvg = static_cast<int>(std::max(0.0, dfGeoLocPixel));
1135 822652 : iXAvg = std::min(iXAvg, psTransform->nGeoLocXSize - 1);
1136 822652 : int iYAvg = static_cast<int>(std::max(0.0, dfGeoLocLine));
1137 822652 : iYAvg = std::min(iYAvg, psTransform->nGeoLocYSize - 1);
1138 822652 : const double dfGLX = pAccessors->geolocXAccessor.Get(iXAvg, iYAvg);
1139 822652 : const double dfGLY = pAccessors->geolocYAccessor.Get(iXAvg, iYAvg);
1140 :
1141 822652 : const unsigned iX = static_cast<unsigned>(dfX);
1142 822652 : const unsigned iY = static_cast<unsigned>(dfY);
1143 1645290 : if (!(psTransform->bHasNoData && dfGLX == psTransform->dfNoDataX) &&
1144 822642 : ((iX >= static_cast<unsigned>(nXSize - 1) ||
1145 808686 : iY >= static_cast<unsigned>(nYSize - 1)) ||
1146 798875 : (fabs(dfGLX - pAccessors->geolocXAccessor.Get(iX, iY)) <=
1147 798875 : 2 * dfPixelXSize &&
1148 756427 : fabs(dfGLY - pAccessors->geolocYAccessor.Get(iX, iY)) <=
1149 756427 : 2 * dfPixelYSize)))
1150 : {
1151 780051 : pAccessors->backMapXAccessor.Set(iBMX, iBMY, fUpdatedBMX);
1152 780051 : pAccessors->backMapYAccessor.Set(iBMX, iBMY, fUpdatedBMY);
1153 780051 : 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 90 : OGRPoint oPoint;
1162 90 : OGRLinearRing oRing;
1163 45 : 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 45 : const double dfStep = 1. / psTransform->dfOversampleFactor;
1174 :
1175 45 : constexpr int TILE_SIZE = GDALGeoLocDatasetAccessors::TILE_SIZE;
1176 45 : const int nYBlocks = DIV_ROUND_UP(nYSize, TILE_SIZE);
1177 45 : 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 90 : std::vector<std::pair<double, double>> yStartEnd(nYBlocks + 1);
1182 90 : std::vector<std::pair<double, double>> xStartEnd(nXBlocks + 1);
1183 :
1184 : {
1185 45 : int iYBlockLast = -1;
1186 45 : double dfY = -dfStep;
1187 2388 : for (; dfY <= static_cast<double>(nYSize) + 2 * dfStep; dfY += dfStep)
1188 : {
1189 2343 : const int iYBlock = static_cast<int>(dfY / TILE_SIZE);
1190 2343 : if (iYBlock != iYBlockLast)
1191 : {
1192 47 : CPLAssert(iYBlock == iYBlockLast + 1);
1193 47 : if (iYBlockLast >= 0)
1194 2 : yStartEnd[iYBlockLast].second = dfY + dfStep / 10;
1195 47 : yStartEnd[iYBlock].first = dfY;
1196 47 : iYBlockLast = iYBlock;
1197 : }
1198 : }
1199 45 : const int iYBlock = static_cast<int>(dfY / TILE_SIZE);
1200 45 : yStartEnd[iYBlock].second = dfY + dfStep / 10;
1201 : }
1202 :
1203 : {
1204 45 : int iXBlockLast = -1;
1205 45 : double dfX = -dfStep;
1206 2766 : for (; dfX <= static_cast<double>(nXSize) + 2 * dfStep; dfX += dfStep)
1207 : {
1208 2721 : const int iXBlock = static_cast<int>(dfX / TILE_SIZE);
1209 2721 : if (iXBlock != iXBlockLast)
1210 : {
1211 47 : CPLAssert(iXBlock == iXBlockLast + 1);
1212 47 : if (iXBlockLast >= 0)
1213 2 : xStartEnd[iXBlockLast].second = dfX + dfStep / 10;
1214 47 : xStartEnd[iXBlock].first = dfX;
1215 47 : iXBlockLast = iXBlock;
1216 : }
1217 : }
1218 45 : const int iXBlock = static_cast<int>(dfX / TILE_SIZE);
1219 45 : xStartEnd[iXBlock].second = dfX + dfStep / 10;
1220 : }
1221 :
1222 92 : for (int iYBlock = 0; iYBlock < nYBlocks; ++iYBlock)
1223 : {
1224 96 : 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 2655 : for (double dfY = yStartEnd[iYBlock].first;
1233 2655 : dfY < yStartEnd[iYBlock].second; dfY += dfStep)
1234 : {
1235 427339 : for (double dfX = xStartEnd[iXBlock].first;
1236 427339 : 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 424733 : if (!PixelLineToXY(psTransform, dfX, dfY, dfGeoLocX,
1243 : dfGeoLocY))
1244 145378 : continue;
1245 :
1246 : // Compute the floating point coordinates in the pixel space
1247 : // of the backmap
1248 418051 : const double dBMX = static_cast<double>(
1249 418051 : (dfGeoLocX - dfMinX) / dfPixelXSize);
1250 :
1251 418051 : const double dBMY = static_cast<double>(
1252 418051 : (dfMaxY - dfGeoLocY) / dfPixelYSize);
1253 :
1254 : // Get top left index by truncation
1255 418051 : const int iBMX = static_cast<int>(std::floor(dBMX));
1256 418051 : const int iBMY = static_cast<int>(std::floor(dBMY));
1257 :
1258 418051 : 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 416300 : double dfGeoX = dfMinX + iBMX * dfPixelXSize;
1264 416300 : const double dfGeoY = dfMaxY - iBMY * dfPixelYSize;
1265 :
1266 416300 : bool bMatchingGeoLocCellFound = false;
1267 :
1268 416300 : const int nOuterIters =
1269 416300 : psTransform->bGeographicSRSWithMinus180Plus180LongRange &&
1270 311962 : fabs(dfGeoX) >= 180
1271 : ? 2
1272 : : 1;
1273 :
1274 832624 : for (int iOuterIter = 0; iOuterIter < nOuterIters;
1275 : ++iOuterIter)
1276 : {
1277 416324 : if (iOuterIter == 1 && dfGeoX >= 180)
1278 0 : dfGeoX -= 360;
1279 416324 : 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 416324 : oPoint.setX(dfGeoX);
1286 416324 : oPoint.setY(dfGeoY);
1287 416324 : const int nX = static_cast<int>(std::floor(dfX));
1288 416324 : const int nY = static_cast<int>(std::floor(dfY));
1289 1196542 : for (int sx = -1;
1290 1196542 : !bMatchingGeoLocCellFound && sx <= 0; sx++)
1291 : {
1292 2290693 : for (int sy = -1;
1293 2290693 : !bMatchingGeoLocCellFound && sy <= 0; sy++)
1294 : {
1295 1511762 : const int pixel = nX + sx;
1296 1511762 : const int line = nY + sy;
1297 : double x0, y0, x1, y1, x2, y2, x3, y3;
1298 1511762 : if (!PixelLineToXY(psTransform, pixel, line,
1299 1511194 : x0, y0) ||
1300 1511194 : !PixelLineToXY(psTransform, pixel + 1,
1301 1510974 : line, x2, y2) ||
1302 1510974 : !PixelLineToXY(psTransform, pixel,
1303 3022950 : line + 1, x1, y1) ||
1304 1510544 : !PixelLineToXY(psTransform, pixel + 1,
1305 : line + 1, x3, y3))
1306 : {
1307 1288 : break;
1308 : }
1309 :
1310 1510474 : int nIters = 1;
1311 1510474 : if (psTransform
1312 1510474 : ->bGeographicSRSWithMinus180Plus180LongRange &&
1313 1178951 : 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 3022430 : for (int iIter = 0; iIter < nIters; ++iIter)
1332 : {
1333 1511964 : if (iIter == 1)
1334 : {
1335 1490 : x0 += 360;
1336 1490 : x1 += 360;
1337 1490 : x2 += 360;
1338 1490 : x3 += 360;
1339 : }
1340 :
1341 1511964 : oRing.setPoint(0, x0, y0);
1342 1511964 : oRing.setPoint(1, x2, y2);
1343 1511964 : oRing.setPoint(2, x3, y3);
1344 1511964 : oRing.setPoint(3, x1, y1);
1345 1511964 : oRing.setPoint(4, x0, y0);
1346 2885890 : if (oRing.isPointInRing(&oPoint) ||
1347 1373929 : oRing.isPointOnRingBoundary(
1348 : &oPoint))
1349 : {
1350 138238 : bMatchingGeoLocCellFound = true;
1351 138238 : double dfBMXValue = pixel;
1352 138238 : double dfBMYValue = line;
1353 138238 : GDALInverseBilinearInterpolation(
1354 : dfGeoX, dfGeoY, x0, y0, x1, y1,
1355 : x2, y2, x3, y3, dfBMXValue,
1356 : dfBMYValue);
1357 :
1358 138238 : dfBMXValue =
1359 138238 : (dfBMXValue +
1360 138238 : dfGeorefConventionOffset) *
1361 138238 : psTransform->dfPIXEL_STEP +
1362 138238 : psTransform->dfPIXEL_OFFSET;
1363 138238 : dfBMYValue =
1364 138238 : (dfBMYValue +
1365 138238 : dfGeorefConventionOffset) *
1366 138238 : psTransform->dfLINE_STEP +
1367 138238 : psTransform->dfLINE_OFFSET;
1368 :
1369 138238 : pAccessors->backMapXAccessor.Set(
1370 : iBMX, iBMY,
1371 : static_cast<float>(dfBMXValue));
1372 138238 : pAccessors->backMapYAccessor.Set(
1373 : iBMX, iBMY,
1374 : static_cast<float>(dfBMYValue));
1375 138238 : pAccessors->backMapWeightAccessor
1376 : .Set(iBMX, iBMY, 1.0f);
1377 : }
1378 : }
1379 : }
1380 : }
1381 : }
1382 416300 : if (bMatchingGeoLocCellFound)
1383 138238 : 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 279813 : if (iBMX < -1 || iBMY < -1 || iBMX > nBMXSize ||
1391 : iBMY > nBMYSize)
1392 458 : continue;
1393 :
1394 279355 : const double fracBMX = dBMX - iBMX;
1395 279355 : const double fracBMY = dBMY - iBMY;
1396 :
1397 : // Check logic for top left pixel
1398 279049 : if ((iBMX >= 0) && (iBMY >= 0) && (iBMX < nBMXSize) &&
1399 558404 : (iBMY < nBMYSize) &&
1400 278062 : pAccessors->backMapWeightAccessor.Get(iBMX, iBMY) !=
1401 : 1.0f)
1402 : {
1403 185849 : const double tempwt = (1.0 - fracBMX) * (1.0 - fracBMY);
1404 185849 : UpdateBackmap(iBMX, iBMY, dfX, dfY, tempwt);
1405 : }
1406 :
1407 : // Check logic for top right pixel
1408 278952 : if ((iBMY >= 0) && (iBMX + 1 < nBMXSize) &&
1409 558307 : (iBMY < nBMYSize) &&
1410 278312 : pAccessors->backMapWeightAccessor.Get(iBMX + 1, iBMY) !=
1411 : 1.0f)
1412 : {
1413 189767 : const double tempwt = fracBMX * (1.0 - fracBMY);
1414 189767 : UpdateBackmap(iBMX + 1, iBMY, dfX, dfY, tempwt);
1415 : }
1416 :
1417 : // Check logic for bottom right pixel
1418 557991 : if ((iBMX + 1 < nBMXSize) && (iBMY + 1 < nBMYSize) &&
1419 278636 : pAccessors->backMapWeightAccessor.Get(iBMX + 1,
1420 278636 : iBMY + 1) != 1.0f)
1421 : {
1422 226559 : const double tempwt = fracBMX * fracBMY;
1423 226559 : UpdateBackmap(iBMX + 1, iBMY + 1, dfX, dfY, tempwt);
1424 : }
1425 :
1426 : // Check logic for bottom left pixel
1427 279049 : if ((iBMX >= 0) && (iBMX < nBMXSize) &&
1428 836802 : (iBMY + 1 < nBMYSize) &&
1429 278398 : pAccessors->backMapWeightAccessor.Get(iBMX, iBMY + 1) !=
1430 : 1.0f)
1431 : {
1432 220480 : const double tempwt = (1.0 - fracBMX) * fracBMY;
1433 220480 : 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 151 : START_ITER_PER_BLOCK(nBMXSize, TILE_SIZE, nBMYSize, TILE_SIZE, (void)0,
1443 : iXStart, iXEnd, iYStart, iYEnd)
1444 : {
1445 2277 : for (int iY = iYStart; iY < iYEnd; ++iY)
1446 : {
1447 328049 : 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 465582 : const auto weight =
1453 325833 : pAccessors->backMapWeightAccessor.Get(iX, iY);
1454 325833 : if (weight > 0)
1455 : {
1456 236318 : pAccessors->backMapXAccessor.Set(
1457 : iX, iY,
1458 170397 : pAccessors->backMapXAccessor.Get(iX, iY) / weight);
1459 236318 : pAccessors->backMapYAccessor.Set(
1460 : iX, iY,
1461 170397 : pAccessors->backMapYAccessor.Get(iX, iY) / weight);
1462 : }
1463 : else
1464 : {
1465 155436 : pAccessors->backMapXAccessor.Set(iX, iY, INVALID_BMXY);
1466 155436 : pAccessors->backMapYAccessor.Set(iX, iY, INVALID_BMXY);
1467 : }
1468 : }
1469 : }
1470 : }
1471 : END_ITER_PER_BLOCK
1472 :
1473 45 : pAccessors->FreeWghtsBackMap();
1474 :
1475 : // Fill holes in backmap
1476 45 : auto poBackmapDS = pAccessors->GetBackmapDataset();
1477 :
1478 45 : 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 45 : constexpr double dfMaxSearchDist = 3.0;
1491 45 : constexpr int nSmoothingIterations = 1;
1492 135 : for (int i = 1; i <= 2; i++)
1493 : {
1494 90 : 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 45 : std::vector<LastValidStruct> lastValid(TILE_SIZE);
1518 225 : const auto reinitLine = [&lastValid]()
1519 : {
1520 45 : const size_t nSize = lastValid.size();
1521 45 : lastValid.clear();
1522 45 : lastValid.resize(nSize);
1523 : };
1524 151 : START_ITER_PER_BLOCK(nBMXSize, TILE_SIZE, nBMYSize, TILE_SIZE, reinitLine(),
1525 : iXStart, iXEnd, iYStart, iYEnd)
1526 : {
1527 61 : const int iYCount = iYEnd - iYStart;
1528 2277 : for (int iYIter = 0; iYIter < iYCount; ++iYIter)
1529 : {
1530 2216 : int iLastValidIX = lastValid[iYIter].iX;
1531 2216 : float bmXLastValid = lastValid[iYIter].bmX;
1532 2216 : const int iBMY = iYStart + iYIter;
1533 328049 : for (int iBMX = iXStart; iBMX < iXEnd; ++iBMX)
1534 : {
1535 325833 : const float bmX = pAccessors->backMapXAccessor.Get(iBMX, iBMY);
1536 325833 : if (bmX == INVALID_BMXY)
1537 127144 : continue;
1538 199241 : 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 198689 : iLastValidIX = iBMX;
1563 198689 : bmXLastValid = bmX;
1564 : }
1565 2216 : lastValid[iYIter].iX = iLastValidIX;
1566 2216 : 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 45 : pAccessors->ReleaseBackmapDataset(poBackmapDS);
1583 45 : CPLDebug("GEOLOC", "Ending backmap generation");
1584 :
1585 45 : 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 10 : CPLStringList GDALCreateGeolocationMetadata(GDALDatasetH hBaseDS,
1650 : const char *pszGeolocationDataset,
1651 : bool bIsSource)
1652 : {
1653 20 : CPLStringList aosMD;
1654 :
1655 : // Open geolocation dataset
1656 : auto poGeolocDS = std::unique_ptr<GDALDataset>(
1657 20 : GDALDataset::Open(pszGeolocationDataset, GDAL_OF_RASTER));
1658 10 : if (poGeolocDS == nullptr)
1659 : {
1660 2 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid dataset: %s",
1661 : pszGeolocationDataset);
1662 2 : return CPLStringList();
1663 : }
1664 8 : const int nGeoLocXSize = poGeolocDS->GetRasterXSize();
1665 8 : const int nGeoLocYSize = poGeolocDS->GetRasterYSize();
1666 8 : 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 8 : auto papszGeolocMDFromGeolocDS = poGeolocDS->GetMetadata("GEOLOCATION");
1676 8 : if (papszGeolocMDFromGeolocDS)
1677 3 : aosMD = CSLDuplicate(papszGeolocMDFromGeolocDS);
1678 :
1679 8 : aosMD.SetNameValue("X_DATASET", pszGeolocationDataset);
1680 8 : aosMD.SetNameValue("Y_DATASET", pszGeolocationDataset);
1681 :
1682 : // Set X_BAND, Y_BAND to 1, 2 if they are not specified in the initial
1683 : // GEOLOC metadata domain.and the geolocation dataset has 2 bands.
1684 13 : if (aosMD.FetchNameValue("X_BAND") == nullptr &&
1685 5 : aosMD.FetchNameValue("Y_BAND") == nullptr)
1686 : {
1687 5 : if (poGeolocDS->GetRasterCount() != 2)
1688 : {
1689 1 : CPLError(CE_Failure, CPLE_AppDefined,
1690 : "Expected 2 bands for %s. Got %d", pszGeolocationDataset,
1691 : poGeolocDS->GetRasterCount());
1692 1 : return CPLStringList();
1693 : }
1694 4 : aosMD.SetNameValue("X_BAND", "1");
1695 4 : aosMD.SetNameValue("Y_BAND", "2");
1696 : }
1697 :
1698 : // Set the geoloc SRS from the geolocation dataset SRS if there's no
1699 : // explicit one in the initial GEOLOC metadata domain.
1700 7 : if (aosMD.FetchNameValue("SRS") == nullptr)
1701 : {
1702 7 : auto poSRS = poGeolocDS->GetSpatialRef();
1703 7 : if (poSRS)
1704 : {
1705 3 : char *pszWKT = nullptr;
1706 3 : poSRS->exportToWkt(&pszWKT);
1707 3 : aosMD.SetNameValue("SRS", pszWKT);
1708 3 : CPLFree(pszWKT);
1709 : }
1710 : }
1711 7 : if (aosMD.FetchNameValue("SRS") == nullptr)
1712 : {
1713 4 : aosMD.SetNameValue("SRS", SRS_WKT_WGS84_LAT_LONG);
1714 : }
1715 :
1716 : // Set default values for PIXEL/LINE_OFFSET/STEP if not present.
1717 7 : if (aosMD.FetchNameValue("PIXEL_OFFSET") == nullptr)
1718 4 : aosMD.SetNameValue("PIXEL_OFFSET", "0");
1719 :
1720 7 : if (aosMD.FetchNameValue("LINE_OFFSET") == nullptr)
1721 4 : aosMD.SetNameValue("LINE_OFFSET", "0");
1722 :
1723 7 : if (aosMD.FetchNameValue("PIXEL_STEP") == nullptr)
1724 : {
1725 : aosMD.SetNameValue(
1726 4 : "PIXEL_STEP", CPLSPrintf("%.17g", static_cast<double>(
1727 4 : GDALGetRasterXSize(hBaseDS)) /
1728 4 : nGeoLocXSize));
1729 : }
1730 :
1731 7 : if (aosMD.FetchNameValue("LINE_STEP") == nullptr)
1732 : {
1733 : aosMD.SetNameValue(
1734 4 : "LINE_STEP", CPLSPrintf("%.17g", static_cast<double>(
1735 4 : GDALGetRasterYSize(hBaseDS)) /
1736 4 : nGeoLocYSize));
1737 : }
1738 :
1739 7 : if (aosMD.FetchNameValue("GEOREFERENCING_CONVENTION") == nullptr)
1740 : {
1741 : const char *pszConvention =
1742 7 : poGeolocDS->GetMetadataItem("GEOREFERENCING_CONVENTION");
1743 7 : if (pszConvention)
1744 2 : aosMD.SetNameValue("GEOREFERENCING_CONVENTION", pszConvention);
1745 : }
1746 :
1747 14 : std::string osDebugMsg;
1748 7 : osDebugMsg = "Synthesized GEOLOCATION metadata for ";
1749 7 : osDebugMsg += bIsSource ? "source" : "target";
1750 7 : osDebugMsg += ":\n";
1751 72 : for (int i = 0; i < aosMD.size(); ++i)
1752 : {
1753 65 : osDebugMsg += " ";
1754 65 : osDebugMsg += aosMD[i];
1755 65 : osDebugMsg += '\n';
1756 : }
1757 7 : CPLDebug("GEOLOC", "%s", osDebugMsg.c_str());
1758 :
1759 7 : return aosMD;
1760 : }
1761 :
1762 : /************************************************************************/
1763 : /* GDALCreateGeoLocTransformer() */
1764 : /************************************************************************/
1765 :
1766 51 : void *GDALCreateGeoLocTransformerEx(GDALDatasetH hBaseDS,
1767 : CSLConstList papszGeolocationInfo,
1768 : int bReversed, const char *pszSourceDataset,
1769 : CSLConstList papszTransformOptions)
1770 :
1771 : {
1772 :
1773 51 : if (CSLFetchNameValue(papszGeolocationInfo, "PIXEL_OFFSET") == nullptr ||
1774 49 : CSLFetchNameValue(papszGeolocationInfo, "LINE_OFFSET") == nullptr ||
1775 49 : CSLFetchNameValue(papszGeolocationInfo, "PIXEL_STEP") == nullptr ||
1776 49 : CSLFetchNameValue(papszGeolocationInfo, "LINE_STEP") == nullptr ||
1777 149 : CSLFetchNameValue(papszGeolocationInfo, "X_BAND") == nullptr ||
1778 49 : CSLFetchNameValue(papszGeolocationInfo, "Y_BAND") == nullptr)
1779 : {
1780 2 : CPLError(CE_Failure, CPLE_AppDefined,
1781 : "Missing some geolocation fields in "
1782 : "GDALCreateGeoLocTransformer()");
1783 2 : return nullptr;
1784 : }
1785 :
1786 : /* -------------------------------------------------------------------- */
1787 : /* Initialize core info. */
1788 : /* -------------------------------------------------------------------- */
1789 : GDALGeoLocTransformInfo *psTransform =
1790 : static_cast<GDALGeoLocTransformInfo *>(
1791 49 : CPLCalloc(sizeof(GDALGeoLocTransformInfo), 1));
1792 :
1793 49 : psTransform->bReversed = CPL_TO_BOOL(bReversed);
1794 49 : psTransform->dfOversampleFactor = std::max(
1795 98 : 0.1,
1796 98 : std::min(2.0,
1797 49 : CPLAtof(CSLFetchNameValueDef(
1798 : papszTransformOptions, "GEOLOC_BACKMAP_OVERSAMPLE_FACTOR",
1799 : CPLGetConfigOption("GDAL_GEOLOC_BACKMAP_OVERSAMPLE_FACTOR",
1800 49 : "1.3")))));
1801 :
1802 49 : memcpy(psTransform->sTI.abySignature, GDAL_GTI2_SIGNATURE,
1803 : strlen(GDAL_GTI2_SIGNATURE));
1804 49 : psTransform->sTI.pszClassName = "GDALGeoLocTransformer";
1805 49 : psTransform->sTI.pfnTransform = GDALGeoLocTransform;
1806 49 : psTransform->sTI.pfnCleanup = GDALDestroyGeoLocTransformer;
1807 49 : psTransform->sTI.pfnSerialize = GDALSerializeGeoLocTransformer;
1808 49 : psTransform->sTI.pfnCreateSimilar = GDALCreateSimilarGeoLocTransformer;
1809 :
1810 49 : psTransform->papszGeolocationInfo = CSLDuplicate(papszGeolocationInfo);
1811 :
1812 49 : psTransform->bGeographicSRSWithMinus180Plus180LongRange =
1813 49 : CPLTestBool(CSLFetchNameValueDef(
1814 : papszTransformOptions,
1815 : "GEOLOC_NORMALIZE_LONGITUDE_MINUS_180_PLUS_180", "NO"));
1816 :
1817 : /* -------------------------------------------------------------------- */
1818 : /* Pull geolocation info from the options/metadata. */
1819 : /* -------------------------------------------------------------------- */
1820 49 : psTransform->dfPIXEL_OFFSET =
1821 49 : CPLAtof(CSLFetchNameValue(papszGeolocationInfo, "PIXEL_OFFSET"));
1822 49 : psTransform->dfLINE_OFFSET =
1823 49 : CPLAtof(CSLFetchNameValue(papszGeolocationInfo, "LINE_OFFSET"));
1824 49 : psTransform->dfPIXEL_STEP =
1825 49 : CPLAtof(CSLFetchNameValue(papszGeolocationInfo, "PIXEL_STEP"));
1826 49 : psTransform->dfLINE_STEP =
1827 49 : CPLAtof(CSLFetchNameValue(papszGeolocationInfo, "LINE_STEP"));
1828 :
1829 49 : psTransform->bOriginIsTopLeftCorner = EQUAL(
1830 : CSLFetchNameValueDef(papszGeolocationInfo, "GEOREFERENCING_CONVENTION",
1831 : "TOP_LEFT_CORNER"),
1832 : "TOP_LEFT_CORNER");
1833 :
1834 : /* -------------------------------------------------------------------- */
1835 : /* Establish access to geolocation dataset(s). */
1836 : /* -------------------------------------------------------------------- */
1837 : const char *pszDSName =
1838 49 : CSLFetchNameValue(papszGeolocationInfo, "X_DATASET");
1839 49 : if (pszDSName != nullptr)
1840 : {
1841 98 : CPLConfigOptionSetter oSetter("CPL_ALLOW_VSISTDIN", "NO", true);
1842 49 : if (CPLTestBool(CSLFetchNameValueDef(
1843 50 : papszGeolocationInfo, "X_DATASET_RELATIVE_TO_SOURCE", "NO")) &&
1844 1 : (hBaseDS != nullptr || pszSourceDataset))
1845 : {
1846 8 : const CPLString osFilename = CPLProjectRelativeFilenameSafe(
1847 7 : CPLGetDirnameSafe(pszSourceDataset
1848 : ? pszSourceDataset
1849 3 : : GDALGetDescription(hBaseDS))
1850 : .c_str(),
1851 4 : pszDSName);
1852 4 : psTransform->hDS_X =
1853 4 : GDALOpenShared(osFilename.c_str(), GA_ReadOnly);
1854 : }
1855 : else
1856 : {
1857 45 : psTransform->hDS_X = GDALOpenShared(pszDSName, GA_ReadOnly);
1858 : }
1859 : }
1860 : else
1861 : {
1862 0 : psTransform->hDS_X = hBaseDS;
1863 0 : if (hBaseDS)
1864 : {
1865 0 : GDALReferenceDataset(psTransform->hDS_X);
1866 0 : psTransform->papszGeolocationInfo =
1867 0 : CSLSetNameValue(psTransform->papszGeolocationInfo, "X_DATASET",
1868 : GDALGetDescription(hBaseDS));
1869 : }
1870 : }
1871 :
1872 49 : pszDSName = CSLFetchNameValue(papszGeolocationInfo, "Y_DATASET");
1873 49 : if (pszDSName != nullptr)
1874 : {
1875 98 : CPLConfigOptionSetter oSetter("CPL_ALLOW_VSISTDIN", "NO", true);
1876 49 : if (CPLTestBool(CSLFetchNameValueDef(
1877 50 : papszGeolocationInfo, "Y_DATASET_RELATIVE_TO_SOURCE", "NO")) &&
1878 1 : (hBaseDS != nullptr || pszSourceDataset))
1879 : {
1880 8 : const CPLString osFilename = CPLProjectRelativeFilenameSafe(
1881 7 : CPLGetDirnameSafe(pszSourceDataset
1882 : ? pszSourceDataset
1883 3 : : GDALGetDescription(hBaseDS))
1884 : .c_str(),
1885 4 : pszDSName);
1886 4 : psTransform->hDS_Y =
1887 4 : GDALOpenShared(osFilename.c_str(), GA_ReadOnly);
1888 : }
1889 : else
1890 : {
1891 45 : psTransform->hDS_Y = GDALOpenShared(pszDSName, GA_ReadOnly);
1892 : }
1893 : }
1894 : else
1895 : {
1896 0 : psTransform->hDS_Y = hBaseDS;
1897 0 : if (hBaseDS)
1898 : {
1899 0 : GDALReferenceDataset(psTransform->hDS_Y);
1900 0 : psTransform->papszGeolocationInfo =
1901 0 : CSLSetNameValue(psTransform->papszGeolocationInfo, "Y_DATASET",
1902 : GDALGetDescription(hBaseDS));
1903 : }
1904 : }
1905 :
1906 49 : if (psTransform->hDS_X == nullptr || psTransform->hDS_Y == nullptr)
1907 : {
1908 0 : GDALDestroyGeoLocTransformer(psTransform);
1909 0 : return nullptr;
1910 : }
1911 :
1912 : /* -------------------------------------------------------------------- */
1913 : /* Get the band handles. */
1914 : /* -------------------------------------------------------------------- */
1915 : const int nXBand =
1916 49 : std::max(1, atoi(CSLFetchNameValue(papszGeolocationInfo, "X_BAND")));
1917 49 : psTransform->hBand_X = GDALGetRasterBand(psTransform->hDS_X, nXBand);
1918 :
1919 49 : psTransform->dfNoDataX = GDALGetRasterNoDataValue(
1920 : psTransform->hBand_X, &(psTransform->bHasNoData));
1921 :
1922 : const int nYBand =
1923 49 : std::max(1, atoi(CSLFetchNameValue(papszGeolocationInfo, "Y_BAND")));
1924 49 : psTransform->hBand_Y = GDALGetRasterBand(psTransform->hDS_Y, nYBand);
1925 :
1926 49 : if (psTransform->hBand_X == nullptr || psTransform->hBand_Y == nullptr)
1927 : {
1928 0 : GDALDestroyGeoLocTransformer(psTransform);
1929 0 : return nullptr;
1930 : }
1931 :
1932 49 : psTransform->bSwapXY = CPLTestBool(
1933 : CSLFetchNameValueDef(papszGeolocationInfo, "SWAP_XY", "NO"));
1934 :
1935 : /* -------------------------------------------------------------------- */
1936 : /* Check that X and Y bands have the same dimensions */
1937 : /* -------------------------------------------------------------------- */
1938 49 : const int nXSize_XBand = GDALGetRasterXSize(psTransform->hDS_X);
1939 49 : const int nYSize_XBand = GDALGetRasterYSize(psTransform->hDS_X);
1940 49 : const int nXSize_YBand = GDALGetRasterXSize(psTransform->hDS_Y);
1941 49 : const int nYSize_YBand = GDALGetRasterYSize(psTransform->hDS_Y);
1942 49 : if (nYSize_XBand == 1 || nYSize_YBand == 1)
1943 : {
1944 15 : if (nYSize_XBand != 1 || nYSize_YBand != 1)
1945 : {
1946 0 : CPLError(CE_Failure, CPLE_AppDefined,
1947 : "X_BAND and Y_BAND should have both nYSize == 1");
1948 0 : GDALDestroyGeoLocTransformer(psTransform);
1949 0 : return nullptr;
1950 : }
1951 : }
1952 34 : else if (nXSize_XBand != nXSize_YBand || nYSize_XBand != nYSize_YBand)
1953 : {
1954 0 : CPLError(CE_Failure, CPLE_AppDefined,
1955 : "X_BAND and Y_BAND do not have the same dimensions");
1956 0 : GDALDestroyGeoLocTransformer(psTransform);
1957 0 : return nullptr;
1958 : }
1959 :
1960 49 : if (nXSize_XBand <= 0 || nYSize_XBand <= 0 || nXSize_YBand <= 0 ||
1961 : nYSize_YBand <= 0)
1962 : {
1963 0 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid X_BAND / Y_BAND size");
1964 0 : GDALDestroyGeoLocTransformer(psTransform);
1965 0 : return nullptr;
1966 : }
1967 :
1968 : // Is it a regular grid ? That is:
1969 : // The XBAND contains the x coordinates for all lines.
1970 : // The YBAND contains the y coordinates for all columns.
1971 49 : const bool bIsRegularGrid = (nYSize_XBand == 1 && nYSize_YBand == 1);
1972 :
1973 49 : const int nXSize = nXSize_XBand;
1974 49 : const int nYSize = bIsRegularGrid ? nXSize_YBand : nYSize_XBand;
1975 :
1976 98 : if (static_cast<size_t>(nXSize) >
1977 49 : std::numeric_limits<size_t>::max() / static_cast<size_t>(nYSize))
1978 : {
1979 0 : CPLError(CE_Failure, CPLE_AppDefined, "Int overflow : %d x %d", nXSize,
1980 : nYSize);
1981 0 : GDALDestroyGeoLocTransformer(psTransform);
1982 0 : return nullptr;
1983 : }
1984 :
1985 49 : psTransform->nGeoLocXSize = nXSize;
1986 49 : psTransform->nGeoLocYSize = nYSize;
1987 :
1988 49 : if (hBaseDS && psTransform->dfPIXEL_OFFSET == 0 &&
1989 45 : psTransform->dfLINE_OFFSET == 0 && psTransform->dfPIXEL_STEP == 1 &&
1990 38 : psTransform->dfLINE_STEP == 1)
1991 : {
1992 75 : if (GDALGetRasterXSize(hBaseDS) > nXSize ||
1993 37 : GDALGetRasterYSize(hBaseDS) > nYSize)
1994 : {
1995 2 : CPLError(CE_Warning, CPLE_AppDefined,
1996 : "Geolocation array is %d x %d large, "
1997 : "whereas dataset is %d x %d large. Result might be "
1998 : "incorrect due to lack of values in geolocation array.",
1999 : nXSize, nYSize, GDALGetRasterXSize(hBaseDS),
2000 : GDALGetRasterYSize(hBaseDS));
2001 : }
2002 : }
2003 :
2004 : /* -------------------------------------------------------------------- */
2005 : /* Load the geolocation array. */
2006 : /* -------------------------------------------------------------------- */
2007 :
2008 : // The quadtree method is experimental. It simplifies the code
2009 : // significantly, but unfortunately burns more RAM and is slower.
2010 : const bool bUseQuadtree =
2011 49 : EQUAL(CPLGetConfigOption("GDAL_GEOLOC_INVERSE_METHOD", "BACKMAP"),
2012 : "QUADTREE");
2013 :
2014 : // Decide if we should C-arrays for geoloc and backmap, or on-disk
2015 : // temporary datasets.
2016 49 : const char *pszUseTempDatasets = CSLFetchNameValueDef(
2017 : papszTransformOptions, "GEOLOC_USE_TEMP_DATASETS",
2018 : CPLGetConfigOption("GDAL_GEOLOC_USE_TEMP_DATASETS", nullptr));
2019 49 : if (pszUseTempDatasets)
2020 4 : psTransform->bUseArray = !CPLTestBool(pszUseTempDatasets);
2021 : else
2022 : {
2023 45 : constexpr int MEGAPIXEL_LIMIT = 24;
2024 45 : psTransform->bUseArray =
2025 45 : nXSize < MEGAPIXEL_LIMIT * 1000 * 1000 / nYSize;
2026 45 : if (!psTransform->bUseArray)
2027 : {
2028 0 : CPLDebug("GEOLOC",
2029 : "Using temporary GTiff backing to store backmap, because "
2030 : "geoloc arrays require %d megapixels, exceeding the %d "
2031 : "megapixels limit. You can set the "
2032 : "GDAL_GEOLOC_USE_TEMP_DATASETS configuration option to "
2033 : "NO to force RAM storage of backmap",
2034 0 : static_cast<int>(static_cast<int64_t>(nXSize) * nYSize /
2035 : (1000 * 1000)),
2036 : MEGAPIXEL_LIMIT);
2037 : }
2038 : }
2039 :
2040 49 : if (psTransform->bUseArray)
2041 : {
2042 47 : auto pAccessors = new GDALGeoLocCArrayAccessors(psTransform);
2043 47 : psTransform->pAccessors = pAccessors;
2044 47 : if (!pAccessors->Load(bIsRegularGrid, bUseQuadtree))
2045 : {
2046 0 : GDALDestroyGeoLocTransformer(psTransform);
2047 0 : return nullptr;
2048 : }
2049 : }
2050 : else
2051 : {
2052 2 : auto pAccessors = new GDALGeoLocDatasetAccessors(psTransform);
2053 2 : psTransform->pAccessors = pAccessors;
2054 2 : if (!pAccessors->Load(bIsRegularGrid, bUseQuadtree))
2055 : {
2056 0 : GDALDestroyGeoLocTransformer(psTransform);
2057 0 : return nullptr;
2058 : }
2059 : }
2060 :
2061 49 : return psTransform;
2062 : }
2063 :
2064 : /** Create GeoLocation transformer */
2065 1 : void *GDALCreateGeoLocTransformer(GDALDatasetH hBaseDS,
2066 : char **papszGeolocationInfo, int bReversed)
2067 :
2068 : {
2069 1 : return GDALCreateGeoLocTransformerEx(hBaseDS, papszGeolocationInfo,
2070 1 : bReversed, nullptr, nullptr);
2071 : }
2072 :
2073 : /************************************************************************/
2074 : /* GDALDestroyGeoLocTransformer() */
2075 : /************************************************************************/
2076 :
2077 : /** Destroy GeoLocation transformer */
2078 49 : void GDALDestroyGeoLocTransformer(void *pTransformAlg)
2079 :
2080 : {
2081 49 : if (pTransformAlg == nullptr)
2082 0 : return;
2083 :
2084 49 : GDALGeoLocTransformInfo *psTransform =
2085 : static_cast<GDALGeoLocTransformInfo *>(pTransformAlg);
2086 :
2087 49 : CSLDestroy(psTransform->papszGeolocationInfo);
2088 :
2089 49 : if (psTransform->bUseArray)
2090 : delete static_cast<GDALGeoLocCArrayAccessors *>(
2091 47 : psTransform->pAccessors);
2092 : else
2093 : delete static_cast<GDALGeoLocDatasetAccessors *>(
2094 2 : psTransform->pAccessors);
2095 :
2096 98 : if (psTransform->hDS_X != nullptr &&
2097 49 : GDALDereferenceDataset(psTransform->hDS_X) == 0)
2098 28 : GDALClose(psTransform->hDS_X);
2099 :
2100 98 : if (psTransform->hDS_Y != nullptr &&
2101 49 : GDALDereferenceDataset(psTransform->hDS_Y) == 0)
2102 41 : GDALClose(psTransform->hDS_Y);
2103 :
2104 49 : if (psTransform->hQuadTree != nullptr)
2105 4 : CPLQuadTreeDestroy(psTransform->hQuadTree);
2106 :
2107 49 : CPLFree(pTransformAlg);
2108 : }
2109 :
2110 : /** Use GeoLocation transformer */
2111 31180 : int GDALGeoLocTransform(void *pTransformArg, int bDstToSrc, int nPointCount,
2112 : double *padfX, double *padfY, double *padfZ,
2113 : int *panSuccess)
2114 : {
2115 31180 : GDALGeoLocTransformInfo *psTransform =
2116 : static_cast<GDALGeoLocTransformInfo *>(pTransformArg);
2117 31180 : if (psTransform->bUseArray)
2118 : {
2119 29059 : return GDALGeoLoc<GDALGeoLocCArrayAccessors>::Transform(
2120 : pTransformArg, bDstToSrc, nPointCount, padfX, padfY, padfZ,
2121 29059 : panSuccess);
2122 : }
2123 : else
2124 : {
2125 2121 : return GDALGeoLoc<GDALGeoLocDatasetAccessors>::Transform(
2126 : pTransformArg, bDstToSrc, nPointCount, padfX, padfY, padfZ,
2127 2121 : panSuccess);
2128 : }
2129 : }
2130 :
2131 : /************************************************************************/
2132 : /* GDALSerializeGeoLocTransformer() */
2133 : /************************************************************************/
2134 :
2135 1 : CPLXMLNode *GDALSerializeGeoLocTransformer(void *pTransformArg)
2136 :
2137 : {
2138 1 : VALIDATE_POINTER1(pTransformArg, "GDALSerializeGeoLocTransformer", nullptr);
2139 :
2140 1 : GDALGeoLocTransformInfo *psInfo =
2141 : static_cast<GDALGeoLocTransformInfo *>(pTransformArg);
2142 :
2143 : CPLXMLNode *psTree =
2144 1 : CPLCreateXMLNode(nullptr, CXT_Element, "GeoLocTransformer");
2145 :
2146 : /* -------------------------------------------------------------------- */
2147 : /* Serialize bReversed. */
2148 : /* -------------------------------------------------------------------- */
2149 1 : CPLCreateXMLElementAndValue(
2150 : psTree, "Reversed",
2151 2 : CPLString().Printf("%d", static_cast<int>(psInfo->bReversed)));
2152 :
2153 : /* -------------------------------------------------------------------- */
2154 : /* geoloc metadata. */
2155 : /* -------------------------------------------------------------------- */
2156 1 : char **papszMD = psInfo->papszGeolocationInfo;
2157 1 : CPLXMLNode *psMD = CPLCreateXMLNode(psTree, CXT_Element, "Metadata");
2158 :
2159 10 : for (int i = 0; papszMD != nullptr && papszMD[i] != nullptr; i++)
2160 : {
2161 9 : char *pszKey = nullptr;
2162 9 : const char *pszRawValue = CPLParseNameValue(papszMD[i], &pszKey);
2163 :
2164 9 : CPLXMLNode *psMDI = CPLCreateXMLNode(psMD, CXT_Element, "MDI");
2165 9 : CPLSetXMLValue(psMDI, "#key", pszKey);
2166 9 : CPLCreateXMLNode(psMDI, CXT_Text, pszRawValue);
2167 :
2168 9 : CPLFree(pszKey);
2169 : }
2170 :
2171 1 : return psTree;
2172 : }
2173 :
2174 : /************************************************************************/
2175 : /* GDALDeserializeGeoLocTransformer() */
2176 : /************************************************************************/
2177 :
2178 1 : void *GDALDeserializeGeoLocTransformer(CPLXMLNode *psTree)
2179 :
2180 : {
2181 : /* -------------------------------------------------------------------- */
2182 : /* Collect metadata. */
2183 : /* -------------------------------------------------------------------- */
2184 1 : CPLXMLNode *psMetadata = CPLGetXMLNode(psTree, "Metadata");
2185 :
2186 1 : if (psMetadata == nullptr || psMetadata->eType != CXT_Element ||
2187 1 : !EQUAL(psMetadata->pszValue, "Metadata"))
2188 0 : return nullptr;
2189 :
2190 1 : char **papszMD = nullptr;
2191 :
2192 12 : for (CPLXMLNode *psMDI = psMetadata->psChild; psMDI != nullptr;
2193 11 : psMDI = psMDI->psNext)
2194 : {
2195 11 : if (!EQUAL(psMDI->pszValue, "MDI") || psMDI->eType != CXT_Element ||
2196 11 : psMDI->psChild == nullptr || psMDI->psChild->psNext == nullptr ||
2197 11 : psMDI->psChild->eType != CXT_Attribute ||
2198 11 : psMDI->psChild->psChild == nullptr)
2199 0 : continue;
2200 :
2201 11 : papszMD = CSLSetNameValue(papszMD, psMDI->psChild->psChild->pszValue,
2202 11 : psMDI->psChild->psNext->pszValue);
2203 : }
2204 :
2205 : /* -------------------------------------------------------------------- */
2206 : /* Get other flags. */
2207 : /* -------------------------------------------------------------------- */
2208 1 : const int bReversed = atoi(CPLGetXMLValue(psTree, "Reversed", "0"));
2209 :
2210 : /* -------------------------------------------------------------------- */
2211 : /* Generate transformation. */
2212 : /* -------------------------------------------------------------------- */
2213 :
2214 : const char *pszSourceDataset =
2215 1 : CPLGetXMLValue(psTree, "SourceDataset", nullptr);
2216 :
2217 1 : void *pResult = GDALCreateGeoLocTransformerEx(nullptr, papszMD, bReversed,
2218 : pszSourceDataset, nullptr);
2219 :
2220 : /* -------------------------------------------------------------------- */
2221 : /* Cleanup GCP copy. */
2222 : /* -------------------------------------------------------------------- */
2223 1 : CSLDestroy(papszMD);
2224 :
2225 1 : return pResult;
2226 : }
|
