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