Coverage Report

Created: 2025-06-09 07:11

/src/gdal/alg/gdalgeolocquadtree.cpp
Line
Count
Source (jump to first uncovered line)
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
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
0
{
35
0
    const size_t nExtendedWidth = psTransform->nGeoLocXSize +
36
0
                                  (psTransform->bOriginIsTopLeftCorner ? 0 : 1);
37
0
    int nX = static_cast<int>(nIdx % nExtendedWidth);
38
0
    int nY = static_cast<int>(nIdx / nExtendedWidth);
39
40
0
    if (!psTransform->bOriginIsTopLeftCorner)
41
0
    {
42
0
        nX--;
43
0
        nY--;
44
0
    }
45
46
0
    return GDALGeoLocExtractSquare(psTransform, static_cast<int>(nX),
47
0
                                   static_cast<int>(nY), x0, y0, x1, y1, x2, y2,
48
0
                                   x3, y3);
49
0
}
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
static void GDALGeoLocQuadTreeGetFeatureBounds(const void *hFeature,
62
                                               void *pUserData,
63
                                               CPLRectObj *pBounds)
64
0
{
65
0
    const GDALGeoLocTransformInfo *psTransform =
66
0
        static_cast<const GDALGeoLocTransformInfo *>(pUserData);
67
0
    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
0
    const bool bXRefAt180 = (nIdx >> BIT_IDX_RANGE_180) != 0;
71
    // Clear that bit.
72
0
    nIdx &= ~BIT_IDX_RANGE_180_SET;
73
74
0
    double x0 = 0, y0 = 0, x1 = 0, y1 = 0, x2 = 0, y2 = 0, x3 = 0, y3 = 0;
75
0
    GDALGeoLocQuadTreeGetFeatureCorners(psTransform, nIdx, x0, y0, x1, y1, x2,
76
0
                                        y2, x3, y3);
77
78
0
    if (psTransform->bGeographicSRSWithMinus180Plus180LongRange &&
79
0
        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
0
    {
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
0
    }
90
0
    pBounds->minx = std::min(std::min(x0, x1), std::min(x2, x3));
91
0
    pBounds->miny = std::min(std::min(y0, y1), std::min(y2, y3));
92
0
    pBounds->maxx = std::max(std::max(x0, x1), std::max(x2, x3));
93
0
    pBounds->maxy = std::max(std::max(y0, y1), std::max(y2, y3));
94
0
}
95
96
/************************************************************************/
97
/*                      GDALGeoLocBuildQuadTree()                       */
98
/************************************************************************/
99
100
bool GDALGeoLocBuildQuadTree(GDALGeoLocTransformInfo *psTransform)
101
0
{
102
    // For the pixel-center convention, insert a "virtual" row and column
103
    // at top and left of the geoloc array.
104
0
    const int nExtraPixel = psTransform->bOriginIsTopLeftCorner ? 0 : 1;
105
106
0
    if (psTransform->nGeoLocXSize > INT_MAX - nExtraPixel ||
107
0
        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
0
        static_cast<size_t>(psTransform->nGeoLocXSize + nExtraPixel) >
112
0
            (std::numeric_limits<size_t>::max() >> 1) /
113
0
                static_cast<size_t>(psTransform->nGeoLocYSize + nExtraPixel))
114
0
    {
115
0
        CPLError(CE_Failure, CPLE_AppDefined, "Too big geolocation array");
116
0
        return false;
117
0
    }
118
119
0
    const int nExtendedWidth = psTransform->nGeoLocXSize + nExtraPixel;
120
0
    const int nExtendedHeight = psTransform->nGeoLocYSize + nExtraPixel;
121
0
    const size_t nExtendedXYCount =
122
0
        static_cast<size_t>(nExtendedWidth) * nExtendedHeight;
123
124
0
    CPLDebug("GEOLOC", "Start quadtree construction");
125
126
0
    CPLRectObj globalBounds;
127
0
    globalBounds.minx = psTransform->dfMinX;
128
0
    globalBounds.miny = psTransform->dfMinY;
129
0
    globalBounds.maxx = psTransform->dfMaxX;
130
0
    globalBounds.maxy = psTransform->dfMaxY;
131
0
    psTransform->hQuadTree = CPLQuadTreeCreateEx(
132
0
        &globalBounds, GDALGeoLocQuadTreeGetFeatureBounds, psTransform);
133
134
0
    CPLQuadTreeForceUseOfSubNodes(psTransform->hQuadTree);
135
136
0
    for (size_t i = 0; i < nExtendedXYCount; i++)
137
0
    {
138
0
        double x0, y0, x1, y1, x2, y2, x3, y3;
139
0
        if (!GDALGeoLocQuadTreeGetFeatureCorners(psTransform, i, x0, y0, x1, y1,
140
0
                                                 x2, y2, x3, y3))
141
0
        {
142
0
            continue;
143
0
        }
144
145
        // Skip too large geometries (typically at very high latitudes)
146
        // that would fill too many nodes in the quadtree
147
0
        if (psTransform->bGeographicSRSWithMinus180Plus180LongRange &&
148
0
            (std::fabs(x0) > 170 || std::fabs(x1) > 170 ||
149
0
             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
0
        {
155
0
            continue;
156
0
        }
157
158
0
        CPLQuadTreeInsert(psTransform->hQuadTree,
159
0
                          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
0
        if (psTransform->bGeographicSRSWithMinus180Plus180LongRange &&
165
0
            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
0
        {
170
0
            CPLQuadTreeInsert(psTransform->hQuadTree,
171
0
                              reinterpret_cast<void *>(static_cast<uintptr_t>(
172
0
                                  i | BIT_IDX_RANGE_180_SET)));
173
0
        }
174
0
    }
175
176
0
    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
0
    return true;
193
0
}
194
195
/************************************************************************/
196
/*                  GDALGeoLocInverseTransformQuadtree()                */
197
/************************************************************************/
198
199
void GDALGeoLocInverseTransformQuadtree(
200
    const GDALGeoLocTransformInfo *psTransform, int nPointCount, double *padfX,
201
    double *padfY, int *panSuccess)
202
0
{
203
    // Keep those objects in this outer scope, so they are re-used, to
204
    // save memory allocations.
205
0
    OGRPoint oPoint;
206
0
    OGRLinearRing oRing;
207
0
    oRing.setNumPoints(5);
208
209
0
    const double dfGeorefConventionOffset =
210
0
        psTransform->bOriginIsTopLeftCorner ? 0 : 0.5;
211
212
0
    for (int i = 0; i < nPointCount; i++)
213
0
    {
214
0
        if (padfX[i] == HUGE_VAL || padfY[i] == HUGE_VAL)
215
0
        {
216
0
            panSuccess[i] = FALSE;
217
0
            continue;
218
0
        }
219
220
0
        if (psTransform->bSwapXY)
221
0
        {
222
0
            std::swap(padfX[i], padfY[i]);
223
0
        }
224
225
0
        const double dfGeoX = padfX[i];
226
0
        const double dfGeoY = padfY[i];
227
228
0
        bool bDone = false;
229
230
0
        CPLRectObj aoi;
231
0
        aoi.minx = dfGeoX;
232
0
        aoi.maxx = dfGeoX;
233
0
        aoi.miny = dfGeoY;
234
0
        aoi.maxy = dfGeoY;
235
0
        int nFeatureCount = 0;
236
0
        void **pahFeatures =
237
0
            CPLQuadTreeSearch(psTransform->hQuadTree, &aoi, &nFeatureCount);
238
0
        if (nFeatureCount != 0)
239
0
        {
240
0
            oPoint.setX(dfGeoX);
241
0
            oPoint.setY(dfGeoY);
242
0
            for (int iFeat = 0; iFeat < nFeatureCount; iFeat++)
243
0
            {
244
0
                size_t nIdx = reinterpret_cast<size_t>(pahFeatures[iFeat]);
245
0
                const bool bXRefAt180 = (nIdx >> BIT_IDX_RANGE_180) != 0;
246
                // Clear that bit.
247
0
                nIdx &= ~BIT_IDX_RANGE_180_SET;
248
249
0
                double x0 = 0, y0 = 0, x1 = 0, y1 = 0, x2 = 0, y2 = 0, x3 = 0,
250
0
                       y3 = 0;
251
0
                GDALGeoLocQuadTreeGetFeatureCorners(psTransform, nIdx, x0, y0,
252
0
                                                    x2, y2, x1, y1, x3, y3);
253
254
0
                if (psTransform->bGeographicSRSWithMinus180Plus180LongRange &&
255
0
                    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
0
                {
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
0
                }
266
267
0
                oRing.setPoint(0, x0, y0);
268
0
                oRing.setPoint(1, x2, y2);
269
0
                oRing.setPoint(2, x3, y3);
270
0
                oRing.setPoint(3, x1, y1);
271
0
                oRing.setPoint(4, x0, y0);
272
273
0
                if (oRing.isPointInRing(&oPoint) ||
274
0
                    oRing.isPointOnRingBoundary(&oPoint))
275
0
                {
276
0
                    const size_t nExtendedWidth =
277
0
                        psTransform->nGeoLocXSize +
278
0
                        (psTransform->bOriginIsTopLeftCorner ? 0 : 1);
279
0
                    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
0
                    const size_t nY = nIdx / nExtendedWidth;
284
0
                    double dfY = static_cast<double>(nY);
285
0
                    if (!psTransform->bOriginIsTopLeftCorner)
286
0
                    {
287
0
                        dfX -= 1.0;
288
0
                        dfY -= 1.0;
289
0
                    }
290
0
                    GDALInverseBilinearInterpolation(dfGeoX, dfGeoY, x0, y0, x1,
291
0
                                                     y1, x2, y2, x3, y3, dfX,
292
0
                                                     dfY);
293
294
0
                    dfX = (dfX + dfGeorefConventionOffset) *
295
0
                              psTransform->dfPIXEL_STEP +
296
0
                          psTransform->dfPIXEL_OFFSET;
297
0
                    dfY = (dfY + dfGeorefConventionOffset) *
298
0
                              psTransform->dfLINE_STEP +
299
0
                          psTransform->dfLINE_OFFSET;
300
301
0
                    bDone = true;
302
0
                    panSuccess[i] = TRUE;
303
0
                    padfX[i] = dfX;
304
0
                    padfY[i] = dfY;
305
0
                    break;
306
0
                }
307
0
            }
308
0
        }
309
0
        CPLFree(pahFeatures);
310
311
0
        if (!bDone)
312
0
        {
313
0
            panSuccess[i] = FALSE;
314
0
            padfX[i] = HUGE_VAL;
315
0
            padfY[i] = HUGE_VAL;
316
0
        }
317
0
    }
318
0
}