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 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include "gdalgeoloc.h"
31 : #include "gdalgeolocquadtree.h"
32 :
33 : #include "cpl_quad_tree.h"
34 :
35 : #include "ogr_geometry.h"
36 :
37 : #include <algorithm>
38 : #include <cstddef>
39 : #include <limits>
40 :
41 : /************************************************************************/
42 : /* GDALGeoLocQuadTreeGetFeatureCorners() */
43 : /************************************************************************/
44 :
45 : static bool
46 3200 : GDALGeoLocQuadTreeGetFeatureCorners(const GDALGeoLocTransformInfo *psTransform,
47 : size_t nIdx, double &x0, double &y0,
48 : double &x1, double &y1, double &x2,
49 : double &y2, double &x3, double &y3)
50 : {
51 6400 : const size_t nExtendedWidth = psTransform->nGeoLocXSize +
52 3200 : (psTransform->bOriginIsTopLeftCorner ? 0 : 1);
53 3200 : int nX = static_cast<int>(nIdx % nExtendedWidth);
54 3200 : int nY = static_cast<int>(nIdx / nExtendedWidth);
55 :
56 3200 : if (!psTransform->bOriginIsTopLeftCorner)
57 : {
58 1681 : nX--;
59 1681 : nY--;
60 : }
61 :
62 3200 : return GDALGeoLocExtractSquare(psTransform, static_cast<int>(nX),
63 : static_cast<int>(nY), x0, y0, x1, y1, x2, y2,
64 3200 : x3, y3);
65 : }
66 :
67 : /************************************************************************/
68 : /* GDALGeoLocQuadTreeGetFeatureBounds() */
69 : /************************************************************************/
70 :
71 : constexpr size_t BIT_IDX_RANGE_180 = 8 * sizeof(size_t) - 1;
72 : constexpr size_t BIT_IDX_RANGE_180_SET = static_cast<size_t>(1)
73 : << BIT_IDX_RANGE_180;
74 :
75 : // Callback used by quadtree to retrieve the bounding box, in georeferenced
76 : // space, of a cell of the geolocation array.
77 2114 : static void GDALGeoLocQuadTreeGetFeatureBounds(const void *hFeature,
78 : void *pUserData,
79 : CPLRectObj *pBounds)
80 : {
81 2114 : const GDALGeoLocTransformInfo *psTransform =
82 : static_cast<const GDALGeoLocTransformInfo *>(pUserData);
83 2114 : size_t nIdx = reinterpret_cast<size_t>(hFeature);
84 : // Most significant bit set means that geometries crossing the antimeridian
85 : // should have their longitudes lower or greater than 180 deg.
86 2114 : const bool bXRefAt180 = (nIdx >> BIT_IDX_RANGE_180) != 0;
87 : // Clear that bit.
88 2114 : nIdx &= ~BIT_IDX_RANGE_180_SET;
89 :
90 2114 : double x0 = 0, y0 = 0, x1 = 0, y1 = 0, x2 = 0, y2 = 0, x3 = 0, y3 = 0;
91 2114 : GDALGeoLocQuadTreeGetFeatureCorners(psTransform, nIdx, x0, y0, x1, y1, x2,
92 : y2, x3, y3);
93 :
94 2114 : if (psTransform->bGeographicSRSWithMinus180Plus180LongRange &&
95 2114 : std::fabs(x0) > 170 && std::fabs(x1) > 170 && std::fabs(x2) > 170 &&
96 0 : std::fabs(x3) > 170 &&
97 0 : (std::fabs(x1 - x0) > 180 || std::fabs(x2 - x0) > 180 ||
98 0 : std::fabs(x3 - x0) > 180))
99 : {
100 0 : const double dfXRef = bXRefAt180 ? 180 : -180;
101 0 : x0 = ShiftGeoX(psTransform, dfXRef, x0);
102 0 : x1 = ShiftGeoX(psTransform, dfXRef, x1);
103 0 : x2 = ShiftGeoX(psTransform, dfXRef, x2);
104 0 : x3 = ShiftGeoX(psTransform, dfXRef, x3);
105 : }
106 2114 : pBounds->minx = std::min(std::min(x0, x1), std::min(x2, x3));
107 2114 : pBounds->miny = std::min(std::min(y0, y1), std::min(y2, y3));
108 2114 : pBounds->maxx = std::max(std::max(x0, x1), std::max(x2, x3));
109 2114 : pBounds->maxy = std::max(std::max(y0, y1), std::max(y2, y3));
110 2114 : }
111 :
112 : /************************************************************************/
113 : /* GDALGeoLocBuildQuadTree() */
114 : /************************************************************************/
115 :
116 4 : bool GDALGeoLocBuildQuadTree(GDALGeoLocTransformInfo *psTransform)
117 : {
118 : // For the pixel-center convention, insert a "virtual" row and column
119 : // at top and left of the geoloc array.
120 4 : const int nExtraPixel = psTransform->bOriginIsTopLeftCorner ? 0 : 1;
121 :
122 12 : if (psTransform->nGeoLocXSize > INT_MAX - nExtraPixel ||
123 8 : psTransform->nGeoLocYSize > INT_MAX - nExtraPixel ||
124 : // The >> 1 shift is because we need to reserve the most-significant-bit
125 : // for the second 'version' of anti-meridian crossing quadrilaterals.
126 : // See below
127 4 : static_cast<size_t>(psTransform->nGeoLocXSize + nExtraPixel) >
128 4 : (std::numeric_limits<size_t>::max() >> 1) /
129 4 : static_cast<size_t>(psTransform->nGeoLocYSize + nExtraPixel))
130 : {
131 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too big geolocation array");
132 0 : return false;
133 : }
134 :
135 4 : const int nExtendedWidth = psTransform->nGeoLocXSize + nExtraPixel;
136 4 : const int nExtendedHeight = psTransform->nGeoLocYSize + nExtraPixel;
137 4 : const size_t nExtendedXYCount =
138 4 : static_cast<size_t>(nExtendedWidth) * nExtendedHeight;
139 :
140 4 : CPLDebug("GEOLOC", "Start quadtree construction");
141 :
142 : CPLRectObj globalBounds;
143 4 : globalBounds.minx = psTransform->dfMinX;
144 4 : globalBounds.miny = psTransform->dfMinY;
145 4 : globalBounds.maxx = psTransform->dfMaxX;
146 4 : globalBounds.maxy = psTransform->dfMaxY;
147 4 : psTransform->hQuadTree = CPLQuadTreeCreateEx(
148 : &globalBounds, GDALGeoLocQuadTreeGetFeatureBounds, psTransform);
149 :
150 4 : CPLQuadTreeForceUseOfSubNodes(psTransform->hQuadTree);
151 :
152 1066 : for (size_t i = 0; i < nExtendedXYCount; i++)
153 : {
154 : double x0, y0, x1, y1, x2, y2, x3, y3;
155 1062 : if (!GDALGeoLocQuadTreeGetFeatureCorners(psTransform, i, x0, y0, x1, y1,
156 : x2, y2, x3, y3))
157 : {
158 0 : continue;
159 : }
160 :
161 : // Skip too large geometries (typically at very high latitudes)
162 : // that would fill too many nodes in the quadtree
163 1062 : if (psTransform->bGeographicSRSWithMinus180Plus180LongRange &&
164 1062 : (std::fabs(x0) > 170 || std::fabs(x1) > 170 ||
165 1062 : std::fabs(x2) > 170 || std::fabs(x3) > 170) &&
166 0 : (std::fabs(x1 - x0) > 180 || std::fabs(x2 - x0) > 180 ||
167 0 : std::fabs(x3 - x0) > 180) &&
168 0 : !(std::fabs(x0) > 170 && std::fabs(x1) > 170 &&
169 0 : std::fabs(x2) > 170 && std::fabs(x3) > 170))
170 : {
171 0 : continue;
172 : }
173 :
174 1062 : CPLQuadTreeInsert(psTransform->hQuadTree,
175 : reinterpret_cast<void *>(static_cast<uintptr_t>(i)));
176 :
177 : // For a geometry crossing the antimeridian, we've insert before
178 : // the "version" around -180 deg. Insert its corresponding version
179 : // around +180 deg.
180 1062 : if (psTransform->bGeographicSRSWithMinus180Plus180LongRange &&
181 1062 : std::fabs(x0) > 170 && std::fabs(x1) > 170 && std::fabs(x2) > 170 &&
182 0 : std::fabs(x3) > 170 &&
183 0 : (std::fabs(x1 - x0) > 180 || std::fabs(x2 - x0) > 180 ||
184 0 : std::fabs(x3 - x0) > 180))
185 : {
186 0 : CPLQuadTreeInsert(psTransform->hQuadTree,
187 : reinterpret_cast<void *>(static_cast<uintptr_t>(
188 0 : i | BIT_IDX_RANGE_180_SET)));
189 : }
190 : }
191 :
192 4 : CPLDebug("GEOLOC", "End of quadtree construction");
193 :
194 : #ifdef DEBUG_GEOLOC
195 : int nFeatureCount = 0;
196 : int nNodeCount = 0;
197 : int nMaxDepth = 0;
198 : int nMaxBucketCapacity = 0;
199 : CPLQuadTreeGetStats(psTransform->hQuadTree, &nFeatureCount, &nNodeCount,
200 : &nMaxDepth, &nMaxBucketCapacity);
201 : CPLDebug("GEOLOC", "Quadtree stats:");
202 : CPLDebug("GEOLOC", " nFeatureCount = %d", nFeatureCount);
203 : CPLDebug("GEOLOC", " nNodeCount = %d", nNodeCount);
204 : CPLDebug("GEOLOC", " nMaxDepth = %d", nMaxDepth);
205 : CPLDebug("GEOLOC", " nMaxBucketCapacity = %d", nMaxBucketCapacity);
206 : #endif
207 :
208 4 : return true;
209 : }
210 :
211 : /************************************************************************/
212 : /* GDALGeoLocInverseTransformQuadtree() */
213 : /************************************************************************/
214 :
215 24 : void GDALGeoLocInverseTransformQuadtree(
216 : const GDALGeoLocTransformInfo *psTransform, int nPointCount, double *padfX,
217 : double *padfY, int *panSuccess)
218 : {
219 : // Keep those objects in this outer scope, so they are re-used, to
220 : // save memory allocations.
221 48 : OGRPoint oPoint;
222 48 : OGRLinearRing oRing;
223 24 : oRing.setNumPoints(5);
224 :
225 24 : const double dfGeorefConventionOffset =
226 24 : psTransform->bOriginIsTopLeftCorner ? 0 : 0.5;
227 :
228 48 : for (int i = 0; i < nPointCount; i++)
229 : {
230 24 : if (padfX[i] == HUGE_VAL || padfY[i] == HUGE_VAL)
231 : {
232 0 : panSuccess[i] = FALSE;
233 0 : continue;
234 : }
235 :
236 24 : if (psTransform->bSwapXY)
237 : {
238 0 : std::swap(padfX[i], padfY[i]);
239 : }
240 :
241 24 : const double dfGeoX = padfX[i];
242 24 : const double dfGeoY = padfY[i];
243 :
244 24 : bool bDone = false;
245 :
246 : CPLRectObj aoi;
247 24 : aoi.minx = dfGeoX;
248 24 : aoi.maxx = dfGeoX;
249 24 : aoi.miny = dfGeoY;
250 24 : aoi.maxy = dfGeoY;
251 24 : int nFeatureCount = 0;
252 : void **pahFeatures =
253 24 : CPLQuadTreeSearch(psTransform->hQuadTree, &aoi, &nFeatureCount);
254 24 : if (nFeatureCount != 0)
255 : {
256 24 : oPoint.setX(dfGeoX);
257 24 : oPoint.setY(dfGeoY);
258 24 : for (int iFeat = 0; iFeat < nFeatureCount; iFeat++)
259 : {
260 24 : size_t nIdx = reinterpret_cast<size_t>(pahFeatures[iFeat]);
261 24 : const bool bXRefAt180 = (nIdx >> BIT_IDX_RANGE_180) != 0;
262 : // Clear that bit.
263 24 : nIdx &= ~BIT_IDX_RANGE_180_SET;
264 :
265 24 : double x0 = 0, y0 = 0, x1 = 0, y1 = 0, x2 = 0, y2 = 0, x3 = 0,
266 24 : y3 = 0;
267 24 : GDALGeoLocQuadTreeGetFeatureCorners(psTransform, nIdx, x0, y0,
268 : x2, y2, x1, y1, x3, y3);
269 :
270 24 : if (psTransform->bGeographicSRSWithMinus180Plus180LongRange &&
271 24 : std::fabs(x0) > 170 && std::fabs(x1) > 170 &&
272 0 : std::fabs(x2) > 170 && std::fabs(x3) > 170 &&
273 0 : (std::fabs(x1 - x0) > 180 || std::fabs(x2 - x0) > 180 ||
274 0 : std::fabs(x3 - x0) > 180))
275 : {
276 0 : const double dfXRef = bXRefAt180 ? 180 : -180;
277 0 : x0 = ShiftGeoX(psTransform, dfXRef, x0);
278 0 : x1 = ShiftGeoX(psTransform, dfXRef, x1);
279 0 : x2 = ShiftGeoX(psTransform, dfXRef, x2);
280 0 : x3 = ShiftGeoX(psTransform, dfXRef, x3);
281 : }
282 :
283 24 : oRing.setPoint(0, x0, y0);
284 24 : oRing.setPoint(1, x2, y2);
285 24 : oRing.setPoint(2, x3, y3);
286 24 : oRing.setPoint(3, x1, y1);
287 24 : oRing.setPoint(4, x0, y0);
288 :
289 32 : if (oRing.isPointInRing(&oPoint) ||
290 8 : oRing.isPointOnRingBoundary(&oPoint))
291 : {
292 24 : const size_t nExtendedWidth =
293 48 : psTransform->nGeoLocXSize +
294 24 : (psTransform->bOriginIsTopLeftCorner ? 0 : 1);
295 24 : double dfX = static_cast<double>(nIdx % nExtendedWidth);
296 : // store the result as int, and then cast to double, to
297 : // avoid Coverity Scan warning about
298 : // UNINTENDED_INTEGER_DIVISION
299 24 : const size_t nY = nIdx / nExtendedWidth;
300 24 : double dfY = static_cast<double>(nY);
301 24 : if (!psTransform->bOriginIsTopLeftCorner)
302 : {
303 12 : dfX -= 1.0;
304 12 : dfY -= 1.0;
305 : }
306 24 : GDALInverseBilinearInterpolation(dfGeoX, dfGeoY, x0, y0, x1,
307 : y1, x2, y2, x3, y3, dfX,
308 : dfY);
309 :
310 24 : dfX = (dfX + dfGeorefConventionOffset) *
311 24 : psTransform->dfPIXEL_STEP +
312 24 : psTransform->dfPIXEL_OFFSET;
313 24 : dfY = (dfY + dfGeorefConventionOffset) *
314 24 : psTransform->dfLINE_STEP +
315 24 : psTransform->dfLINE_OFFSET;
316 :
317 24 : bDone = true;
318 24 : panSuccess[i] = TRUE;
319 24 : padfX[i] = dfX;
320 24 : padfY[i] = dfY;
321 24 : break;
322 : }
323 : }
324 : }
325 24 : CPLFree(pahFeatures);
326 :
327 24 : if (!bDone)
328 : {
329 0 : panSuccess[i] = FALSE;
330 0 : padfX[i] = HUGE_VAL;
331 0 : padfY[i] = HUGE_VAL;
332 : }
333 : }
334 24 : }
|
