Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: Implements Geolocation array based transformer, using a quadtree
5 : * for inverse
6 : * Author: Even Rouault, <even.rouault at spatialys.com>
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2022, Planet Labs
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "gdalgeoloc.h"
15 : #include "gdalgeolocquadtree.h"
16 :
17 : #include "cpl_quad_tree.h"
18 :
19 : #include "ogr_geometry.h"
20 :
21 : #include <algorithm>
22 : #include <cstddef>
23 : #include <limits>
24 :
25 : /************************************************************************/
26 : /* GDALGeoLocQuadTreeGetFeatureCorners() */
27 : /************************************************************************/
28 :
29 : static bool
30 3200 : GDALGeoLocQuadTreeGetFeatureCorners(const GDALGeoLocTransformInfo *psTransform,
31 : size_t nIdx, double &x0, double &y0,
32 : double &x1, double &y1, double &x2,
33 : double &y2, double &x3, double &y3)
34 : {
35 6400 : const size_t nExtendedWidth = psTransform->nGeoLocXSize +
36 3200 : (psTransform->bOriginIsTopLeftCorner ? 0 : 1);
37 3200 : int nX = static_cast<int>(nIdx % nExtendedWidth);
38 3200 : int nY = static_cast<int>(nIdx / nExtendedWidth);
39 :
40 3200 : if (!psTransform->bOriginIsTopLeftCorner)
41 : {
42 1681 : nX--;
43 1681 : nY--;
44 : }
45 :
46 3200 : return GDALGeoLocExtractSquare(psTransform, static_cast<int>(nX),
47 : static_cast<int>(nY), x0, y0, x1, y1, x2, y2,
48 3200 : x3, y3);
49 : }
50 :
51 : /************************************************************************/
52 : /* GDALGeoLocQuadTreeGetFeatureBounds() */
53 : /************************************************************************/
54 :
55 : constexpr size_t BIT_IDX_RANGE_180 = 8 * sizeof(size_t) - 1;
56 : constexpr size_t BIT_IDX_RANGE_180_SET = static_cast<size_t>(1)
57 : << BIT_IDX_RANGE_180;
58 :
59 : // Callback used by quadtree to retrieve the bounding box, in georeferenced
60 : // space, of a cell of the geolocation array.
61 2114 : static void GDALGeoLocQuadTreeGetFeatureBounds(const void *hFeature,
62 : void *pUserData,
63 : CPLRectObj *pBounds)
64 : {
65 2114 : const GDALGeoLocTransformInfo *psTransform =
66 : static_cast<const GDALGeoLocTransformInfo *>(pUserData);
67 2114 : size_t nIdx = reinterpret_cast<size_t>(hFeature);
68 : // Most significant bit set means that geometries crossing the antimeridian
69 : // should have their longitudes lower or greater than 180 deg.
70 2114 : const bool bXRefAt180 = (nIdx >> BIT_IDX_RANGE_180) != 0;
71 : // Clear that bit.
72 2114 : nIdx &= ~BIT_IDX_RANGE_180_SET;
73 :
74 2114 : double x0 = 0, y0 = 0, x1 = 0, y1 = 0, x2 = 0, y2 = 0, x3 = 0, y3 = 0;
75 2114 : GDALGeoLocQuadTreeGetFeatureCorners(psTransform, nIdx, x0, y0, x1, y1, x2,
76 : y2, x3, y3);
77 :
78 2114 : if (psTransform->bGeographicSRSWithMinus180Plus180LongRange &&
79 2114 : std::fabs(x0) > 170 && std::fabs(x1) > 170 && std::fabs(x2) > 170 &&
80 0 : std::fabs(x3) > 170 &&
81 0 : (std::fabs(x1 - x0) > 180 || std::fabs(x2 - x0) > 180 ||
82 0 : std::fabs(x3 - x0) > 180))
83 : {
84 0 : const double dfXRef = bXRefAt180 ? 180 : -180;
85 0 : x0 = ShiftGeoX(psTransform, dfXRef, x0);
86 0 : x1 = ShiftGeoX(psTransform, dfXRef, x1);
87 0 : x2 = ShiftGeoX(psTransform, dfXRef, x2);
88 0 : x3 = ShiftGeoX(psTransform, dfXRef, x3);
89 : }
90 2114 : pBounds->minx = std::min(std::min(x0, x1), std::min(x2, x3));
91 2114 : pBounds->miny = std::min(std::min(y0, y1), std::min(y2, y3));
92 2114 : pBounds->maxx = std::max(std::max(x0, x1), std::max(x2, x3));
93 2114 : pBounds->maxy = std::max(std::max(y0, y1), std::max(y2, y3));
94 2114 : }
95 :
96 : /************************************************************************/
97 : /* GDALGeoLocBuildQuadTree() */
98 : /************************************************************************/
99 :
100 4 : bool GDALGeoLocBuildQuadTree(GDALGeoLocTransformInfo *psTransform)
101 : {
102 : // For the pixel-center convention, insert a "virtual" row and column
103 : // at top and left of the geoloc array.
104 4 : const int nExtraPixel = psTransform->bOriginIsTopLeftCorner ? 0 : 1;
105 :
106 12 : if (psTransform->nGeoLocXSize > INT_MAX - nExtraPixel ||
107 8 : psTransform->nGeoLocYSize > INT_MAX - nExtraPixel ||
108 : // The >> 1 shift is because we need to reserve the most-significant-bit
109 : // for the second 'version' of anti-meridian crossing quadrilaterals.
110 : // See below
111 4 : static_cast<size_t>(psTransform->nGeoLocXSize + nExtraPixel) >
112 4 : (std::numeric_limits<size_t>::max() >> 1) /
113 4 : static_cast<size_t>(psTransform->nGeoLocYSize + nExtraPixel))
114 : {
115 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too big geolocation array");
116 0 : return false;
117 : }
118 :
119 4 : const int nExtendedWidth = psTransform->nGeoLocXSize + nExtraPixel;
120 4 : const int nExtendedHeight = psTransform->nGeoLocYSize + nExtraPixel;
121 4 : const size_t nExtendedXYCount =
122 4 : static_cast<size_t>(nExtendedWidth) * nExtendedHeight;
123 :
124 4 : CPLDebug("GEOLOC", "Start quadtree construction");
125 :
126 : CPLRectObj globalBounds;
127 4 : globalBounds.minx = psTransform->dfMinX;
128 4 : globalBounds.miny = psTransform->dfMinY;
129 4 : globalBounds.maxx = psTransform->dfMaxX;
130 4 : globalBounds.maxy = psTransform->dfMaxY;
131 4 : psTransform->hQuadTree = CPLQuadTreeCreateEx(
132 : &globalBounds, GDALGeoLocQuadTreeGetFeatureBounds, psTransform);
133 :
134 4 : CPLQuadTreeForceUseOfSubNodes(psTransform->hQuadTree);
135 :
136 1066 : for (size_t i = 0; i < nExtendedXYCount; i++)
137 : {
138 : double x0, y0, x1, y1, x2, y2, x3, y3;
139 1062 : if (!GDALGeoLocQuadTreeGetFeatureCorners(psTransform, i, x0, y0, x1, y1,
140 : x2, y2, x3, y3))
141 : {
142 0 : continue;
143 : }
144 :
145 : // Skip too large geometries (typically at very high latitudes)
146 : // that would fill too many nodes in the quadtree
147 1062 : if (psTransform->bGeographicSRSWithMinus180Plus180LongRange &&
148 1062 : (std::fabs(x0) > 170 || std::fabs(x1) > 170 ||
149 1062 : std::fabs(x2) > 170 || std::fabs(x3) > 170) &&
150 0 : (std::fabs(x1 - x0) > 180 || std::fabs(x2 - x0) > 180 ||
151 0 : std::fabs(x3 - x0) > 180) &&
152 0 : !(std::fabs(x0) > 170 && std::fabs(x1) > 170 &&
153 0 : std::fabs(x2) > 170 && std::fabs(x3) > 170))
154 : {
155 0 : continue;
156 : }
157 :
158 1062 : CPLQuadTreeInsert(psTransform->hQuadTree,
159 : reinterpret_cast<void *>(static_cast<uintptr_t>(i)));
160 :
161 : // For a geometry crossing the antimeridian, we've insert before
162 : // the "version" around -180 deg. Insert its corresponding version
163 : // around +180 deg.
164 1062 : if (psTransform->bGeographicSRSWithMinus180Plus180LongRange &&
165 1062 : std::fabs(x0) > 170 && std::fabs(x1) > 170 && std::fabs(x2) > 170 &&
166 0 : std::fabs(x3) > 170 &&
167 0 : (std::fabs(x1 - x0) > 180 || std::fabs(x2 - x0) > 180 ||
168 0 : std::fabs(x3 - x0) > 180))
169 : {
170 0 : CPLQuadTreeInsert(psTransform->hQuadTree,
171 : reinterpret_cast<void *>(static_cast<uintptr_t>(
172 0 : i | BIT_IDX_RANGE_180_SET)));
173 : }
174 : }
175 :
176 4 : CPLDebug("GEOLOC", "End of quadtree construction");
177 :
178 : #ifdef DEBUG_GEOLOC
179 : int nFeatureCount = 0;
180 : int nNodeCount = 0;
181 : int nMaxDepth = 0;
182 : int nMaxBucketCapacity = 0;
183 : CPLQuadTreeGetStats(psTransform->hQuadTree, &nFeatureCount, &nNodeCount,
184 : &nMaxDepth, &nMaxBucketCapacity);
185 : CPLDebug("GEOLOC", "Quadtree stats:");
186 : CPLDebug("GEOLOC", " nFeatureCount = %d", nFeatureCount);
187 : CPLDebug("GEOLOC", " nNodeCount = %d", nNodeCount);
188 : CPLDebug("GEOLOC", " nMaxDepth = %d", nMaxDepth);
189 : CPLDebug("GEOLOC", " nMaxBucketCapacity = %d", nMaxBucketCapacity);
190 : #endif
191 :
192 4 : return true;
193 : }
194 :
195 : /************************************************************************/
196 : /* GDALGeoLocInverseTransformQuadtree() */
197 : /************************************************************************/
198 :
199 24 : void GDALGeoLocInverseTransformQuadtree(
200 : const GDALGeoLocTransformInfo *psTransform, int nPointCount, double *padfX,
201 : double *padfY, int *panSuccess)
202 : {
203 : // Keep those objects in this outer scope, so they are re-used, to
204 : // save memory allocations.
205 48 : OGRPoint oPoint;
206 48 : OGRLinearRing oRing;
207 24 : oRing.setNumPoints(5);
208 :
209 24 : const double dfGeorefConventionOffset =
210 24 : psTransform->bOriginIsTopLeftCorner ? 0 : 0.5;
211 :
212 48 : for (int i = 0; i < nPointCount; i++)
213 : {
214 24 : if (padfX[i] == HUGE_VAL || padfY[i] == HUGE_VAL)
215 : {
216 0 : panSuccess[i] = FALSE;
217 0 : continue;
218 : }
219 :
220 24 : if (psTransform->bSwapXY)
221 : {
222 0 : std::swap(padfX[i], padfY[i]);
223 : }
224 :
225 24 : const double dfGeoX = padfX[i];
226 24 : const double dfGeoY = padfY[i];
227 :
228 24 : bool bDone = false;
229 :
230 : CPLRectObj aoi;
231 24 : aoi.minx = dfGeoX;
232 24 : aoi.maxx = dfGeoX;
233 24 : aoi.miny = dfGeoY;
234 24 : aoi.maxy = dfGeoY;
235 24 : int nFeatureCount = 0;
236 : void **pahFeatures =
237 24 : CPLQuadTreeSearch(psTransform->hQuadTree, &aoi, &nFeatureCount);
238 24 : if (nFeatureCount != 0)
239 : {
240 24 : oPoint.setX(dfGeoX);
241 24 : oPoint.setY(dfGeoY);
242 24 : for (int iFeat = 0; iFeat < nFeatureCount; iFeat++)
243 : {
244 24 : size_t nIdx = reinterpret_cast<size_t>(pahFeatures[iFeat]);
245 24 : const bool bXRefAt180 = (nIdx >> BIT_IDX_RANGE_180) != 0;
246 : // Clear that bit.
247 24 : nIdx &= ~BIT_IDX_RANGE_180_SET;
248 :
249 24 : double x0 = 0, y0 = 0, x1 = 0, y1 = 0, x2 = 0, y2 = 0, x3 = 0,
250 24 : y3 = 0;
251 24 : GDALGeoLocQuadTreeGetFeatureCorners(psTransform, nIdx, x0, y0,
252 : x2, y2, x1, y1, x3, y3);
253 :
254 24 : if (psTransform->bGeographicSRSWithMinus180Plus180LongRange &&
255 24 : std::fabs(x0) > 170 && std::fabs(x1) > 170 &&
256 0 : std::fabs(x2) > 170 && std::fabs(x3) > 170 &&
257 0 : (std::fabs(x1 - x0) > 180 || std::fabs(x2 - x0) > 180 ||
258 0 : std::fabs(x3 - x0) > 180))
259 : {
260 0 : const double dfXRef = bXRefAt180 ? 180 : -180;
261 0 : x0 = ShiftGeoX(psTransform, dfXRef, x0);
262 0 : x1 = ShiftGeoX(psTransform, dfXRef, x1);
263 0 : x2 = ShiftGeoX(psTransform, dfXRef, x2);
264 0 : x3 = ShiftGeoX(psTransform, dfXRef, x3);
265 : }
266 :
267 24 : oRing.setPoint(0, x0, y0);
268 24 : oRing.setPoint(1, x2, y2);
269 24 : oRing.setPoint(2, x3, y3);
270 24 : oRing.setPoint(3, x1, y1);
271 24 : oRing.setPoint(4, x0, y0);
272 :
273 32 : if (oRing.isPointInRing(&oPoint) ||
274 8 : oRing.isPointOnRingBoundary(&oPoint))
275 : {
276 24 : const size_t nExtendedWidth =
277 48 : psTransform->nGeoLocXSize +
278 24 : (psTransform->bOriginIsTopLeftCorner ? 0 : 1);
279 24 : double dfX = static_cast<double>(nIdx % nExtendedWidth);
280 : // store the result as int, and then cast to double, to
281 : // avoid Coverity Scan warning about
282 : // UNINTENDED_INTEGER_DIVISION
283 24 : const size_t nY = nIdx / nExtendedWidth;
284 24 : double dfY = static_cast<double>(nY);
285 24 : if (!psTransform->bOriginIsTopLeftCorner)
286 : {
287 12 : dfX -= 1.0;
288 12 : dfY -= 1.0;
289 : }
290 24 : GDALInverseBilinearInterpolation(dfGeoX, dfGeoY, x0, y0, x1,
291 : y1, x2, y2, x3, y3, dfX,
292 : dfY);
293 :
294 24 : dfX = (dfX + dfGeorefConventionOffset) *
295 24 : psTransform->dfPIXEL_STEP +
296 24 : psTransform->dfPIXEL_OFFSET;
297 24 : dfY = (dfY + dfGeorefConventionOffset) *
298 24 : psTransform->dfLINE_STEP +
299 24 : psTransform->dfLINE_OFFSET;
300 :
301 24 : bDone = true;
302 24 : panSuccess[i] = TRUE;
303 24 : padfX[i] = dfX;
304 24 : padfY[i] = dfY;
305 24 : break;
306 : }
307 : }
308 : }
309 24 : CPLFree(pahFeatures);
310 :
311 24 : if (!bDone)
312 : {
313 0 : panSuccess[i] = FALSE;
314 0 : padfX[i] = HUGE_VAL;
315 0 : padfY[i] = HUGE_VAL;
316 : }
317 : }
318 24 : }
|
