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