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 10501780 : bool GDALGeoLoc<Accessors>::PixelLineToXY(
524 : const GDALGeoLocTransformInfo *psTransform, const int nGeoLocPixel,
525 : const int nGeoLocLine, double &dfX, double &dfY)
526 : {
527 10501780 : if (nGeoLocPixel >= 0 && nGeoLocPixel < psTransform->nGeoLocXSize &&
528 9826260 : nGeoLocLine >= 0 && nGeoLocLine < psTransform->nGeoLocYSize)
529 : {
530 9711660 : auto pAccessors = static_cast<Accessors *>(psTransform->pAccessors);
531 9711660 : const double dfGLX =
532 : pAccessors->geolocXAccessor.Get(nGeoLocPixel, nGeoLocLine);
533 9711660 : const double dfGLY =
534 : pAccessors->geolocYAccessor.Get(nGeoLocPixel, nGeoLocLine);
535 9711660 : if (psTransform->bHasNoData && dfGLX == psTransform->dfNoDataX)
536 : {
537 46946 : return false;
538 : }
539 9664710 : dfX = dfGLX;
540 9664710 : dfY = dfGLY;
541 9664710 : 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 31164 : int GDALGeoLoc<Accessors>::Transform(void *pTransformArg, int bDstToSrc,
588 : int nPointCount, double *padfX,
589 : double *padfY, double * /* padfZ */,
590 : int *panSuccess)
591 : {
592 31164 : int bSuccess = TRUE;
593 31164 : GDALGeoLocTransformInfo *psTransform =
594 : static_cast<GDALGeoLocTransformInfo *>(pTransformArg);
595 :
596 31164 : if (psTransform->bReversed)
597 0 : bDstToSrc = !bDstToSrc;
598 :
599 31164 : const double dfGeorefConventionOffset =
600 31164 : psTransform->bOriginIsTopLeftCorner ? 0 : 0.5;
601 :
602 : /* -------------------------------------------------------------------- */
603 : /* Do original pixel line to target geox/geoy. */
604 : /* -------------------------------------------------------------------- */
605 31164 : if (!bDstToSrc)
606 : {
607 70066 : for (int i = 0; i < nPointCount; i++)
608 : {
609 45165 : if (padfX[i] == HUGE_VAL || padfY[i] == HUGE_VAL)
610 : {
611 2745 : bSuccess = FALSE;
612 2745 : panSuccess[i] = FALSE;
613 2745 : continue;
614 : }
615 :
616 42420 : const double dfGeoLocPixel =
617 42420 : (padfX[i] - psTransform->dfPIXEL_OFFSET) /
618 42420 : psTransform->dfPIXEL_STEP -
619 : dfGeorefConventionOffset;
620 42420 : const double dfGeoLocLine =
621 42420 : (padfY[i] - psTransform->dfLINE_OFFSET) /
622 42420 : psTransform->dfLINE_STEP -
623 : dfGeorefConventionOffset;
624 :
625 42420 : if (!PixelLineToXY(psTransform, dfGeoLocPixel, dfGeoLocLine,
626 42420 : padfX[i], padfY[i]))
627 : {
628 2826 : bSuccess = FALSE;
629 2826 : panSuccess[i] = FALSE;
630 2826 : padfX[i] = HUGE_VAL;
631 2826 : padfY[i] = HUGE_VAL;
632 2826 : continue;
633 : }
634 :
635 39594 : if (psTransform->bSwapXY)
636 : {
637 568 : std::swap(padfX[i], padfY[i]);
638 : }
639 :
640 39594 : panSuccess[i] = TRUE;
641 : }
642 : }
643 :
644 : /* -------------------------------------------------------------------- */
645 : /* geox/geoy to pixel/line using backmap. */
646 : /* -------------------------------------------------------------------- */
647 : else
648 : {
649 6263 : if (psTransform->hQuadTree)
650 : {
651 24 : GDALGeoLocInverseTransformQuadtree(psTransform, nPointCount, padfX,
652 : padfY, panSuccess);
653 24 : return TRUE;
654 : }
655 :
656 6239 : const bool bGeolocMaxAccuracy = CPLTestBool(
657 : CPLGetConfigOption("GDAL_GEOLOC_USE_MAX_ACCURACY", "YES"));
658 :
659 : // Keep those objects in this outer scope, so they are re-used, to
660 : // save memory allocations.
661 12478 : OGRPoint oPoint;
662 12478 : OGRLinearRing oRing;
663 6239 : oRing.setNumPoints(5);
664 :
665 6239 : auto pAccessors = static_cast<Accessors *>(psTransform->pAccessors);
666 :
667 216480 : for (int i = 0; i < nPointCount; i++)
668 : {
669 210241 : if (padfX[i] == HUGE_VAL || padfY[i] == HUGE_VAL)
670 : {
671 0 : bSuccess = FALSE;
672 0 : panSuccess[i] = FALSE;
673 0 : continue;
674 : }
675 :
676 210241 : if (psTransform->bSwapXY)
677 : {
678 836 : std::swap(padfX[i], padfY[i]);
679 : }
680 :
681 210241 : const double dfGeoX = padfX[i];
682 210241 : const double dfGeoY = padfY[i];
683 :
684 210241 : const double dfBMX =
685 210241 : ((padfX[i] - psTransform->adfBackMapGeoTransform[0]) /
686 : psTransform->adfBackMapGeoTransform[1]);
687 210241 : const double dfBMY =
688 210241 : ((padfY[i] - psTransform->adfBackMapGeoTransform[3]) /
689 : psTransform->adfBackMapGeoTransform[5]);
690 :
691 210241 : if (!(dfBMX >= 0 && dfBMY >= 0 &&
692 210111 : dfBMX + 1 < psTransform->nBackMapWidth &&
693 208952 : dfBMY + 1 < psTransform->nBackMapHeight))
694 : {
695 2068 : bSuccess = FALSE;
696 2068 : panSuccess[i] = FALSE;
697 2068 : padfX[i] = HUGE_VAL;
698 2068 : padfY[i] = HUGE_VAL;
699 2068 : continue;
700 : }
701 :
702 208173 : const int iBMX = static_cast<int>(dfBMX);
703 208173 : const int iBMY = static_cast<int>(dfBMY);
704 :
705 208173 : const auto fBMX_0_0 = pAccessors->backMapXAccessor.Get(iBMX, iBMY);
706 208173 : const auto fBMY_0_0 = pAccessors->backMapYAccessor.Get(iBMX, iBMY);
707 208173 : if (fBMX_0_0 == INVALID_BMXY)
708 : {
709 86072 : bSuccess = FALSE;
710 86072 : panSuccess[i] = FALSE;
711 86072 : padfX[i] = HUGE_VAL;
712 86072 : padfY[i] = HUGE_VAL;
713 86072 : continue;
714 : }
715 :
716 122101 : const auto fBMX_1_0 =
717 : pAccessors->backMapXAccessor.Get(iBMX + 1, iBMY);
718 122101 : const auto fBMY_1_0 =
719 : pAccessors->backMapYAccessor.Get(iBMX + 1, iBMY);
720 122101 : const auto fBMX_0_1 =
721 : pAccessors->backMapXAccessor.Get(iBMX, iBMY + 1);
722 122101 : const auto fBMY_0_1 =
723 : pAccessors->backMapYAccessor.Get(iBMX, iBMY + 1);
724 122101 : const auto fBMX_1_1 =
725 : pAccessors->backMapXAccessor.Get(iBMX + 1, iBMY + 1);
726 122101 : const auto fBMY_1_1 =
727 : pAccessors->backMapYAccessor.Get(iBMX + 1, iBMY + 1);
728 122101 : if (fBMX_1_0 != INVALID_BMXY && fBMX_0_1 != INVALID_BMXY &&
729 : fBMX_1_1 != INVALID_BMXY)
730 : {
731 119870 : padfX[i] =
732 119870 : (1 - (dfBMY - iBMY)) *
733 119870 : (fBMX_0_0 + (dfBMX - iBMX) * (fBMX_1_0 - fBMX_0_0)) +
734 119870 : (dfBMY - iBMY) *
735 119870 : (fBMX_0_1 + (dfBMX - iBMX) * (fBMX_1_1 - fBMX_0_1));
736 119870 : padfY[i] =
737 119870 : (1 - (dfBMY - iBMY)) *
738 119870 : (fBMY_0_0 + (dfBMX - iBMX) * (fBMY_1_0 - fBMY_0_0)) +
739 119870 : (dfBMY - iBMY) *
740 119870 : (fBMY_0_1 + (dfBMX - iBMX) * (fBMY_1_1 - fBMY_0_1));
741 : }
742 2231 : else if (fBMX_1_0 != INVALID_BMXY)
743 : {
744 1686 : padfX[i] = fBMX_0_0 + (dfBMX - iBMX) * (fBMX_1_0 - fBMX_0_0);
745 1686 : padfY[i] = fBMY_0_0 + (dfBMX - iBMX) * (fBMY_1_0 - fBMY_0_0);
746 : }
747 545 : else if (fBMX_0_1 != INVALID_BMXY)
748 : {
749 410 : padfX[i] = fBMX_0_0 + (dfBMY - iBMY) * (fBMX_0_1 - fBMX_0_0);
750 410 : padfY[i] = fBMY_0_0 + (dfBMY - iBMY) * (fBMY_0_1 - fBMY_0_0);
751 : }
752 : else
753 : {
754 135 : padfX[i] = fBMX_0_0;
755 135 : padfY[i] = fBMY_0_0;
756 : }
757 :
758 122101 : const double dfGeoLocPixel =
759 122101 : (padfX[i] - psTransform->dfPIXEL_OFFSET) /
760 122101 : psTransform->dfPIXEL_STEP -
761 : dfGeorefConventionOffset;
762 122101 : const double dfGeoLocLine =
763 122101 : (padfY[i] - psTransform->dfLINE_OFFSET) /
764 122101 : psTransform->dfLINE_STEP -
765 : dfGeorefConventionOffset;
766 : #if 0
767 : CPLDebug("GEOLOC", "%f %f %f %f", padfX[i], padfY[i], dfGeoLocPixel, dfGeoLocLine);
768 : if( !psTransform->bOriginIsTopLeftCorner )
769 : {
770 : if( dfGeoLocPixel + dfGeorefConventionOffset > psTransform->nGeoLocXSize-1 ||
771 : dfGeoLocLine + dfGeorefConventionOffset > psTransform->nGeoLocYSize-1 )
772 : {
773 : panSuccess[i] = FALSE;
774 : padfX[i] = HUGE_VAL;
775 : padfY[i] = HUGE_VAL;
776 : continue;
777 : }
778 : }
779 : #endif
780 122101 : if (!bGeolocMaxAccuracy)
781 : {
782 0 : panSuccess[i] = TRUE;
783 0 : continue;
784 : }
785 :
786 : // Now that we have an approximate solution, identify a matching
787 : // cell in the geolocation array, where we can use inverse bilinear
788 : // interpolation to find the exact solution.
789 :
790 : // NOTE: if the geolocation array is an affine transformation,
791 : // the approximate solution should match the exact one, if the
792 : // backmap has correctly been built.
793 :
794 122101 : oPoint.setX(dfGeoX);
795 122101 : oPoint.setY(dfGeoY);
796 : // The thresholds and radius are rather empirical and have been
797 : // tuned on the product
798 : // S5P_TEST_L2__NO2____20190509T220707_20190509T234837_08137_01_010400_20200220T091343.nc
799 : // that includes the north pole.
800 : // Amended with the test case of
801 : // https://github.com/OSGeo/gdal/issues/5823
802 122101 : const int nSearchRadius =
803 122101 : psTransform->bGeographicSRSWithMinus180Plus180LongRange &&
804 72344 : fabs(dfGeoY) >= 85
805 : ? 5
806 : : 3;
807 122101 : const int nGeoLocPixel =
808 122101 : static_cast<int>(std::floor(dfGeoLocPixel));
809 122101 : const int nGeoLocLine = static_cast<int>(std::floor(dfGeoLocLine));
810 :
811 122101 : bool bDone = false;
812 : // Using the above approximate nGeoLocPixel, nGeoLocLine, try to
813 : // find a forward cell that includes (dfGeoX, dfGeoY), with an
814 : // increasing search radius, up to nSearchRadius.
815 329732 : for (int r = 0; !bDone && r <= nSearchRadius; r++)
816 : {
817 1672277 : for (int iter = 0; !bDone && iter < (r == 0 ? 1 : 8 * r);
818 : ++iter)
819 : {
820 : // For r=1, the below formulas will give the following
821 : // offsets:
822 : // (-1,1), (0,1), (1,1), (1,0), (1,-1), (0,-1), (1,-1)
823 2807190 : const int sx = (r == 0) ? 0
824 1342545 : : (iter < 2 * r) ? -r + iter
825 985782 : : (iter < 4 * r) ? r
826 646809 : : (iter < 6 * r) ? r - (iter - 4 * r)
827 : : -r;
828 2807190 : const int sy = (r == 0) ? 0
829 1342545 : : (iter < 2 * r) ? r
830 985782 : : (iter < 4 * r) ? r - (iter - 2 * r)
831 646809 : : (iter < 6 * r) ? -r
832 317228 : : -r + (iter - 6 * r);
833 1783986 : if (nGeoLocPixel >=
834 1464646 : static_cast<int>(psTransform->nGeoLocXSize) - sx ||
835 : nGeoLocLine >=
836 1233276 : static_cast<int>(psTransform->nGeoLocYSize) - sy)
837 : {
838 319344 : continue;
839 : }
840 1145302 : const int iX = nGeoLocPixel + sx;
841 1145302 : const int iY = nGeoLocLine + sy;
842 1145302 : if (iX >= -1 || iY >= -1)
843 : {
844 : double x0, y0, x1, y1, x2, y2, x3, y3;
845 :
846 1145118 : if (!PixelLineToXY(psTransform, iX, iY, x0, y0) ||
847 1105472 : !PixelLineToXY(psTransform, iX + 1, iY, x2, y2) ||
848 3347450 : !PixelLineToXY(psTransform, iX, iY + 1, x1, y1) ||
849 1096860 : !PixelLineToXY(psTransform, iX + 1, iY + 1, x3, y3))
850 : {
851 48686 : continue;
852 : }
853 :
854 1096432 : int nIters = 1;
855 : // For a bounding box crossing the anti-meridian, check
856 : // both around -180 and +180 deg.
857 1096432 : if (psTransform
858 1096432 : ->bGeographicSRSWithMinus180Plus180LongRange &&
859 950423 : std::fabs(x0) > 170 && std::fabs(x1) > 170 &&
860 50764 : std::fabs(x2) > 170 && std::fabs(x3) > 170 &&
861 37304 : (std::fabs(x1 - x0) > 180 ||
862 37012 : std::fabs(x2 - x0) > 180 ||
863 30492 : std::fabs(x3 - x0) > 180))
864 : {
865 7104 : nIters = 2;
866 7104 : if (x0 > 0)
867 0 : x0 -= 360;
868 7104 : if (x1 > 0)
869 292 : x1 -= 360;
870 7104 : if (x2 > 0)
871 6812 : x2 -= 360;
872 7104 : if (x3 > 0)
873 7104 : x3 -= 360;
874 : }
875 2199849 : for (int iIter = 0; !bDone && iIter < nIters; ++iIter)
876 : {
877 1103418 : if (iIter == 1)
878 : {
879 6986 : x0 += 360;
880 6986 : x1 += 360;
881 6986 : x2 += 360;
882 6986 : x3 += 360;
883 : }
884 1103418 : oRing.setPoint(0, x0, y0);
885 1103418 : oRing.setPoint(1, x2, y2);
886 1103418 : oRing.setPoint(2, x3, y3);
887 1103418 : oRing.setPoint(3, x1, y1);
888 1103418 : oRing.setPoint(4, x0, y0);
889 2103446 : if (oRing.isPointInRing(&oPoint) ||
890 1000033 : oRing.isPointOnRingBoundary(&oPoint))
891 : {
892 104799 : double dfX = static_cast<double>(iX);
893 104799 : double dfY = static_cast<double>(iY);
894 104799 : GDALInverseBilinearInterpolation(
895 : dfGeoX, dfGeoY, x0, y0, x1, y1, x2, y2, x3,
896 : y3, dfX, dfY);
897 :
898 104799 : dfX = (dfX + dfGeorefConventionOffset) *
899 104799 : psTransform->dfPIXEL_STEP +
900 104799 : psTransform->dfPIXEL_OFFSET;
901 104799 : dfY = (dfY + dfGeorefConventionOffset) *
902 104799 : psTransform->dfLINE_STEP +
903 104799 : psTransform->dfLINE_OFFSET;
904 :
905 : #ifdef DEBUG_GEOLOC_REALLY_VERBOSE
906 : CPLDebug("GEOLOC",
907 : "value before adjustment: %f %f, "
908 : "after adjustment: %f %f",
909 : padfX[i], padfY[i], dfX, dfY);
910 : #endif
911 :
912 104799 : padfX[i] = dfX;
913 104799 : padfY[i] = dfY;
914 :
915 104799 : bDone = true;
916 : }
917 : }
918 : }
919 : }
920 : }
921 122101 : if (!bDone)
922 : {
923 17302 : bSuccess = FALSE;
924 17302 : panSuccess[i] = FALSE;
925 17302 : padfX[i] = HUGE_VAL;
926 17302 : padfY[i] = HUGE_VAL;
927 17302 : continue;
928 : }
929 :
930 104799 : panSuccess[i] = TRUE;
931 : }
932 : }
933 :
934 31140 : return bSuccess;
935 : }
936 :
937 : /*! @endcond */
938 :
939 : /************************************************************************/
940 : /* GDALInverseBilinearInterpolation() */
941 : /************************************************************************/
942 :
943 : // (i,j) before the call should correspond to the input coordinates that give
944 : // (x0,y0) as output of the forward interpolation
945 : // After the call it will be updated to the input coordinates that give (x,y)
946 : // This assumes that (x,y) is within the polygon formed by
947 : // (x0, y0), (x2, y2), (x3, y3), (x1, y1), (x0, y0)
948 242996 : void GDALInverseBilinearInterpolation(const double x, const double y,
949 : const double x0, const double y0,
950 : const double x1, const double y1,
951 : const double x2, const double y2,
952 : const double x3, const double y3,
953 : double &i, double &j)
954 : {
955 : // Exact inverse bilinear interpolation method.
956 : // Maths from https://stackoverflow.com/a/812077
957 :
958 242996 : const double A = (x0 - x) * (y0 - y2) - (y0 - y) * (x0 - x2);
959 242996 : const double B = (((x0 - x) * (y1 - y3) - (y0 - y) * (x1 - x3)) +
960 242996 : ((x1 - x) * (y0 - y2) - (y1 - y) * (x0 - x2))) /
961 : 2;
962 242996 : const double C = (x1 - x) * (y1 - y3) - (y1 - y) * (x1 - x3);
963 242996 : const double denom = A - 2 * B + C;
964 : double s;
965 242996 : const double magnitudeOfValues = fabs(A) + fabs(B) + fabs(C);
966 242996 : if (fabs(denom) <= 1e-12 * magnitudeOfValues)
967 : {
968 : // Happens typically when the x_i,y_i points form a rectangle
969 : // Can also happen when they form a triangle.
970 130853 : s = A / (A - C);
971 : }
972 : else
973 : {
974 112143 : const double sqrtTerm = sqrt(B * B - A * C);
975 112143 : const double s1 = ((A - B) + sqrtTerm) / denom;
976 112143 : const double s2 = ((A - B) - sqrtTerm) / denom;
977 112143 : if (s1 < 0 || s1 > 1)
978 95694 : s = s2;
979 : else
980 16449 : s = s1;
981 : }
982 :
983 242996 : const double t_denom_x = (1 - s) * (x0 - x2) + s * (x1 - x3);
984 242996 : if (fabs(t_denom_x) > 1e-12 * magnitudeOfValues)
985 : {
986 241486 : i += ((1 - s) * (x0 - x) + s * (x1 - x)) / t_denom_x;
987 : }
988 : else
989 : {
990 1510 : const double t_denom_y = (1 - s) * (y0 - y2) + s * (y1 - y3);
991 1510 : if (fabs(t_denom_y) > 1e-12 * magnitudeOfValues)
992 : {
993 1510 : i += ((1 - s) * (y0 - y) + s * (y1 - y)) / t_denom_y;
994 : }
995 : }
996 :
997 242996 : j += s;
998 242996 : }
999 :
1000 : /************************************************************************/
1001 : /* GeoLocGenerateBackMap() */
1002 : /************************************************************************/
1003 :
1004 : /*! @cond Doxygen_Suppress */
1005 :
1006 : template <class Accessors>
1007 44 : bool GDALGeoLoc<Accessors>::GenerateBackMap(
1008 : GDALGeoLocTransformInfo *psTransform)
1009 :
1010 : {
1011 44 : CPLDebug("GEOLOC", "Starting backmap generation");
1012 44 : const int nXSize = psTransform->nGeoLocXSize;
1013 44 : const int nYSize = psTransform->nGeoLocYSize;
1014 :
1015 : /* -------------------------------------------------------------------- */
1016 : /* Decide on resolution for backmap. We aim for slightly */
1017 : /* higher resolution than the source but we can't easily */
1018 : /* establish how much dead space there is in the backmap, so it */
1019 : /* is approximate. */
1020 : /* -------------------------------------------------------------------- */
1021 44 : const double dfTargetPixels =
1022 44 : static_cast<double>(nXSize) * nYSize * psTransform->dfOversampleFactor;
1023 : const double dfPixelSizeSquare =
1024 44 : sqrt((psTransform->dfMaxX - psTransform->dfMinX) *
1025 44 : (psTransform->dfMaxY - psTransform->dfMinY) / dfTargetPixels);
1026 44 : if (dfPixelSizeSquare == 0.0)
1027 : {
1028 0 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid pixel size for backmap");
1029 0 : return false;
1030 : }
1031 :
1032 44 : const double dfMinX = psTransform->dfMinX - dfPixelSizeSquare / 2.0;
1033 44 : const double dfMaxX = psTransform->dfMaxX + dfPixelSizeSquare / 2.0;
1034 44 : const double dfMaxY = psTransform->dfMaxY + dfPixelSizeSquare / 2.0;
1035 44 : const double dfMinY = psTransform->dfMinY - dfPixelSizeSquare / 2.0;
1036 44 : const double dfBMXSize = std::ceil((dfMaxX - dfMinX) / dfPixelSizeSquare);
1037 44 : const double dfBMYSize = std::ceil((dfMaxY - dfMinY) / dfPixelSizeSquare);
1038 :
1039 : // +2 : +1 due to afterwards nBMXSize++, and another +1 as security margin
1040 : // for other computations.
1041 44 : if (!(dfBMXSize > 0 && dfBMXSize + 2 < INT_MAX) ||
1042 44 : !(dfBMYSize > 0 && dfBMYSize + 2 < INT_MAX))
1043 : {
1044 0 : CPLError(CE_Failure, CPLE_AppDefined, "Int overflow : %f x %f",
1045 : dfBMXSize, dfBMYSize);
1046 0 : return false;
1047 : }
1048 :
1049 44 : int nBMXSize = static_cast<int>(dfBMXSize);
1050 44 : int nBMYSize = static_cast<int>(dfBMYSize);
1051 :
1052 88 : if (static_cast<size_t>(1 + nBMYSize) >
1053 44 : std::numeric_limits<size_t>::max() / static_cast<size_t>(1 + nBMXSize))
1054 : {
1055 0 : CPLError(CE_Failure, CPLE_AppDefined, "Int overflow : %f x %f",
1056 : dfBMXSize, dfBMYSize);
1057 0 : return false;
1058 : }
1059 :
1060 44 : const double dfPixelXSize = (dfMaxX - dfMinX) / nBMXSize;
1061 44 : const double dfPixelYSize = (dfMaxY - dfMinY) / nBMYSize;
1062 :
1063 : // Extra pixel for right-edge and bottom-edge extensions in TOP_LEFT_CORNER
1064 : // convention.
1065 44 : nBMXSize++;
1066 44 : nBMYSize++;
1067 44 : psTransform->nBackMapWidth = nBMXSize;
1068 44 : psTransform->nBackMapHeight = nBMYSize;
1069 :
1070 44 : psTransform->adfBackMapGeoTransform[0] = dfMinX;
1071 44 : psTransform->adfBackMapGeoTransform[1] = dfPixelXSize;
1072 44 : psTransform->adfBackMapGeoTransform[2] = 0.0;
1073 44 : psTransform->adfBackMapGeoTransform[3] = dfMaxY;
1074 44 : psTransform->adfBackMapGeoTransform[4] = 0.0;
1075 44 : psTransform->adfBackMapGeoTransform[5] = -dfPixelYSize;
1076 :
1077 : /* -------------------------------------------------------------------- */
1078 : /* Allocate backmap. */
1079 : /* -------------------------------------------------------------------- */
1080 44 : auto pAccessors = static_cast<Accessors *>(psTransform->pAccessors);
1081 44 : if (!pAccessors->AllocateBackMap())
1082 0 : return false;
1083 :
1084 44 : const double dfGeorefConventionOffset =
1085 44 : psTransform->bOriginIsTopLeftCorner ? 0 : 0.5;
1086 :
1087 44 : const auto UpdateBackmap =
1088 8957995 : [&](int iBMX, int iBMY, double dfX, double dfY, double tempwt)
1089 : {
1090 822374 : const auto fBMX = pAccessors->backMapXAccessor.Get(iBMX, iBMY);
1091 822374 : const auto fBMY = pAccessors->backMapYAccessor.Get(iBMX, iBMY);
1092 822374 : const float fUpdatedBMX =
1093 : fBMX +
1094 822374 : static_cast<float>(tempwt * ((dfX + dfGeorefConventionOffset) *
1095 822374 : psTransform->dfPIXEL_STEP +
1096 822374 : psTransform->dfPIXEL_OFFSET));
1097 822374 : const float fUpdatedBMY =
1098 : fBMY +
1099 822374 : static_cast<float>(tempwt * ((dfY + dfGeorefConventionOffset) *
1100 822374 : psTransform->dfLINE_STEP +
1101 822374 : psTransform->dfLINE_OFFSET));
1102 822374 : const float fUpdatedWeight =
1103 822374 : pAccessors->backMapWeightAccessor.Get(iBMX, iBMY) +
1104 822374 : static_cast<float>(tempwt);
1105 :
1106 : // Only update the backmap if the updated averaged value results in a
1107 : // geoloc position that isn't too different from the original one.
1108 : // (there's no guarantee that if padfGeoLocX[i] ~= padfGeoLoc[j],
1109 : // padfGeoLoc[alpha * i + (1 - alpha) * j] ~= padfGeoLoc[i] )
1110 822374 : if (fUpdatedWeight > 0)
1111 : {
1112 822371 : const float fX = fUpdatedBMX / fUpdatedWeight;
1113 822371 : const float fY = fUpdatedBMY / fUpdatedWeight;
1114 822371 : const double dfGeoLocPixel =
1115 822371 : (fX - psTransform->dfPIXEL_OFFSET) / psTransform->dfPIXEL_STEP -
1116 : dfGeorefConventionOffset;
1117 822371 : const double dfGeoLocLine =
1118 822371 : (fY - psTransform->dfLINE_OFFSET) / psTransform->dfLINE_STEP -
1119 : dfGeorefConventionOffset;
1120 822371 : int iXAvg = static_cast<int>(std::max(0.0, dfGeoLocPixel));
1121 822371 : iXAvg = std::min(iXAvg, psTransform->nGeoLocXSize - 1);
1122 822371 : int iYAvg = static_cast<int>(std::max(0.0, dfGeoLocLine));
1123 822371 : iYAvg = std::min(iYAvg, psTransform->nGeoLocYSize - 1);
1124 822371 : const double dfGLX = pAccessors->geolocXAccessor.Get(iXAvg, iYAvg);
1125 822371 : const double dfGLY = pAccessors->geolocYAccessor.Get(iXAvg, iYAvg);
1126 :
1127 822371 : const unsigned iX = static_cast<unsigned>(dfX);
1128 822371 : const unsigned iY = static_cast<unsigned>(dfY);
1129 1644730 : if (!(psTransform->bHasNoData && dfGLX == psTransform->dfNoDataX) &&
1130 822361 : ((iX >= static_cast<unsigned>(nXSize - 1) ||
1131 808482 : iY >= static_cast<unsigned>(nYSize - 1)) ||
1132 798718 : (fabs(dfGLX - pAccessors->geolocXAccessor.Get(iX, iY)) <=
1133 798718 : 2 * dfPixelXSize &&
1134 756277 : fabs(dfGLY - pAccessors->geolocYAccessor.Get(iX, iY)) <=
1135 756277 : 2 * dfPixelYSize)))
1136 : {
1137 779782 : pAccessors->backMapXAccessor.Set(iBMX, iBMY, fUpdatedBMX);
1138 779782 : pAccessors->backMapYAccessor.Set(iBMX, iBMY, fUpdatedBMY);
1139 779782 : pAccessors->backMapWeightAccessor.Set(iBMX, iBMY,
1140 : fUpdatedWeight);
1141 : }
1142 : }
1143 : };
1144 :
1145 : // Keep those objects in this outer scope, so they are re-used, to
1146 : // save memory allocations.
1147 88 : OGRPoint oPoint;
1148 88 : OGRLinearRing oRing;
1149 44 : oRing.setNumPoints(5);
1150 :
1151 : /* -------------------------------------------------------------------- */
1152 : /* Run through the whole geoloc array forward projecting and */
1153 : /* pushing into the backmap. */
1154 : /* -------------------------------------------------------------------- */
1155 :
1156 : // Iterate over the (i,j) pixel space of the geolocation array, in a
1157 : // sufficiently dense way that if the geolocation array expressed an affine
1158 : // transformation, we would hit every node of the backmap.
1159 44 : const double dfStep = 1. / psTransform->dfOversampleFactor;
1160 :
1161 44 : constexpr int TILE_SIZE = GDALGeoLocDatasetAccessors::TILE_SIZE;
1162 44 : const int nYBlocks = DIV_ROUND_UP(nYSize, TILE_SIZE);
1163 44 : const int nXBlocks = DIV_ROUND_UP(nXSize, TILE_SIZE);
1164 :
1165 : // First compute for each block the start end ending floating-point
1166 : // pixel/line values
1167 88 : std::vector<std::pair<double, double>> yStartEnd(nYBlocks + 1);
1168 88 : std::vector<std::pair<double, double>> xStartEnd(nXBlocks + 1);
1169 :
1170 : {
1171 44 : int iYBlockLast = -1;
1172 44 : double dfY = -dfStep;
1173 2371 : for (; dfY <= static_cast<double>(nYSize) + 2 * dfStep; dfY += dfStep)
1174 : {
1175 2327 : const int iYBlock = static_cast<int>(dfY / TILE_SIZE);
1176 2327 : if (iYBlock != iYBlockLast)
1177 : {
1178 46 : CPLAssert(iYBlock == iYBlockLast + 1);
1179 46 : if (iYBlockLast >= 0)
1180 2 : yStartEnd[iYBlockLast].second = dfY + dfStep / 10;
1181 46 : yStartEnd[iYBlock].first = dfY;
1182 46 : iYBlockLast = iYBlock;
1183 : }
1184 : }
1185 44 : const int iYBlock = static_cast<int>(dfY / TILE_SIZE);
1186 44 : yStartEnd[iYBlock].second = dfY + dfStep / 10;
1187 : }
1188 :
1189 : {
1190 44 : int iXBlockLast = -1;
1191 44 : double dfX = -dfStep;
1192 2755 : for (; dfX <= static_cast<double>(nXSize) + 2 * dfStep; dfX += dfStep)
1193 : {
1194 2711 : const int iXBlock = static_cast<int>(dfX / TILE_SIZE);
1195 2711 : if (iXBlock != iXBlockLast)
1196 : {
1197 46 : CPLAssert(iXBlock == iXBlockLast + 1);
1198 46 : if (iXBlockLast >= 0)
1199 2 : xStartEnd[iXBlockLast].second = dfX + dfStep / 10;
1200 46 : xStartEnd[iXBlock].first = dfX;
1201 46 : iXBlockLast = iXBlock;
1202 : }
1203 : }
1204 44 : const int iXBlock = static_cast<int>(dfX / TILE_SIZE);
1205 44 : xStartEnd[iXBlock].second = dfX + dfStep / 10;
1206 : }
1207 :
1208 90 : for (int iYBlock = 0; iYBlock < nYBlocks; ++iYBlock)
1209 : {
1210 94 : for (int iXBlock = 0; iXBlock < nXBlocks; ++iXBlock)
1211 : {
1212 : #if 0
1213 : CPLDebug("Process geoloc block (y=%d,x=%d) for y in [%f, %f] and x in [%f, %f]",
1214 : iYBlock, iXBlock,
1215 : yStartEnd[iYBlock].first, yStartEnd[iYBlock].second,
1216 : xStartEnd[iXBlock].first, xStartEnd[iXBlock].second);
1217 : #endif
1218 2637 : for (double dfY = yStartEnd[iYBlock].first;
1219 2637 : dfY < yStartEnd[iYBlock].second; dfY += dfStep)
1220 : {
1221 427135 : for (double dfX = xStartEnd[iXBlock].first;
1222 427135 : dfX < xStartEnd[iXBlock].second; dfX += dfStep)
1223 : {
1224 : // Use forward geolocation array interpolation to compute
1225 : // the georeferenced position corresponding to (dfX, dfY)
1226 : double dfGeoLocX;
1227 : double dfGeoLocY;
1228 424546 : if (!PixelLineToXY(psTransform, dfX, dfY, dfGeoLocX,
1229 : dfGeoLocY))
1230 145313 : continue;
1231 :
1232 : // Compute the floating point coordinates in the pixel space
1233 : // of the backmap
1234 417864 : const double dBMX = static_cast<double>(
1235 417864 : (dfGeoLocX - dfMinX) / dfPixelXSize);
1236 :
1237 417864 : const double dBMY = static_cast<double>(
1238 417864 : (dfMaxY - dfGeoLocY) / dfPixelYSize);
1239 :
1240 : // Get top left index by truncation
1241 417864 : const int iBMX = static_cast<int>(std::floor(dBMX));
1242 417864 : const int iBMY = static_cast<int>(std::floor(dBMY));
1243 :
1244 417864 : if (iBMX >= 0 && iBMX < nBMXSize && iBMY >= 0 &&
1245 : iBMY < nBMYSize)
1246 : {
1247 : // Compute the georeferenced position of the top-left
1248 : // index of the backmap
1249 416115 : double dfGeoX = dfMinX + iBMX * dfPixelXSize;
1250 416115 : const double dfGeoY = dfMaxY - iBMY * dfPixelYSize;
1251 :
1252 416115 : bool bMatchingGeoLocCellFound = false;
1253 :
1254 416115 : const int nOuterIters =
1255 416115 : psTransform->bGeographicSRSWithMinus180Plus180LongRange &&
1256 311777 : fabs(dfGeoX) >= 180
1257 : ? 2
1258 : : 1;
1259 :
1260 832254 : for (int iOuterIter = 0; iOuterIter < nOuterIters;
1261 : ++iOuterIter)
1262 : {
1263 416139 : if (iOuterIter == 1 && dfGeoX >= 180)
1264 0 : dfGeoX -= 360;
1265 416139 : else if (iOuterIter == 1 && dfGeoX <= -180)
1266 24 : dfGeoX += 360;
1267 :
1268 : // Identify a cell (quadrilateral in georeferenced
1269 : // space) in the geolocation array in which dfGeoX,
1270 : // dfGeoY falls into.
1271 416139 : oPoint.setX(dfGeoX);
1272 416139 : oPoint.setY(dfGeoY);
1273 416139 : const int nX = static_cast<int>(std::floor(dfX));
1274 416139 : const int nY = static_cast<int>(std::floor(dfY));
1275 1196012 : for (int sx = -1;
1276 1196012 : !bMatchingGeoLocCellFound && sx <= 0; sx++)
1277 : {
1278 2289683 : for (int sy = -1;
1279 2289683 : !bMatchingGeoLocCellFound && sy <= 0; sy++)
1280 : {
1281 1511101 : const int pixel = nX + sx;
1282 1511101 : const int line = nY + sy;
1283 : double x0, y0, x1, y1, x2, y2, x3, y3;
1284 1511101 : if (!PixelLineToXY(psTransform, pixel, line,
1285 1510533 : x0, y0) ||
1286 1510533 : !PixelLineToXY(psTransform, pixel + 1,
1287 1510313 : line, x2, y2) ||
1288 1510313 : !PixelLineToXY(psTransform, pixel,
1289 3021630 : line + 1, x1, y1) ||
1290 1509883 : !PixelLineToXY(psTransform, pixel + 1,
1291 : line + 1, x3, y3))
1292 : {
1293 1288 : break;
1294 : }
1295 :
1296 1509813 : int nIters = 1;
1297 1509813 : if (psTransform
1298 1509813 : ->bGeographicSRSWithMinus180Plus180LongRange &&
1299 1178290 : std::fabs(x0) > 170 &&
1300 13548 : std::fabs(x1) > 170 &&
1301 13310 : std::fabs(x2) > 170 &&
1302 12092 : std::fabs(x3) > 170 &&
1303 11914 : (std::fabs(x1 - x0) > 180 ||
1304 11858 : std::fabs(x2 - x0) > 180 ||
1305 10496 : std::fabs(x3 - x0) > 180))
1306 : {
1307 1486 : nIters = 2;
1308 1486 : if (x0 > 0)
1309 0 : x0 -= 360;
1310 1486 : if (x1 > 0)
1311 56 : x1 -= 360;
1312 1486 : if (x2 > 0)
1313 1418 : x2 -= 360;
1314 1486 : if (x3 > 0)
1315 1486 : x3 -= 360;
1316 : }
1317 3021110 : for (int iIter = 0; iIter < nIters; ++iIter)
1318 : {
1319 1511299 : if (iIter == 1)
1320 : {
1321 1486 : x0 += 360;
1322 1486 : x1 += 360;
1323 1486 : x2 += 360;
1324 1486 : x3 += 360;
1325 : }
1326 :
1327 1511299 : oRing.setPoint(0, x0, y0);
1328 1511299 : oRing.setPoint(1, x2, y2);
1329 1511299 : oRing.setPoint(2, x3, y3);
1330 1511299 : oRing.setPoint(3, x1, y1);
1331 1511299 : oRing.setPoint(4, x0, y0);
1332 2884630 : if (oRing.isPointInRing(&oPoint) ||
1333 1373329 : oRing.isPointOnRingBoundary(
1334 : &oPoint))
1335 : {
1336 138173 : bMatchingGeoLocCellFound = true;
1337 138173 : double dfBMXValue = pixel;
1338 138173 : double dfBMYValue = line;
1339 138173 : GDALInverseBilinearInterpolation(
1340 : dfGeoX, dfGeoY, x0, y0, x1, y1,
1341 : x2, y2, x3, y3, dfBMXValue,
1342 : dfBMYValue);
1343 :
1344 138173 : dfBMXValue =
1345 138173 : (dfBMXValue +
1346 138173 : dfGeorefConventionOffset) *
1347 138173 : psTransform->dfPIXEL_STEP +
1348 138173 : psTransform->dfPIXEL_OFFSET;
1349 138173 : dfBMYValue =
1350 138173 : (dfBMYValue +
1351 138173 : dfGeorefConventionOffset) *
1352 138173 : psTransform->dfLINE_STEP +
1353 138173 : psTransform->dfLINE_OFFSET;
1354 :
1355 138173 : pAccessors->backMapXAccessor.Set(
1356 : iBMX, iBMY,
1357 : static_cast<float>(dfBMXValue));
1358 138173 : pAccessors->backMapYAccessor.Set(
1359 : iBMX, iBMY,
1360 : static_cast<float>(dfBMYValue));
1361 138173 : pAccessors->backMapWeightAccessor
1362 : .Set(iBMX, iBMY, 1.0f);
1363 : }
1364 : }
1365 : }
1366 : }
1367 : }
1368 416115 : if (bMatchingGeoLocCellFound)
1369 138173 : continue;
1370 : }
1371 :
1372 : // We will end up here in non-nominal cases, with nodata,
1373 : // holes, etc.
1374 :
1375 : // Check if the center is in range
1376 279691 : if (iBMX < -1 || iBMY < -1 || iBMX > nBMXSize ||
1377 : iBMY > nBMYSize)
1378 458 : continue;
1379 :
1380 279233 : const double fracBMX = dBMX - iBMX;
1381 279233 : const double fracBMY = dBMY - iBMY;
1382 :
1383 : // Check logic for top left pixel
1384 278929 : if ((iBMX >= 0) && (iBMY >= 0) && (iBMX < nBMXSize) &&
1385 558162 : (iBMY < nBMYSize) &&
1386 277942 : pAccessors->backMapWeightAccessor.Get(iBMX, iBMY) !=
1387 : 1.0f)
1388 : {
1389 185775 : const double tempwt = (1.0 - fracBMX) * (1.0 - fracBMY);
1390 185775 : UpdateBackmap(iBMX, iBMY, dfX, dfY, tempwt);
1391 : }
1392 :
1393 : // Check logic for top right pixel
1394 278830 : if ((iBMY >= 0) && (iBMX + 1 < nBMXSize) &&
1395 558063 : (iBMY < nBMYSize) &&
1396 278190 : pAccessors->backMapWeightAccessor.Get(iBMX + 1, iBMY) !=
1397 : 1.0f)
1398 : {
1399 189705 : const double tempwt = fracBMX * (1.0 - fracBMY);
1400 189705 : UpdateBackmap(iBMX + 1, iBMY, dfX, dfY, tempwt);
1401 : }
1402 :
1403 : // Check logic for bottom right pixel
1404 557747 : if ((iBMX + 1 < nBMXSize) && (iBMY + 1 < nBMYSize) &&
1405 278514 : pAccessors->backMapWeightAccessor.Get(iBMX + 1,
1406 278514 : iBMY + 1) != 1.0f)
1407 : {
1408 226485 : const double tempwt = fracBMX * fracBMY;
1409 226485 : UpdateBackmap(iBMX + 1, iBMY + 1, dfX, dfY, tempwt);
1410 : }
1411 :
1412 : // Check logic for bottom left pixel
1413 278929 : if ((iBMX >= 0) && (iBMX < nBMXSize) &&
1414 836440 : (iBMY + 1 < nBMYSize) &&
1415 278278 : pAccessors->backMapWeightAccessor.Get(iBMX, iBMY + 1) !=
1416 : 1.0f)
1417 : {
1418 220409 : const double tempwt = (1.0 - fracBMX) * fracBMY;
1419 220409 : UpdateBackmap(iBMX, iBMY + 1, dfX, dfY, tempwt);
1420 : }
1421 : }
1422 : }
1423 : }
1424 : }
1425 :
1426 : // Each pixel in the backmap may have multiple entries.
1427 : // We now go in average it out using the weights
1428 148 : START_ITER_PER_BLOCK(nBMXSize, TILE_SIZE, nBMYSize, TILE_SIZE, (void)0,
1429 : iXStart, iXEnd, iYStart, iYEnd)
1430 : {
1431 2268 : for (int iY = iYStart; iY < iYEnd; ++iY)
1432 : {
1433 327929 : for (int iX = iXStart; iX < iXEnd; ++iX)
1434 : {
1435 : // Check if pixel was only touch during neighbor scan
1436 : // But no real weight was added as source point matched
1437 : // backmap grid node
1438 465470 : const auto weight =
1439 325721 : pAccessors->backMapWeightAccessor.Get(iX, iY);
1440 325721 : if (weight > 0)
1441 : {
1442 236251 : pAccessors->backMapXAccessor.Set(
1443 : iX, iY,
1444 170330 : pAccessors->backMapXAccessor.Get(iX, iY) / weight);
1445 236251 : pAccessors->backMapYAccessor.Set(
1446 : iX, iY,
1447 170330 : pAccessors->backMapYAccessor.Get(iX, iY) / weight);
1448 : }
1449 : else
1450 : {
1451 155391 : pAccessors->backMapXAccessor.Set(iX, iY, INVALID_BMXY);
1452 155391 : pAccessors->backMapYAccessor.Set(iX, iY, INVALID_BMXY);
1453 : }
1454 : }
1455 : }
1456 : }
1457 : END_ITER_PER_BLOCK
1458 :
1459 44 : pAccessors->FreeWghtsBackMap();
1460 :
1461 : // Fill holes in backmap
1462 44 : auto poBackmapDS = pAccessors->GetBackmapDataset();
1463 :
1464 44 : pAccessors->FlushBackmapCaches();
1465 :
1466 : #ifdef DEBUG_GEOLOC
1467 : if (CPLTestBool(CPLGetConfigOption("GEOLOC_DUMP", "NO")))
1468 : {
1469 : poBackmapDS->SetGeoTransform(psTransform->adfBackMapGeoTransform);
1470 : GDALClose(GDALCreateCopy(GDALGetDriverByName("GTiff"),
1471 : "/tmp/geoloc_before_fill.tif", poBackmapDS,
1472 : false, nullptr, nullptr, nullptr));
1473 : }
1474 : #endif
1475 :
1476 44 : constexpr double dfMaxSearchDist = 3.0;
1477 44 : constexpr int nSmoothingIterations = 1;
1478 132 : for (int i = 1; i <= 2; i++)
1479 : {
1480 88 : GDALFillNodata(GDALRasterBand::ToHandle(poBackmapDS->GetRasterBand(i)),
1481 : nullptr, dfMaxSearchDist,
1482 : 0, // unused parameter
1483 : nSmoothingIterations, nullptr, nullptr, nullptr);
1484 : }
1485 :
1486 : #ifdef DEBUG_GEOLOC
1487 : if (CPLTestBool(CPLGetConfigOption("GEOLOC_DUMP", "NO")))
1488 : {
1489 : GDALClose(GDALCreateCopy(GDALGetDriverByName("GTiff"),
1490 : "/tmp/geoloc_after_fill.tif", poBackmapDS,
1491 : false, nullptr, nullptr, nullptr));
1492 : }
1493 : #endif
1494 :
1495 : // A final hole filling logic, proceeding line by line, and filling
1496 : // holes when the backmap values surrounding the hole are close enough.
1497 : struct LastValidStruct
1498 : {
1499 : int iX = -1;
1500 : float bmX = 0;
1501 : };
1502 :
1503 44 : std::vector<LastValidStruct> lastValid(TILE_SIZE);
1504 220 : const auto reinitLine = [&lastValid]()
1505 : {
1506 44 : const size_t nSize = lastValid.size();
1507 44 : lastValid.clear();
1508 44 : lastValid.resize(nSize);
1509 : };
1510 148 : START_ITER_PER_BLOCK(nBMXSize, TILE_SIZE, nBMYSize, TILE_SIZE, reinitLine(),
1511 : iXStart, iXEnd, iYStart, iYEnd)
1512 : {
1513 60 : const int iYCount = iYEnd - iYStart;
1514 2268 : for (int iYIter = 0; iYIter < iYCount; ++iYIter)
1515 : {
1516 2208 : int iLastValidIX = lastValid[iYIter].iX;
1517 2208 : float bmXLastValid = lastValid[iYIter].bmX;
1518 2208 : const int iBMY = iYStart + iYIter;
1519 327929 : for (int iBMX = iXStart; iBMX < iXEnd; ++iBMX)
1520 : {
1521 325721 : const float bmX = pAccessors->backMapXAccessor.Get(iBMX, iBMY);
1522 325721 : if (bmX == INVALID_BMXY)
1523 127144 : continue;
1524 199129 : if (iLastValidIX != -1 && iBMX > iLastValidIX + 1 &&
1525 552 : fabs(bmX - bmXLastValid) <= 2)
1526 : {
1527 354 : const float bmY =
1528 239 : pAccessors->backMapYAccessor.Get(iBMX, iBMY);
1529 354 : const float bmYLastValid =
1530 239 : pAccessors->backMapYAccessor.Get(iLastValidIX, iBMY);
1531 239 : if (fabs(bmY - bmYLastValid) <= 2)
1532 : {
1533 266 : for (int iBMXInner = iLastValidIX + 1; iBMXInner < iBMX;
1534 : ++iBMXInner)
1535 : {
1536 145 : const float alpha =
1537 145 : static_cast<float>(iBMXInner - iLastValidIX) /
1538 145 : (iBMX - iLastValidIX);
1539 145 : pAccessors->backMapXAccessor.Set(
1540 : iBMXInner, iBMY,
1541 145 : (1.0f - alpha) * bmXLastValid + alpha * bmX);
1542 145 : pAccessors->backMapYAccessor.Set(
1543 : iBMXInner, iBMY,
1544 145 : (1.0f - alpha) * bmYLastValid + alpha * bmY);
1545 : }
1546 : }
1547 : }
1548 198577 : iLastValidIX = iBMX;
1549 198577 : bmXLastValid = bmX;
1550 : }
1551 2208 : lastValid[iYIter].iX = iLastValidIX;
1552 2208 : lastValid[iYIter].bmX = bmXLastValid;
1553 : }
1554 : }
1555 : END_ITER_PER_BLOCK
1556 :
1557 : #ifdef DEBUG_GEOLOC
1558 : if (CPLTestBool(CPLGetConfigOption("GEOLOC_DUMP", "NO")))
1559 : {
1560 : pAccessors->FlushBackmapCaches();
1561 :
1562 : GDALClose(GDALCreateCopy(GDALGetDriverByName("GTiff"),
1563 : "/tmp/geoloc_after_line_fill.tif", poBackmapDS,
1564 : false, nullptr, nullptr, nullptr));
1565 : }
1566 : #endif
1567 :
1568 44 : pAccessors->ReleaseBackmapDataset(poBackmapDS);
1569 44 : CPLDebug("GEOLOC", "Ending backmap generation");
1570 :
1571 44 : return true;
1572 : }
1573 :
1574 : /*! @endcond */
1575 :
1576 : /************************************************************************/
1577 : /* GDALGeoLocRescale() */
1578 : /************************************************************************/
1579 :
1580 4 : static void GDALGeoLocRescale(char **&papszMD, const char *pszItem,
1581 : double dfRatio, double dfDefaultVal)
1582 : {
1583 : const double dfVal =
1584 4 : dfRatio * CPLAtofM(CSLFetchNameValueDef(
1585 4 : papszMD, pszItem, CPLSPrintf("%.17g", dfDefaultVal)));
1586 :
1587 4 : papszMD = CSLSetNameValue(papszMD, pszItem, CPLSPrintf("%.17g", dfVal));
1588 4 : }
1589 :
1590 : /************************************************************************/
1591 : /* GDALCreateSimilarGeoLocTransformer() */
1592 : /************************************************************************/
1593 :
1594 1 : static void *GDALCreateSimilarGeoLocTransformer(void *hTransformArg,
1595 : double dfRatioX,
1596 : double dfRatioY)
1597 : {
1598 1 : VALIDATE_POINTER1(hTransformArg, "GDALCreateSimilarGeoLocTransformer",
1599 : nullptr);
1600 :
1601 1 : GDALGeoLocTransformInfo *psInfo =
1602 : static_cast<GDALGeoLocTransformInfo *>(hTransformArg);
1603 :
1604 1 : char **papszGeolocationInfo = CSLDuplicate(psInfo->papszGeolocationInfo);
1605 :
1606 1 : if (dfRatioX != 1.0 || dfRatioY != 1.0)
1607 : {
1608 1 : GDALGeoLocRescale(papszGeolocationInfo, "PIXEL_OFFSET", dfRatioX, 0.0);
1609 1 : GDALGeoLocRescale(papszGeolocationInfo, "LINE_OFFSET", dfRatioY, 0.0);
1610 1 : GDALGeoLocRescale(papszGeolocationInfo, "PIXEL_STEP", 1.0 / dfRatioX,
1611 : 1.0);
1612 1 : GDALGeoLocRescale(papszGeolocationInfo, "LINE_STEP", 1.0 / dfRatioY,
1613 : 1.0);
1614 : }
1615 :
1616 : auto psInfoNew =
1617 2 : static_cast<GDALGeoLocTransformInfo *>(GDALCreateGeoLocTransformer(
1618 1 : nullptr, papszGeolocationInfo, psInfo->bReversed));
1619 1 : psInfoNew->dfOversampleFactor = psInfo->dfOversampleFactor;
1620 :
1621 1 : CSLDestroy(papszGeolocationInfo);
1622 :
1623 1 : return psInfoNew;
1624 : }
1625 :
1626 : /************************************************************************/
1627 : /* GDALCreateGeolocationMetadata() */
1628 : /************************************************************************/
1629 :
1630 : /** Synthetize the content of a GEOLOCATION metadata domain from a
1631 : * geolocation dataset.
1632 : *
1633 : * This is used when doing gdalwarp -to GEOLOC_ARRAY=some.tif
1634 : */
1635 10 : CPLStringList GDALCreateGeolocationMetadata(GDALDatasetH hBaseDS,
1636 : const char *pszGeolocationDataset,
1637 : bool bIsSource)
1638 : {
1639 20 : CPLStringList aosMD;
1640 :
1641 : // Open geolocation dataset
1642 : auto poGeolocDS = std::unique_ptr<GDALDataset>(
1643 20 : GDALDataset::Open(pszGeolocationDataset, GDAL_OF_RASTER));
1644 10 : if (poGeolocDS == nullptr)
1645 : {
1646 2 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid dataset: %s",
1647 : pszGeolocationDataset);
1648 2 : return CPLStringList();
1649 : }
1650 8 : const int nGeoLocXSize = poGeolocDS->GetRasterXSize();
1651 8 : const int nGeoLocYSize = poGeolocDS->GetRasterYSize();
1652 8 : if (nGeoLocXSize == 0 || nGeoLocYSize == 0)
1653 : {
1654 0 : CPLError(CE_Failure, CPLE_AppDefined,
1655 : "Invalid dataset dimension for %s: %dx%d",
1656 : pszGeolocationDataset, nGeoLocXSize, nGeoLocYSize);
1657 0 : return CPLStringList();
1658 : }
1659 :
1660 : // Import the GEOLOCATION metadata from the geolocation dataset, if existing
1661 8 : auto papszGeolocMDFromGeolocDS = poGeolocDS->GetMetadata("GEOLOCATION");
1662 8 : if (papszGeolocMDFromGeolocDS)
1663 3 : aosMD = CSLDuplicate(papszGeolocMDFromGeolocDS);
1664 :
1665 8 : aosMD.SetNameValue("X_DATASET", pszGeolocationDataset);
1666 8 : aosMD.SetNameValue("Y_DATASET", pszGeolocationDataset);
1667 :
1668 : // Set X_BAND, Y_BAND to 1, 2 if they are not specified in the initial
1669 : // GEOLOC metadata domain.and the geolocation dataset has 2 bands.
1670 13 : if (aosMD.FetchNameValue("X_BAND") == nullptr &&
1671 5 : aosMD.FetchNameValue("Y_BAND") == nullptr)
1672 : {
1673 5 : if (poGeolocDS->GetRasterCount() != 2)
1674 : {
1675 1 : CPLError(CE_Failure, CPLE_AppDefined,
1676 : "Expected 2 bands for %s. Got %d", pszGeolocationDataset,
1677 : poGeolocDS->GetRasterCount());
1678 1 : return CPLStringList();
1679 : }
1680 4 : aosMD.SetNameValue("X_BAND", "1");
1681 4 : aosMD.SetNameValue("Y_BAND", "2");
1682 : }
1683 :
1684 : // Set the geoloc SRS from the geolocation dataset SRS if there's no
1685 : // explicit one in the initial GEOLOC metadata domain.
1686 7 : if (aosMD.FetchNameValue("SRS") == nullptr)
1687 : {
1688 7 : auto poSRS = poGeolocDS->GetSpatialRef();
1689 7 : if (poSRS)
1690 : {
1691 3 : char *pszWKT = nullptr;
1692 3 : poSRS->exportToWkt(&pszWKT);
1693 3 : aosMD.SetNameValue("SRS", pszWKT);
1694 3 : CPLFree(pszWKT);
1695 : }
1696 : }
1697 7 : if (aosMD.FetchNameValue("SRS") == nullptr)
1698 : {
1699 4 : aosMD.SetNameValue("SRS", SRS_WKT_WGS84_LAT_LONG);
1700 : }
1701 :
1702 : // Set default values for PIXEL/LINE_OFFSET/STEP if not present.
1703 7 : if (aosMD.FetchNameValue("PIXEL_OFFSET") == nullptr)
1704 4 : aosMD.SetNameValue("PIXEL_OFFSET", "0");
1705 :
1706 7 : if (aosMD.FetchNameValue("LINE_OFFSET") == nullptr)
1707 4 : aosMD.SetNameValue("LINE_OFFSET", "0");
1708 :
1709 7 : if (aosMD.FetchNameValue("PIXEL_STEP") == nullptr)
1710 : {
1711 : aosMD.SetNameValue(
1712 4 : "PIXEL_STEP", CPLSPrintf("%.17g", static_cast<double>(
1713 4 : GDALGetRasterXSize(hBaseDS)) /
1714 4 : nGeoLocXSize));
1715 : }
1716 :
1717 7 : if (aosMD.FetchNameValue("LINE_STEP") == nullptr)
1718 : {
1719 : aosMD.SetNameValue(
1720 4 : "LINE_STEP", CPLSPrintf("%.17g", static_cast<double>(
1721 4 : GDALGetRasterYSize(hBaseDS)) /
1722 4 : nGeoLocYSize));
1723 : }
1724 :
1725 7 : if (aosMD.FetchNameValue("GEOREFERENCING_CONVENTION") == nullptr)
1726 : {
1727 : const char *pszConvention =
1728 7 : poGeolocDS->GetMetadataItem("GEOREFERENCING_CONVENTION");
1729 7 : if (pszConvention)
1730 2 : aosMD.SetNameValue("GEOREFERENCING_CONVENTION", pszConvention);
1731 : }
1732 :
1733 14 : std::string osDebugMsg;
1734 7 : osDebugMsg = "Synthetized GEOLOCATION metadata for ";
1735 7 : osDebugMsg += bIsSource ? "source" : "target";
1736 7 : osDebugMsg += ":\n";
1737 72 : for (int i = 0; i < aosMD.size(); ++i)
1738 : {
1739 65 : osDebugMsg += " ";
1740 65 : osDebugMsg += aosMD[i];
1741 65 : osDebugMsg += '\n';
1742 : }
1743 7 : CPLDebug("GEOLOC", "%s", osDebugMsg.c_str());
1744 :
1745 7 : return aosMD;
1746 : }
1747 :
1748 : /************************************************************************/
1749 : /* GDALCreateGeoLocTransformer() */
1750 : /************************************************************************/
1751 :
1752 50 : void *GDALCreateGeoLocTransformerEx(GDALDatasetH hBaseDS,
1753 : CSLConstList papszGeolocationInfo,
1754 : int bReversed, const char *pszSourceDataset,
1755 : CSLConstList papszTransformOptions)
1756 :
1757 : {
1758 :
1759 50 : if (CSLFetchNameValue(papszGeolocationInfo, "PIXEL_OFFSET") == nullptr ||
1760 48 : CSLFetchNameValue(papszGeolocationInfo, "LINE_OFFSET") == nullptr ||
1761 48 : CSLFetchNameValue(papszGeolocationInfo, "PIXEL_STEP") == nullptr ||
1762 48 : CSLFetchNameValue(papszGeolocationInfo, "LINE_STEP") == nullptr ||
1763 146 : CSLFetchNameValue(papszGeolocationInfo, "X_BAND") == nullptr ||
1764 48 : CSLFetchNameValue(papszGeolocationInfo, "Y_BAND") == nullptr)
1765 : {
1766 2 : CPLError(CE_Failure, CPLE_AppDefined,
1767 : "Missing some geolocation fields in "
1768 : "GDALCreateGeoLocTransformer()");
1769 2 : return nullptr;
1770 : }
1771 :
1772 : /* -------------------------------------------------------------------- */
1773 : /* Initialize core info. */
1774 : /* -------------------------------------------------------------------- */
1775 : GDALGeoLocTransformInfo *psTransform =
1776 : static_cast<GDALGeoLocTransformInfo *>(
1777 48 : CPLCalloc(sizeof(GDALGeoLocTransformInfo), 1));
1778 :
1779 48 : psTransform->bReversed = CPL_TO_BOOL(bReversed);
1780 48 : psTransform->dfOversampleFactor = std::max(
1781 96 : 0.1,
1782 96 : std::min(2.0,
1783 48 : CPLAtof(CSLFetchNameValueDef(
1784 : papszTransformOptions, "GEOLOC_BACKMAP_OVERSAMPLE_FACTOR",
1785 : CPLGetConfigOption("GDAL_GEOLOC_BACKMAP_OVERSAMPLE_FACTOR",
1786 48 : "1.3")))));
1787 :
1788 48 : memcpy(psTransform->sTI.abySignature, GDAL_GTI2_SIGNATURE,
1789 : strlen(GDAL_GTI2_SIGNATURE));
1790 48 : psTransform->sTI.pszClassName = "GDALGeoLocTransformer";
1791 48 : psTransform->sTI.pfnTransform = GDALGeoLocTransform;
1792 48 : psTransform->sTI.pfnCleanup = GDALDestroyGeoLocTransformer;
1793 48 : psTransform->sTI.pfnSerialize = GDALSerializeGeoLocTransformer;
1794 48 : psTransform->sTI.pfnCreateSimilar = GDALCreateSimilarGeoLocTransformer;
1795 :
1796 48 : psTransform->papszGeolocationInfo = CSLDuplicate(papszGeolocationInfo);
1797 :
1798 : /* -------------------------------------------------------------------- */
1799 : /* Pull geolocation info from the options/metadata. */
1800 : /* -------------------------------------------------------------------- */
1801 48 : psTransform->dfPIXEL_OFFSET =
1802 48 : CPLAtof(CSLFetchNameValue(papszGeolocationInfo, "PIXEL_OFFSET"));
1803 48 : psTransform->dfLINE_OFFSET =
1804 48 : CPLAtof(CSLFetchNameValue(papszGeolocationInfo, "LINE_OFFSET"));
1805 48 : psTransform->dfPIXEL_STEP =
1806 48 : CPLAtof(CSLFetchNameValue(papszGeolocationInfo, "PIXEL_STEP"));
1807 48 : psTransform->dfLINE_STEP =
1808 48 : CPLAtof(CSLFetchNameValue(papszGeolocationInfo, "LINE_STEP"));
1809 :
1810 48 : psTransform->bOriginIsTopLeftCorner = EQUAL(
1811 : CSLFetchNameValueDef(papszGeolocationInfo, "GEOREFERENCING_CONVENTION",
1812 : "TOP_LEFT_CORNER"),
1813 : "TOP_LEFT_CORNER");
1814 :
1815 : /* -------------------------------------------------------------------- */
1816 : /* Establish access to geolocation dataset(s). */
1817 : /* -------------------------------------------------------------------- */
1818 : const char *pszDSName =
1819 48 : CSLFetchNameValue(papszGeolocationInfo, "X_DATASET");
1820 48 : if (pszDSName != nullptr)
1821 : {
1822 96 : CPLConfigOptionSetter oSetter("CPL_ALLOW_VSISTDIN", "NO", true);
1823 48 : if (CPLTestBool(CSLFetchNameValueDef(
1824 49 : papszGeolocationInfo, "X_DATASET_RELATIVE_TO_SOURCE", "NO")) &&
1825 1 : (hBaseDS != nullptr || pszSourceDataset))
1826 : {
1827 8 : const CPLString osFilename = CPLProjectRelativeFilenameSafe(
1828 7 : CPLGetDirnameSafe(pszSourceDataset
1829 : ? pszSourceDataset
1830 3 : : GDALGetDescription(hBaseDS))
1831 : .c_str(),
1832 4 : pszDSName);
1833 4 : psTransform->hDS_X =
1834 4 : GDALOpenShared(osFilename.c_str(), GA_ReadOnly);
1835 : }
1836 : else
1837 : {
1838 44 : psTransform->hDS_X = GDALOpenShared(pszDSName, GA_ReadOnly);
1839 : }
1840 : }
1841 : else
1842 : {
1843 0 : psTransform->hDS_X = hBaseDS;
1844 0 : if (hBaseDS)
1845 : {
1846 0 : GDALReferenceDataset(psTransform->hDS_X);
1847 0 : psTransform->papszGeolocationInfo =
1848 0 : CSLSetNameValue(psTransform->papszGeolocationInfo, "X_DATASET",
1849 : GDALGetDescription(hBaseDS));
1850 : }
1851 : }
1852 :
1853 48 : pszDSName = CSLFetchNameValue(papszGeolocationInfo, "Y_DATASET");
1854 48 : if (pszDSName != nullptr)
1855 : {
1856 96 : CPLConfigOptionSetter oSetter("CPL_ALLOW_VSISTDIN", "NO", true);
1857 48 : if (CPLTestBool(CSLFetchNameValueDef(
1858 49 : papszGeolocationInfo, "Y_DATASET_RELATIVE_TO_SOURCE", "NO")) &&
1859 1 : (hBaseDS != nullptr || pszSourceDataset))
1860 : {
1861 8 : const CPLString osFilename = CPLProjectRelativeFilenameSafe(
1862 7 : CPLGetDirnameSafe(pszSourceDataset
1863 : ? pszSourceDataset
1864 3 : : GDALGetDescription(hBaseDS))
1865 : .c_str(),
1866 4 : pszDSName);
1867 4 : psTransform->hDS_Y =
1868 4 : GDALOpenShared(osFilename.c_str(), GA_ReadOnly);
1869 : }
1870 : else
1871 : {
1872 44 : psTransform->hDS_Y = GDALOpenShared(pszDSName, GA_ReadOnly);
1873 : }
1874 : }
1875 : else
1876 : {
1877 0 : psTransform->hDS_Y = hBaseDS;
1878 0 : if (hBaseDS)
1879 : {
1880 0 : GDALReferenceDataset(psTransform->hDS_Y);
1881 0 : psTransform->papszGeolocationInfo =
1882 0 : CSLSetNameValue(psTransform->papszGeolocationInfo, "Y_DATASET",
1883 : GDALGetDescription(hBaseDS));
1884 : }
1885 : }
1886 :
1887 48 : if (psTransform->hDS_X == nullptr || psTransform->hDS_Y == nullptr)
1888 : {
1889 0 : GDALDestroyGeoLocTransformer(psTransform);
1890 0 : return nullptr;
1891 : }
1892 :
1893 : /* -------------------------------------------------------------------- */
1894 : /* Get the band handles. */
1895 : /* -------------------------------------------------------------------- */
1896 : const int nXBand =
1897 48 : std::max(1, atoi(CSLFetchNameValue(papszGeolocationInfo, "X_BAND")));
1898 48 : psTransform->hBand_X = GDALGetRasterBand(psTransform->hDS_X, nXBand);
1899 :
1900 48 : psTransform->dfNoDataX = GDALGetRasterNoDataValue(
1901 : psTransform->hBand_X, &(psTransform->bHasNoData));
1902 :
1903 : const int nYBand =
1904 48 : std::max(1, atoi(CSLFetchNameValue(papszGeolocationInfo, "Y_BAND")));
1905 48 : psTransform->hBand_Y = GDALGetRasterBand(psTransform->hDS_Y, nYBand);
1906 :
1907 48 : if (psTransform->hBand_X == nullptr || psTransform->hBand_Y == nullptr)
1908 : {
1909 0 : GDALDestroyGeoLocTransformer(psTransform);
1910 0 : return nullptr;
1911 : }
1912 :
1913 48 : psTransform->bSwapXY = CPLTestBool(
1914 : CSLFetchNameValueDef(papszGeolocationInfo, "SWAP_XY", "NO"));
1915 :
1916 : /* -------------------------------------------------------------------- */
1917 : /* Check that X and Y bands have the same dimensions */
1918 : /* -------------------------------------------------------------------- */
1919 48 : const int nXSize_XBand = GDALGetRasterXSize(psTransform->hDS_X);
1920 48 : const int nYSize_XBand = GDALGetRasterYSize(psTransform->hDS_X);
1921 48 : const int nXSize_YBand = GDALGetRasterXSize(psTransform->hDS_Y);
1922 48 : const int nYSize_YBand = GDALGetRasterYSize(psTransform->hDS_Y);
1923 48 : if (nYSize_XBand == 1 || nYSize_YBand == 1)
1924 : {
1925 15 : if (nYSize_XBand != 1 || nYSize_YBand != 1)
1926 : {
1927 0 : CPLError(CE_Failure, CPLE_AppDefined,
1928 : "X_BAND and Y_BAND should have both nYSize == 1");
1929 0 : GDALDestroyGeoLocTransformer(psTransform);
1930 0 : return nullptr;
1931 : }
1932 : }
1933 33 : else if (nXSize_XBand != nXSize_YBand || nYSize_XBand != nYSize_YBand)
1934 : {
1935 0 : CPLError(CE_Failure, CPLE_AppDefined,
1936 : "X_BAND and Y_BAND do not have the same dimensions");
1937 0 : GDALDestroyGeoLocTransformer(psTransform);
1938 0 : return nullptr;
1939 : }
1940 :
1941 48 : if (nXSize_XBand <= 0 || nYSize_XBand <= 0 || nXSize_YBand <= 0 ||
1942 : nYSize_YBand <= 0)
1943 : {
1944 0 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid X_BAND / Y_BAND size");
1945 0 : GDALDestroyGeoLocTransformer(psTransform);
1946 0 : return nullptr;
1947 : }
1948 :
1949 : // Is it a regular grid ? That is:
1950 : // The XBAND contains the x coordinates for all lines.
1951 : // The YBAND contains the y coordinates for all columns.
1952 48 : const bool bIsRegularGrid = (nYSize_XBand == 1 && nYSize_YBand == 1);
1953 :
1954 48 : const int nXSize = nXSize_XBand;
1955 48 : const int nYSize = bIsRegularGrid ? nXSize_YBand : nYSize_XBand;
1956 :
1957 96 : if (static_cast<size_t>(nXSize) >
1958 48 : std::numeric_limits<size_t>::max() / static_cast<size_t>(nYSize))
1959 : {
1960 0 : CPLError(CE_Failure, CPLE_AppDefined, "Int overflow : %d x %d", nXSize,
1961 : nYSize);
1962 0 : GDALDestroyGeoLocTransformer(psTransform);
1963 0 : return nullptr;
1964 : }
1965 :
1966 48 : psTransform->nGeoLocXSize = nXSize;
1967 48 : psTransform->nGeoLocYSize = nYSize;
1968 :
1969 48 : if (hBaseDS && psTransform->dfPIXEL_OFFSET == 0 &&
1970 44 : psTransform->dfLINE_OFFSET == 0 && psTransform->dfPIXEL_STEP == 1 &&
1971 37 : psTransform->dfLINE_STEP == 1)
1972 : {
1973 73 : if (GDALGetRasterXSize(hBaseDS) > nXSize ||
1974 36 : GDALGetRasterYSize(hBaseDS) > nYSize)
1975 : {
1976 2 : CPLError(CE_Warning, CPLE_AppDefined,
1977 : "Geolocation array is %d x %d large, "
1978 : "whereas dataset is %d x %d large. Result might be "
1979 : "incorrect due to lack of values in geolocation array.",
1980 : nXSize, nYSize, GDALGetRasterXSize(hBaseDS),
1981 : GDALGetRasterYSize(hBaseDS));
1982 : }
1983 : }
1984 :
1985 : /* -------------------------------------------------------------------- */
1986 : /* Load the geolocation array. */
1987 : /* -------------------------------------------------------------------- */
1988 :
1989 : // The quadtree method is experimental. It simplifies the code
1990 : // significantly, but unfortunately burns more RAM and is slower.
1991 : const bool bUseQuadtree =
1992 48 : EQUAL(CPLGetConfigOption("GDAL_GEOLOC_INVERSE_METHOD", "BACKMAP"),
1993 : "QUADTREE");
1994 :
1995 : // Decide if we should C-arrays for geoloc and backmap, or on-disk
1996 : // temporary datasets.
1997 48 : const char *pszUseTempDatasets = CSLFetchNameValueDef(
1998 : papszTransformOptions, "GEOLOC_USE_TEMP_DATASETS",
1999 : CPLGetConfigOption("GDAL_GEOLOC_USE_TEMP_DATASETS", nullptr));
2000 48 : if (pszUseTempDatasets)
2001 4 : psTransform->bUseArray = !CPLTestBool(pszUseTempDatasets);
2002 : else
2003 : {
2004 44 : constexpr int MEGAPIXEL_LIMIT = 24;
2005 44 : psTransform->bUseArray =
2006 44 : nXSize < MEGAPIXEL_LIMIT * 1000 * 1000 / nYSize;
2007 44 : if (!psTransform->bUseArray)
2008 : {
2009 0 : CPLDebug("GEOLOC",
2010 : "Using temporary GTiff backing to store backmap, because "
2011 : "geoloc arrays require %d megapixels, exceeding the %d "
2012 : "megapixels limit. You can set the "
2013 : "GDAL_GEOLOC_USE_TEMP_DATASETS configuration option to "
2014 : "NO to force RAM storage of backmap",
2015 0 : static_cast<int>(static_cast<int64_t>(nXSize) * nYSize /
2016 : (1000 * 1000)),
2017 : MEGAPIXEL_LIMIT);
2018 : }
2019 : }
2020 :
2021 48 : if (psTransform->bUseArray)
2022 : {
2023 46 : auto pAccessors = new GDALGeoLocCArrayAccessors(psTransform);
2024 46 : psTransform->pAccessors = pAccessors;
2025 46 : if (!pAccessors->Load(bIsRegularGrid, bUseQuadtree))
2026 : {
2027 0 : GDALDestroyGeoLocTransformer(psTransform);
2028 0 : return nullptr;
2029 : }
2030 : }
2031 : else
2032 : {
2033 2 : auto pAccessors = new GDALGeoLocDatasetAccessors(psTransform);
2034 2 : psTransform->pAccessors = pAccessors;
2035 2 : if (!pAccessors->Load(bIsRegularGrid, bUseQuadtree))
2036 : {
2037 0 : GDALDestroyGeoLocTransformer(psTransform);
2038 0 : return nullptr;
2039 : }
2040 : }
2041 :
2042 48 : return psTransform;
2043 : }
2044 :
2045 : /** Create GeoLocation transformer */
2046 1 : void *GDALCreateGeoLocTransformer(GDALDatasetH hBaseDS,
2047 : char **papszGeolocationInfo, int bReversed)
2048 :
2049 : {
2050 1 : return GDALCreateGeoLocTransformerEx(hBaseDS, papszGeolocationInfo,
2051 1 : bReversed, nullptr, nullptr);
2052 : }
2053 :
2054 : /************************************************************************/
2055 : /* GDALDestroyGeoLocTransformer() */
2056 : /************************************************************************/
2057 :
2058 : /** Destroy GeoLocation transformer */
2059 48 : void GDALDestroyGeoLocTransformer(void *pTransformAlg)
2060 :
2061 : {
2062 48 : if (pTransformAlg == nullptr)
2063 0 : return;
2064 :
2065 48 : GDALGeoLocTransformInfo *psTransform =
2066 : static_cast<GDALGeoLocTransformInfo *>(pTransformAlg);
2067 :
2068 48 : CSLDestroy(psTransform->papszGeolocationInfo);
2069 :
2070 48 : if (psTransform->bUseArray)
2071 : delete static_cast<GDALGeoLocCArrayAccessors *>(
2072 46 : psTransform->pAccessors);
2073 : else
2074 : delete static_cast<GDALGeoLocDatasetAccessors *>(
2075 2 : psTransform->pAccessors);
2076 :
2077 96 : if (psTransform->hDS_X != nullptr &&
2078 48 : GDALDereferenceDataset(psTransform->hDS_X) == 0)
2079 27 : GDALClose(psTransform->hDS_X);
2080 :
2081 96 : if (psTransform->hDS_Y != nullptr &&
2082 48 : GDALDereferenceDataset(psTransform->hDS_Y) == 0)
2083 40 : GDALClose(psTransform->hDS_Y);
2084 :
2085 48 : if (psTransform->hQuadTree != nullptr)
2086 4 : CPLQuadTreeDestroy(psTransform->hQuadTree);
2087 :
2088 48 : CPLFree(pTransformAlg);
2089 : }
2090 :
2091 : /** Use GeoLocation transformer */
2092 31164 : int GDALGeoLocTransform(void *pTransformArg, int bDstToSrc, int nPointCount,
2093 : double *padfX, double *padfY, double *padfZ,
2094 : int *panSuccess)
2095 : {
2096 31164 : GDALGeoLocTransformInfo *psTransform =
2097 : static_cast<GDALGeoLocTransformInfo *>(pTransformArg);
2098 31164 : if (psTransform->bUseArray)
2099 : {
2100 29043 : return GDALGeoLoc<GDALGeoLocCArrayAccessors>::Transform(
2101 : pTransformArg, bDstToSrc, nPointCount, padfX, padfY, padfZ,
2102 29043 : panSuccess);
2103 : }
2104 : else
2105 : {
2106 2121 : return GDALGeoLoc<GDALGeoLocDatasetAccessors>::Transform(
2107 : pTransformArg, bDstToSrc, nPointCount, padfX, padfY, padfZ,
2108 2121 : panSuccess);
2109 : }
2110 : }
2111 :
2112 : /************************************************************************/
2113 : /* GDALSerializeGeoLocTransformer() */
2114 : /************************************************************************/
2115 :
2116 1 : CPLXMLNode *GDALSerializeGeoLocTransformer(void *pTransformArg)
2117 :
2118 : {
2119 1 : VALIDATE_POINTER1(pTransformArg, "GDALSerializeGeoLocTransformer", nullptr);
2120 :
2121 1 : GDALGeoLocTransformInfo *psInfo =
2122 : static_cast<GDALGeoLocTransformInfo *>(pTransformArg);
2123 :
2124 : CPLXMLNode *psTree =
2125 1 : CPLCreateXMLNode(nullptr, CXT_Element, "GeoLocTransformer");
2126 :
2127 : /* -------------------------------------------------------------------- */
2128 : /* Serialize bReversed. */
2129 : /* -------------------------------------------------------------------- */
2130 1 : CPLCreateXMLElementAndValue(
2131 : psTree, "Reversed",
2132 2 : CPLString().Printf("%d", static_cast<int>(psInfo->bReversed)));
2133 :
2134 : /* -------------------------------------------------------------------- */
2135 : /* geoloc metadata. */
2136 : /* -------------------------------------------------------------------- */
2137 1 : char **papszMD = psInfo->papszGeolocationInfo;
2138 1 : CPLXMLNode *psMD = CPLCreateXMLNode(psTree, CXT_Element, "Metadata");
2139 :
2140 10 : for (int i = 0; papszMD != nullptr && papszMD[i] != nullptr; i++)
2141 : {
2142 9 : char *pszKey = nullptr;
2143 9 : const char *pszRawValue = CPLParseNameValue(papszMD[i], &pszKey);
2144 :
2145 9 : CPLXMLNode *psMDI = CPLCreateXMLNode(psMD, CXT_Element, "MDI");
2146 9 : CPLSetXMLValue(psMDI, "#key", pszKey);
2147 9 : CPLCreateXMLNode(psMDI, CXT_Text, pszRawValue);
2148 :
2149 9 : CPLFree(pszKey);
2150 : }
2151 :
2152 1 : return psTree;
2153 : }
2154 :
2155 : /************************************************************************/
2156 : /* GDALDeserializeGeoLocTransformer() */
2157 : /************************************************************************/
2158 :
2159 1 : void *GDALDeserializeGeoLocTransformer(CPLXMLNode *psTree)
2160 :
2161 : {
2162 : /* -------------------------------------------------------------------- */
2163 : /* Collect metadata. */
2164 : /* -------------------------------------------------------------------- */
2165 1 : CPLXMLNode *psMetadata = CPLGetXMLNode(psTree, "Metadata");
2166 :
2167 1 : if (psMetadata == nullptr || psMetadata->eType != CXT_Element ||
2168 1 : !EQUAL(psMetadata->pszValue, "Metadata"))
2169 0 : return nullptr;
2170 :
2171 1 : char **papszMD = nullptr;
2172 :
2173 12 : for (CPLXMLNode *psMDI = psMetadata->psChild; psMDI != nullptr;
2174 11 : psMDI = psMDI->psNext)
2175 : {
2176 11 : if (!EQUAL(psMDI->pszValue, "MDI") || psMDI->eType != CXT_Element ||
2177 11 : psMDI->psChild == nullptr || psMDI->psChild->psNext == nullptr ||
2178 11 : psMDI->psChild->eType != CXT_Attribute ||
2179 11 : psMDI->psChild->psChild == nullptr)
2180 0 : continue;
2181 :
2182 11 : papszMD = CSLSetNameValue(papszMD, psMDI->psChild->psChild->pszValue,
2183 11 : psMDI->psChild->psNext->pszValue);
2184 : }
2185 :
2186 : /* -------------------------------------------------------------------- */
2187 : /* Get other flags. */
2188 : /* -------------------------------------------------------------------- */
2189 1 : const int bReversed = atoi(CPLGetXMLValue(psTree, "Reversed", "0"));
2190 :
2191 : /* -------------------------------------------------------------------- */
2192 : /* Generate transformation. */
2193 : /* -------------------------------------------------------------------- */
2194 :
2195 : const char *pszSourceDataset =
2196 1 : CPLGetXMLValue(psTree, "SourceDataset", nullptr);
2197 :
2198 1 : void *pResult = GDALCreateGeoLocTransformerEx(nullptr, papszMD, bReversed,
2199 : pszSourceDataset, nullptr);
2200 :
2201 : /* -------------------------------------------------------------------- */
2202 : /* Cleanup GCP copy. */
2203 : /* -------------------------------------------------------------------- */
2204 1 : CSLDestroy(papszMD);
2205 :
2206 1 : return pResult;
2207 : }
|
