/src/gdal/alg/gdalgeoloc.cpp
Line | Count | Source (jump to first uncovered line) |
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 | | static double UnshiftGeoX(const GDALGeoLocTransformInfo *psTransform, |
73 | | double dfX) |
74 | 0 | { |
75 | 0 | if (!psTransform->bGeographicSRSWithMinus180Plus180LongRange) |
76 | 0 | return dfX; |
77 | 0 | if (dfX > 180) |
78 | 0 | return dfX - 360; |
79 | 0 | if (dfX < -180) |
80 | 0 | return dfX + 360; |
81 | 0 | return dfX; |
82 | 0 | } |
83 | | |
84 | | /************************************************************************/ |
85 | | /* UpdateMinMax() */ |
86 | | /************************************************************************/ |
87 | | |
88 | | inline void UpdateMinMax(GDALGeoLocTransformInfo *psTransform, double dfGeoLocX, |
89 | | double dfGeoLocY) |
90 | 0 | { |
91 | 0 | if (dfGeoLocX < psTransform->dfMinX) |
92 | 0 | { |
93 | 0 | psTransform->dfMinX = dfGeoLocX; |
94 | 0 | psTransform->dfYAtMinX = dfGeoLocY; |
95 | 0 | } |
96 | 0 | if (dfGeoLocX > psTransform->dfMaxX) |
97 | 0 | { |
98 | 0 | psTransform->dfMaxX = dfGeoLocX; |
99 | 0 | psTransform->dfYAtMaxX = dfGeoLocY; |
100 | 0 | } |
101 | 0 | if (dfGeoLocY < psTransform->dfMinY) |
102 | 0 | { |
103 | 0 | psTransform->dfMinY = dfGeoLocY; |
104 | 0 | psTransform->dfXAtMinY = dfGeoLocX; |
105 | 0 | } |
106 | 0 | if (dfGeoLocY > psTransform->dfMaxY) |
107 | 0 | { |
108 | 0 | psTransform->dfMaxY = dfGeoLocY; |
109 | 0 | psTransform->dfXAtMaxY = dfGeoLocX; |
110 | 0 | } |
111 | 0 | } |
112 | | |
113 | | /************************************************************************/ |
114 | | /* Clamp() */ |
115 | | /************************************************************************/ |
116 | | |
117 | | inline double Clamp(double v, double minV, double maxV) |
118 | 0 | { |
119 | 0 | return std::min(std::max(v, minV), maxV); |
120 | 0 | } |
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 | 0 | { \ |
130 | 0 | const int _nYBlocks = DIV_ROUND_UP(_rasterYSize, _tileYSize); \ |
131 | 0 | const int _nXBlocks = DIV_ROUND_UP(_rasterXSize, _tileXSize); \ |
132 | 0 | for (int _iYBlock = 0; _iYBlock < _nYBlocks; ++_iYBlock) \ |
133 | 0 | { \ |
134 | 0 | const int _iYStart = _iYBlock * _tileYSize; \ |
135 | 0 | const int _iYEnd = _iYBlock == _nYBlocks - 1 \ |
136 | 0 | ? _rasterYSize \ |
137 | 0 | : _iYStart + _tileYSize; \ |
138 | 0 | INIT_YBLOCK; \ |
139 | 0 | for (int _iXBlock = 0; _iXBlock < _nXBlocks; ++_iXBlock) \ |
140 | 0 | { \ |
141 | 0 | const int _iXStart = _iXBlock * _tileXSize; \ |
142 | 0 | const int _iXEnd = _iXBlock == _nXBlocks - 1 \ |
143 | 0 | ? _rasterXSize \ |
144 | 0 | : _iXStart + _tileXSize; |
145 | | |
146 | | #define END_ITER_PER_BLOCK \ |
147 | 0 | } \ |
148 | 0 | } \ |
149 | 0 | } |
150 | | |
151 | | /************************************************************************/ |
152 | | /* GDALGeoLoc::LoadGeolocFinish() */ |
153 | | /************************************************************************/ |
154 | | |
155 | | /*! @cond Doxygen_Suppress */ |
156 | | |
157 | | template <class Accessors> |
158 | | void GDALGeoLoc<Accessors>::LoadGeolocFinish( |
159 | | GDALGeoLocTransformInfo *psTransform) |
160 | 0 | { |
161 | 0 | auto pAccessors = static_cast<Accessors *>(psTransform->pAccessors); |
162 | 0 | CSLConstList papszGeolocationInfo = psTransform->papszGeolocationInfo; |
163 | | |
164 | | /* -------------------------------------------------------------------- */ |
165 | | /* Scan forward map for lat/long extents. */ |
166 | | /* -------------------------------------------------------------------- */ |
167 | 0 | psTransform->dfMinX = std::numeric_limits<double>::max(); |
168 | 0 | psTransform->dfMaxX = -std::numeric_limits<double>::max(); |
169 | 0 | psTransform->dfMinY = std::numeric_limits<double>::max(); |
170 | 0 | psTransform->dfMaxY = -std::numeric_limits<double>::max(); |
171 | |
|
172 | 0 | constexpr int TILE_SIZE = GDALGeoLocDatasetAccessors::TILE_SIZE; |
173 | 0 | START_ITER_PER_BLOCK(psTransform->nGeoLocXSize, TILE_SIZE, |
174 | 0 | psTransform->nGeoLocYSize, TILE_SIZE, (void)0, iXStart, |
175 | 0 | iXEnd, iYStart, iYEnd) |
176 | 0 | { |
177 | 0 | for (int iY = iYStart; iY < iYEnd; ++iY) |
178 | 0 | { |
179 | 0 | for (int iX = iXStart; iX < iXEnd; ++iX) |
180 | 0 | { |
181 | 0 | const auto dfX = pAccessors->geolocXAccessor.Get(iX, iY); |
182 | 0 | if (!psTransform->bHasNoData || dfX != psTransform->dfNoDataX) |
183 | 0 | { |
184 | 0 | UpdateMinMax(psTransform, dfX, |
185 | 0 | pAccessors->geolocYAccessor.Get(iX, iY)); |
186 | 0 | } |
187 | 0 | } |
188 | 0 | } |
189 | 0 | } |
190 | 0 | END_ITER_PER_BLOCK |
191 | | |
192 | | // Check if the SRS is geographic and the geoloc longitudes are in |
193 | | // [-180,180] |
194 | 0 | psTransform->bGeographicSRSWithMinus180Plus180LongRange = false; |
195 | 0 | const char *pszSRS = CSLFetchNameValue(papszGeolocationInfo, "SRS"); |
196 | 0 | if (pszSRS && psTransform->dfMinX >= -180.0 && |
197 | 0 | psTransform->dfMaxX <= 180.0 && !psTransform->bSwapXY) |
198 | 0 | { |
199 | 0 | OGRSpatialReference oSRS; |
200 | 0 | psTransform->bGeographicSRSWithMinus180Plus180LongRange = |
201 | 0 | oSRS.importFromWkt(pszSRS) == OGRERR_NONE && |
202 | 0 | CPL_TO_BOOL(oSRS.IsGeographic()); |
203 | 0 | } |
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 | 0 | if (psTransform->bOriginIsTopLeftCorner) |
279 | 0 | { |
280 | | // Add "virtual" edge at Y=nGeoLocYSize |
281 | 0 | for (int iX = 0; iX <= psTransform->nGeoLocXSize; iX++) |
282 | 0 | { |
283 | 0 | double dfGeoLocX; |
284 | 0 | double dfGeoLocY; |
285 | 0 | if (!PixelLineToXY(psTransform, static_cast<double>(iX), |
286 | 0 | static_cast<double>(psTransform->nGeoLocYSize), |
287 | 0 | dfGeoLocX, dfGeoLocY)) |
288 | 0 | continue; |
289 | 0 | if (psTransform->bGeographicSRSWithMinus180Plus180LongRange) |
290 | 0 | dfGeoLocX = Clamp(dfGeoLocX, -180.0, 180.0); |
291 | 0 | UpdateMinMax(psTransform, dfGeoLocX, dfGeoLocY); |
292 | 0 | } |
293 | | |
294 | | // Add "virtual" edge at X=nGeoLocXSize |
295 | 0 | for (int iY = 0; iY <= psTransform->nGeoLocYSize; iY++) |
296 | 0 | { |
297 | 0 | double dfGeoLocX; |
298 | 0 | double dfGeoLocY; |
299 | 0 | if (!PixelLineToXY(psTransform, |
300 | 0 | static_cast<double>(psTransform->nGeoLocXSize), |
301 | 0 | static_cast<double>(iY), dfGeoLocX, dfGeoLocY)) |
302 | 0 | continue; |
303 | 0 | if (psTransform->bGeographicSRSWithMinus180Plus180LongRange) |
304 | 0 | dfGeoLocX = Clamp(dfGeoLocX, -180.0, 180.0); |
305 | 0 | UpdateMinMax(psTransform, dfGeoLocX, dfGeoLocY); |
306 | 0 | } |
307 | 0 | } |
308 | 0 | else |
309 | 0 | { |
310 | | // Extend by half-pixel on 4 edges for pixel-center convention |
311 | |
|
312 | 0 | for (int iX = 0; iX <= psTransform->nGeoLocXSize; iX++) |
313 | 0 | { |
314 | 0 | double dfGeoLocX; |
315 | 0 | double dfGeoLocY; |
316 | 0 | if (!PixelLineToXY(psTransform, static_cast<double>(iX), -0.5, |
317 | 0 | dfGeoLocX, dfGeoLocY)) |
318 | 0 | continue; |
319 | 0 | if (psTransform->bGeographicSRSWithMinus180Plus180LongRange) |
320 | 0 | dfGeoLocX = Clamp(dfGeoLocX, -180.0, 180.0); |
321 | 0 | UpdateMinMax(psTransform, dfGeoLocX, dfGeoLocY); |
322 | 0 | } |
323 | |
|
324 | 0 | for (int iX = 0; iX <= psTransform->nGeoLocXSize; iX++) |
325 | 0 | { |
326 | 0 | double dfGeoLocX; |
327 | 0 | double dfGeoLocY; |
328 | 0 | if (!PixelLineToXY( |
329 | 0 | psTransform, static_cast<double>(iX), |
330 | 0 | static_cast<double>(psTransform->nGeoLocYSize - 1 + 0.5), |
331 | 0 | dfGeoLocX, dfGeoLocY)) |
332 | 0 | continue; |
333 | 0 | if (psTransform->bGeographicSRSWithMinus180Plus180LongRange) |
334 | 0 | dfGeoLocX = Clamp(dfGeoLocX, -180.0, 180.0); |
335 | 0 | UpdateMinMax(psTransform, dfGeoLocX, dfGeoLocY); |
336 | 0 | } |
337 | |
|
338 | 0 | for (int iY = 0; iY <= psTransform->nGeoLocYSize; iY++) |
339 | 0 | { |
340 | 0 | double dfGeoLocX; |
341 | 0 | double dfGeoLocY; |
342 | 0 | if (!PixelLineToXY(psTransform, -0.5, static_cast<double>(iY), |
343 | 0 | dfGeoLocX, dfGeoLocY)) |
344 | 0 | continue; |
345 | 0 | if (psTransform->bGeographicSRSWithMinus180Plus180LongRange) |
346 | 0 | dfGeoLocX = Clamp(dfGeoLocX, -180.0, 180.0); |
347 | 0 | UpdateMinMax(psTransform, dfGeoLocX, dfGeoLocY); |
348 | 0 | } |
349 | |
|
350 | 0 | for (int iY = 0; iY <= psTransform->nGeoLocYSize; iY++) |
351 | 0 | { |
352 | 0 | double dfGeoLocX; |
353 | 0 | double dfGeoLocY; |
354 | 0 | if (!PixelLineToXY(psTransform, psTransform->nGeoLocXSize - 1 + 0.5, |
355 | 0 | static_cast<double>(iY), dfGeoLocX, dfGeoLocY)) |
356 | 0 | continue; |
357 | 0 | if (psTransform->bGeographicSRSWithMinus180Plus180LongRange) |
358 | 0 | dfGeoLocX = Clamp(dfGeoLocX, -180.0, 180.0); |
359 | 0 | UpdateMinMax(psTransform, dfGeoLocX, dfGeoLocY); |
360 | 0 | } |
361 | 0 | } |
362 | 0 | } Unexecuted instantiation: GDALGeoLoc<GDALGeoLocCArrayAccessors>::LoadGeolocFinish(GDALGeoLocTransformInfo*) Unexecuted instantiation: GDALGeoLoc<GDALGeoLocDatasetAccessors>::LoadGeolocFinish(GDALGeoLocTransformInfo*) |
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 | | bool GDALGeoLoc<Accessors>::PixelLineToXY( |
387 | | const GDALGeoLocTransformInfo *psTransform, const double dfGeoLocPixel, |
388 | | const double dfGeoLocLine, double &dfX, double &dfY) |
389 | 0 | { |
390 | 0 | int iX = static_cast<int>( |
391 | 0 | std::min(std::max(0.0, dfGeoLocPixel), |
392 | 0 | static_cast<double>(psTransform->nGeoLocXSize - 1))); |
393 | 0 | int iY = static_cast<int>( |
394 | 0 | std::min(std::max(0.0, dfGeoLocLine), |
395 | 0 | static_cast<double>(psTransform->nGeoLocYSize - 1))); |
396 | |
|
397 | 0 | auto pAccessors = static_cast<Accessors *>(psTransform->pAccessors); |
398 | |
|
399 | 0 | for (int iAttempt = 0; iAttempt < 2; ++iAttempt) |
400 | 0 | { |
401 | 0 | const double dfGLX_0_0 = pAccessors->geolocXAccessor.Get(iX, iY); |
402 | 0 | const double dfGLY_0_0 = pAccessors->geolocYAccessor.Get(iX, iY); |
403 | 0 | if (psTransform->bHasNoData && dfGLX_0_0 == psTransform->dfNoDataX) |
404 | 0 | { |
405 | 0 | return false; |
406 | 0 | } |
407 | | |
408 | | // This assumes infinite extension beyond borders of available |
409 | | // data based on closest grid square. |
410 | 0 | if (iX + 1 < psTransform->nGeoLocXSize && |
411 | 0 | iY + 1 < psTransform->nGeoLocYSize) |
412 | 0 | { |
413 | 0 | const double dfGLX_1_0 = |
414 | 0 | pAccessors->geolocXAccessor.Get(iX + 1, iY); |
415 | 0 | const double dfGLY_1_0 = |
416 | 0 | pAccessors->geolocYAccessor.Get(iX + 1, iY); |
417 | 0 | const double dfGLX_0_1 = |
418 | 0 | pAccessors->geolocXAccessor.Get(iX, iY + 1); |
419 | 0 | const double dfGLY_0_1 = |
420 | 0 | pAccessors->geolocYAccessor.Get(iX, iY + 1); |
421 | 0 | const double dfGLX_1_1 = |
422 | 0 | pAccessors->geolocXAccessor.Get(iX + 1, iY + 1); |
423 | 0 | const double dfGLY_1_1 = |
424 | 0 | pAccessors->geolocYAccessor.Get(iX + 1, iY + 1); |
425 | 0 | if (!psTransform->bHasNoData || |
426 | 0 | (dfGLX_1_0 != psTransform->dfNoDataX && |
427 | 0 | dfGLX_0_1 != psTransform->dfNoDataX && |
428 | 0 | dfGLX_1_1 != psTransform->dfNoDataX)) |
429 | 0 | { |
430 | 0 | const double dfGLX_1_0_adjusted = |
431 | 0 | ShiftGeoX(psTransform, dfGLX_0_0, dfGLX_1_0); |
432 | 0 | const double dfGLX_0_1_adjusted = |
433 | 0 | ShiftGeoX(psTransform, dfGLX_0_0, dfGLX_0_1); |
434 | 0 | const double dfGLX_1_1_adjusted = |
435 | 0 | ShiftGeoX(psTransform, dfGLX_0_0, dfGLX_1_1); |
436 | 0 | dfX = (1 - (dfGeoLocLine - iY)) * |
437 | 0 | (dfGLX_0_0 + (dfGeoLocPixel - iX) * |
438 | 0 | (dfGLX_1_0_adjusted - dfGLX_0_0)) + |
439 | 0 | (dfGeoLocLine - iY) * |
440 | 0 | (dfGLX_0_1_adjusted + |
441 | 0 | (dfGeoLocPixel - iX) * |
442 | 0 | (dfGLX_1_1_adjusted - dfGLX_0_1_adjusted)); |
443 | 0 | dfX = UnshiftGeoX(psTransform, dfX); |
444 | |
|
445 | 0 | dfY = (1 - (dfGeoLocLine - iY)) * |
446 | 0 | (dfGLY_0_0 + |
447 | 0 | (dfGeoLocPixel - iX) * (dfGLY_1_0 - dfGLY_0_0)) + |
448 | 0 | (dfGeoLocLine - iY) * |
449 | 0 | (dfGLY_0_1 + |
450 | 0 | (dfGeoLocPixel - iX) * (dfGLY_1_1 - dfGLY_0_1)); |
451 | 0 | break; |
452 | 0 | } |
453 | 0 | } |
454 | | |
455 | 0 | if (iX == psTransform->nGeoLocXSize - 1 && iX >= 1 && |
456 | 0 | iY + 1 < psTransform->nGeoLocYSize) |
457 | 0 | { |
458 | | // If we are after the right edge, then go one pixel left |
459 | | // and retry |
460 | 0 | iX--; |
461 | 0 | continue; |
462 | 0 | } |
463 | 0 | else if (iY == psTransform->nGeoLocYSize - 1 && iY >= 1 && |
464 | 0 | iX + 1 < psTransform->nGeoLocXSize) |
465 | 0 | { |
466 | | // If we are after the bottom edge, then go one pixel up |
467 | | // and retry |
468 | 0 | iY--; |
469 | 0 | continue; |
470 | 0 | } |
471 | 0 | else if (iX == psTransform->nGeoLocXSize - 1 && |
472 | 0 | iY == psTransform->nGeoLocYSize - 1 && iX >= 1 && iY >= 1) |
473 | 0 | { |
474 | | // If we are after the right and bottom edge, then go one pixel left |
475 | | // and up and retry |
476 | 0 | iX--; |
477 | 0 | iY--; |
478 | 0 | continue; |
479 | 0 | } |
480 | 0 | else if (iX + 1 < psTransform->nGeoLocXSize && |
481 | 0 | (!psTransform->bHasNoData || |
482 | 0 | pAccessors->geolocXAccessor.Get(iX + 1, iY) != |
483 | 0 | psTransform->dfNoDataX)) |
484 | 0 | { |
485 | 0 | const double dfGLX_1_0 = |
486 | 0 | pAccessors->geolocXAccessor.Get(iX + 1, iY); |
487 | 0 | const double dfGLY_1_0 = |
488 | 0 | pAccessors->geolocYAccessor.Get(iX + 1, iY); |
489 | 0 | dfX = |
490 | 0 | dfGLX_0_0 + |
491 | 0 | (dfGeoLocPixel - iX) * |
492 | 0 | (ShiftGeoX(psTransform, dfGLX_0_0, dfGLX_1_0) - dfGLX_0_0); |
493 | 0 | dfX = UnshiftGeoX(psTransform, dfX); |
494 | 0 | dfY = dfGLY_0_0 + (dfGeoLocPixel - iX) * (dfGLY_1_0 - dfGLY_0_0); |
495 | 0 | } |
496 | 0 | else if (iY + 1 < psTransform->nGeoLocYSize && |
497 | 0 | (!psTransform->bHasNoData || |
498 | 0 | pAccessors->geolocXAccessor.Get(iX, iY + 1) != |
499 | 0 | psTransform->dfNoDataX)) |
500 | 0 | { |
501 | 0 | const double dfGLX_0_1 = |
502 | 0 | pAccessors->geolocXAccessor.Get(iX, iY + 1); |
503 | 0 | const double dfGLY_0_1 = |
504 | 0 | pAccessors->geolocYAccessor.Get(iX, iY + 1); |
505 | 0 | dfX = |
506 | 0 | dfGLX_0_0 + |
507 | 0 | (dfGeoLocLine - iY) * |
508 | 0 | (ShiftGeoX(psTransform, dfGLX_0_0, dfGLX_0_1) - dfGLX_0_0); |
509 | 0 | dfX = UnshiftGeoX(psTransform, dfX); |
510 | 0 | dfY = dfGLY_0_0 + (dfGeoLocLine - iY) * (dfGLY_0_1 - dfGLY_0_0); |
511 | 0 | } |
512 | 0 | else |
513 | 0 | { |
514 | 0 | dfX = dfGLX_0_0; |
515 | 0 | dfY = dfGLY_0_0; |
516 | 0 | } |
517 | 0 | break; |
518 | 0 | } |
519 | 0 | return true; |
520 | 0 | } Unexecuted instantiation: GDALGeoLoc<GDALGeoLocCArrayAccessors>::PixelLineToXY(GDALGeoLocTransformInfo const*, double, double, double&, double&) Unexecuted instantiation: GDALGeoLoc<GDALGeoLocDatasetAccessors>::PixelLineToXY(GDALGeoLocTransformInfo const*, double, double, double&, double&) |
521 | | |
522 | | template <class Accessors> |
523 | | bool GDALGeoLoc<Accessors>::PixelLineToXY( |
524 | | const GDALGeoLocTransformInfo *psTransform, const int nGeoLocPixel, |
525 | | const int nGeoLocLine, double &dfX, double &dfY) |
526 | 0 | { |
527 | 0 | if (nGeoLocPixel >= 0 && nGeoLocPixel < psTransform->nGeoLocXSize && |
528 | 0 | nGeoLocLine >= 0 && nGeoLocLine < psTransform->nGeoLocYSize) |
529 | 0 | { |
530 | 0 | auto pAccessors = static_cast<Accessors *>(psTransform->pAccessors); |
531 | 0 | const double dfGLX = |
532 | 0 | pAccessors->geolocXAccessor.Get(nGeoLocPixel, nGeoLocLine); |
533 | 0 | const double dfGLY = |
534 | 0 | pAccessors->geolocYAccessor.Get(nGeoLocPixel, nGeoLocLine); |
535 | 0 | if (psTransform->bHasNoData && dfGLX == psTransform->dfNoDataX) |
536 | 0 | { |
537 | 0 | return false; |
538 | 0 | } |
539 | 0 | dfX = dfGLX; |
540 | 0 | dfY = dfGLY; |
541 | 0 | return true; |
542 | 0 | } |
543 | 0 | return PixelLineToXY(psTransform, static_cast<double>(nGeoLocPixel), |
544 | 0 | static_cast<double>(nGeoLocLine), dfX, dfY); |
545 | 0 | } Unexecuted instantiation: GDALGeoLoc<GDALGeoLocCArrayAccessors>::PixelLineToXY(GDALGeoLocTransformInfo const*, int, int, double&, double&) Unexecuted instantiation: GDALGeoLoc<GDALGeoLocDatasetAccessors>::PixelLineToXY(GDALGeoLocTransformInfo const*, int, int, double&, double&) |
546 | | |
547 | | /************************************************************************/ |
548 | | /* GDALGeoLoc::ExtractSquare() */ |
549 | | /************************************************************************/ |
550 | | |
551 | | template <class Accessors> |
552 | | 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 | 0 | { |
557 | 0 | return PixelLineToXY(psTransform, nX, nY, dfX_0_0, dfY_0_0) && |
558 | 0 | PixelLineToXY(psTransform, nX + 1, nY, dfX_1_0, dfY_1_0) && |
559 | 0 | PixelLineToXY(psTransform, nX, nY + 1, dfX_0_1, dfY_0_1) && |
560 | 0 | PixelLineToXY(psTransform, nX + 1, nY + 1, dfX_1_1, dfY_1_1); |
561 | 0 | } Unexecuted instantiation: GDALGeoLoc<GDALGeoLocCArrayAccessors>::ExtractSquare(GDALGeoLocTransformInfo const*, int, int, double&, double&, double&, double&, double&, double&, double&, double&) Unexecuted instantiation: GDALGeoLoc<GDALGeoLocDatasetAccessors>::ExtractSquare(GDALGeoLocTransformInfo const*, int, int, double&, double&, double&, double&, double&, double&, double&, double&) |
562 | | |
563 | | 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 | 0 | { |
568 | 0 | if (psTransform->bUseArray) |
569 | 0 | { |
570 | 0 | return GDALGeoLoc<GDALGeoLocCArrayAccessors>::ExtractSquare( |
571 | 0 | psTransform, nX, nY, dfX_0_0, dfY_0_0, dfX_1_0, dfY_1_0, dfX_0_1, |
572 | 0 | dfY_0_1, dfX_1_1, dfY_1_1); |
573 | 0 | } |
574 | 0 | else |
575 | 0 | { |
576 | 0 | return GDALGeoLoc<GDALGeoLocDatasetAccessors>::ExtractSquare( |
577 | 0 | 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 | 0 | } |
580 | 0 | } |
581 | | |
582 | | /************************************************************************/ |
583 | | /* GDALGeoLocTransform() */ |
584 | | /************************************************************************/ |
585 | | |
586 | | template <class Accessors> |
587 | | int GDALGeoLoc<Accessors>::Transform(void *pTransformArg, int bDstToSrc, |
588 | | int nPointCount, double *padfX, |
589 | | double *padfY, double * /* padfZ */, |
590 | | int *panSuccess) |
591 | 0 | { |
592 | 0 | int bSuccess = TRUE; |
593 | 0 | GDALGeoLocTransformInfo *psTransform = |
594 | 0 | static_cast<GDALGeoLocTransformInfo *>(pTransformArg); |
595 | |
|
596 | 0 | if (psTransform->bReversed) |
597 | 0 | bDstToSrc = !bDstToSrc; |
598 | |
|
599 | 0 | const double dfGeorefConventionOffset = |
600 | 0 | psTransform->bOriginIsTopLeftCorner ? 0 : 0.5; |
601 | | |
602 | | /* -------------------------------------------------------------------- */ |
603 | | /* Do original pixel line to target geox/geoy. */ |
604 | | /* -------------------------------------------------------------------- */ |
605 | 0 | if (!bDstToSrc) |
606 | 0 | { |
607 | 0 | for (int i = 0; i < nPointCount; i++) |
608 | 0 | { |
609 | 0 | if (padfX[i] == HUGE_VAL || padfY[i] == HUGE_VAL) |
610 | 0 | { |
611 | 0 | bSuccess = FALSE; |
612 | 0 | panSuccess[i] = FALSE; |
613 | 0 | continue; |
614 | 0 | } |
615 | | |
616 | 0 | const double dfGeoLocPixel = |
617 | 0 | (padfX[i] - psTransform->dfPIXEL_OFFSET) / |
618 | 0 | psTransform->dfPIXEL_STEP - |
619 | 0 | dfGeorefConventionOffset; |
620 | 0 | const double dfGeoLocLine = |
621 | 0 | (padfY[i] - psTransform->dfLINE_OFFSET) / |
622 | 0 | psTransform->dfLINE_STEP - |
623 | 0 | dfGeorefConventionOffset; |
624 | |
|
625 | 0 | if (!PixelLineToXY(psTransform, dfGeoLocPixel, dfGeoLocLine, |
626 | 0 | padfX[i], padfY[i])) |
627 | 0 | { |
628 | 0 | bSuccess = FALSE; |
629 | 0 | panSuccess[i] = FALSE; |
630 | 0 | padfX[i] = HUGE_VAL; |
631 | 0 | padfY[i] = HUGE_VAL; |
632 | 0 | continue; |
633 | 0 | } |
634 | | |
635 | 0 | if (psTransform->bSwapXY) |
636 | 0 | { |
637 | 0 | std::swap(padfX[i], padfY[i]); |
638 | 0 | } |
639 | |
|
640 | 0 | panSuccess[i] = TRUE; |
641 | 0 | } |
642 | 0 | } |
643 | | |
644 | | /* -------------------------------------------------------------------- */ |
645 | | /* geox/geoy to pixel/line using backmap. */ |
646 | | /* -------------------------------------------------------------------- */ |
647 | 0 | else |
648 | 0 | { |
649 | 0 | if (psTransform->hQuadTree) |
650 | 0 | { |
651 | 0 | GDALGeoLocInverseTransformQuadtree(psTransform, nPointCount, padfX, |
652 | 0 | padfY, panSuccess); |
653 | 0 | return TRUE; |
654 | 0 | } |
655 | | |
656 | 0 | const bool bGeolocMaxAccuracy = CPLTestBool( |
657 | 0 | 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 | 0 | OGRPoint oPoint; |
662 | 0 | OGRLinearRing oRing; |
663 | 0 | oRing.setNumPoints(5); |
664 | |
|
665 | 0 | auto pAccessors = static_cast<Accessors *>(psTransform->pAccessors); |
666 | |
|
667 | 0 | for (int i = 0; i < nPointCount; i++) |
668 | 0 | { |
669 | 0 | if (padfX[i] == HUGE_VAL || padfY[i] == HUGE_VAL) |
670 | 0 | { |
671 | 0 | bSuccess = FALSE; |
672 | 0 | panSuccess[i] = FALSE; |
673 | 0 | continue; |
674 | 0 | } |
675 | | |
676 | 0 | if (psTransform->bSwapXY) |
677 | 0 | { |
678 | 0 | std::swap(padfX[i], padfY[i]); |
679 | 0 | } |
680 | |
|
681 | 0 | const double dfGeoX = padfX[i]; |
682 | 0 | const double dfGeoY = padfY[i]; |
683 | |
|
684 | 0 | const double dfBMX = |
685 | 0 | ((padfX[i] - psTransform->adfBackMapGeoTransform[0]) / |
686 | 0 | psTransform->adfBackMapGeoTransform[1]); |
687 | 0 | const double dfBMY = |
688 | 0 | ((padfY[i] - psTransform->adfBackMapGeoTransform[3]) / |
689 | 0 | psTransform->adfBackMapGeoTransform[5]); |
690 | |
|
691 | 0 | if (!(dfBMX >= 0 && dfBMY >= 0 && |
692 | 0 | dfBMX + 1 < psTransform->nBackMapWidth && |
693 | 0 | dfBMY + 1 < psTransform->nBackMapHeight)) |
694 | 0 | { |
695 | 0 | bSuccess = FALSE; |
696 | 0 | panSuccess[i] = FALSE; |
697 | 0 | padfX[i] = HUGE_VAL; |
698 | 0 | padfY[i] = HUGE_VAL; |
699 | 0 | continue; |
700 | 0 | } |
701 | | |
702 | 0 | const int iBMX = static_cast<int>(dfBMX); |
703 | 0 | const int iBMY = static_cast<int>(dfBMY); |
704 | |
|
705 | 0 | const auto fBMX_0_0 = pAccessors->backMapXAccessor.Get(iBMX, iBMY); |
706 | 0 | const auto fBMY_0_0 = pAccessors->backMapYAccessor.Get(iBMX, iBMY); |
707 | 0 | if (fBMX_0_0 == INVALID_BMXY) |
708 | 0 | { |
709 | 0 | bSuccess = FALSE; |
710 | 0 | panSuccess[i] = FALSE; |
711 | 0 | padfX[i] = HUGE_VAL; |
712 | 0 | padfY[i] = HUGE_VAL; |
713 | 0 | continue; |
714 | 0 | } |
715 | | |
716 | 0 | const auto fBMX_1_0 = |
717 | 0 | pAccessors->backMapXAccessor.Get(iBMX + 1, iBMY); |
718 | 0 | const auto fBMY_1_0 = |
719 | 0 | pAccessors->backMapYAccessor.Get(iBMX + 1, iBMY); |
720 | 0 | const auto fBMX_0_1 = |
721 | 0 | pAccessors->backMapXAccessor.Get(iBMX, iBMY + 1); |
722 | 0 | const auto fBMY_0_1 = |
723 | 0 | pAccessors->backMapYAccessor.Get(iBMX, iBMY + 1); |
724 | 0 | const auto fBMX_1_1 = |
725 | 0 | pAccessors->backMapXAccessor.Get(iBMX + 1, iBMY + 1); |
726 | 0 | const auto fBMY_1_1 = |
727 | 0 | pAccessors->backMapYAccessor.Get(iBMX + 1, iBMY + 1); |
728 | 0 | if (fBMX_1_0 != INVALID_BMXY && fBMX_0_1 != INVALID_BMXY && |
729 | 0 | fBMX_1_1 != INVALID_BMXY) |
730 | 0 | { |
731 | 0 | padfX[i] = |
732 | 0 | (1 - (dfBMY - iBMY)) * |
733 | 0 | (fBMX_0_0 + (dfBMX - iBMX) * (fBMX_1_0 - fBMX_0_0)) + |
734 | 0 | (dfBMY - iBMY) * |
735 | 0 | (fBMX_0_1 + (dfBMX - iBMX) * (fBMX_1_1 - fBMX_0_1)); |
736 | 0 | padfY[i] = |
737 | 0 | (1 - (dfBMY - iBMY)) * |
738 | 0 | (fBMY_0_0 + (dfBMX - iBMX) * (fBMY_1_0 - fBMY_0_0)) + |
739 | 0 | (dfBMY - iBMY) * |
740 | 0 | (fBMY_0_1 + (dfBMX - iBMX) * (fBMY_1_1 - fBMY_0_1)); |
741 | 0 | } |
742 | 0 | else if (fBMX_1_0 != INVALID_BMXY) |
743 | 0 | { |
744 | 0 | padfX[i] = fBMX_0_0 + (dfBMX - iBMX) * (fBMX_1_0 - fBMX_0_0); |
745 | 0 | padfY[i] = fBMY_0_0 + (dfBMX - iBMX) * (fBMY_1_0 - fBMY_0_0); |
746 | 0 | } |
747 | 0 | else if (fBMX_0_1 != INVALID_BMXY) |
748 | 0 | { |
749 | 0 | padfX[i] = fBMX_0_0 + (dfBMY - iBMY) * (fBMX_0_1 - fBMX_0_0); |
750 | 0 | padfY[i] = fBMY_0_0 + (dfBMY - iBMY) * (fBMY_0_1 - fBMY_0_0); |
751 | 0 | } |
752 | 0 | else |
753 | 0 | { |
754 | 0 | padfX[i] = fBMX_0_0; |
755 | 0 | padfY[i] = fBMY_0_0; |
756 | 0 | } |
757 | |
|
758 | 0 | const double dfGeoLocPixel = |
759 | 0 | (padfX[i] - psTransform->dfPIXEL_OFFSET) / |
760 | 0 | psTransform->dfPIXEL_STEP - |
761 | 0 | dfGeorefConventionOffset; |
762 | 0 | const double dfGeoLocLine = |
763 | 0 | (padfY[i] - psTransform->dfLINE_OFFSET) / |
764 | 0 | psTransform->dfLINE_STEP - |
765 | 0 | 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 | 0 | if (!bGeolocMaxAccuracy) |
781 | 0 | { |
782 | 0 | panSuccess[i] = TRUE; |
783 | 0 | continue; |
784 | 0 | } |
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 | 0 | oPoint.setX(dfGeoX); |
795 | 0 | 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 | 0 | const int nSearchRadius = |
803 | 0 | psTransform->bGeographicSRSWithMinus180Plus180LongRange && |
804 | 0 | fabs(dfGeoY) >= 85 |
805 | 0 | ? 5 |
806 | 0 | : 3; |
807 | 0 | const int nGeoLocPixel = |
808 | 0 | static_cast<int>(std::floor(dfGeoLocPixel)); |
809 | 0 | const int nGeoLocLine = static_cast<int>(std::floor(dfGeoLocLine)); |
810 | |
|
811 | 0 | 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 | 0 | for (int r = 0; !bDone && r <= nSearchRadius; r++) |
816 | 0 | { |
817 | 0 | for (int iter = 0; !bDone && iter < (r == 0 ? 1 : 8 * r); |
818 | 0 | ++iter) |
819 | 0 | { |
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 | 0 | const int sx = (r == 0) ? 0 |
824 | 0 | : (iter < 2 * r) ? -r + iter |
825 | 0 | : (iter < 4 * r) ? r |
826 | 0 | : (iter < 6 * r) ? r - (iter - 4 * r) |
827 | 0 | : -r; |
828 | 0 | const int sy = (r == 0) ? 0 |
829 | 0 | : (iter < 2 * r) ? r |
830 | 0 | : (iter < 4 * r) ? r - (iter - 2 * r) |
831 | 0 | : (iter < 6 * r) ? -r |
832 | 0 | : -r + (iter - 6 * r); |
833 | 0 | if (nGeoLocPixel >= |
834 | 0 | static_cast<int>(psTransform->nGeoLocXSize) - sx || |
835 | 0 | nGeoLocLine >= |
836 | 0 | static_cast<int>(psTransform->nGeoLocYSize) - sy) |
837 | 0 | { |
838 | 0 | continue; |
839 | 0 | } |
840 | 0 | const int iX = nGeoLocPixel + sx; |
841 | 0 | const int iY = nGeoLocLine + sy; |
842 | 0 | if (iX >= -1 || iY >= -1) |
843 | 0 | { |
844 | 0 | double x0, y0, x1, y1, x2, y2, x3, y3; |
845 | |
|
846 | 0 | if (!PixelLineToXY(psTransform, iX, iY, x0, y0) || |
847 | 0 | !PixelLineToXY(psTransform, iX + 1, iY, x2, y2) || |
848 | 0 | !PixelLineToXY(psTransform, iX, iY + 1, x1, y1) || |
849 | 0 | !PixelLineToXY(psTransform, iX + 1, iY + 1, x3, y3)) |
850 | 0 | { |
851 | 0 | continue; |
852 | 0 | } |
853 | | |
854 | 0 | int nIters = 1; |
855 | | // For a bounding box crossing the anti-meridian, check |
856 | | // both around -180 and +180 deg. |
857 | 0 | if (psTransform |
858 | 0 | ->bGeographicSRSWithMinus180Plus180LongRange && |
859 | 0 | std::fabs(x0) > 170 && std::fabs(x1) > 170 && |
860 | 0 | std::fabs(x2) > 170 && std::fabs(x3) > 170 && |
861 | 0 | (std::fabs(x1 - x0) > 180 || |
862 | 0 | std::fabs(x2 - x0) > 180 || |
863 | 0 | std::fabs(x3 - x0) > 180)) |
864 | 0 | { |
865 | 0 | nIters = 2; |
866 | 0 | if (x0 > 0) |
867 | 0 | x0 -= 360; |
868 | 0 | if (x1 > 0) |
869 | 0 | x1 -= 360; |
870 | 0 | if (x2 > 0) |
871 | 0 | x2 -= 360; |
872 | 0 | if (x3 > 0) |
873 | 0 | x3 -= 360; |
874 | 0 | } |
875 | 0 | for (int iIter = 0; !bDone && iIter < nIters; ++iIter) |
876 | 0 | { |
877 | 0 | if (iIter == 1) |
878 | 0 | { |
879 | 0 | x0 += 360; |
880 | 0 | x1 += 360; |
881 | 0 | x2 += 360; |
882 | 0 | x3 += 360; |
883 | 0 | } |
884 | 0 | oRing.setPoint(0, x0, y0); |
885 | 0 | oRing.setPoint(1, x2, y2); |
886 | 0 | oRing.setPoint(2, x3, y3); |
887 | 0 | oRing.setPoint(3, x1, y1); |
888 | 0 | oRing.setPoint(4, x0, y0); |
889 | 0 | if (oRing.isPointInRing(&oPoint) || |
890 | 0 | oRing.isPointOnRingBoundary(&oPoint)) |
891 | 0 | { |
892 | 0 | double dfX = static_cast<double>(iX); |
893 | 0 | double dfY = static_cast<double>(iY); |
894 | 0 | GDALInverseBilinearInterpolation( |
895 | 0 | dfGeoX, dfGeoY, x0, y0, x1, y1, x2, y2, x3, |
896 | 0 | y3, dfX, dfY); |
897 | |
|
898 | 0 | dfX = (dfX + dfGeorefConventionOffset) * |
899 | 0 | psTransform->dfPIXEL_STEP + |
900 | 0 | psTransform->dfPIXEL_OFFSET; |
901 | 0 | dfY = (dfY + dfGeorefConventionOffset) * |
902 | 0 | psTransform->dfLINE_STEP + |
903 | 0 | 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 | 0 | padfX[i] = dfX; |
913 | 0 | padfY[i] = dfY; |
914 | |
|
915 | 0 | bDone = true; |
916 | 0 | } |
917 | 0 | } |
918 | 0 | } |
919 | 0 | } |
920 | 0 | } |
921 | 0 | if (!bDone) |
922 | 0 | { |
923 | 0 | bSuccess = FALSE; |
924 | 0 | panSuccess[i] = FALSE; |
925 | 0 | padfX[i] = HUGE_VAL; |
926 | 0 | padfY[i] = HUGE_VAL; |
927 | 0 | continue; |
928 | 0 | } |
929 | | |
930 | 0 | panSuccess[i] = TRUE; |
931 | 0 | } |
932 | 0 | } |
933 | | |
934 | 0 | return bSuccess; |
935 | 0 | } Unexecuted instantiation: GDALGeoLoc<GDALGeoLocCArrayAccessors>::Transform(void*, int, int, double*, double*, double*, int*) Unexecuted instantiation: GDALGeoLoc<GDALGeoLocDatasetAccessors>::Transform(void*, int, int, double*, double*, double*, int*) |
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 | | 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 | 0 | { |
955 | | // Exact inverse bilinear interpolation method. |
956 | | // Maths from https://stackoverflow.com/a/812077 |
957 | |
|
958 | 0 | const double A = (x0 - x) * (y0 - y2) - (y0 - y) * (x0 - x2); |
959 | 0 | const double B = (((x0 - x) * (y1 - y3) - (y0 - y) * (x1 - x3)) + |
960 | 0 | ((x1 - x) * (y0 - y2) - (y1 - y) * (x0 - x2))) / |
961 | 0 | 2; |
962 | 0 | const double C = (x1 - x) * (y1 - y3) - (y1 - y) * (x1 - x3); |
963 | 0 | const double denom = A - 2 * B + C; |
964 | 0 | double s; |
965 | 0 | const double magnitudeOfValues = fabs(A) + fabs(B) + fabs(C); |
966 | 0 | if (fabs(denom) <= 1e-12 * magnitudeOfValues) |
967 | 0 | { |
968 | | // Happens typically when the x_i,y_i points form a rectangle |
969 | | // Can also happen when they form a triangle. |
970 | 0 | s = A / (A - C); |
971 | 0 | } |
972 | 0 | else |
973 | 0 | { |
974 | 0 | const double sqrtTerm = sqrt(B * B - A * C); |
975 | 0 | const double s1 = ((A - B) + sqrtTerm) / denom; |
976 | 0 | const double s2 = ((A - B) - sqrtTerm) / denom; |
977 | 0 | if (s1 < 0 || s1 > 1) |
978 | 0 | s = s2; |
979 | 0 | else |
980 | 0 | s = s1; |
981 | 0 | } |
982 | |
|
983 | 0 | const double t_denom_x = (1 - s) * (x0 - x2) + s * (x1 - x3); |
984 | 0 | if (fabs(t_denom_x) > 1e-12 * magnitudeOfValues) |
985 | 0 | { |
986 | 0 | i += ((1 - s) * (x0 - x) + s * (x1 - x)) / t_denom_x; |
987 | 0 | } |
988 | 0 | else |
989 | 0 | { |
990 | 0 | const double t_denom_y = (1 - s) * (y0 - y2) + s * (y1 - y3); |
991 | 0 | if (fabs(t_denom_y) > 1e-12 * magnitudeOfValues) |
992 | 0 | { |
993 | 0 | i += ((1 - s) * (y0 - y) + s * (y1 - y)) / t_denom_y; |
994 | 0 | } |
995 | 0 | } |
996 | |
|
997 | 0 | j += s; |
998 | 0 | } |
999 | | |
1000 | | /************************************************************************/ |
1001 | | /* GeoLocGenerateBackMap() */ |
1002 | | /************************************************************************/ |
1003 | | |
1004 | | /*! @cond Doxygen_Suppress */ |
1005 | | |
1006 | | template <class Accessors> |
1007 | | bool GDALGeoLoc<Accessors>::GenerateBackMap( |
1008 | | GDALGeoLocTransformInfo *psTransform) |
1009 | | |
1010 | 0 | { |
1011 | 0 | CPLDebug("GEOLOC", "Starting backmap generation"); |
1012 | 0 | const int nXSize = psTransform->nGeoLocXSize; |
1013 | 0 | 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 | 0 | const double dfTargetPixels = |
1022 | 0 | static_cast<double>(nXSize) * nYSize * psTransform->dfOversampleFactor; |
1023 | 0 | const double dfPixelSizeSquare = |
1024 | 0 | sqrt((psTransform->dfMaxX - psTransform->dfMinX) * |
1025 | 0 | (psTransform->dfMaxY - psTransform->dfMinY) / dfTargetPixels); |
1026 | 0 | if (dfPixelSizeSquare == 0.0) |
1027 | 0 | { |
1028 | 0 | CPLError(CE_Failure, CPLE_AppDefined, "Invalid pixel size for backmap"); |
1029 | 0 | return false; |
1030 | 0 | } |
1031 | | |
1032 | 0 | const double dfMinX = psTransform->dfMinX - dfPixelSizeSquare / 2.0; |
1033 | 0 | const double dfMaxX = psTransform->dfMaxX + dfPixelSizeSquare / 2.0; |
1034 | 0 | const double dfMaxY = psTransform->dfMaxY + dfPixelSizeSquare / 2.0; |
1035 | 0 | const double dfMinY = psTransform->dfMinY - dfPixelSizeSquare / 2.0; |
1036 | 0 | const double dfBMXSize = std::ceil((dfMaxX - dfMinX) / dfPixelSizeSquare); |
1037 | 0 | 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 | 0 | if (!(dfBMXSize > 0 && dfBMXSize + 2 < INT_MAX) || |
1042 | 0 | !(dfBMYSize > 0 && dfBMYSize + 2 < INT_MAX)) |
1043 | 0 | { |
1044 | 0 | CPLError(CE_Failure, CPLE_AppDefined, "Int overflow : %f x %f", |
1045 | 0 | dfBMXSize, dfBMYSize); |
1046 | 0 | return false; |
1047 | 0 | } |
1048 | | |
1049 | 0 | int nBMXSize = static_cast<int>(dfBMXSize); |
1050 | 0 | int nBMYSize = static_cast<int>(dfBMYSize); |
1051 | |
|
1052 | 0 | if (static_cast<size_t>(1 + nBMYSize) > |
1053 | 0 | std::numeric_limits<size_t>::max() / static_cast<size_t>(1 + nBMXSize)) |
1054 | 0 | { |
1055 | 0 | CPLError(CE_Failure, CPLE_AppDefined, "Int overflow : %f x %f", |
1056 | 0 | dfBMXSize, dfBMYSize); |
1057 | 0 | return false; |
1058 | 0 | } |
1059 | | |
1060 | 0 | const double dfPixelXSize = (dfMaxX - dfMinX) / nBMXSize; |
1061 | 0 | const double dfPixelYSize = (dfMaxY - dfMinY) / nBMYSize; |
1062 | | |
1063 | | // Extra pixel for right-edge and bottom-edge extensions in TOP_LEFT_CORNER |
1064 | | // convention. |
1065 | 0 | nBMXSize++; |
1066 | 0 | nBMYSize++; |
1067 | 0 | psTransform->nBackMapWidth = nBMXSize; |
1068 | 0 | psTransform->nBackMapHeight = nBMYSize; |
1069 | |
|
1070 | 0 | psTransform->adfBackMapGeoTransform[0] = dfMinX; |
1071 | 0 | psTransform->adfBackMapGeoTransform[1] = dfPixelXSize; |
1072 | 0 | psTransform->adfBackMapGeoTransform[2] = 0.0; |
1073 | 0 | psTransform->adfBackMapGeoTransform[3] = dfMaxY; |
1074 | 0 | psTransform->adfBackMapGeoTransform[4] = 0.0; |
1075 | 0 | psTransform->adfBackMapGeoTransform[5] = -dfPixelYSize; |
1076 | | |
1077 | | /* -------------------------------------------------------------------- */ |
1078 | | /* Allocate backmap. */ |
1079 | | /* -------------------------------------------------------------------- */ |
1080 | 0 | auto pAccessors = static_cast<Accessors *>(psTransform->pAccessors); |
1081 | 0 | if (!pAccessors->AllocateBackMap()) |
1082 | 0 | return false; |
1083 | | |
1084 | 0 | const double dfGeorefConventionOffset = |
1085 | 0 | psTransform->bOriginIsTopLeftCorner ? 0 : 0.5; |
1086 | |
|
1087 | 0 | const auto UpdateBackmap = |
1088 | 0 | [&](int iBMX, int iBMY, double dfX, double dfY, double tempwt) |
1089 | 0 | { |
1090 | 0 | const auto fBMX = pAccessors->backMapXAccessor.Get(iBMX, iBMY); |
1091 | 0 | const auto fBMY = pAccessors->backMapYAccessor.Get(iBMX, iBMY); |
1092 | 0 | const float fUpdatedBMX = |
1093 | 0 | fBMX + |
1094 | 0 | static_cast<float>(tempwt * ((dfX + dfGeorefConventionOffset) * |
1095 | 0 | psTransform->dfPIXEL_STEP + |
1096 | 0 | psTransform->dfPIXEL_OFFSET)); |
1097 | 0 | const float fUpdatedBMY = |
1098 | 0 | fBMY + |
1099 | 0 | static_cast<float>(tempwt * ((dfY + dfGeorefConventionOffset) * |
1100 | 0 | psTransform->dfLINE_STEP + |
1101 | 0 | psTransform->dfLINE_OFFSET)); |
1102 | 0 | const float fUpdatedWeight = |
1103 | 0 | pAccessors->backMapWeightAccessor.Get(iBMX, iBMY) + |
1104 | 0 | 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 | 0 | if (fUpdatedWeight > 0) |
1111 | 0 | { |
1112 | 0 | const float fX = fUpdatedBMX / fUpdatedWeight; |
1113 | 0 | const float fY = fUpdatedBMY / fUpdatedWeight; |
1114 | 0 | const double dfGeoLocPixel = |
1115 | 0 | (fX - psTransform->dfPIXEL_OFFSET) / psTransform->dfPIXEL_STEP - |
1116 | 0 | dfGeorefConventionOffset; |
1117 | 0 | const double dfGeoLocLine = |
1118 | 0 | (fY - psTransform->dfLINE_OFFSET) / psTransform->dfLINE_STEP - |
1119 | 0 | dfGeorefConventionOffset; |
1120 | 0 | int iXAvg = static_cast<int>(std::max(0.0, dfGeoLocPixel)); |
1121 | 0 | iXAvg = std::min(iXAvg, psTransform->nGeoLocXSize - 1); |
1122 | 0 | int iYAvg = static_cast<int>(std::max(0.0, dfGeoLocLine)); |
1123 | 0 | iYAvg = std::min(iYAvg, psTransform->nGeoLocYSize - 1); |
1124 | 0 | const double dfGLX = pAccessors->geolocXAccessor.Get(iXAvg, iYAvg); |
1125 | 0 | const double dfGLY = pAccessors->geolocYAccessor.Get(iXAvg, iYAvg); |
1126 | |
|
1127 | 0 | const unsigned iX = static_cast<unsigned>(dfX); |
1128 | 0 | const unsigned iY = static_cast<unsigned>(dfY); |
1129 | 0 | if (!(psTransform->bHasNoData && dfGLX == psTransform->dfNoDataX) && |
1130 | 0 | ((iX >= static_cast<unsigned>(nXSize - 1) || |
1131 | 0 | iY >= static_cast<unsigned>(nYSize - 1)) || |
1132 | 0 | (fabs(dfGLX - pAccessors->geolocXAccessor.Get(iX, iY)) <= |
1133 | 0 | 2 * dfPixelXSize && |
1134 | 0 | fabs(dfGLY - pAccessors->geolocYAccessor.Get(iX, iY)) <= |
1135 | 0 | 2 * dfPixelYSize))) |
1136 | 0 | { |
1137 | 0 | pAccessors->backMapXAccessor.Set(iBMX, iBMY, fUpdatedBMX); |
1138 | 0 | pAccessors->backMapYAccessor.Set(iBMX, iBMY, fUpdatedBMY); |
1139 | 0 | pAccessors->backMapWeightAccessor.Set(iBMX, iBMY, |
1140 | 0 | fUpdatedWeight); |
1141 | 0 | } |
1142 | 0 | } |
1143 | 0 | }; Unexecuted instantiation: GDALGeoLoc<GDALGeoLocCArrayAccessors>::GenerateBackMap(GDALGeoLocTransformInfo*)::{lambda(int, int, double, double, double)#1}::operator()(int, int, double, double, double) constUnexecuted instantiation: GDALGeoLoc<GDALGeoLocDatasetAccessors>::GenerateBackMap(GDALGeoLocTransformInfo*)::{lambda(int, int, double, double, double)#1}::operator()(int, int, double, double, double) const |
1144 | | |
1145 | | // Keep those objects in this outer scope, so they are re-used, to |
1146 | | // save memory allocations. |
1147 | 0 | OGRPoint oPoint; |
1148 | 0 | OGRLinearRing oRing; |
1149 | 0 | 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 | 0 | const double dfStep = 1. / psTransform->dfOversampleFactor; |
1160 | |
|
1161 | 0 | constexpr int TILE_SIZE = GDALGeoLocDatasetAccessors::TILE_SIZE; |
1162 | 0 | const int nYBlocks = DIV_ROUND_UP(nYSize, TILE_SIZE); |
1163 | 0 | 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 | 0 | std::vector<std::pair<double, double>> yStartEnd(nYBlocks + 1); |
1168 | 0 | std::vector<std::pair<double, double>> xStartEnd(nXBlocks + 1); |
1169 | |
|
1170 | 0 | { |
1171 | 0 | int iYBlockLast = -1; |
1172 | 0 | double dfY = -dfStep; |
1173 | 0 | for (; dfY <= static_cast<double>(nYSize) + 2 * dfStep; dfY += dfStep) |
1174 | 0 | { |
1175 | 0 | const int iYBlock = static_cast<int>(dfY / TILE_SIZE); |
1176 | 0 | if (iYBlock != iYBlockLast) |
1177 | 0 | { |
1178 | 0 | CPLAssert(iYBlock == iYBlockLast + 1); |
1179 | 0 | if (iYBlockLast >= 0) |
1180 | 0 | yStartEnd[iYBlockLast].second = dfY + dfStep / 10; |
1181 | 0 | yStartEnd[iYBlock].first = dfY; |
1182 | 0 | iYBlockLast = iYBlock; |
1183 | 0 | } |
1184 | 0 | } |
1185 | 0 | const int iYBlock = static_cast<int>(dfY / TILE_SIZE); |
1186 | 0 | yStartEnd[iYBlock].second = dfY + dfStep / 10; |
1187 | 0 | } |
1188 | |
|
1189 | 0 | { |
1190 | 0 | int iXBlockLast = -1; |
1191 | 0 | double dfX = -dfStep; |
1192 | 0 | for (; dfX <= static_cast<double>(nXSize) + 2 * dfStep; dfX += dfStep) |
1193 | 0 | { |
1194 | 0 | const int iXBlock = static_cast<int>(dfX / TILE_SIZE); |
1195 | 0 | if (iXBlock != iXBlockLast) |
1196 | 0 | { |
1197 | 0 | CPLAssert(iXBlock == iXBlockLast + 1); |
1198 | 0 | if (iXBlockLast >= 0) |
1199 | 0 | xStartEnd[iXBlockLast].second = dfX + dfStep / 10; |
1200 | 0 | xStartEnd[iXBlock].first = dfX; |
1201 | 0 | iXBlockLast = iXBlock; |
1202 | 0 | } |
1203 | 0 | } |
1204 | 0 | const int iXBlock = static_cast<int>(dfX / TILE_SIZE); |
1205 | 0 | xStartEnd[iXBlock].second = dfX + dfStep / 10; |
1206 | 0 | } |
1207 | |
|
1208 | 0 | for (int iYBlock = 0; iYBlock < nYBlocks; ++iYBlock) |
1209 | 0 | { |
1210 | 0 | for (int iXBlock = 0; iXBlock < nXBlocks; ++iXBlock) |
1211 | 0 | { |
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 | 0 | for (double dfY = yStartEnd[iYBlock].first; |
1219 | 0 | dfY < yStartEnd[iYBlock].second; dfY += dfStep) |
1220 | 0 | { |
1221 | 0 | for (double dfX = xStartEnd[iXBlock].first; |
1222 | 0 | dfX < xStartEnd[iXBlock].second; dfX += dfStep) |
1223 | 0 | { |
1224 | | // Use forward geolocation array interpolation to compute |
1225 | | // the georeferenced position corresponding to (dfX, dfY) |
1226 | 0 | double dfGeoLocX; |
1227 | 0 | double dfGeoLocY; |
1228 | 0 | if (!PixelLineToXY(psTransform, dfX, dfY, dfGeoLocX, |
1229 | 0 | dfGeoLocY)) |
1230 | 0 | continue; |
1231 | | |
1232 | | // Compute the floating point coordinates in the pixel space |
1233 | | // of the backmap |
1234 | 0 | const double dBMX = static_cast<double>( |
1235 | 0 | (dfGeoLocX - dfMinX) / dfPixelXSize); |
1236 | |
|
1237 | 0 | const double dBMY = static_cast<double>( |
1238 | 0 | (dfMaxY - dfGeoLocY) / dfPixelYSize); |
1239 | | |
1240 | | // Get top left index by truncation |
1241 | 0 | const int iBMX = static_cast<int>(std::floor(dBMX)); |
1242 | 0 | const int iBMY = static_cast<int>(std::floor(dBMY)); |
1243 | |
|
1244 | 0 | if (iBMX >= 0 && iBMX < nBMXSize && iBMY >= 0 && |
1245 | 0 | iBMY < nBMYSize) |
1246 | 0 | { |
1247 | | // Compute the georeferenced position of the top-left |
1248 | | // index of the backmap |
1249 | 0 | double dfGeoX = dfMinX + iBMX * dfPixelXSize; |
1250 | 0 | const double dfGeoY = dfMaxY - iBMY * dfPixelYSize; |
1251 | |
|
1252 | 0 | bool bMatchingGeoLocCellFound = false; |
1253 | |
|
1254 | 0 | const int nOuterIters = |
1255 | 0 | psTransform->bGeographicSRSWithMinus180Plus180LongRange && |
1256 | 0 | fabs(dfGeoX) >= 180 |
1257 | 0 | ? 2 |
1258 | 0 | : 1; |
1259 | |
|
1260 | 0 | for (int iOuterIter = 0; iOuterIter < nOuterIters; |
1261 | 0 | ++iOuterIter) |
1262 | 0 | { |
1263 | 0 | if (iOuterIter == 1 && dfGeoX >= 180) |
1264 | 0 | dfGeoX -= 360; |
1265 | 0 | else if (iOuterIter == 1 && dfGeoX <= -180) |
1266 | 0 | dfGeoX += 360; |
1267 | | |
1268 | | // Identify a cell (quadrilateral in georeferenced |
1269 | | // space) in the geolocation array in which dfGeoX, |
1270 | | // dfGeoY falls into. |
1271 | 0 | oPoint.setX(dfGeoX); |
1272 | 0 | oPoint.setY(dfGeoY); |
1273 | 0 | const int nX = static_cast<int>(std::floor(dfX)); |
1274 | 0 | const int nY = static_cast<int>(std::floor(dfY)); |
1275 | 0 | for (int sx = -1; |
1276 | 0 | !bMatchingGeoLocCellFound && sx <= 0; sx++) |
1277 | 0 | { |
1278 | 0 | for (int sy = -1; |
1279 | 0 | !bMatchingGeoLocCellFound && sy <= 0; sy++) |
1280 | 0 | { |
1281 | 0 | const int pixel = nX + sx; |
1282 | 0 | const int line = nY + sy; |
1283 | 0 | double x0, y0, x1, y1, x2, y2, x3, y3; |
1284 | 0 | if (!PixelLineToXY(psTransform, pixel, line, |
1285 | 0 | x0, y0) || |
1286 | 0 | !PixelLineToXY(psTransform, pixel + 1, |
1287 | 0 | line, x2, y2) || |
1288 | 0 | !PixelLineToXY(psTransform, pixel, |
1289 | 0 | line + 1, x1, y1) || |
1290 | 0 | !PixelLineToXY(psTransform, pixel + 1, |
1291 | 0 | line + 1, x3, y3)) |
1292 | 0 | { |
1293 | 0 | break; |
1294 | 0 | } |
1295 | | |
1296 | 0 | int nIters = 1; |
1297 | 0 | if (psTransform |
1298 | 0 | ->bGeographicSRSWithMinus180Plus180LongRange && |
1299 | 0 | std::fabs(x0) > 170 && |
1300 | 0 | std::fabs(x1) > 170 && |
1301 | 0 | std::fabs(x2) > 170 && |
1302 | 0 | std::fabs(x3) > 170 && |
1303 | 0 | (std::fabs(x1 - x0) > 180 || |
1304 | 0 | std::fabs(x2 - x0) > 180 || |
1305 | 0 | std::fabs(x3 - x0) > 180)) |
1306 | 0 | { |
1307 | 0 | nIters = 2; |
1308 | 0 | if (x0 > 0) |
1309 | 0 | x0 -= 360; |
1310 | 0 | if (x1 > 0) |
1311 | 0 | x1 -= 360; |
1312 | 0 | if (x2 > 0) |
1313 | 0 | x2 -= 360; |
1314 | 0 | if (x3 > 0) |
1315 | 0 | x3 -= 360; |
1316 | 0 | } |
1317 | 0 | for (int iIter = 0; iIter < nIters; ++iIter) |
1318 | 0 | { |
1319 | 0 | if (iIter == 1) |
1320 | 0 | { |
1321 | 0 | x0 += 360; |
1322 | 0 | x1 += 360; |
1323 | 0 | x2 += 360; |
1324 | 0 | x3 += 360; |
1325 | 0 | } |
1326 | |
|
1327 | 0 | oRing.setPoint(0, x0, y0); |
1328 | 0 | oRing.setPoint(1, x2, y2); |
1329 | 0 | oRing.setPoint(2, x3, y3); |
1330 | 0 | oRing.setPoint(3, x1, y1); |
1331 | 0 | oRing.setPoint(4, x0, y0); |
1332 | 0 | if (oRing.isPointInRing(&oPoint) || |
1333 | 0 | oRing.isPointOnRingBoundary( |
1334 | 0 | &oPoint)) |
1335 | 0 | { |
1336 | 0 | bMatchingGeoLocCellFound = true; |
1337 | 0 | double dfBMXValue = pixel; |
1338 | 0 | double dfBMYValue = line; |
1339 | 0 | GDALInverseBilinearInterpolation( |
1340 | 0 | dfGeoX, dfGeoY, x0, y0, x1, y1, |
1341 | 0 | x2, y2, x3, y3, dfBMXValue, |
1342 | 0 | dfBMYValue); |
1343 | |
|
1344 | 0 | dfBMXValue = |
1345 | 0 | (dfBMXValue + |
1346 | 0 | dfGeorefConventionOffset) * |
1347 | 0 | psTransform->dfPIXEL_STEP + |
1348 | 0 | psTransform->dfPIXEL_OFFSET; |
1349 | 0 | dfBMYValue = |
1350 | 0 | (dfBMYValue + |
1351 | 0 | dfGeorefConventionOffset) * |
1352 | 0 | psTransform->dfLINE_STEP + |
1353 | 0 | psTransform->dfLINE_OFFSET; |
1354 | |
|
1355 | 0 | pAccessors->backMapXAccessor.Set( |
1356 | 0 | iBMX, iBMY, |
1357 | 0 | static_cast<float>(dfBMXValue)); |
1358 | 0 | pAccessors->backMapYAccessor.Set( |
1359 | 0 | iBMX, iBMY, |
1360 | 0 | static_cast<float>(dfBMYValue)); |
1361 | 0 | pAccessors->backMapWeightAccessor |
1362 | 0 | .Set(iBMX, iBMY, 1.0f); |
1363 | 0 | } |
1364 | 0 | } |
1365 | 0 | } |
1366 | 0 | } |
1367 | 0 | } |
1368 | 0 | if (bMatchingGeoLocCellFound) |
1369 | 0 | continue; |
1370 | 0 | } |
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 | 0 | if (iBMX < -1 || iBMY < -1 || iBMX > nBMXSize || |
1377 | 0 | iBMY > nBMYSize) |
1378 | 0 | continue; |
1379 | | |
1380 | 0 | const double fracBMX = dBMX - iBMX; |
1381 | 0 | const double fracBMY = dBMY - iBMY; |
1382 | | |
1383 | | // Check logic for top left pixel |
1384 | 0 | if ((iBMX >= 0) && (iBMY >= 0) && (iBMX < nBMXSize) && |
1385 | 0 | (iBMY < nBMYSize) && |
1386 | 0 | pAccessors->backMapWeightAccessor.Get(iBMX, iBMY) != |
1387 | 0 | 1.0f) |
1388 | 0 | { |
1389 | 0 | const double tempwt = (1.0 - fracBMX) * (1.0 - fracBMY); |
1390 | 0 | UpdateBackmap(iBMX, iBMY, dfX, dfY, tempwt); |
1391 | 0 | } |
1392 | | |
1393 | | // Check logic for top right pixel |
1394 | 0 | if ((iBMY >= 0) && (iBMX + 1 < nBMXSize) && |
1395 | 0 | (iBMY < nBMYSize) && |
1396 | 0 | pAccessors->backMapWeightAccessor.Get(iBMX + 1, iBMY) != |
1397 | 0 | 1.0f) |
1398 | 0 | { |
1399 | 0 | const double tempwt = fracBMX * (1.0 - fracBMY); |
1400 | 0 | UpdateBackmap(iBMX + 1, iBMY, dfX, dfY, tempwt); |
1401 | 0 | } |
1402 | | |
1403 | | // Check logic for bottom right pixel |
1404 | 0 | if ((iBMX + 1 < nBMXSize) && (iBMY + 1 < nBMYSize) && |
1405 | 0 | pAccessors->backMapWeightAccessor.Get(iBMX + 1, |
1406 | 0 | iBMY + 1) != 1.0f) |
1407 | 0 | { |
1408 | 0 | const double tempwt = fracBMX * fracBMY; |
1409 | 0 | UpdateBackmap(iBMX + 1, iBMY + 1, dfX, dfY, tempwt); |
1410 | 0 | } |
1411 | | |
1412 | | // Check logic for bottom left pixel |
1413 | 0 | if ((iBMX >= 0) && (iBMX < nBMXSize) && |
1414 | 0 | (iBMY + 1 < nBMYSize) && |
1415 | 0 | pAccessors->backMapWeightAccessor.Get(iBMX, iBMY + 1) != |
1416 | 0 | 1.0f) |
1417 | 0 | { |
1418 | 0 | const double tempwt = (1.0 - fracBMX) * fracBMY; |
1419 | 0 | UpdateBackmap(iBMX, iBMY + 1, dfX, dfY, tempwt); |
1420 | 0 | } |
1421 | 0 | } |
1422 | 0 | } |
1423 | 0 | } |
1424 | 0 | } |
1425 | | |
1426 | | // Each pixel in the backmap may have multiple entries. |
1427 | | // We now go in average it out using the weights |
1428 | 0 | START_ITER_PER_BLOCK(nBMXSize, TILE_SIZE, nBMYSize, TILE_SIZE, (void)0, |
1429 | 0 | iXStart, iXEnd, iYStart, iYEnd) |
1430 | 0 | { |
1431 | 0 | for (int iY = iYStart; iY < iYEnd; ++iY) |
1432 | 0 | { |
1433 | 0 | for (int iX = iXStart; iX < iXEnd; ++iX) |
1434 | 0 | { |
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 | 0 | const auto weight = |
1439 | 0 | pAccessors->backMapWeightAccessor.Get(iX, iY); |
1440 | 0 | if (weight > 0) |
1441 | 0 | { |
1442 | 0 | pAccessors->backMapXAccessor.Set( |
1443 | 0 | iX, iY, |
1444 | 0 | pAccessors->backMapXAccessor.Get(iX, iY) / weight); |
1445 | 0 | pAccessors->backMapYAccessor.Set( |
1446 | 0 | iX, iY, |
1447 | 0 | pAccessors->backMapYAccessor.Get(iX, iY) / weight); |
1448 | 0 | } |
1449 | 0 | else |
1450 | 0 | { |
1451 | 0 | pAccessors->backMapXAccessor.Set(iX, iY, INVALID_BMXY); |
1452 | 0 | pAccessors->backMapYAccessor.Set(iX, iY, INVALID_BMXY); |
1453 | 0 | } |
1454 | 0 | } |
1455 | 0 | } |
1456 | 0 | } |
1457 | 0 | END_ITER_PER_BLOCK |
1458 | |
|
1459 | 0 | pAccessors->FreeWghtsBackMap(); |
1460 | | |
1461 | | // Fill holes in backmap |
1462 | 0 | auto poBackmapDS = pAccessors->GetBackmapDataset(); |
1463 | |
|
1464 | 0 | 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 | 0 | constexpr double dfMaxSearchDist = 3.0; |
1477 | 0 | constexpr int nSmoothingIterations = 1; |
1478 | 0 | for (int i = 1; i <= 2; i++) |
1479 | 0 | { |
1480 | 0 | GDALFillNodata(GDALRasterBand::ToHandle(poBackmapDS->GetRasterBand(i)), |
1481 | 0 | nullptr, dfMaxSearchDist, |
1482 | 0 | 0, // unused parameter |
1483 | 0 | nSmoothingIterations, nullptr, nullptr, nullptr); |
1484 | 0 | } |
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 | 0 | struct LastValidStruct |
1498 | 0 | { |
1499 | 0 | int iX = -1; |
1500 | 0 | float bmX = 0; |
1501 | 0 | }; |
1502 | |
|
1503 | 0 | std::vector<LastValidStruct> lastValid(TILE_SIZE); |
1504 | 0 | const auto reinitLine = [&lastValid]() |
1505 | 0 | { |
1506 | 0 | const size_t nSize = lastValid.size(); |
1507 | 0 | lastValid.clear(); |
1508 | 0 | lastValid.resize(nSize); |
1509 | 0 | }; Unexecuted instantiation: GDALGeoLoc<GDALGeoLocCArrayAccessors>::GenerateBackMap(GDALGeoLocTransformInfo*)::{lambda()#1}::operator()() constUnexecuted instantiation: GDALGeoLoc<GDALGeoLocDatasetAccessors>::GenerateBackMap(GDALGeoLocTransformInfo*)::{lambda()#1}::operator()() const |
1510 | 0 | START_ITER_PER_BLOCK(nBMXSize, TILE_SIZE, nBMYSize, TILE_SIZE, reinitLine(), |
1511 | 0 | iXStart, iXEnd, iYStart, iYEnd) |
1512 | 0 | { |
1513 | 0 | const int iYCount = iYEnd - iYStart; |
1514 | 0 | for (int iYIter = 0; iYIter < iYCount; ++iYIter) |
1515 | 0 | { |
1516 | 0 | int iLastValidIX = lastValid[iYIter].iX; |
1517 | 0 | float bmXLastValid = lastValid[iYIter].bmX; |
1518 | 0 | const int iBMY = iYStart + iYIter; |
1519 | 0 | for (int iBMX = iXStart; iBMX < iXEnd; ++iBMX) |
1520 | 0 | { |
1521 | 0 | const float bmX = pAccessors->backMapXAccessor.Get(iBMX, iBMY); |
1522 | 0 | if (bmX == INVALID_BMXY) |
1523 | 0 | continue; |
1524 | 0 | if (iLastValidIX != -1 && iBMX > iLastValidIX + 1 && |
1525 | 0 | fabs(bmX - bmXLastValid) <= 2) |
1526 | 0 | { |
1527 | 0 | const float bmY = |
1528 | 0 | pAccessors->backMapYAccessor.Get(iBMX, iBMY); |
1529 | 0 | const float bmYLastValid = |
1530 | 0 | pAccessors->backMapYAccessor.Get(iLastValidIX, iBMY); |
1531 | 0 | if (fabs(bmY - bmYLastValid) <= 2) |
1532 | 0 | { |
1533 | 0 | for (int iBMXInner = iLastValidIX + 1; iBMXInner < iBMX; |
1534 | 0 | ++iBMXInner) |
1535 | 0 | { |
1536 | 0 | const float alpha = |
1537 | 0 | static_cast<float>(iBMXInner - iLastValidIX) / |
1538 | 0 | (iBMX - iLastValidIX); |
1539 | 0 | pAccessors->backMapXAccessor.Set( |
1540 | 0 | iBMXInner, iBMY, |
1541 | 0 | (1.0f - alpha) * bmXLastValid + alpha * bmX); |
1542 | 0 | pAccessors->backMapYAccessor.Set( |
1543 | 0 | iBMXInner, iBMY, |
1544 | 0 | (1.0f - alpha) * bmYLastValid + alpha * bmY); |
1545 | 0 | } |
1546 | 0 | } |
1547 | 0 | } |
1548 | 0 | iLastValidIX = iBMX; |
1549 | 0 | bmXLastValid = bmX; |
1550 | 0 | } |
1551 | 0 | lastValid[iYIter].iX = iLastValidIX; |
1552 | 0 | lastValid[iYIter].bmX = bmXLastValid; |
1553 | 0 | } |
1554 | 0 | } |
1555 | 0 | 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 | 0 | pAccessors->ReleaseBackmapDataset(poBackmapDS); |
1569 | 0 | CPLDebug("GEOLOC", "Ending backmap generation"); |
1570 | |
|
1571 | 0 | return true; |
1572 | 0 | } Unexecuted instantiation: GDALGeoLoc<GDALGeoLocCArrayAccessors>::GenerateBackMap(GDALGeoLocTransformInfo*) Unexecuted instantiation: GDALGeoLoc<GDALGeoLocDatasetAccessors>::GenerateBackMap(GDALGeoLocTransformInfo*) |
1573 | | |
1574 | | /*! @endcond */ |
1575 | | |
1576 | | /************************************************************************/ |
1577 | | /* GDALGeoLocRescale() */ |
1578 | | /************************************************************************/ |
1579 | | |
1580 | | static void GDALGeoLocRescale(char **&papszMD, const char *pszItem, |
1581 | | double dfRatio, double dfDefaultVal) |
1582 | 0 | { |
1583 | 0 | const double dfVal = |
1584 | 0 | dfRatio * CPLAtofM(CSLFetchNameValueDef( |
1585 | 0 | papszMD, pszItem, CPLSPrintf("%.17g", dfDefaultVal))); |
1586 | |
|
1587 | 0 | papszMD = CSLSetNameValue(papszMD, pszItem, CPLSPrintf("%.17g", dfVal)); |
1588 | 0 | } |
1589 | | |
1590 | | /************************************************************************/ |
1591 | | /* GDALCreateSimilarGeoLocTransformer() */ |
1592 | | /************************************************************************/ |
1593 | | |
1594 | | static void *GDALCreateSimilarGeoLocTransformer(void *hTransformArg, |
1595 | | double dfRatioX, |
1596 | | double dfRatioY) |
1597 | 0 | { |
1598 | 0 | VALIDATE_POINTER1(hTransformArg, "GDALCreateSimilarGeoLocTransformer", |
1599 | 0 | nullptr); |
1600 | | |
1601 | 0 | GDALGeoLocTransformInfo *psInfo = |
1602 | 0 | static_cast<GDALGeoLocTransformInfo *>(hTransformArg); |
1603 | |
|
1604 | 0 | char **papszGeolocationInfo = CSLDuplicate(psInfo->papszGeolocationInfo); |
1605 | |
|
1606 | 0 | if (dfRatioX != 1.0 || dfRatioY != 1.0) |
1607 | 0 | { |
1608 | 0 | GDALGeoLocRescale(papszGeolocationInfo, "PIXEL_OFFSET", dfRatioX, 0.0); |
1609 | 0 | GDALGeoLocRescale(papszGeolocationInfo, "LINE_OFFSET", dfRatioY, 0.0); |
1610 | 0 | GDALGeoLocRescale(papszGeolocationInfo, "PIXEL_STEP", 1.0 / dfRatioX, |
1611 | 0 | 1.0); |
1612 | 0 | GDALGeoLocRescale(papszGeolocationInfo, "LINE_STEP", 1.0 / dfRatioY, |
1613 | 0 | 1.0); |
1614 | 0 | } |
1615 | |
|
1616 | 0 | auto psInfoNew = |
1617 | 0 | static_cast<GDALGeoLocTransformInfo *>(GDALCreateGeoLocTransformer( |
1618 | 0 | nullptr, papszGeolocationInfo, psInfo->bReversed)); |
1619 | 0 | psInfoNew->dfOversampleFactor = psInfo->dfOversampleFactor; |
1620 | |
|
1621 | 0 | CSLDestroy(papszGeolocationInfo); |
1622 | |
|
1623 | 0 | return psInfoNew; |
1624 | 0 | } |
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 | | CPLStringList GDALCreateGeolocationMetadata(GDALDatasetH hBaseDS, |
1636 | | const char *pszGeolocationDataset, |
1637 | | bool bIsSource) |
1638 | 0 | { |
1639 | 0 | CPLStringList aosMD; |
1640 | | |
1641 | | // Open geolocation dataset |
1642 | 0 | auto poGeolocDS = std::unique_ptr<GDALDataset>( |
1643 | 0 | GDALDataset::Open(pszGeolocationDataset, GDAL_OF_RASTER)); |
1644 | 0 | if (poGeolocDS == nullptr) |
1645 | 0 | { |
1646 | 0 | CPLError(CE_Failure, CPLE_AppDefined, "Invalid dataset: %s", |
1647 | 0 | pszGeolocationDataset); |
1648 | 0 | return CPLStringList(); |
1649 | 0 | } |
1650 | 0 | const int nGeoLocXSize = poGeolocDS->GetRasterXSize(); |
1651 | 0 | const int nGeoLocYSize = poGeolocDS->GetRasterYSize(); |
1652 | 0 | if (nGeoLocXSize == 0 || nGeoLocYSize == 0) |
1653 | 0 | { |
1654 | 0 | CPLError(CE_Failure, CPLE_AppDefined, |
1655 | 0 | "Invalid dataset dimension for %s: %dx%d", |
1656 | 0 | pszGeolocationDataset, nGeoLocXSize, nGeoLocYSize); |
1657 | 0 | return CPLStringList(); |
1658 | 0 | } |
1659 | | |
1660 | | // Import the GEOLOCATION metadata from the geolocation dataset, if existing |
1661 | 0 | auto papszGeolocMDFromGeolocDS = poGeolocDS->GetMetadata("GEOLOCATION"); |
1662 | 0 | if (papszGeolocMDFromGeolocDS) |
1663 | 0 | aosMD = CSLDuplicate(papszGeolocMDFromGeolocDS); |
1664 | |
|
1665 | 0 | aosMD.SetNameValue("X_DATASET", pszGeolocationDataset); |
1666 | 0 | 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 | 0 | if (aosMD.FetchNameValue("X_BAND") == nullptr && |
1671 | 0 | aosMD.FetchNameValue("Y_BAND") == nullptr) |
1672 | 0 | { |
1673 | 0 | if (poGeolocDS->GetRasterCount() != 2) |
1674 | 0 | { |
1675 | 0 | CPLError(CE_Failure, CPLE_AppDefined, |
1676 | 0 | "Expected 2 bands for %s. Got %d", pszGeolocationDataset, |
1677 | 0 | poGeolocDS->GetRasterCount()); |
1678 | 0 | return CPLStringList(); |
1679 | 0 | } |
1680 | 0 | aosMD.SetNameValue("X_BAND", "1"); |
1681 | 0 | aosMD.SetNameValue("Y_BAND", "2"); |
1682 | 0 | } |
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 | 0 | if (aosMD.FetchNameValue("SRS") == nullptr) |
1687 | 0 | { |
1688 | 0 | auto poSRS = poGeolocDS->GetSpatialRef(); |
1689 | 0 | if (poSRS) |
1690 | 0 | { |
1691 | 0 | char *pszWKT = nullptr; |
1692 | 0 | poSRS->exportToWkt(&pszWKT); |
1693 | 0 | aosMD.SetNameValue("SRS", pszWKT); |
1694 | 0 | CPLFree(pszWKT); |
1695 | 0 | } |
1696 | 0 | } |
1697 | 0 | if (aosMD.FetchNameValue("SRS") == nullptr) |
1698 | 0 | { |
1699 | 0 | aosMD.SetNameValue("SRS", SRS_WKT_WGS84_LAT_LONG); |
1700 | 0 | } |
1701 | | |
1702 | | // Set default values for PIXEL/LINE_OFFSET/STEP if not present. |
1703 | 0 | if (aosMD.FetchNameValue("PIXEL_OFFSET") == nullptr) |
1704 | 0 | aosMD.SetNameValue("PIXEL_OFFSET", "0"); |
1705 | |
|
1706 | 0 | if (aosMD.FetchNameValue("LINE_OFFSET") == nullptr) |
1707 | 0 | aosMD.SetNameValue("LINE_OFFSET", "0"); |
1708 | |
|
1709 | 0 | if (aosMD.FetchNameValue("PIXEL_STEP") == nullptr) |
1710 | 0 | { |
1711 | 0 | aosMD.SetNameValue( |
1712 | 0 | "PIXEL_STEP", CPLSPrintf("%.17g", static_cast<double>( |
1713 | 0 | GDALGetRasterXSize(hBaseDS)) / |
1714 | 0 | nGeoLocXSize)); |
1715 | 0 | } |
1716 | |
|
1717 | 0 | if (aosMD.FetchNameValue("LINE_STEP") == nullptr) |
1718 | 0 | { |
1719 | 0 | aosMD.SetNameValue( |
1720 | 0 | "LINE_STEP", CPLSPrintf("%.17g", static_cast<double>( |
1721 | 0 | GDALGetRasterYSize(hBaseDS)) / |
1722 | 0 | nGeoLocYSize)); |
1723 | 0 | } |
1724 | |
|
1725 | 0 | if (aosMD.FetchNameValue("GEOREFERENCING_CONVENTION") == nullptr) |
1726 | 0 | { |
1727 | 0 | const char *pszConvention = |
1728 | 0 | poGeolocDS->GetMetadataItem("GEOREFERENCING_CONVENTION"); |
1729 | 0 | if (pszConvention) |
1730 | 0 | aosMD.SetNameValue("GEOREFERENCING_CONVENTION", pszConvention); |
1731 | 0 | } |
1732 | |
|
1733 | 0 | std::string osDebugMsg; |
1734 | 0 | osDebugMsg = "Synthetized GEOLOCATION metadata for "; |
1735 | 0 | osDebugMsg += bIsSource ? "source" : "target"; |
1736 | 0 | osDebugMsg += ":\n"; |
1737 | 0 | for (int i = 0; i < aosMD.size(); ++i) |
1738 | 0 | { |
1739 | 0 | osDebugMsg += " "; |
1740 | 0 | osDebugMsg += aosMD[i]; |
1741 | 0 | osDebugMsg += '\n'; |
1742 | 0 | } |
1743 | 0 | CPLDebug("GEOLOC", "%s", osDebugMsg.c_str()); |
1744 | |
|
1745 | 0 | return aosMD; |
1746 | 0 | } |
1747 | | |
1748 | | /************************************************************************/ |
1749 | | /* GDALCreateGeoLocTransformer() */ |
1750 | | /************************************************************************/ |
1751 | | |
1752 | | void *GDALCreateGeoLocTransformerEx(GDALDatasetH hBaseDS, |
1753 | | CSLConstList papszGeolocationInfo, |
1754 | | int bReversed, const char *pszSourceDataset, |
1755 | | CSLConstList papszTransformOptions) |
1756 | | |
1757 | 10 | { |
1758 | | |
1759 | 10 | if (CSLFetchNameValue(papszGeolocationInfo, "PIXEL_OFFSET") == nullptr || |
1760 | 10 | CSLFetchNameValue(papszGeolocationInfo, "LINE_OFFSET") == nullptr || |
1761 | 10 | CSLFetchNameValue(papszGeolocationInfo, "PIXEL_STEP") == nullptr || |
1762 | 10 | CSLFetchNameValue(papszGeolocationInfo, "LINE_STEP") == nullptr || |
1763 | 10 | CSLFetchNameValue(papszGeolocationInfo, "X_BAND") == nullptr || |
1764 | 10 | CSLFetchNameValue(papszGeolocationInfo, "Y_BAND") == nullptr) |
1765 | 6 | { |
1766 | 6 | CPLError(CE_Failure, CPLE_AppDefined, |
1767 | 6 | "Missing some geolocation fields in " |
1768 | 6 | "GDALCreateGeoLocTransformer()"); |
1769 | 6 | return nullptr; |
1770 | 6 | } |
1771 | | |
1772 | | /* -------------------------------------------------------------------- */ |
1773 | | /* Initialize core info. */ |
1774 | | /* -------------------------------------------------------------------- */ |
1775 | 4 | GDALGeoLocTransformInfo *psTransform = |
1776 | 4 | static_cast<GDALGeoLocTransformInfo *>( |
1777 | 4 | CPLCalloc(sizeof(GDALGeoLocTransformInfo), 1)); |
1778 | | |
1779 | 4 | psTransform->bReversed = CPL_TO_BOOL(bReversed); |
1780 | 4 | psTransform->dfOversampleFactor = std::max( |
1781 | 4 | 0.1, |
1782 | 4 | std::min(2.0, |
1783 | 4 | CPLAtof(CSLFetchNameValueDef( |
1784 | 4 | papszTransformOptions, "GEOLOC_BACKMAP_OVERSAMPLE_FACTOR", |
1785 | 4 | CPLGetConfigOption("GDAL_GEOLOC_BACKMAP_OVERSAMPLE_FACTOR", |
1786 | 4 | "1.3"))))); |
1787 | | |
1788 | 4 | memcpy(psTransform->sTI.abySignature, GDAL_GTI2_SIGNATURE, |
1789 | 4 | strlen(GDAL_GTI2_SIGNATURE)); |
1790 | 4 | psTransform->sTI.pszClassName = "GDALGeoLocTransformer"; |
1791 | 4 | psTransform->sTI.pfnTransform = GDALGeoLocTransform; |
1792 | 4 | psTransform->sTI.pfnCleanup = GDALDestroyGeoLocTransformer; |
1793 | 4 | psTransform->sTI.pfnSerialize = GDALSerializeGeoLocTransformer; |
1794 | 4 | psTransform->sTI.pfnCreateSimilar = GDALCreateSimilarGeoLocTransformer; |
1795 | | |
1796 | 4 | psTransform->papszGeolocationInfo = CSLDuplicate(papszGeolocationInfo); |
1797 | | |
1798 | | /* -------------------------------------------------------------------- */ |
1799 | | /* Pull geolocation info from the options/metadata. */ |
1800 | | /* -------------------------------------------------------------------- */ |
1801 | 4 | psTransform->dfPIXEL_OFFSET = |
1802 | 4 | CPLAtof(CSLFetchNameValue(papszGeolocationInfo, "PIXEL_OFFSET")); |
1803 | 4 | psTransform->dfLINE_OFFSET = |
1804 | 4 | CPLAtof(CSLFetchNameValue(papszGeolocationInfo, "LINE_OFFSET")); |
1805 | 4 | psTransform->dfPIXEL_STEP = |
1806 | 4 | CPLAtof(CSLFetchNameValue(papszGeolocationInfo, "PIXEL_STEP")); |
1807 | 4 | psTransform->dfLINE_STEP = |
1808 | 4 | CPLAtof(CSLFetchNameValue(papszGeolocationInfo, "LINE_STEP")); |
1809 | | |
1810 | 4 | psTransform->bOriginIsTopLeftCorner = EQUAL( |
1811 | 4 | CSLFetchNameValueDef(papszGeolocationInfo, "GEOREFERENCING_CONVENTION", |
1812 | 4 | "TOP_LEFT_CORNER"), |
1813 | 4 | "TOP_LEFT_CORNER"); |
1814 | | |
1815 | | /* -------------------------------------------------------------------- */ |
1816 | | /* Establish access to geolocation dataset(s). */ |
1817 | | /* -------------------------------------------------------------------- */ |
1818 | 4 | const char *pszDSName = |
1819 | 4 | CSLFetchNameValue(papszGeolocationInfo, "X_DATASET"); |
1820 | 4 | if (pszDSName != nullptr) |
1821 | 2 | { |
1822 | 2 | CPLConfigOptionSetter oSetter("CPL_ALLOW_VSISTDIN", "NO", true); |
1823 | 2 | if (CPLTestBool(CSLFetchNameValueDef( |
1824 | 2 | papszGeolocationInfo, "X_DATASET_RELATIVE_TO_SOURCE", "NO")) && |
1825 | 2 | (hBaseDS != nullptr || pszSourceDataset)) |
1826 | 2 | { |
1827 | 2 | const CPLString osFilename = CPLProjectRelativeFilenameSafe( |
1828 | 2 | CPLGetDirnameSafe(pszSourceDataset |
1829 | 2 | ? pszSourceDataset |
1830 | 2 | : GDALGetDescription(hBaseDS)) |
1831 | 2 | .c_str(), |
1832 | 2 | pszDSName); |
1833 | 2 | psTransform->hDS_X = |
1834 | 2 | GDALOpenShared(osFilename.c_str(), GA_ReadOnly); |
1835 | 2 | } |
1836 | 0 | else |
1837 | 0 | { |
1838 | 0 | psTransform->hDS_X = GDALOpenShared(pszDSName, GA_ReadOnly); |
1839 | 0 | } |
1840 | 2 | } |
1841 | 2 | else |
1842 | 2 | { |
1843 | 2 | psTransform->hDS_X = hBaseDS; |
1844 | 2 | if (hBaseDS) |
1845 | 0 | { |
1846 | 0 | GDALReferenceDataset(psTransform->hDS_X); |
1847 | 0 | psTransform->papszGeolocationInfo = |
1848 | 0 | CSLSetNameValue(psTransform->papszGeolocationInfo, "X_DATASET", |
1849 | 0 | GDALGetDescription(hBaseDS)); |
1850 | 0 | } |
1851 | 2 | } |
1852 | | |
1853 | 4 | pszDSName = CSLFetchNameValue(papszGeolocationInfo, "Y_DATASET"); |
1854 | 4 | if (pszDSName != nullptr) |
1855 | 4 | { |
1856 | 4 | CPLConfigOptionSetter oSetter("CPL_ALLOW_VSISTDIN", "NO", true); |
1857 | 4 | if (CPLTestBool(CSLFetchNameValueDef( |
1858 | 4 | papszGeolocationInfo, "Y_DATASET_RELATIVE_TO_SOURCE", "NO")) && |
1859 | 4 | (hBaseDS != nullptr || pszSourceDataset)) |
1860 | 4 | { |
1861 | 4 | const CPLString osFilename = CPLProjectRelativeFilenameSafe( |
1862 | 4 | CPLGetDirnameSafe(pszSourceDataset |
1863 | 4 | ? pszSourceDataset |
1864 | 4 | : GDALGetDescription(hBaseDS)) |
1865 | 4 | .c_str(), |
1866 | 4 | pszDSName); |
1867 | 4 | psTransform->hDS_Y = |
1868 | 4 | GDALOpenShared(osFilename.c_str(), GA_ReadOnly); |
1869 | 4 | } |
1870 | 0 | else |
1871 | 0 | { |
1872 | 0 | psTransform->hDS_Y = GDALOpenShared(pszDSName, GA_ReadOnly); |
1873 | 0 | } |
1874 | 4 | } |
1875 | 0 | else |
1876 | 0 | { |
1877 | 0 | psTransform->hDS_Y = hBaseDS; |
1878 | 0 | if (hBaseDS) |
1879 | 0 | { |
1880 | 0 | GDALReferenceDataset(psTransform->hDS_Y); |
1881 | 0 | psTransform->papszGeolocationInfo = |
1882 | 0 | CSLSetNameValue(psTransform->papszGeolocationInfo, "Y_DATASET", |
1883 | 0 | GDALGetDescription(hBaseDS)); |
1884 | 0 | } |
1885 | 0 | } |
1886 | | |
1887 | 4 | if (psTransform->hDS_X == nullptr || psTransform->hDS_Y == nullptr) |
1888 | 4 | { |
1889 | 4 | GDALDestroyGeoLocTransformer(psTransform); |
1890 | 4 | return nullptr; |
1891 | 4 | } |
1892 | | |
1893 | | /* -------------------------------------------------------------------- */ |
1894 | | /* Get the band handles. */ |
1895 | | /* -------------------------------------------------------------------- */ |
1896 | 0 | const int nXBand = |
1897 | 0 | std::max(1, atoi(CSLFetchNameValue(papszGeolocationInfo, "X_BAND"))); |
1898 | 0 | psTransform->hBand_X = GDALGetRasterBand(psTransform->hDS_X, nXBand); |
1899 | |
|
1900 | 0 | psTransform->dfNoDataX = GDALGetRasterNoDataValue( |
1901 | 0 | psTransform->hBand_X, &(psTransform->bHasNoData)); |
1902 | |
|
1903 | 0 | const int nYBand = |
1904 | 0 | std::max(1, atoi(CSLFetchNameValue(papszGeolocationInfo, "Y_BAND"))); |
1905 | 0 | psTransform->hBand_Y = GDALGetRasterBand(psTransform->hDS_Y, nYBand); |
1906 | |
|
1907 | 0 | if (psTransform->hBand_X == nullptr || psTransform->hBand_Y == nullptr) |
1908 | 0 | { |
1909 | 0 | GDALDestroyGeoLocTransformer(psTransform); |
1910 | 0 | return nullptr; |
1911 | 0 | } |
1912 | | |
1913 | 0 | psTransform->bSwapXY = CPLTestBool( |
1914 | 0 | CSLFetchNameValueDef(papszGeolocationInfo, "SWAP_XY", "NO")); |
1915 | | |
1916 | | /* -------------------------------------------------------------------- */ |
1917 | | /* Check that X and Y bands have the same dimensions */ |
1918 | | /* -------------------------------------------------------------------- */ |
1919 | 0 | const int nXSize_XBand = GDALGetRasterXSize(psTransform->hDS_X); |
1920 | 0 | const int nYSize_XBand = GDALGetRasterYSize(psTransform->hDS_X); |
1921 | 0 | const int nXSize_YBand = GDALGetRasterXSize(psTransform->hDS_Y); |
1922 | 0 | const int nYSize_YBand = GDALGetRasterYSize(psTransform->hDS_Y); |
1923 | 0 | if (nYSize_XBand == 1 || nYSize_YBand == 1) |
1924 | 0 | { |
1925 | 0 | if (nYSize_XBand != 1 || nYSize_YBand != 1) |
1926 | 0 | { |
1927 | 0 | CPLError(CE_Failure, CPLE_AppDefined, |
1928 | 0 | "X_BAND and Y_BAND should have both nYSize == 1"); |
1929 | 0 | GDALDestroyGeoLocTransformer(psTransform); |
1930 | 0 | return nullptr; |
1931 | 0 | } |
1932 | 0 | } |
1933 | 0 | else if (nXSize_XBand != nXSize_YBand || nYSize_XBand != nYSize_YBand) |
1934 | 0 | { |
1935 | 0 | CPLError(CE_Failure, CPLE_AppDefined, |
1936 | 0 | "X_BAND and Y_BAND do not have the same dimensions"); |
1937 | 0 | GDALDestroyGeoLocTransformer(psTransform); |
1938 | 0 | return nullptr; |
1939 | 0 | } |
1940 | | |
1941 | 0 | if (nXSize_XBand <= 0 || nYSize_XBand <= 0 || nXSize_YBand <= 0 || |
1942 | 0 | nYSize_YBand <= 0) |
1943 | 0 | { |
1944 | 0 | CPLError(CE_Failure, CPLE_AppDefined, "Invalid X_BAND / Y_BAND size"); |
1945 | 0 | GDALDestroyGeoLocTransformer(psTransform); |
1946 | 0 | return nullptr; |
1947 | 0 | } |
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 | 0 | const bool bIsRegularGrid = (nYSize_XBand == 1 && nYSize_YBand == 1); |
1953 | |
|
1954 | 0 | const int nXSize = nXSize_XBand; |
1955 | 0 | const int nYSize = bIsRegularGrid ? nXSize_YBand : nYSize_XBand; |
1956 | |
|
1957 | 0 | if (static_cast<size_t>(nXSize) > |
1958 | 0 | std::numeric_limits<size_t>::max() / static_cast<size_t>(nYSize)) |
1959 | 0 | { |
1960 | 0 | CPLError(CE_Failure, CPLE_AppDefined, "Int overflow : %d x %d", nXSize, |
1961 | 0 | nYSize); |
1962 | 0 | GDALDestroyGeoLocTransformer(psTransform); |
1963 | 0 | return nullptr; |
1964 | 0 | } |
1965 | | |
1966 | 0 | psTransform->nGeoLocXSize = nXSize; |
1967 | 0 | psTransform->nGeoLocYSize = nYSize; |
1968 | |
|
1969 | 0 | if (hBaseDS && psTransform->dfPIXEL_OFFSET == 0 && |
1970 | 0 | psTransform->dfLINE_OFFSET == 0 && psTransform->dfPIXEL_STEP == 1 && |
1971 | 0 | psTransform->dfLINE_STEP == 1) |
1972 | 0 | { |
1973 | 0 | if (GDALGetRasterXSize(hBaseDS) > nXSize || |
1974 | 0 | GDALGetRasterYSize(hBaseDS) > nYSize) |
1975 | 0 | { |
1976 | 0 | CPLError(CE_Warning, CPLE_AppDefined, |
1977 | 0 | "Geolocation array is %d x %d large, " |
1978 | 0 | "whereas dataset is %d x %d large. Result might be " |
1979 | 0 | "incorrect due to lack of values in geolocation array.", |
1980 | 0 | nXSize, nYSize, GDALGetRasterXSize(hBaseDS), |
1981 | 0 | GDALGetRasterYSize(hBaseDS)); |
1982 | 0 | } |
1983 | 0 | } |
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 | 0 | const bool bUseQuadtree = |
1992 | 0 | EQUAL(CPLGetConfigOption("GDAL_GEOLOC_INVERSE_METHOD", "BACKMAP"), |
1993 | 0 | "QUADTREE"); |
1994 | | |
1995 | | // Decide if we should C-arrays for geoloc and backmap, or on-disk |
1996 | | // temporary datasets. |
1997 | 0 | const char *pszUseTempDatasets = CSLFetchNameValueDef( |
1998 | 0 | papszTransformOptions, "GEOLOC_USE_TEMP_DATASETS", |
1999 | 0 | CPLGetConfigOption("GDAL_GEOLOC_USE_TEMP_DATASETS", nullptr)); |
2000 | 0 | if (pszUseTempDatasets) |
2001 | 0 | psTransform->bUseArray = !CPLTestBool(pszUseTempDatasets); |
2002 | 0 | else |
2003 | 0 | { |
2004 | 0 | constexpr int MEGAPIXEL_LIMIT = 24; |
2005 | 0 | psTransform->bUseArray = |
2006 | 0 | nXSize < MEGAPIXEL_LIMIT * 1000 * 1000 / nYSize; |
2007 | 0 | if (!psTransform->bUseArray) |
2008 | 0 | { |
2009 | 0 | CPLDebug("GEOLOC", |
2010 | 0 | "Using temporary GTiff backing to store backmap, because " |
2011 | 0 | "geoloc arrays require %d megapixels, exceeding the %d " |
2012 | 0 | "megapixels limit. You can set the " |
2013 | 0 | "GDAL_GEOLOC_USE_TEMP_DATASETS configuration option to " |
2014 | 0 | "NO to force RAM storage of backmap", |
2015 | 0 | static_cast<int>(static_cast<int64_t>(nXSize) * nYSize / |
2016 | 0 | (1000 * 1000)), |
2017 | 0 | MEGAPIXEL_LIMIT); |
2018 | 0 | } |
2019 | 0 | } |
2020 | |
|
2021 | 0 | if (psTransform->bUseArray) |
2022 | 0 | { |
2023 | 0 | auto pAccessors = new GDALGeoLocCArrayAccessors(psTransform); |
2024 | 0 | psTransform->pAccessors = pAccessors; |
2025 | 0 | if (!pAccessors->Load(bIsRegularGrid, bUseQuadtree)) |
2026 | 0 | { |
2027 | 0 | GDALDestroyGeoLocTransformer(psTransform); |
2028 | 0 | return nullptr; |
2029 | 0 | } |
2030 | 0 | } |
2031 | 0 | else |
2032 | 0 | { |
2033 | 0 | auto pAccessors = new GDALGeoLocDatasetAccessors(psTransform); |
2034 | 0 | psTransform->pAccessors = pAccessors; |
2035 | 0 | if (!pAccessors->Load(bIsRegularGrid, bUseQuadtree)) |
2036 | 0 | { |
2037 | 0 | GDALDestroyGeoLocTransformer(psTransform); |
2038 | 0 | return nullptr; |
2039 | 0 | } |
2040 | 0 | } |
2041 | | |
2042 | 0 | return psTransform; |
2043 | 0 | } |
2044 | | |
2045 | | /** Create GeoLocation transformer */ |
2046 | | void *GDALCreateGeoLocTransformer(GDALDatasetH hBaseDS, |
2047 | | char **papszGeolocationInfo, int bReversed) |
2048 | | |
2049 | 0 | { |
2050 | 0 | return GDALCreateGeoLocTransformerEx(hBaseDS, papszGeolocationInfo, |
2051 | 0 | bReversed, nullptr, nullptr); |
2052 | 0 | } |
2053 | | |
2054 | | /************************************************************************/ |
2055 | | /* GDALDestroyGeoLocTransformer() */ |
2056 | | /************************************************************************/ |
2057 | | |
2058 | | /** Destroy GeoLocation transformer */ |
2059 | | void GDALDestroyGeoLocTransformer(void *pTransformAlg) |
2060 | | |
2061 | 4 | { |
2062 | 4 | if (pTransformAlg == nullptr) |
2063 | 0 | return; |
2064 | | |
2065 | 4 | GDALGeoLocTransformInfo *psTransform = |
2066 | 4 | static_cast<GDALGeoLocTransformInfo *>(pTransformAlg); |
2067 | | |
2068 | 4 | CSLDestroy(psTransform->papszGeolocationInfo); |
2069 | | |
2070 | 4 | if (psTransform->bUseArray) |
2071 | 0 | delete static_cast<GDALGeoLocCArrayAccessors *>( |
2072 | 0 | psTransform->pAccessors); |
2073 | 4 | else |
2074 | 4 | delete static_cast<GDALGeoLocDatasetAccessors *>( |
2075 | 4 | psTransform->pAccessors); |
2076 | | |
2077 | 4 | if (psTransform->hDS_X != nullptr && |
2078 | 4 | GDALDereferenceDataset(psTransform->hDS_X) == 0) |
2079 | 0 | GDALClose(psTransform->hDS_X); |
2080 | | |
2081 | 4 | if (psTransform->hDS_Y != nullptr && |
2082 | 4 | GDALDereferenceDataset(psTransform->hDS_Y) == 0) |
2083 | 0 | GDALClose(psTransform->hDS_Y); |
2084 | | |
2085 | 4 | if (psTransform->hQuadTree != nullptr) |
2086 | 0 | CPLQuadTreeDestroy(psTransform->hQuadTree); |
2087 | | |
2088 | 4 | CPLFree(pTransformAlg); |
2089 | 4 | } |
2090 | | |
2091 | | /** Use GeoLocation transformer */ |
2092 | | int GDALGeoLocTransform(void *pTransformArg, int bDstToSrc, int nPointCount, |
2093 | | double *padfX, double *padfY, double *padfZ, |
2094 | | int *panSuccess) |
2095 | 0 | { |
2096 | 0 | GDALGeoLocTransformInfo *psTransform = |
2097 | 0 | static_cast<GDALGeoLocTransformInfo *>(pTransformArg); |
2098 | 0 | if (psTransform->bUseArray) |
2099 | 0 | { |
2100 | 0 | return GDALGeoLoc<GDALGeoLocCArrayAccessors>::Transform( |
2101 | 0 | pTransformArg, bDstToSrc, nPointCount, padfX, padfY, padfZ, |
2102 | 0 | panSuccess); |
2103 | 0 | } |
2104 | 0 | else |
2105 | 0 | { |
2106 | 0 | return GDALGeoLoc<GDALGeoLocDatasetAccessors>::Transform( |
2107 | 0 | pTransformArg, bDstToSrc, nPointCount, padfX, padfY, padfZ, |
2108 | 0 | panSuccess); |
2109 | 0 | } |
2110 | 0 | } |
2111 | | |
2112 | | /************************************************************************/ |
2113 | | /* GDALSerializeGeoLocTransformer() */ |
2114 | | /************************************************************************/ |
2115 | | |
2116 | | CPLXMLNode *GDALSerializeGeoLocTransformer(void *pTransformArg) |
2117 | | |
2118 | 0 | { |
2119 | 0 | VALIDATE_POINTER1(pTransformArg, "GDALSerializeGeoLocTransformer", nullptr); |
2120 | | |
2121 | 0 | GDALGeoLocTransformInfo *psInfo = |
2122 | 0 | static_cast<GDALGeoLocTransformInfo *>(pTransformArg); |
2123 | |
|
2124 | 0 | CPLXMLNode *psTree = |
2125 | 0 | CPLCreateXMLNode(nullptr, CXT_Element, "GeoLocTransformer"); |
2126 | | |
2127 | | /* -------------------------------------------------------------------- */ |
2128 | | /* Serialize bReversed. */ |
2129 | | /* -------------------------------------------------------------------- */ |
2130 | 0 | CPLCreateXMLElementAndValue( |
2131 | 0 | psTree, "Reversed", |
2132 | 0 | CPLString().Printf("%d", static_cast<int>(psInfo->bReversed))); |
2133 | | |
2134 | | /* -------------------------------------------------------------------- */ |
2135 | | /* geoloc metadata. */ |
2136 | | /* -------------------------------------------------------------------- */ |
2137 | 0 | char **papszMD = psInfo->papszGeolocationInfo; |
2138 | 0 | CPLXMLNode *psMD = CPLCreateXMLNode(psTree, CXT_Element, "Metadata"); |
2139 | |
|
2140 | 0 | for (int i = 0; papszMD != nullptr && papszMD[i] != nullptr; i++) |
2141 | 0 | { |
2142 | 0 | char *pszKey = nullptr; |
2143 | 0 | const char *pszRawValue = CPLParseNameValue(papszMD[i], &pszKey); |
2144 | |
|
2145 | 0 | CPLXMLNode *psMDI = CPLCreateXMLNode(psMD, CXT_Element, "MDI"); |
2146 | 0 | CPLSetXMLValue(psMDI, "#key", pszKey); |
2147 | 0 | CPLCreateXMLNode(psMDI, CXT_Text, pszRawValue); |
2148 | |
|
2149 | 0 | CPLFree(pszKey); |
2150 | 0 | } |
2151 | |
|
2152 | 0 | return psTree; |
2153 | 0 | } |
2154 | | |
2155 | | /************************************************************************/ |
2156 | | /* GDALDeserializeGeoLocTransformer() */ |
2157 | | /************************************************************************/ |
2158 | | |
2159 | | void *GDALDeserializeGeoLocTransformer(CPLXMLNode *psTree) |
2160 | | |
2161 | 10 | { |
2162 | | /* -------------------------------------------------------------------- */ |
2163 | | /* Collect metadata. */ |
2164 | | /* -------------------------------------------------------------------- */ |
2165 | 10 | CPLXMLNode *psMetadata = CPLGetXMLNode(psTree, "Metadata"); |
2166 | | |
2167 | 10 | if (psMetadata == nullptr || psMetadata->eType != CXT_Element || |
2168 | 10 | !EQUAL(psMetadata->pszValue, "Metadata")) |
2169 | 0 | return nullptr; |
2170 | | |
2171 | 10 | char **papszMD = nullptr; |
2172 | | |
2173 | 60 | for (CPLXMLNode *psMDI = psMetadata->psChild; psMDI != nullptr; |
2174 | 50 | psMDI = psMDI->psNext) |
2175 | 50 | { |
2176 | 50 | if (!EQUAL(psMDI->pszValue, "MDI") || psMDI->eType != CXT_Element || |
2177 | 50 | psMDI->psChild == nullptr || psMDI->psChild->psNext == nullptr || |
2178 | 50 | psMDI->psChild->eType != CXT_Attribute || |
2179 | 50 | psMDI->psChild->psChild == nullptr) |
2180 | 0 | continue; |
2181 | | |
2182 | 50 | papszMD = CSLSetNameValue(papszMD, psMDI->psChild->psChild->pszValue, |
2183 | 50 | psMDI->psChild->psNext->pszValue); |
2184 | 50 | } |
2185 | | |
2186 | | /* -------------------------------------------------------------------- */ |
2187 | | /* Get other flags. */ |
2188 | | /* -------------------------------------------------------------------- */ |
2189 | 10 | const int bReversed = atoi(CPLGetXMLValue(psTree, "Reversed", "0")); |
2190 | | |
2191 | | /* -------------------------------------------------------------------- */ |
2192 | | /* Generate transformation. */ |
2193 | | /* -------------------------------------------------------------------- */ |
2194 | | |
2195 | 10 | const char *pszSourceDataset = |
2196 | 10 | CPLGetXMLValue(psTree, "SourceDataset", nullptr); |
2197 | | |
2198 | 10 | void *pResult = GDALCreateGeoLocTransformerEx(nullptr, papszMD, bReversed, |
2199 | 10 | pszSourceDataset, nullptr); |
2200 | | |
2201 | | /* -------------------------------------------------------------------- */ |
2202 | | /* Cleanup GCP copy. */ |
2203 | | /* -------------------------------------------------------------------- */ |
2204 | 10 | CSLDestroy(papszMD); |
2205 | | |
2206 | 10 | return pResult; |
2207 | 10 | } |