/src/gdal/alg/gdalgeolocquadtree.cpp

Source (jump to first uncovered line)
/******************************************************************************
 *
 * Project:  GDAL
 * Purpose:  Implements Geolocation array based transformer, using a quadtree
 *           for inverse
 * Author:   Even Rouault, <even.rouault at spatialys.com>
 *
 ******************************************************************************
 * Copyright (c) 2022, Planet Labs
 *
 * SPDX-License-Identifier: MIT
 ****************************************************************************/

#include "gdalgeoloc.h"
#include "gdalgeolocquadtree.h"

#include "cpl_quad_tree.h"

#include "ogr_geometry.h"

#include <algorithm>
#include <cstddef>
#include <limits>

/************************************************************************/
/*               GDALGeoLocQuadTreeGetFeatureCorners()                  */
/************************************************************************/

static bool
GDALGeoLocQuadTreeGetFeatureCorners(const GDALGeoLocTransformInfo *psTransform,
                                    size_t nIdx, double &x0, double &y0,
                                    double &x1, double &y1, double &x2,
                                    double &y2, double &x3, double &y3)
{
    const size_t nExtendedWidth = psTransform->nGeoLocXSize +
                                  (psTransform->bOriginIsTopLeftCorner ? 0 : 1);
    int nX = static_cast<int>(nIdx % nExtendedWidth);
    int nY = static_cast<int>(nIdx / nExtendedWidth);

    if (!psTransform->bOriginIsTopLeftCorner)
    {
        nX--;
        nY--;
    }

    return GDALGeoLocExtractSquare(psTransform, static_cast<int>(nX),
                                   static_cast<int>(nY), x0, y0, x1, y1, x2, y2,
                                   x3, y3);
}

/************************************************************************/
/*               GDALGeoLocQuadTreeGetFeatureBounds()                   */
/************************************************************************/

constexpr size_t BIT_IDX_RANGE_180 = 8 * sizeof(size_t) - 1;
constexpr size_t BIT_IDX_RANGE_180_SET = static_cast<size_t>(1)
                                         << BIT_IDX_RANGE_180;

// Callback used by quadtree to retrieve the bounding box, in georeferenced
// space, of a cell of the geolocation array.
static void GDALGeoLocQuadTreeGetFeatureBounds(const void *hFeature,
                                               void *pUserData,
                                               CPLRectObj *pBounds)
{
    const GDALGeoLocTransformInfo *psTransform =
        static_cast<const GDALGeoLocTransformInfo *>(pUserData);
    size_t nIdx = reinterpret_cast<size_t>(hFeature);
    // Most significant bit set means that geometries crossing the antimeridian
    // should have their longitudes lower or greater than 180 deg.
    const bool bXRefAt180 = (nIdx >> BIT_IDX_RANGE_180) != 0;
    // Clear that bit.
    nIdx &= ~BIT_IDX_RANGE_180_SET;

    double x0 = 0, y0 = 0, x1 = 0, y1 = 0, x2 = 0, y2 = 0, x3 = 0, y3 = 0;
    GDALGeoLocQuadTreeGetFeatureCorners(psTransform, nIdx, x0, y0, x1, y1, x2,
                                        y2, x3, y3);

    if (psTransform->bGeographicSRSWithMinus180Plus180LongRange &&
        std::fabs(x0) > 170 && std::fabs(x1) > 170 && std::fabs(x2) > 170 &&
        std::fabs(x3) > 170 &&
        (std::fabs(x1 - x0) > 180 || std::fabs(x2 - x0) > 180 ||
         std::fabs(x3 - x0) > 180))
    {
        const double dfXRef = bXRefAt180 ? 180 : -180;
        x0 = ShiftGeoX(psTransform, dfXRef, x0);
        x1 = ShiftGeoX(psTransform, dfXRef, x1);
        x2 = ShiftGeoX(psTransform, dfXRef, x2);
        x3 = ShiftGeoX(psTransform, dfXRef, x3);
    }
    pBounds->minx = std::min(std::min(x0, x1), std::min(x2, x3));
    pBounds->miny = std::min(std::min(y0, y1), std::min(y2, y3));
    pBounds->maxx = std::max(std::max(x0, x1), std::max(x2, x3));
    pBounds->maxy = std::max(std::max(y0, y1), std::max(y2, y3));
}

/************************************************************************/
/*                      GDALGeoLocBuildQuadTree()                       */
/************************************************************************/

bool GDALGeoLocBuildQuadTree(GDALGeoLocTransformInfo *psTransform)
{
    // For the pixel-center convention, insert a "virtual" row and column
    // at top and left of the geoloc array.
    const int nExtraPixel = psTransform->bOriginIsTopLeftCorner ? 0 : 1;

    if (psTransform->nGeoLocXSize > INT_MAX - nExtraPixel ||
        psTransform->nGeoLocYSize > INT_MAX - nExtraPixel ||
        // The >> 1 shift is because we need to reserve the most-significant-bit
        // for the second 'version' of anti-meridian crossing quadrilaterals.
        // See below
        static_cast<size_t>(psTransform->nGeoLocXSize + nExtraPixel) >
            (std::numeric_limits<size_t>::max() >> 1) /
                static_cast<size_t>(psTransform->nGeoLocYSize + nExtraPixel))
    {
        CPLError(CE_Failure, CPLE_AppDefined, "Too big geolocation array");
        return false;
    }

    const int nExtendedWidth = psTransform->nGeoLocXSize + nExtraPixel;
    const int nExtendedHeight = psTransform->nGeoLocYSize + nExtraPixel;
    const size_t nExtendedXYCount =
        static_cast<size_t>(nExtendedWidth) * nExtendedHeight;

    CPLDebug("GEOLOC", "Start quadtree construction");

    CPLRectObj globalBounds;
    globalBounds.minx = psTransform->dfMinX;
    globalBounds.miny = psTransform->dfMinY;
    globalBounds.maxx = psTransform->dfMaxX;
    globalBounds.maxy = psTransform->dfMaxY;
    psTransform->hQuadTree = CPLQuadTreeCreateEx(
        &globalBounds, GDALGeoLocQuadTreeGetFeatureBounds, psTransform);

    CPLQuadTreeForceUseOfSubNodes(psTransform->hQuadTree);

    for (size_t i = 0; i < nExtendedXYCount; i++)
    {
        double x0, y0, x1, y1, x2, y2, x3, y3;
        if (!GDALGeoLocQuadTreeGetFeatureCorners(psTransform, i, x0, y0, x1, y1,
                                                 x2, y2, x3, y3))
        {
            continue;
        }

        // Skip too large geometries (typically at very high latitudes)
        // that would fill too many nodes in the quadtree
        if (psTransform->bGeographicSRSWithMinus180Plus180LongRange &&
            (std::fabs(x0) > 170 || std::fabs(x1) > 170 ||
             std::fabs(x2) > 170 || std::fabs(x3) > 170) &&
            (std::fabs(x1 - x0) > 180 || std::fabs(x2 - x0) > 180 ||
             std::fabs(x3 - x0) > 180) &&
            !(std::fabs(x0) > 170 && std::fabs(x1) > 170 &&
              std::fabs(x2) > 170 && std::fabs(x3) > 170))
        {
            continue;
        }

        CPLQuadTreeInsert(psTransform->hQuadTree,
                          reinterpret_cast<void *>(static_cast<uintptr_t>(i)));

        // For a geometry crossing the antimeridian, we've insert before
        // the "version" around -180 deg. Insert its corresponding version
        // around +180 deg.
        if (psTransform->bGeographicSRSWithMinus180Plus180LongRange &&
            std::fabs(x0) > 170 && std::fabs(x1) > 170 && std::fabs(x2) > 170 &&
            std::fabs(x3) > 170 &&
            (std::fabs(x1 - x0) > 180 || std::fabs(x2 - x0) > 180 ||
             std::fabs(x3 - x0) > 180))
        {
            CPLQuadTreeInsert(psTransform->hQuadTree,
                              reinterpret_cast<void *>(static_cast<uintptr_t>(
                                  i | BIT_IDX_RANGE_180_SET)));
        }
    }

    CPLDebug("GEOLOC", "End of quadtree construction");

#ifdef DEBUG_GEOLOC
    int nFeatureCount = 0;
    int nNodeCount = 0;
    int nMaxDepth = 0;
    int nMaxBucketCapacity = 0;
    CPLQuadTreeGetStats(psTransform->hQuadTree, &nFeatureCount, &nNodeCount,
                        &nMaxDepth, &nMaxBucketCapacity);
    CPLDebug("GEOLOC", "Quadtree stats:");
    CPLDebug("GEOLOC", "  nFeatureCount = %d", nFeatureCount);
    CPLDebug("GEOLOC", "  nNodeCount = %d", nNodeCount);
    CPLDebug("GEOLOC", "  nMaxDepth = %d", nMaxDepth);
    CPLDebug("GEOLOC", "  nMaxBucketCapacity = %d", nMaxBucketCapacity);
#endif

    return true;
}

/************************************************************************/
/*                  GDALGeoLocInverseTransformQuadtree()                */
/************************************************************************/

void GDALGeoLocInverseTransformQuadtree(
    const GDALGeoLocTransformInfo *psTransform, int nPointCount, double *padfX,
    double *padfY, int *panSuccess)
{
    // Keep those objects in this outer scope, so they are re-used, to
    // save memory allocations.
    OGRPoint oPoint;
    OGRLinearRing oRing;
    oRing.setNumPoints(5);

    const double dfGeorefConventionOffset =
        psTransform->bOriginIsTopLeftCorner ? 0 : 0.5;

    for (int i = 0; i < nPointCount; i++)
    {
        if (padfX[i] == HUGE_VAL || padfY[i] == HUGE_VAL)
        {
            panSuccess[i] = FALSE;
            continue;
        }

        if (psTransform->bSwapXY)
        {
            std::swap(padfX[i], padfY[i]);
        }

        const double dfGeoX = padfX[i];
        const double dfGeoY = padfY[i];

        bool bDone = false;

        CPLRectObj aoi;
        aoi.minx = dfGeoX;
        aoi.maxx = dfGeoX;
        aoi.miny = dfGeoY;
        aoi.maxy = dfGeoY;
        int nFeatureCount = 0;
        void **pahFeatures =
            CPLQuadTreeSearch(psTransform->hQuadTree, &aoi, &nFeatureCount);
        if (nFeatureCount != 0)
        {
            oPoint.setX(dfGeoX);
            oPoint.setY(dfGeoY);
            for (int iFeat = 0; iFeat < nFeatureCount; iFeat++)
            {
                size_t nIdx = reinterpret_cast<size_t>(pahFeatures[iFeat]);
                const bool bXRefAt180 = (nIdx >> BIT_IDX_RANGE_180) != 0;
                // Clear that bit.
                nIdx &= ~BIT_IDX_RANGE_180_SET;

                double x0 = 0, y0 = 0, x1 = 0, y1 = 0, x2 = 0, y2 = 0, x3 = 0,
                       y3 = 0;
                GDALGeoLocQuadTreeGetFeatureCorners(psTransform, nIdx, x0, y0,
                                                    x2, y2, x1, y1, x3, y3);

                if (psTransform->bGeographicSRSWithMinus180Plus180LongRange &&
                    std::fabs(x0) > 170 && std::fabs(x1) > 170 &&
                    std::fabs(x2) > 170 && std::fabs(x3) > 170 &&
                    (std::fabs(x1 - x0) > 180 || std::fabs(x2 - x0) > 180 ||
                     std::fabs(x3 - x0) > 180))
                {
                    const double dfXRef = bXRefAt180 ? 180 : -180;
                    x0 = ShiftGeoX(psTransform, dfXRef, x0);
                    x1 = ShiftGeoX(psTransform, dfXRef, x1);
                    x2 = ShiftGeoX(psTransform, dfXRef, x2);
                    x3 = ShiftGeoX(psTransform, dfXRef, x3);
                }

                oRing.setPoint(0, x0, y0);
                oRing.setPoint(1, x2, y2);
                oRing.setPoint(2, x3, y3);
                oRing.setPoint(3, x1, y1);
                oRing.setPoint(4, x0, y0);

                if (oRing.isPointInRing(&oPoint) ||
                    oRing.isPointOnRingBoundary(&oPoint))
                {
                    const size_t nExtendedWidth =
                        psTransform->nGeoLocXSize +
                        (psTransform->bOriginIsTopLeftCorner ? 0 : 1);
                    double dfX = static_cast<double>(nIdx % nExtendedWidth);
                    // store the result as int, and then cast to double, to
                    // avoid Coverity Scan warning about
                    // UNINTENDED_INTEGER_DIVISION
                    const size_t nY = nIdx / nExtendedWidth;
                    double dfY = static_cast<double>(nY);
                    if (!psTransform->bOriginIsTopLeftCorner)
                    {
                        dfX -= 1.0;
                        dfY -= 1.0;
                    }
                    GDALInverseBilinearInterpolation(dfGeoX, dfGeoY, x0, y0, x1,
                                                     y1, x2, y2, x3, y3, dfX,
                                                     dfY);

                    dfX = (dfX + dfGeorefConventionOffset) *
                              psTransform->dfPIXEL_STEP +
                          psTransform->dfPIXEL_OFFSET;
                    dfY = (dfY + dfGeorefConventionOffset) *
                              psTransform->dfLINE_STEP +
                          psTransform->dfLINE_OFFSET;

                    bDone = true;
                    panSuccess[i] = TRUE;
                    padfX[i] = dfX;
                    padfY[i] = dfY;
                    break;
                }
            }
        }
        CPLFree(pahFeatures);

        if (!bDone)
        {
            panSuccess[i] = FALSE;
            padfX[i] = HUGE_VAL;
            padfY[i] = HUGE_VAL;
        }
    }
}

Coverage Report

Created: 2025-06-13 06:29

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	}