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