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 11409700 : static double UnshiftGeoX(const GDALGeoLocTransformInfo *psTransform,
73 : double dfX)
74 : {
75 11409700 : if (!psTransform->bGeographicSRSWithMinus180Plus180LongRange)
76 2267220 : return dfX;
77 9142510 : if (dfX > 180 || dfX < -180)
78 : {
79 35493 : dfX = fmod(dfX + 180.0, 360.0);
80 35493 : if (dfX < 0)
81 372 : dfX += 180.0;
82 : else
83 35121 : dfX -= 180.0;
84 35493 : return dfX;
85 : }
86 9107010 : return dfX;
87 : }
88 :
89 : /************************************************************************/
90 : /* UpdateMinMax() */
91 : /************************************************************************/
92 :
93 242727 : inline void UpdateMinMax(GDALGeoLocTransformInfo *psTransform, double dfGeoLocX,
94 : double dfGeoLocY)
95 : {
96 242727 : if (dfGeoLocX < psTransform->dfMinX)
97 : {
98 1225 : psTransform->dfMinX = dfGeoLocX;
99 1225 : psTransform->dfYAtMinX = dfGeoLocY;
100 : }
101 242727 : if (dfGeoLocX > psTransform->dfMaxX)
102 : {
103 871 : psTransform->dfMaxX = dfGeoLocX;
104 871 : psTransform->dfYAtMaxX = dfGeoLocY;
105 : }
106 242727 : if (dfGeoLocY < psTransform->dfMinY)
107 : {
108 1329 : psTransform->dfMinY = dfGeoLocY;
109 1329 : psTransform->dfXAtMinY = dfGeoLocX;
110 : }
111 242727 : if (dfGeoLocY > psTransform->dfMaxY)
112 : {
113 3413 : psTransform->dfMaxY = dfGeoLocY;
114 3413 : psTransform->dfXAtMaxY = dfGeoLocX;
115 : }
116 242727 : }
117 :
118 : /************************************************************************/
119 : /* Clamp() */
120 : /************************************************************************/
121 :
122 3100 : inline double Clamp(double v, double minV, double maxV)
123 : {
124 3100 : 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 50 : void GDALGeoLoc<Accessors>::LoadGeolocFinish(
164 : GDALGeoLocTransformInfo *psTransform)
165 : {
166 50 : auto pAccessors = static_cast<Accessors *>(psTransform->pAccessors);
167 50 : CSLConstList papszGeolocationInfo = psTransform->papszGeolocationInfo;
168 :
169 : /* -------------------------------------------------------------------- */
170 : /* Scan forward map for lat/long extents. */
171 : /* -------------------------------------------------------------------- */
172 50 : psTransform->dfMinX = std::numeric_limits<double>::max();
173 50 : psTransform->dfMaxX = -std::numeric_limits<double>::max();
174 50 : psTransform->dfMinY = std::numeric_limits<double>::max();
175 50 : psTransform->dfMaxY = -std::numeric_limits<double>::max();
176 :
177 50 : constexpr int TILE_SIZE = GDALGeoLocDatasetAccessors::TILE_SIZE;
178 156 : START_ITER_PER_BLOCK(psTransform->nGeoLocXSize, TILE_SIZE,
179 : psTransform->nGeoLocYSize, TILE_SIZE, (void)0, iXStart,
180 : iXEnd, iYStart, iYEnd)
181 : {
182 1984 : for (int iY = iYStart; iY < iYEnd; ++iY)
183 : {
184 243439 : for (int iX = iXStart; iX < iXEnd; ++iX)
185 : {
186 241509 : double dfX = pAccessors->geolocXAccessor.Get(iX, iY);
187 241509 : if (!psTransform->bHasNoData || dfX != psTransform->dfNoDataX)
188 : {
189 238567 : dfX = UnshiftGeoX(psTransform, dfX);
190 238567 : 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 50 : const char *pszSRS = CSLFetchNameValue(papszGeolocationInfo, "SRS");
201 50 : if (!psTransform->bGeographicSRSWithMinus180Plus180LongRange && pszSRS &&
202 47 : psTransform->dfMinX >= -180.0 && psTransform->dfMaxX <= 180.0 &&
203 46 : !psTransform->bSwapXY)
204 : {
205 45 : OGRSpatialReference oSRS;
206 45 : psTransform->bGeographicSRSWithMinus180Plus180LongRange =
207 90 : oSRS.importFromWkt(pszSRS) == OGRERR_NONE &&
208 45 : 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 50 : if (psTransform->bOriginIsTopLeftCorner)
285 : {
286 : // Add "virtual" edge at Y=nGeoLocYSize
287 1881 : for (int iX = 0; iX <= psTransform->nGeoLocXSize; iX++)
288 : {
289 : double dfGeoLocX;
290 : double dfGeoLocY;
291 1853 : if (!PixelLineToXY(psTransform, static_cast<double>(iX),
292 1853 : static_cast<double>(psTransform->nGeoLocYSize),
293 : dfGeoLocX, dfGeoLocY))
294 0 : continue;
295 1853 : if (psTransform->bGeographicSRSWithMinus180Plus180LongRange)
296 1089 : dfGeoLocX = Clamp(dfGeoLocX, -180.0, 180.0);
297 1853 : UpdateMinMax(psTransform, dfGeoLocX, dfGeoLocY);
298 : }
299 :
300 : // Add "virtual" edge at X=nGeoLocXSize
301 1453 : for (int iY = 0; iY <= psTransform->nGeoLocYSize; iY++)
302 : {
303 : double dfGeoLocX;
304 : double dfGeoLocY;
305 1425 : if (!PixelLineToXY(psTransform,
306 1425 : static_cast<double>(psTransform->nGeoLocXSize),
307 : static_cast<double>(iY), dfGeoLocX, dfGeoLocY))
308 0 : continue;
309 1425 : if (psTransform->bGeographicSRSWithMinus180Plus180LongRange)
310 1221 : dfGeoLocX = Clamp(dfGeoLocX, -180.0, 180.0);
311 1425 : 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 50 : }
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 1323437 : bool GDALGeoLoc<Accessors>::PixelLineToXY(
393 : const GDALGeoLocTransformInfo *psTransform, const double dfGeoLocPixel,
394 : const double dfGeoLocLine, double &dfX, double &dfY)
395 : {
396 1323437 : int iX = static_cast<int>(
397 2646878 : std::min(std::max(0.0, dfGeoLocPixel),
398 1323437 : static_cast<double>(psTransform->nGeoLocXSize - 1)));
399 1323437 : int iY = static_cast<int>(
400 2646878 : std::min(std::max(0.0, dfGeoLocLine),
401 1323437 : static_cast<double>(psTransform->nGeoLocYSize - 1)));
402 :
403 1323437 : auto pAccessors = static_cast<Accessors *>(psTransform->pAccessors);
404 :
405 1663901 : for (int iAttempt = 0; iAttempt < 2; ++iAttempt)
406 : {
407 1663901 : const double dfGLX_0_0 = pAccessors->geolocXAccessor.Get(iX, iY);
408 1663901 : const double dfGLY_0_0 = pAccessors->geolocYAccessor.Get(iX, iY);
409 1663901 : if (psTransform->bHasNoData && dfGLX_0_0 == psTransform->dfNoDataX)
410 : {
411 13386 : return false;
412 : }
413 :
414 : // This assumes infinite extension beyond borders of available
415 : // data based on closest grid square.
416 1650511 : if (iX + 1 < psTransform->nGeoLocXSize &&
417 1452173 : iY + 1 < psTransform->nGeoLocYSize)
418 : {
419 1310051 : const double dfGLX_1_0 =
420 : pAccessors->geolocXAccessor.Get(iX + 1, iY);
421 1310051 : const double dfGLY_1_0 =
422 : pAccessors->geolocYAccessor.Get(iX + 1, iY);
423 1310051 : const double dfGLX_0_1 =
424 : pAccessors->geolocXAccessor.Get(iX, iY + 1);
425 1310051 : const double dfGLY_0_1 =
426 : pAccessors->geolocYAccessor.Get(iX, iY + 1);
427 1310051 : const double dfGLX_1_1 =
428 : pAccessors->geolocXAccessor.Get(iX + 1, iY + 1);
429 1310051 : const double dfGLY_1_1 =
430 : pAccessors->geolocYAccessor.Get(iX + 1, iY + 1);
431 1310051 : if (!psTransform->bHasNoData ||
432 132991 : (dfGLX_1_0 != psTransform->dfNoDataX &&
433 131453 : dfGLX_0_1 != psTransform->dfNoDataX &&
434 130303 : dfGLX_1_1 != psTransform->dfNoDataX))
435 : {
436 : const double dfGLX_1_0_adjusted =
437 1307171 : ShiftGeoX(psTransform, dfGLX_0_0, dfGLX_1_0);
438 : const double dfGLX_0_1_adjusted =
439 1307171 : ShiftGeoX(psTransform, dfGLX_0_0, dfGLX_0_1);
440 : const double dfGLX_1_1_adjusted =
441 1307171 : ShiftGeoX(psTransform, dfGLX_0_0, dfGLX_1_1);
442 1307171 : dfX = (1 - (dfGeoLocLine - iY)) *
443 1307171 : (dfGLX_0_0 + (dfGeoLocPixel - iX) *
444 1307171 : (dfGLX_1_0_adjusted - dfGLX_0_0)) +
445 1307171 : (dfGeoLocLine - iY) *
446 1307171 : (dfGLX_0_1_adjusted +
447 1307171 : (dfGeoLocPixel - iX) *
448 1307171 : (dfGLX_1_1_adjusted - dfGLX_0_1_adjusted));
449 1307171 : dfX = UnshiftGeoX(psTransform, dfX);
450 :
451 1307171 : dfY = (1 - (dfGeoLocLine - iY)) *
452 1307171 : (dfGLY_0_0 +
453 1307171 : (dfGeoLocPixel - iX) * (dfGLY_1_0 - dfGLY_0_0)) +
454 1307171 : (dfGeoLocLine - iY) *
455 1307171 : (dfGLY_0_1 +
456 1307171 : (dfGeoLocPixel - iX) * (dfGLY_1_1 - dfGLY_0_1));
457 1307171 : break;
458 : }
459 : }
460 :
461 343341 : if (iX == psTransform->nGeoLocXSize - 1 && iX >= 1 &&
462 198339 : iY + 1 < psTransform->nGeoLocYSize)
463 : {
464 : // If we are after the right edge, then go one pixel left
465 : // and retry
466 183917 : iX--;
467 183917 : continue;
468 : }
469 159424 : else if (iY == psTransform->nGeoLocYSize - 1 && iY >= 1 &&
470 156544 : iX + 1 < psTransform->nGeoLocXSize)
471 : {
472 : // If we are after the bottom edge, then go one pixel up
473 : // and retry
474 142122 : iY--;
475 142122 : continue;
476 : }
477 17302 : else if (iX == psTransform->nGeoLocXSize - 1 &&
478 14422 : 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 14422 : iX--;
483 14422 : iY--;
484 14422 : continue;
485 : }
486 5760 : else if (iX + 1 < psTransform->nGeoLocXSize &&
487 2880 : (!psTransform->bHasNoData ||
488 2880 : pAccessors->geolocXAccessor.Get(iX + 1, iY) !=
489 2880 : psTransform->dfNoDataX))
490 : {
491 1342 : const double dfGLX_1_0 =
492 : pAccessors->geolocXAccessor.Get(iX + 1, iY);
493 1342 : const double dfGLY_1_0 =
494 : pAccessors->geolocYAccessor.Get(iX + 1, iY);
495 1342 : dfX =
496 1342 : dfGLX_0_0 +
497 2684 : (dfGeoLocPixel - iX) *
498 1342 : (ShiftGeoX(psTransform, dfGLX_0_0, dfGLX_1_0) - dfGLX_0_0);
499 1342 : dfX = UnshiftGeoX(psTransform, dfX);
500 1342 : dfY = dfGLY_0_0 + (dfGeoLocPixel - iX) * (dfGLY_1_0 - dfGLY_0_0);
501 : }
502 3076 : else if (iY + 1 < psTransform->nGeoLocYSize &&
503 1538 : (!psTransform->bHasNoData ||
504 1538 : pAccessors->geolocXAccessor.Get(iX, iY + 1) !=
505 1538 : psTransform->dfNoDataX))
506 : {
507 1486 : const double dfGLX_0_1 =
508 : pAccessors->geolocXAccessor.Get(iX, iY + 1);
509 1486 : const double dfGLY_0_1 =
510 : pAccessors->geolocYAccessor.Get(iX, iY + 1);
511 1486 : dfX =
512 1486 : dfGLX_0_0 +
513 2972 : (dfGeoLocLine - iY) *
514 1486 : (ShiftGeoX(psTransform, dfGLX_0_0, dfGLX_0_1) - dfGLX_0_0);
515 1486 : dfX = UnshiftGeoX(psTransform, dfX);
516 1486 : dfY = dfGLY_0_0 + (dfGeoLocLine - iY) * (dfGLY_0_1 - dfGLY_0_0);
517 : }
518 : else
519 : {
520 52 : dfX = UnshiftGeoX(psTransform, dfGLX_0_0);
521 52 : dfY = dfGLY_0_0;
522 : }
523 2880 : break;
524 : }
525 1310051 : return true;
526 : }
527 :
528 : template <class Accessors>
529 10756830 : bool GDALGeoLoc<Accessors>::PixelLineToXY(
530 : const GDALGeoLocTransformInfo *psTransform, const int nGeoLocPixel,
531 : const int nGeoLocLine, double &dfX, double &dfY)
532 : {
533 10756830 : if (nGeoLocPixel >= 0 && nGeoLocPixel < psTransform->nGeoLocXSize &&
534 10034460 : nGeoLocLine >= 0 && nGeoLocLine < psTransform->nGeoLocYSize)
535 : {
536 9911850 : auto pAccessors = static_cast<Accessors *>(psTransform->pAccessors);
537 9911850 : const double dfGLX =
538 : pAccessors->geolocXAccessor.Get(nGeoLocPixel, nGeoLocLine);
539 9911850 : const double dfGLY =
540 : pAccessors->geolocYAccessor.Get(nGeoLocPixel, nGeoLocLine);
541 9911850 : if (psTransform->bHasNoData && dfGLX == psTransform->dfNoDataX)
542 : {
543 50742 : return false;
544 : }
545 9861110 : dfX = UnshiftGeoX(psTransform, dfGLX);
546 9861110 : dfY = dfGLY;
547 9861110 : return true;
548 : }
549 844978 : return PixelLineToXY(psTransform, static_cast<double>(nGeoLocPixel),
550 844978 : 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 31483 : int GDALGeoLoc<Accessors>::Transform(void *pTransformArg, int bDstToSrc,
594 : int nPointCount, double *padfX,
595 : double *padfY, double * /* padfZ */,
596 : int *panSuccess)
597 : {
598 31483 : int bSuccess = TRUE;
599 31483 : GDALGeoLocTransformInfo *psTransform =
600 : static_cast<GDALGeoLocTransformInfo *>(pTransformArg);
601 :
602 31483 : if (psTransform->bReversed)
603 0 : bDstToSrc = !bDstToSrc;
604 :
605 31483 : const double dfGeorefConventionOffset =
606 31483 : psTransform->bOriginIsTopLeftCorner ? 0 : 0.5;
607 :
608 : /* -------------------------------------------------------------------- */
609 : /* Do original pixel line to target geox/geoy. */
610 : /* -------------------------------------------------------------------- */
611 31483 : if (!bDstToSrc)
612 : {
613 72672 : for (int i = 0; i < nPointCount; i++)
614 : {
615 47754 : if (padfX[i] == HUGE_VAL || padfY[i] == HUGE_VAL)
616 : {
617 2786 : bSuccess = FALSE;
618 2786 : panSuccess[i] = FALSE;
619 2786 : continue;
620 : }
621 :
622 44968 : const double dfGeoLocPixel =
623 44968 : (padfX[i] - psTransform->dfPIXEL_OFFSET) /
624 44968 : psTransform->dfPIXEL_STEP -
625 : dfGeorefConventionOffset;
626 44968 : const double dfGeoLocLine =
627 44968 : (padfY[i] - psTransform->dfLINE_OFFSET) /
628 44968 : psTransform->dfLINE_STEP -
629 : dfGeorefConventionOffset;
630 :
631 44968 : if (!PixelLineToXY(psTransform, dfGeoLocPixel, dfGeoLocLine,
632 44968 : 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 42142 : if (psTransform->bSwapXY)
642 : {
643 568 : std::swap(padfX[i], padfY[i]);
644 : }
645 :
646 42142 : panSuccess[i] = TRUE;
647 : }
648 : }
649 :
650 : /* -------------------------------------------------------------------- */
651 : /* geox/geoy to pixel/line using backmap. */
652 : /* -------------------------------------------------------------------- */
653 : else
654 : {
655 6565 : if (psTransform->hQuadTree)
656 : {
657 24 : GDALGeoLocInverseTransformQuadtree(psTransform, nPointCount, padfX,
658 : padfY, panSuccess);
659 24 : return TRUE;
660 : }
661 :
662 6541 : 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 13082 : OGRPoint oPoint;
668 13082 : OGRLinearRing oRing;
669 6541 : oRing.setNumPoints(5);
670 :
671 6541 : auto pAccessors = static_cast<Accessors *>(psTransform->pAccessors);
672 :
673 224618 : for (int i = 0; i < nPointCount; i++)
674 : {
675 218077 : 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 218077 : if (psTransform->bSwapXY)
683 : {
684 964 : std::swap(padfX[i], padfY[i]);
685 : }
686 :
687 218077 : const double dfGeoX = padfX[i];
688 218077 : const double dfGeoY = padfY[i];
689 :
690 218077 : const double dfBMX =
691 218077 : ((padfX[i] - psTransform->adfBackMapGeoTransform[0]) /
692 : psTransform->adfBackMapGeoTransform[1]);
693 218077 : const double dfBMY =
694 218077 : ((padfY[i] - psTransform->adfBackMapGeoTransform[3]) /
695 : psTransform->adfBackMapGeoTransform[5]);
696 :
697 218077 : if (!(dfBMX >= 0 && dfBMY >= 0 &&
698 217912 : dfBMX + 1 < psTransform->nBackMapWidth &&
699 216613 : dfBMY + 1 < psTransform->nBackMapHeight))
700 : {
701 2378 : bSuccess = FALSE;
702 2378 : panSuccess[i] = FALSE;
703 2378 : padfX[i] = HUGE_VAL;
704 2378 : padfY[i] = HUGE_VAL;
705 2378 : continue;
706 : }
707 :
708 215699 : const int iBMX = static_cast<int>(dfBMX);
709 215699 : const int iBMY = static_cast<int>(dfBMY);
710 :
711 215699 : const auto fBMX_0_0 = pAccessors->backMapXAccessor.Get(iBMX, iBMY);
712 215699 : const auto fBMY_0_0 = pAccessors->backMapYAccessor.Get(iBMX, iBMY);
713 215699 : if (fBMX_0_0 == INVALID_BMXY)
714 : {
715 87041 : bSuccess = FALSE;
716 87041 : panSuccess[i] = FALSE;
717 87041 : padfX[i] = HUGE_VAL;
718 87041 : padfY[i] = HUGE_VAL;
719 87041 : continue;
720 : }
721 :
722 128658 : const auto fBMX_1_0 =
723 : pAccessors->backMapXAccessor.Get(iBMX + 1, iBMY);
724 128658 : const auto fBMY_1_0 =
725 : pAccessors->backMapYAccessor.Get(iBMX + 1, iBMY);
726 128658 : const auto fBMX_0_1 =
727 : pAccessors->backMapXAccessor.Get(iBMX, iBMY + 1);
728 128658 : const auto fBMY_0_1 =
729 : pAccessors->backMapYAccessor.Get(iBMX, iBMY + 1);
730 128658 : const auto fBMX_1_1 =
731 : pAccessors->backMapXAccessor.Get(iBMX + 1, iBMY + 1);
732 128658 : const auto fBMY_1_1 =
733 : pAccessors->backMapYAccessor.Get(iBMX + 1, iBMY + 1);
734 128658 : if (fBMX_1_0 != INVALID_BMXY && fBMX_0_1 != INVALID_BMXY &&
735 : fBMX_1_1 != INVALID_BMXY)
736 : {
737 126373 : padfX[i] = (1 - (dfBMY - iBMY)) *
738 126373 : (double(fBMX_0_0) +
739 126373 : (dfBMX - iBMX) * double(fBMX_1_0 - fBMX_0_0)) +
740 126373 : (dfBMY - iBMY) *
741 126373 : (double(fBMX_0_1) +
742 126373 : (dfBMX - iBMX) * double(fBMX_1_1 - fBMX_0_1));
743 126373 : padfY[i] = (1 - (dfBMY - iBMY)) *
744 126373 : (double(fBMY_0_0) +
745 126373 : (dfBMX - iBMX) * double(fBMY_1_0 - fBMY_0_0)) +
746 126373 : (dfBMY - iBMY) *
747 126373 : (double(fBMY_0_1) +
748 126373 : (dfBMX - iBMX) * double(fBMY_1_1 - fBMY_0_1));
749 : }
750 2285 : else if (fBMX_1_0 != INVALID_BMXY)
751 : {
752 1726 : padfX[i] = double(fBMX_0_0) +
753 1726 : (dfBMX - iBMX) * double(fBMX_1_0 - fBMX_0_0);
754 1726 : padfY[i] = double(fBMY_0_0) +
755 1726 : (dfBMX - iBMX) * double(fBMY_1_0 - fBMY_0_0);
756 : }
757 559 : else if (fBMX_0_1 != INVALID_BMXY)
758 : {
759 424 : padfX[i] = double(fBMX_0_0) +
760 424 : (dfBMY - iBMY) * double(fBMX_0_1 - fBMX_0_0);
761 424 : padfY[i] = double(fBMY_0_0) +
762 424 : (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 128658 : const double dfGeoLocPixel =
771 128658 : (padfX[i] - psTransform->dfPIXEL_OFFSET) /
772 128658 : psTransform->dfPIXEL_STEP -
773 : dfGeorefConventionOffset;
774 128658 : const double dfGeoLocLine =
775 128658 : (padfY[i] - psTransform->dfLINE_OFFSET) /
776 128658 : 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 128658 : 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 128658 : oPoint.setX(dfGeoX);
807 128658 : 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 128658 : const int nSearchRadius =
815 128658 : psTransform->bGeographicSRSWithMinus180Plus180LongRange &&
816 78128 : fabs(dfGeoY) >= 85
817 : ? 5
818 : : 3;
819 128658 : const int nGeoLocPixel =
820 128658 : static_cast<int>(std::floor(dfGeoLocPixel));
821 128658 : const int nGeoLocLine = static_cast<int>(std::floor(dfGeoLocLine));
822 :
823 128658 : 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 347357 : for (int r = 0; !bDone && r <= nSearchRadius; r++)
828 : {
829 1756583 : 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 2947100 : const int sx = (r == 0) ? 0
836 1409225 : : (iter < 2 * r) ? -r + iter
837 1034853 : : (iter < 4 * r) ? r
838 679080 : : (iter < 6 * r) ? r - (iter - 4 * r)
839 : : -r;
840 2947100 : const int sy = (r == 0) ? 0
841 1409225 : : (iter < 2 * r) ? r
842 1034853 : : (iter < 4 * r) ? r - (iter - 2 * r)
843 679080 : : (iter < 6 * r) ? -r
844 333054 : : -r + (iter - 6 * r);
845 1877407 : if (nGeoLocPixel >=
846 1537883 : static_cast<int>(psTransform->nGeoLocXSize) - sx ||
847 : nGeoLocLine >=
848 1293608 : static_cast<int>(psTransform->nGeoLocYSize) - sy)
849 : {
850 339528 : continue;
851 : }
852 1198355 : const int iX = nGeoLocPixel + sx;
853 1198355 : const int iY = nGeoLocLine + sy;
854 1198355 : if (iX >= -1 || iY >= -1)
855 : {
856 : double x0, y0, x1, y1, x2, y2, x3, y3;
857 :
858 1197996 : if (!PixelLineToXY(psTransform, iX, iY, x0, y0) ||
859 1154940 : !PixelLineToXY(psTransform, iX + 1, iY, x2, y2) ||
860 3498660 : !PixelLineToXY(psTransform, iX, iY + 1, x1, y1) ||
861 1145720 : !PixelLineToXY(psTransform, iX + 1, iY + 1, x3, y3))
862 : {
863 52754 : continue;
864 : }
865 :
866 1145242 : int nIters = 1;
867 : // For a bounding box crossing the anti-meridian, check
868 : // both around -180 and +180 deg.
869 1145242 : if (psTransform
870 1145242 : ->bGeographicSRSWithMinus180Plus180LongRange &&
871 997727 : std::fabs(x0) > 170 && std::fabs(x1) > 170 &&
872 53449 : std::fabs(x2) > 170 && std::fabs(x3) > 170 &&
873 38624 : (std::fabs(x1 - x0) > 180 ||
874 38293 : std::fabs(x2 - x0) > 180 ||
875 31416 : std::fabs(x3 - x0) > 180))
876 : {
877 7510 : nIters = 2;
878 7510 : if (x0 > 0)
879 58 : x0 -= 360;
880 7510 : if (x1 > 0)
881 331 : x1 -= 360;
882 7510 : if (x2 > 0)
883 7179 : x2 -= 360;
884 7510 : if (x3 > 0)
885 7510 : x3 -= 360;
886 : }
887 2297855 : for (int iIter = 0; !bDone && iIter < nIters; ++iIter)
888 : {
889 1152614 : if (iIter == 1)
890 : {
891 7372 : x0 += 360;
892 7372 : x1 += 360;
893 7372 : x2 += 360;
894 7372 : x3 += 360;
895 : }
896 1152614 : oRing.setPoint(0, x0, y0);
897 1152614 : oRing.setPoint(1, x2, y2);
898 1152614 : oRing.setPoint(2, x3, y3);
899 1152614 : oRing.setPoint(3, x1, y1);
900 1152614 : oRing.setPoint(4, x0, y0);
901 2196600 : if (oRing.isPointInRing(&oPoint) ||
902 1043982 : oRing.isPointOnRingBoundary(&oPoint))
903 : {
904 110351 : double dfX = static_cast<double>(iX);
905 110351 : double dfY = static_cast<double>(iY);
906 110351 : GDALInverseBilinearInterpolation(
907 : dfGeoX, dfGeoY, x0, y0, x1, y1, x2, y2, x3,
908 : y3, dfX, dfY);
909 :
910 110351 : dfX = (dfX + dfGeorefConventionOffset) *
911 110351 : psTransform->dfPIXEL_STEP +
912 110351 : psTransform->dfPIXEL_OFFSET;
913 110351 : dfY = (dfY + dfGeorefConventionOffset) *
914 110351 : psTransform->dfLINE_STEP +
915 110351 : 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 110351 : padfX[i] = dfX;
925 110351 : padfY[i] = dfY;
926 :
927 110351 : bDone = true;
928 : }
929 : }
930 : }
931 : }
932 : }
933 128658 : if (!bDone)
934 : {
935 18307 : bSuccess = FALSE;
936 18307 : panSuccess[i] = FALSE;
937 18307 : padfX[i] = HUGE_VAL;
938 18307 : padfY[i] = HUGE_VAL;
939 18307 : continue;
940 : }
941 :
942 110351 : panSuccess[i] = TRUE;
943 : }
944 : }
945 :
946 31459 : 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 251034 : 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 251034 : const double A = (x0 - x) * (y0 - y2) - (y0 - y) * (x0 - x2);
971 251034 : const double B = (((x0 - x) * (y1 - y3) - (y0 - y) * (x1 - x3)) +
972 251034 : ((x1 - x) * (y0 - y2) - (y1 - y) * (x0 - x2))) /
973 : 2;
974 251034 : const double C = (x1 - x) * (y1 - y3) - (y1 - y) * (x1 - x3);
975 251034 : const double denom = A - 2 * B + C;
976 : double s;
977 251034 : const double magnitudeOfValues = fabs(A) + fabs(B) + fabs(C);
978 251034 : 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 133449 : s = A / (A - C);
983 : }
984 : else
985 : {
986 117585 : const double sqrtTerm = sqrt(B * B - A * C);
987 117585 : const double s1 = ((A - B) + sqrtTerm) / denom;
988 117585 : const double s2 = ((A - B) - sqrtTerm) / denom;
989 117585 : if (s1 < 0 || s1 > 1)
990 100530 : s = s2;
991 : else
992 17055 : s = s1;
993 : }
994 :
995 251034 : const double t_denom_x = (1 - s) * (x0 - x2) + s * (x1 - x3);
996 251034 : if (fabs(t_denom_x) > 1e-12 * magnitudeOfValues)
997 : {
998 249275 : i += ((1 - s) * (x0 - x) + s * (x1 - x)) / t_denom_x;
999 : }
1000 : else
1001 : {
1002 1759 : const double t_denom_y = (1 - s) * (y0 - y2) + s * (y1 - y3);
1003 1759 : if (fabs(t_denom_y) > 1e-12 * magnitudeOfValues)
1004 : {
1005 1759 : i += ((1 - s) * (y0 - y) + s * (y1 - y)) / t_denom_y;
1006 : }
1007 : }
1008 :
1009 251034 : j += s;
1010 251034 : }
1011 :
1012 : /************************************************************************/
1013 : /* GeoLocGenerateBackMap() */
1014 : /************************************************************************/
1015 :
1016 : /*! @cond Doxygen_Suppress */
1017 :
1018 : template <class Accessors>
1019 46 : bool GDALGeoLoc<Accessors>::GenerateBackMap(
1020 : GDALGeoLocTransformInfo *psTransform)
1021 :
1022 : {
1023 46 : CPLDebug("GEOLOC", "Starting backmap generation");
1024 46 : const int nXSize = psTransform->nGeoLocXSize;
1025 46 : 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 46 : const double dfTargetPixels =
1034 46 : static_cast<double>(nXSize) * nYSize * psTransform->dfOversampleFactor;
1035 : const double dfPixelSizeSquare =
1036 46 : sqrt((psTransform->dfMaxX - psTransform->dfMinX) *
1037 46 : (psTransform->dfMaxY - psTransform->dfMinY) / dfTargetPixels);
1038 46 : if (dfPixelSizeSquare == 0.0)
1039 : {
1040 0 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid pixel size for backmap");
1041 0 : return false;
1042 : }
1043 :
1044 46 : const double dfMinX = psTransform->dfMinX - dfPixelSizeSquare / 2.0;
1045 46 : const double dfMaxX = psTransform->dfMaxX + dfPixelSizeSquare / 2.0;
1046 46 : const double dfMaxY = psTransform->dfMaxY + dfPixelSizeSquare / 2.0;
1047 46 : const double dfMinY = psTransform->dfMinY - dfPixelSizeSquare / 2.0;
1048 46 : const double dfBMXSize = std::ceil((dfMaxX - dfMinX) / dfPixelSizeSquare);
1049 46 : 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 46 : if (!(dfBMXSize > 0 && dfBMXSize + 2 < INT_MAX) ||
1054 46 : !(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 46 : int nBMXSize = static_cast<int>(dfBMXSize);
1062 46 : int nBMYSize = static_cast<int>(dfBMYSize);
1063 :
1064 92 : if (static_cast<size_t>(1 + nBMYSize) >
1065 46 : 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 46 : const double dfPixelXSize = (dfMaxX - dfMinX) / nBMXSize;
1073 46 : const double dfPixelYSize = (dfMaxY - dfMinY) / nBMYSize;
1074 :
1075 : // Extra pixel for right-edge and bottom-edge extensions in TOP_LEFT_CORNER
1076 : // convention.
1077 46 : nBMXSize++;
1078 46 : nBMYSize++;
1079 46 : psTransform->nBackMapWidth = nBMXSize;
1080 46 : psTransform->nBackMapHeight = nBMYSize;
1081 :
1082 46 : psTransform->adfBackMapGeoTransform[0] = dfMinX;
1083 46 : psTransform->adfBackMapGeoTransform[1] = dfPixelXSize;
1084 46 : psTransform->adfBackMapGeoTransform[2] = 0.0;
1085 46 : psTransform->adfBackMapGeoTransform[3] = dfMaxY;
1086 46 : psTransform->adfBackMapGeoTransform[4] = 0.0;
1087 46 : psTransform->adfBackMapGeoTransform[5] = -dfPixelYSize;
1088 :
1089 : /* -------------------------------------------------------------------- */
1090 : /* Allocate backmap. */
1091 : /* -------------------------------------------------------------------- */
1092 46 : auto pAccessors = static_cast<Accessors *>(psTransform->pAccessors);
1093 46 : if (!pAccessors->AllocateBackMap())
1094 0 : return false;
1095 :
1096 46 : const double dfGeorefConventionOffset =
1097 46 : psTransform->bOriginIsTopLeftCorner ? 0 : 0.5;
1098 :
1099 46 : const auto UpdateBackmap =
1100 9050905 : [&](int iBMX, int iBMY, double dfX, double dfY, double tempwt)
1101 : {
1102 827142 : const auto fBMX = pAccessors->backMapXAccessor.Get(iBMX, iBMY);
1103 827142 : const auto fBMY = pAccessors->backMapYAccessor.Get(iBMX, iBMY);
1104 827142 : const float fUpdatedBMX =
1105 : fBMX +
1106 827142 : static_cast<float>(tempwt * ((dfX + dfGeorefConventionOffset) *
1107 827142 : psTransform->dfPIXEL_STEP +
1108 827142 : psTransform->dfPIXEL_OFFSET));
1109 827142 : const float fUpdatedBMY =
1110 : fBMY +
1111 827142 : static_cast<float>(tempwt * ((dfY + dfGeorefConventionOffset) *
1112 827142 : psTransform->dfLINE_STEP +
1113 827142 : psTransform->dfLINE_OFFSET));
1114 827142 : const float fUpdatedWeight =
1115 827142 : pAccessors->backMapWeightAccessor.Get(iBMX, iBMY) +
1116 827142 : 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 827142 : if (fUpdatedWeight > 0)
1123 : {
1124 827139 : const float fX = fUpdatedBMX / fUpdatedWeight;
1125 827139 : const float fY = fUpdatedBMY / fUpdatedWeight;
1126 827139 : const double dfGeoLocPixel =
1127 827139 : (double(fX) - psTransform->dfPIXEL_OFFSET) /
1128 827139 : psTransform->dfPIXEL_STEP -
1129 : dfGeorefConventionOffset;
1130 827139 : const double dfGeoLocLine =
1131 827139 : (double(fY) - psTransform->dfLINE_OFFSET) /
1132 827139 : psTransform->dfLINE_STEP -
1133 : dfGeorefConventionOffset;
1134 827139 : int iXAvg = static_cast<int>(std::max(0.0, dfGeoLocPixel));
1135 827139 : iXAvg = std::min(iXAvg, psTransform->nGeoLocXSize - 1);
1136 827139 : int iYAvg = static_cast<int>(std::max(0.0, dfGeoLocLine));
1137 827139 : iYAvg = std::min(iYAvg, psTransform->nGeoLocYSize - 1);
1138 827139 : const double dfGLX = pAccessors->geolocXAccessor.Get(iXAvg, iYAvg);
1139 827139 : const double dfGLY = pAccessors->geolocYAccessor.Get(iXAvg, iYAvg);
1140 :
1141 827139 : const unsigned iX = static_cast<unsigned>(dfX);
1142 827139 : const unsigned iY = static_cast<unsigned>(dfY);
1143 1654270 : if (!(psTransform->bHasNoData && dfGLX == psTransform->dfNoDataX) &&
1144 827129 : ((iX >= static_cast<unsigned>(nXSize - 1) ||
1145 813043 : iY >= static_cast<unsigned>(nYSize - 1)) ||
1146 802939 : (fabs(dfGLX - pAccessors->geolocXAccessor.Get(iX, iY)) <=
1147 802939 : 2 * dfPixelXSize &&
1148 760491 : fabs(dfGLY - pAccessors->geolocYAccessor.Get(iX, iY)) <=
1149 760491 : 2 * dfPixelYSize)))
1150 : {
1151 784536 : pAccessors->backMapXAccessor.Set(iBMX, iBMY, fUpdatedBMX);
1152 784536 : pAccessors->backMapYAccessor.Set(iBMX, iBMY, fUpdatedBMY);
1153 784536 : 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 92 : OGRPoint oPoint;
1162 92 : OGRLinearRing oRing;
1163 46 : 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 46 : const double dfStep = 1. / psTransform->dfOversampleFactor;
1174 :
1175 46 : constexpr int TILE_SIZE = GDALGeoLocDatasetAccessors::TILE_SIZE;
1176 46 : const int nYBlocks = DIV_ROUND_UP(nYSize, TILE_SIZE);
1177 46 : 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 92 : std::vector<std::pair<double, double>> yStartEnd(nYBlocks + 1);
1182 92 : std::vector<std::pair<double, double>> xStartEnd(nXBlocks + 1);
1183 :
1184 : {
1185 46 : int iYBlockLast = -1;
1186 46 : double dfY = -dfStep;
1187 2418 : for (; dfY <= static_cast<double>(nYSize) + 2 * dfStep; dfY += dfStep)
1188 : {
1189 2372 : const int iYBlock = static_cast<int>(dfY / TILE_SIZE);
1190 2372 : if (iYBlock != iYBlockLast)
1191 : {
1192 48 : CPLAssert(iYBlock == iYBlockLast + 1);
1193 48 : if (iYBlockLast >= 0)
1194 2 : yStartEnd[iYBlockLast].second = dfY + dfStep / 10;
1195 48 : yStartEnd[iYBlock].first = dfY;
1196 48 : iYBlockLast = iYBlock;
1197 : }
1198 : }
1199 46 : const int iYBlock = static_cast<int>(dfY / TILE_SIZE);
1200 46 : yStartEnd[iYBlock].second = dfY + dfStep / 10;
1201 : }
1202 :
1203 : {
1204 46 : int iXBlockLast = -1;
1205 46 : double dfX = -dfStep;
1206 2900 : for (; dfX <= static_cast<double>(nXSize) + 2 * dfStep; dfX += dfStep)
1207 : {
1208 2854 : const int iXBlock = static_cast<int>(dfX / TILE_SIZE);
1209 2854 : if (iXBlock != iXBlockLast)
1210 : {
1211 48 : CPLAssert(iXBlock == iXBlockLast + 1);
1212 48 : if (iXBlockLast >= 0)
1213 2 : xStartEnd[iXBlockLast].second = dfX + dfStep / 10;
1214 48 : xStartEnd[iXBlock].first = dfX;
1215 48 : iXBlockLast = iXBlock;
1216 : }
1217 : }
1218 46 : const int iXBlock = static_cast<int>(dfX / TILE_SIZE);
1219 46 : xStartEnd[iXBlock].second = dfX + dfStep / 10;
1220 : }
1221 :
1222 94 : for (int iYBlock = 0; iYBlock < nYBlocks; ++iYBlock)
1223 : {
1224 98 : 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 2686 : for (double dfY = yStartEnd[iYBlock].first;
1233 2686 : dfY < yStartEnd[iYBlock].second; dfY += dfStep)
1234 : {
1235 431389 : for (double dfX = xStartEnd[iXBlock].first;
1236 431389 : 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 428753 : if (!PixelLineToXY(psTransform, dfX, dfY, dfGeoLocX,
1243 : dfGeoLocY))
1244 147817 : continue;
1245 :
1246 : // Compute the floating point coordinates in the pixel space
1247 : // of the backmap
1248 422071 : const double dBMX = static_cast<double>(
1249 422071 : (dfGeoLocX - dfMinX) / dfPixelXSize);
1250 :
1251 422071 : const double dBMY = static_cast<double>(
1252 422071 : (dfMaxY - dfGeoLocY) / dfPixelYSize);
1253 :
1254 : // Get top left index by truncation
1255 422071 : const int iBMX = static_cast<int>(std::floor(dBMX));
1256 422071 : const int iBMY = static_cast<int>(std::floor(dBMY));
1257 :
1258 422071 : 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 420174 : double dfGeoX = dfMinX + iBMX * dfPixelXSize;
1264 420174 : const double dfGeoY = dfMaxY - iBMY * dfPixelYSize;
1265 :
1266 420174 : bool bMatchingGeoLocCellFound = false;
1267 :
1268 420174 : const int nOuterIters =
1269 420174 : psTransform->bGeographicSRSWithMinus180Plus180LongRange &&
1270 315836 : fabs(dfGeoX) >= 180
1271 : ? 2
1272 : : 1;
1273 :
1274 840372 : for (int iOuterIter = 0; iOuterIter < nOuterIters;
1275 : ++iOuterIter)
1276 : {
1277 420198 : if (iOuterIter == 1 && dfGeoX >= 180)
1278 0 : dfGeoX -= 360;
1279 420198 : 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 420198 : oPoint.setX(dfGeoX);
1286 420198 : oPoint.setY(dfGeoY);
1287 420198 : const int nX = static_cast<int>(std::floor(dfX));
1288 420198 : const int nY = static_cast<int>(std::floor(dfY));
1289 1206926 : for (int sx = -1;
1290 1206926 : !bMatchingGeoLocCellFound && sx <= 0; sx++)
1291 : {
1292 2310213 : for (int sy = -1;
1293 2310213 : !bMatchingGeoLocCellFound && sy <= 0; sy++)
1294 : {
1295 1524776 : const int pixel = nX + sx;
1296 1524776 : const int line = nY + sy;
1297 : double x0, y0, x1, y1, x2, y2, x3, y3;
1298 1524776 : if (!PixelLineToXY(psTransform, pixel, line,
1299 1524208 : x0, y0) ||
1300 1524208 : !PixelLineToXY(psTransform, pixel + 1,
1301 1523988 : line, x2, y2) ||
1302 1523988 : !PixelLineToXY(psTransform, pixel,
1303 3048980 : line + 1, x1, y1) ||
1304 1523558 : !PixelLineToXY(psTransform, pixel + 1,
1305 : line + 1, x3, y3))
1306 : {
1307 1288 : break;
1308 : }
1309 :
1310 1523488 : int nIters = 1;
1311 1523488 : if (psTransform
1312 1523488 : ->bGeographicSRSWithMinus180Plus180LongRange &&
1313 1191965 : 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 3048460 : for (int iIter = 0; iIter < nIters; ++iIter)
1332 : {
1333 1524978 : if (iIter == 1)
1334 : {
1335 1490 : x0 += 360;
1336 1490 : x1 += 360;
1337 1490 : x2 += 360;
1338 1490 : x3 += 360;
1339 : }
1340 :
1341 1524978 : oRing.setPoint(0, x0, y0);
1342 1524978 : oRing.setPoint(1, x2, y2);
1343 1524978 : oRing.setPoint(2, x3, y3);
1344 1524978 : oRing.setPoint(3, x1, y1);
1345 1524978 : oRing.setPoint(4, x0, y0);
1346 2909500 : if (oRing.isPointInRing(&oPoint) ||
1347 1384522 : oRing.isPointOnRingBoundary(
1348 : &oPoint))
1349 : {
1350 140659 : bMatchingGeoLocCellFound = true;
1351 140659 : double dfBMXValue = pixel;
1352 140659 : double dfBMYValue = line;
1353 140659 : GDALInverseBilinearInterpolation(
1354 : dfGeoX, dfGeoY, x0, y0, x1, y1,
1355 : x2, y2, x3, y3, dfBMXValue,
1356 : dfBMYValue);
1357 :
1358 140659 : dfBMXValue =
1359 140659 : (dfBMXValue +
1360 140659 : dfGeorefConventionOffset) *
1361 140659 : psTransform->dfPIXEL_STEP +
1362 140659 : psTransform->dfPIXEL_OFFSET;
1363 140659 : dfBMYValue =
1364 140659 : (dfBMYValue +
1365 140659 : dfGeorefConventionOffset) *
1366 140659 : psTransform->dfLINE_STEP +
1367 140659 : psTransform->dfLINE_OFFSET;
1368 :
1369 140659 : pAccessors->backMapXAccessor.Set(
1370 : iBMX, iBMY,
1371 : static_cast<float>(dfBMXValue));
1372 140659 : pAccessors->backMapYAccessor.Set(
1373 : iBMX, iBMY,
1374 : static_cast<float>(dfBMYValue));
1375 140659 : pAccessors->backMapWeightAccessor
1376 : .Set(iBMX, iBMY, 1.0f);
1377 : }
1378 : }
1379 : }
1380 : }
1381 : }
1382 420174 : if (bMatchingGeoLocCellFound)
1383 140659 : 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 281412 : if (iBMX < -1 || iBMY < -1 || iBMX > nBMXSize ||
1391 : iBMY > nBMYSize)
1392 476 : continue;
1393 :
1394 280936 : const double fracBMX = dBMX - iBMX;
1395 280936 : const double fracBMY = dBMY - iBMY;
1396 :
1397 : // Check logic for top left pixel
1398 280511 : if ((iBMX >= 0) && (iBMY >= 0) && (iBMX < nBMXSize) &&
1399 561447 : (iBMY < nBMYSize) &&
1400 279515 : pAccessors->backMapWeightAccessor.Get(iBMX, iBMY) !=
1401 : 1.0f)
1402 : {
1403 187302 : const double tempwt = (1.0 - fracBMX) * (1.0 - fracBMY);
1404 187302 : UpdateBackmap(iBMX, iBMY, dfX, dfY, tempwt);
1405 : }
1406 :
1407 : // Check logic for top right pixel
1408 280528 : if ((iBMY >= 0) && (iBMX + 1 < nBMXSize) &&
1409 561464 : (iBMY < nBMYSize) &&
1410 279884 : pAccessors->backMapWeightAccessor.Get(iBMX + 1, iBMY) !=
1411 : 1.0f)
1412 : {
1413 190561 : const double tempwt = fracBMX * (1.0 - fracBMY);
1414 190561 : UpdateBackmap(iBMX + 1, iBMY, dfX, dfY, tempwt);
1415 : }
1416 :
1417 : // Check logic for bottom right pixel
1418 561142 : if ((iBMX + 1 < nBMXSize) && (iBMY + 1 < nBMYSize) &&
1419 280206 : pAccessors->backMapWeightAccessor.Get(iBMX + 1,
1420 280206 : iBMY + 1) != 1.0f)
1421 : {
1422 227350 : const double tempwt = fracBMX * fracBMY;
1423 227350 : UpdateBackmap(iBMX + 1, iBMY + 1, dfX, dfY, tempwt);
1424 : }
1425 :
1426 : // Check logic for bottom left pixel
1427 280511 : if ((iBMX >= 0) && (iBMX < nBMXSize) &&
1428 841296 : (iBMY + 1 < nBMYSize) &&
1429 279849 : pAccessors->backMapWeightAccessor.Get(iBMX, iBMY + 1) !=
1430 : 1.0f)
1431 : {
1432 221929 : const double tempwt = (1.0 - fracBMX) * fracBMY;
1433 221929 : 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 154 : START_ITER_PER_BLOCK(nBMXSize, TILE_SIZE, nBMYSize, TILE_SIZE, (void)0,
1443 : iXStart, iXEnd, iYStart, iYEnd)
1444 : {
1445 2389 : for (int iY = iYStart; iY < iYEnd; ++iY)
1446 : {
1447 331046 : 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 468468 : const auto weight =
1453 328719 : pAccessors->backMapWeightAccessor.Get(iX, iY);
1454 328719 : if (weight > 0)
1455 : {
1456 239100 : pAccessors->backMapXAccessor.Set(
1457 : iX, iY,
1458 173179 : pAccessors->backMapXAccessor.Get(iX, iY) / weight);
1459 239100 : pAccessors->backMapYAccessor.Set(
1460 : iX, iY,
1461 173179 : pAccessors->backMapYAccessor.Get(iX, iY) / weight);
1462 : }
1463 : else
1464 : {
1465 155540 : pAccessors->backMapXAccessor.Set(iX, iY, INVALID_BMXY);
1466 155540 : pAccessors->backMapYAccessor.Set(iX, iY, INVALID_BMXY);
1467 : }
1468 : }
1469 : }
1470 : }
1471 : END_ITER_PER_BLOCK
1472 :
1473 46 : pAccessors->FreeWghtsBackMap();
1474 :
1475 : // Fill holes in backmap
1476 46 : auto poBackmapDS = pAccessors->GetBackmapDataset();
1477 :
1478 46 : 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 46 : constexpr double dfMaxSearchDist = 3.0;
1491 46 : constexpr int nSmoothingIterations = 1;
1492 138 : for (int i = 1; i <= 2; i++)
1493 : {
1494 92 : 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 46 : std::vector<LastValidStruct> lastValid(TILE_SIZE);
1518 230 : const auto reinitLine = [&lastValid]()
1519 : {
1520 46 : const size_t nSize = lastValid.size();
1521 46 : lastValid.clear();
1522 46 : lastValid.resize(nSize);
1523 : };
1524 154 : START_ITER_PER_BLOCK(nBMXSize, TILE_SIZE, nBMYSize, TILE_SIZE, reinitLine(),
1525 : iXStart, iXEnd, iYStart, iYEnd)
1526 : {
1527 62 : const int iYCount = iYEnd - iYStart;
1528 2389 : for (int iYIter = 0; iYIter < iYCount; ++iYIter)
1529 : {
1530 2327 : int iLastValidIX = lastValid[iYIter].iX;
1531 2327 : float bmXLastValid = lastValid[iYIter].bmX;
1532 2327 : const int iBMY = iYStart + iYIter;
1533 331046 : for (int iBMX = iXStart; iBMX < iXEnd; ++iBMX)
1534 : {
1535 328719 : const float bmX = pAccessors->backMapXAccessor.Get(iBMX, iBMY);
1536 328719 : if (bmX == INVALID_BMXY)
1537 127144 : continue;
1538 202127 : 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 201575 : iLastValidIX = iBMX;
1563 201575 : bmXLastValid = bmX;
1564 : }
1565 2327 : lastValid[iYIter].iX = iLastValidIX;
1566 2327 : 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 46 : pAccessors->ReleaseBackmapDataset(poBackmapDS);
1583 46 : CPLDebug("GEOLOC", "Ending backmap generation");
1584 :
1585 46 : return true;
1586 : }
1587 :
1588 : /*! @endcond */
1589 :
1590 : /************************************************************************/
1591 : /* GDALGeoLocRescale() */
1592 : /************************************************************************/
1593 :
1594 4 : static void GDALGeoLocRescale(char **&papszMD, const char *pszItem,
1595 : double dfRatio, double dfDefaultVal)
1596 : {
1597 : const double dfVal =
1598 4 : dfRatio * CPLAtofM(CSLFetchNameValueDef(
1599 4 : papszMD, pszItem, CPLSPrintf("%.17g", dfDefaultVal)));
1600 :
1601 4 : papszMD = CSLSetNameValue(papszMD, pszItem, CPLSPrintf("%.17g", dfVal));
1602 4 : }
1603 :
1604 : /************************************************************************/
1605 : /* GDALCreateSimilarGeoLocTransformer() */
1606 : /************************************************************************/
1607 :
1608 1 : static void *GDALCreateSimilarGeoLocTransformer(void *hTransformArg,
1609 : double dfRatioX,
1610 : double dfRatioY)
1611 : {
1612 1 : VALIDATE_POINTER1(hTransformArg, "GDALCreateSimilarGeoLocTransformer",
1613 : nullptr);
1614 :
1615 1 : GDALGeoLocTransformInfo *psInfo =
1616 : static_cast<GDALGeoLocTransformInfo *>(hTransformArg);
1617 :
1618 1 : char **papszGeolocationInfo = CSLDuplicate(psInfo->papszGeolocationInfo);
1619 :
1620 1 : if (dfRatioX != 1.0 || dfRatioY != 1.0)
1621 : {
1622 1 : GDALGeoLocRescale(papszGeolocationInfo, "PIXEL_OFFSET", dfRatioX, 0.0);
1623 1 : GDALGeoLocRescale(papszGeolocationInfo, "LINE_OFFSET", dfRatioY, 0.0);
1624 1 : GDALGeoLocRescale(papszGeolocationInfo, "PIXEL_STEP", 1.0 / dfRatioX,
1625 : 1.0);
1626 1 : GDALGeoLocRescale(papszGeolocationInfo, "LINE_STEP", 1.0 / dfRatioY,
1627 : 1.0);
1628 : }
1629 :
1630 : auto psInfoNew =
1631 2 : static_cast<GDALGeoLocTransformInfo *>(GDALCreateGeoLocTransformer(
1632 1 : nullptr, papszGeolocationInfo, psInfo->bReversed));
1633 1 : psInfoNew->dfOversampleFactor = psInfo->dfOversampleFactor;
1634 :
1635 1 : CSLDestroy(papszGeolocationInfo);
1636 :
1637 1 : return psInfoNew;
1638 : }
1639 :
1640 : /************************************************************************/
1641 : /* GDALCreateGeolocationMetadata() */
1642 : /************************************************************************/
1643 :
1644 : /** Synthesize the content of a GEOLOCATION metadata domain from a
1645 : * geolocation dataset.
1646 : *
1647 : * This is used when doing gdalwarp -to GEOLOC_ARRAY=some.tif
1648 : */
1649 13 : CPLStringList GDALCreateGeolocationMetadata(GDALDatasetH hBaseDS,
1650 : const char *pszGeolocationDataset,
1651 : bool bIsSource)
1652 : {
1653 26 : CPLStringList aosMD;
1654 :
1655 : // Open geolocation dataset
1656 : auto poGeolocDS = std::unique_ptr<GDALDataset>(
1657 26 : GDALDataset::Open(pszGeolocationDataset, GDAL_OF_RASTER));
1658 13 : if (poGeolocDS == nullptr)
1659 : {
1660 2 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid dataset: %s",
1661 : pszGeolocationDataset);
1662 2 : return CPLStringList();
1663 : }
1664 11 : const int nGeoLocXSize = poGeolocDS->GetRasterXSize();
1665 11 : const int nGeoLocYSize = poGeolocDS->GetRasterYSize();
1666 11 : if (nGeoLocXSize == 0 || nGeoLocYSize == 0)
1667 : {
1668 0 : CPLError(CE_Failure, CPLE_AppDefined,
1669 : "Invalid dataset dimension for %s: %dx%d",
1670 : pszGeolocationDataset, nGeoLocXSize, nGeoLocYSize);
1671 0 : return CPLStringList();
1672 : }
1673 :
1674 : // Import the GEOLOCATION metadata from the geolocation dataset, if existing
1675 11 : auto papszGeolocMDFromGeolocDS = poGeolocDS->GetMetadata("GEOLOCATION");
1676 11 : if (papszGeolocMDFromGeolocDS)
1677 3 : aosMD = CSLDuplicate(papszGeolocMDFromGeolocDS);
1678 :
1679 11 : aosMD.SetNameValue("X_DATASET", pszGeolocationDataset);
1680 11 : 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 19 : if (aosMD.FetchNameValue("X_BAND") == nullptr &&
1685 8 : aosMD.FetchNameValue("Y_BAND") == nullptr)
1686 : {
1687 8 : 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 7 : aosMD.SetNameValue("X_BAND", "1");
1695 7 : 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 10 : if (aosMD.FetchNameValue("SRS") == nullptr)
1701 : {
1702 10 : auto poSRS = poGeolocDS->GetSpatialRef();
1703 10 : 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 10 : if (aosMD.FetchNameValue("SRS") == nullptr)
1712 : {
1713 7 : aosMD.SetNameValue("SRS", SRS_WKT_WGS84_LAT_LONG);
1714 : }
1715 :
1716 : // Set default values for PIXEL/LINE_OFFSET/STEP if not present.
1717 10 : if (aosMD.FetchNameValue("PIXEL_OFFSET") == nullptr)
1718 7 : aosMD.SetNameValue("PIXEL_OFFSET", "0");
1719 :
1720 10 : if (aosMD.FetchNameValue("LINE_OFFSET") == nullptr)
1721 7 : aosMD.SetNameValue("LINE_OFFSET", "0");
1722 :
1723 10 : if (aosMD.FetchNameValue("PIXEL_STEP") == nullptr)
1724 : {
1725 : aosMD.SetNameValue(
1726 7 : "PIXEL_STEP", CPLSPrintf("%.17g", static_cast<double>(
1727 7 : GDALGetRasterXSize(hBaseDS)) /
1728 7 : nGeoLocXSize));
1729 : }
1730 :
1731 10 : if (aosMD.FetchNameValue("LINE_STEP") == nullptr)
1732 : {
1733 : aosMD.SetNameValue(
1734 7 : "LINE_STEP", CPLSPrintf("%.17g", static_cast<double>(
1735 7 : GDALGetRasterYSize(hBaseDS)) /
1736 7 : nGeoLocYSize));
1737 : }
1738 :
1739 10 : if (aosMD.FetchNameValue("GEOREFERENCING_CONVENTION") == nullptr)
1740 : {
1741 : const char *pszConvention =
1742 10 : poGeolocDS->GetMetadataItem("GEOREFERENCING_CONVENTION");
1743 10 : if (pszConvention)
1744 2 : aosMD.SetNameValue("GEOREFERENCING_CONVENTION", pszConvention);
1745 : }
1746 :
1747 20 : std::string osDebugMsg;
1748 10 : osDebugMsg = "Synthesized GEOLOCATION metadata for ";
1749 10 : osDebugMsg += bIsSource ? "source" : "target";
1750 10 : osDebugMsg += ":\n";
1751 102 : for (int i = 0; i < aosMD.size(); ++i)
1752 : {
1753 92 : osDebugMsg += " ";
1754 92 : osDebugMsg += aosMD[i];
1755 92 : osDebugMsg += '\n';
1756 : }
1757 10 : CPLDebug("GEOLOC", "%s", osDebugMsg.c_str());
1758 :
1759 10 : return aosMD;
1760 : }
1761 :
1762 : /************************************************************************/
1763 : /* GDALCreateGeoLocTransformer() */
1764 : /************************************************************************/
1765 :
1766 52 : void *GDALCreateGeoLocTransformerEx(GDALDatasetH hBaseDS,
1767 : CSLConstList papszGeolocationInfo,
1768 : int bReversed, const char *pszSourceDataset,
1769 : CSLConstList papszTransformOptions)
1770 :
1771 : {
1772 :
1773 52 : if (CSLFetchNameValue(papszGeolocationInfo, "PIXEL_OFFSET") == nullptr ||
1774 50 : CSLFetchNameValue(papszGeolocationInfo, "LINE_OFFSET") == nullptr ||
1775 50 : CSLFetchNameValue(papszGeolocationInfo, "PIXEL_STEP") == nullptr ||
1776 50 : CSLFetchNameValue(papszGeolocationInfo, "LINE_STEP") == nullptr ||
1777 152 : CSLFetchNameValue(papszGeolocationInfo, "X_BAND") == nullptr ||
1778 50 : 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 50 : CPLCalloc(sizeof(GDALGeoLocTransformInfo), 1));
1792 :
1793 50 : psTransform->bReversed = CPL_TO_BOOL(bReversed);
1794 50 : psTransform->dfOversampleFactor = std::max(
1795 100 : 0.1,
1796 100 : std::min(2.0,
1797 50 : CPLAtof(CSLFetchNameValueDef(
1798 : papszTransformOptions, "GEOLOC_BACKMAP_OVERSAMPLE_FACTOR",
1799 : CPLGetConfigOption("GDAL_GEOLOC_BACKMAP_OVERSAMPLE_FACTOR",
1800 50 : "1.3")))));
1801 :
1802 50 : memcpy(psTransform->sTI.abySignature, GDAL_GTI2_SIGNATURE,
1803 : strlen(GDAL_GTI2_SIGNATURE));
1804 50 : psTransform->sTI.pszClassName = "GDALGeoLocTransformer";
1805 50 : psTransform->sTI.pfnTransform = GDALGeoLocTransform;
1806 50 : psTransform->sTI.pfnCleanup = GDALDestroyGeoLocTransformer;
1807 50 : psTransform->sTI.pfnSerialize = GDALSerializeGeoLocTransformer;
1808 50 : psTransform->sTI.pfnCreateSimilar = GDALCreateSimilarGeoLocTransformer;
1809 :
1810 50 : psTransform->papszGeolocationInfo = CSLDuplicate(papszGeolocationInfo);
1811 :
1812 50 : psTransform->bGeographicSRSWithMinus180Plus180LongRange =
1813 50 : 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 50 : psTransform->dfPIXEL_OFFSET =
1821 50 : CPLAtof(CSLFetchNameValue(papszGeolocationInfo, "PIXEL_OFFSET"));
1822 50 : psTransform->dfLINE_OFFSET =
1823 50 : CPLAtof(CSLFetchNameValue(papszGeolocationInfo, "LINE_OFFSET"));
1824 50 : psTransform->dfPIXEL_STEP =
1825 50 : CPLAtof(CSLFetchNameValue(papszGeolocationInfo, "PIXEL_STEP"));
1826 50 : psTransform->dfLINE_STEP =
1827 50 : CPLAtof(CSLFetchNameValue(papszGeolocationInfo, "LINE_STEP"));
1828 :
1829 50 : 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 50 : CSLFetchNameValue(papszGeolocationInfo, "X_DATASET");
1839 50 : if (pszDSName != nullptr)
1840 : {
1841 100 : CPLConfigOptionSetter oSetter("CPL_ALLOW_VSISTDIN", "NO", true);
1842 50 : if (CPLTestBool(CSLFetchNameValueDef(
1843 51 : 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 46 : 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 50 : pszDSName = CSLFetchNameValue(papszGeolocationInfo, "Y_DATASET");
1873 50 : if (pszDSName != nullptr)
1874 : {
1875 100 : CPLConfigOptionSetter oSetter("CPL_ALLOW_VSISTDIN", "NO", true);
1876 50 : if (CPLTestBool(CSLFetchNameValueDef(
1877 51 : 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 46 : 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 50 : 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 50 : std::max(1, atoi(CSLFetchNameValue(papszGeolocationInfo, "X_BAND")));
1917 50 : psTransform->hBand_X = GDALGetRasterBand(psTransform->hDS_X, nXBand);
1918 :
1919 50 : psTransform->dfNoDataX = GDALGetRasterNoDataValue(
1920 : psTransform->hBand_X, &(psTransform->bHasNoData));
1921 :
1922 : const int nYBand =
1923 50 : std::max(1, atoi(CSLFetchNameValue(papszGeolocationInfo, "Y_BAND")));
1924 50 : psTransform->hBand_Y = GDALGetRasterBand(psTransform->hDS_Y, nYBand);
1925 :
1926 50 : if (psTransform->hBand_X == nullptr || psTransform->hBand_Y == nullptr)
1927 : {
1928 0 : GDALDestroyGeoLocTransformer(psTransform);
1929 0 : return nullptr;
1930 : }
1931 :
1932 50 : psTransform->bSwapXY = CPLTestBool(
1933 : CSLFetchNameValueDef(papszGeolocationInfo, "SWAP_XY", "NO"));
1934 :
1935 : /* -------------------------------------------------------------------- */
1936 : /* Check that X and Y bands have the same dimensions */
1937 : /* -------------------------------------------------------------------- */
1938 50 : const int nXSize_XBand = GDALGetRasterXSize(psTransform->hDS_X);
1939 50 : const int nYSize_XBand = GDALGetRasterYSize(psTransform->hDS_X);
1940 50 : const int nXSize_YBand = GDALGetRasterXSize(psTransform->hDS_Y);
1941 50 : const int nYSize_YBand = GDALGetRasterYSize(psTransform->hDS_Y);
1942 50 : 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 35 : 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 50 : 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 50 : const bool bIsRegularGrid = (nYSize_XBand == 1 && nYSize_YBand == 1);
1972 :
1973 50 : const int nXSize = nXSize_XBand;
1974 50 : const int nYSize = bIsRegularGrid ? nXSize_YBand : nYSize_XBand;
1975 :
1976 100 : if (static_cast<size_t>(nXSize) >
1977 50 : 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 50 : psTransform->nGeoLocXSize = nXSize;
1986 50 : psTransform->nGeoLocYSize = nYSize;
1987 :
1988 50 : if (hBaseDS && psTransform->dfPIXEL_OFFSET == 0 &&
1989 46 : 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 50 : 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 50 : const char *pszUseTempDatasets = CSLFetchNameValueDef(
2017 : papszTransformOptions, "GEOLOC_USE_TEMP_DATASETS",
2018 : CPLGetConfigOption("GDAL_GEOLOC_USE_TEMP_DATASETS", nullptr));
2019 50 : if (pszUseTempDatasets)
2020 4 : psTransform->bUseArray = !CPLTestBool(pszUseTempDatasets);
2021 : else
2022 : {
2023 46 : constexpr int MEGAPIXEL_LIMIT = 24;
2024 46 : psTransform->bUseArray =
2025 46 : nXSize < MEGAPIXEL_LIMIT * 1000 * 1000 / nYSize;
2026 46 : 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 50 : if (psTransform->bUseArray)
2041 : {
2042 48 : auto pAccessors = new GDALGeoLocCArrayAccessors(psTransform);
2043 48 : psTransform->pAccessors = pAccessors;
2044 48 : 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 50 : 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 50 : void GDALDestroyGeoLocTransformer(void *pTransformAlg)
2079 :
2080 : {
2081 50 : if (pTransformAlg == nullptr)
2082 0 : return;
2083 :
2084 50 : GDALGeoLocTransformInfo *psTransform =
2085 : static_cast<GDALGeoLocTransformInfo *>(pTransformAlg);
2086 :
2087 50 : CSLDestroy(psTransform->papszGeolocationInfo);
2088 :
2089 50 : if (psTransform->bUseArray)
2090 : delete static_cast<GDALGeoLocCArrayAccessors *>(
2091 48 : psTransform->pAccessors);
2092 : else
2093 : delete static_cast<GDALGeoLocDatasetAccessors *>(
2094 2 : psTransform->pAccessors);
2095 :
2096 100 : if (psTransform->hDS_X != nullptr &&
2097 50 : GDALDereferenceDataset(psTransform->hDS_X) == 0)
2098 28 : GDALClose(psTransform->hDS_X);
2099 :
2100 100 : if (psTransform->hDS_Y != nullptr &&
2101 50 : GDALDereferenceDataset(psTransform->hDS_Y) == 0)
2102 42 : GDALClose(psTransform->hDS_Y);
2103 :
2104 50 : if (psTransform->hQuadTree != nullptr)
2105 4 : CPLQuadTreeDestroy(psTransform->hQuadTree);
2106 :
2107 50 : CPLFree(pTransformAlg);
2108 : }
2109 :
2110 : /** Use GeoLocation transformer */
2111 31483 : int GDALGeoLocTransform(void *pTransformArg, int bDstToSrc, int nPointCount,
2112 : double *padfX, double *padfY, double *padfZ,
2113 : int *panSuccess)
2114 : {
2115 31483 : GDALGeoLocTransformInfo *psTransform =
2116 : static_cast<GDALGeoLocTransformInfo *>(pTransformArg);
2117 31483 : if (psTransform->bUseArray)
2118 : {
2119 29360 : return GDALGeoLoc<GDALGeoLocCArrayAccessors>::Transform(
2120 : pTransformArg, bDstToSrc, nPointCount, padfX, padfY, padfZ,
2121 29360 : panSuccess);
2122 : }
2123 : else
2124 : {
2125 2123 : return GDALGeoLoc<GDALGeoLocDatasetAccessors>::Transform(
2126 : pTransformArg, bDstToSrc, nPointCount, padfX, padfY, padfZ,
2127 2123 : 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 : }
|
