/******************************************************************************

 *

 * Project:  GDAL Gridding API.

 * Purpose:  Implementation of GDAL scattered data gridder.

 * Author:   Andrey Kiselev, dron@ak4719.spb.edu

 *

 ******************************************************************************

 * Copyright (c) 2007, Andrey Kiselev <dron@ak4719.spb.edu>

 * Copyright (c) 2009-2013, Even Rouault <even dot rouault at spatialys.com>

 *

 * SPDX-License-Identifier: MIT

 ****************************************************************************/

#include "cpl_port.h"

#include "gdalgrid.h"

#include "gdalgrid_priv.h"

#include <cfloat>

#include <climits>

#include <cmath>

#include <cstdlib>

#include <cstring>

#include <limits>

#include <map>

#include <utility>

#include <algorithm>

#include "cpl_conv.h"

#include "cpl_cpu_features.h"

#include "cpl_error.h"

#include "cpl_multiproc.h"

#include "cpl_progress.h"

#include "cpl_quad_tree.h"

#include "cpl_string.h"

#include "cpl_vsi.h"

#include "cpl_worker_thread_pool.h"

#include "gdal.h"

#include "gdal_thread_pool.h"

constexpr double TO_RADIANS = M_PI / 180.0;

/************************************************************************/

/*                       GDALGridGetPointBounds()                       */

/************************************************************************/

static void GDALGridGetPointBounds(const void *hFeature, CPLRectObj *pBounds)

{

    const GDALGridPoint *psPoint = static_cast<const GDALGridPoint *>(hFeature);

    GDALGridXYArrays *psXYArrays = psPoint->psXYArrays;

    const int i = psPoint->i;

    const double dfX = psXYArrays->padfX[i];

    const double dfY = psXYArrays->padfY[i];

    pBounds->minx = dfX;

    pBounds->miny = dfY;

    pBounds->maxx = dfX;

    pBounds->maxy = dfY;

}

/************************************************************************/

/*                  GDALGridInverseDistanceToAPower()                   */

/************************************************************************/

/**

 * Inverse distance to a power.

 *

 * The Inverse Distance to a Power gridding method is a weighted average

 * interpolator. You should supply the input arrays with the scattered data

 * values including coordinates of every data point and output grid geometry.

 * The function will compute interpolated value for the given position in

 * output grid.

 *

 * For every grid node the resulting value \f$Z\f$ will be calculated using

 * formula:

 *

 * \f[

 *      Z=\frac{\sum_{i=1}^n{\frac{Z_i}{r_i^p}}}{\sum_{i=1}^n{\frac{1}{r_i^p}}}

 * \f]

 *

 *  where

 *  <ul>

 *      <li> \f$Z_i\f$ is a known value at point \f$i\f$,

 *      <li> \f$r_i\f$ is an Euclidean distance from the grid node

 *           to point \f$i\f$,

 *      <li> \f$p\f$ is a weighting power,

 *      <li> \f$n\f$ is a total number of points in search ellipse.

 *  </ul>

 *

 *  In this method the weighting factor \f$w\f$ is

 *

 *  \f[

 *      w=\frac{1}{r^p}

 *  \f]

 *

 * @param poOptionsIn Algorithm parameters. This should point to

 * GDALGridInverseDistanceToAPowerOptions object.

 * @param nPoints Number of elements in input arrays.

 * @param padfX Input array of X coordinates.

 * @param padfY Input array of Y coordinates.

 * @param padfZ Input array of Z values.

 * @param dfXPoint X coordinate of the point to compute.

 * @param dfYPoint Y coordinate of the point to compute.

 * @param pdfValue Pointer to variable where the computed grid node value

 * will be returned.

 * @param hExtraParamsIn extra parameters (unused)

 *

 * @return CE_None on success or CE_Failure if something goes wrong.

 */

CPLErr GDALGridInverseDistanceToAPower(const void *poOptionsIn, GUInt32 nPoints,

                                       const double *padfX, const double *padfY,

                                       const double *padfZ, double dfXPoint,

                                       double dfYPoint, double *pdfValue,

                                       CPL_UNUSED void *hExtraParamsIn)

{

    // TODO: For optimization purposes pre-computed parameters should be moved

    // out of this routine to the calling function.

    const GDALGridInverseDistanceToAPowerOptions *const poOptions =

        static_cast<const GDALGridInverseDistanceToAPowerOptions *>(

            poOptionsIn);

    // Pre-compute search ellipse parameters.

    const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;

    const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;

    const double dfR12Square = dfRadius1Square * dfRadius2Square;

    // Compute coefficients for coordinate system rotation.

    const double dfAngle = TO_RADIANS * poOptions->dfAngle;

    const bool bRotated = dfAngle != 0.0;

    const double dfCoeff1 = bRotated ? cos(dfAngle) : 0.0;

    const double dfCoeff2 = bRotated ? sin(dfAngle) : 0.0;

    const double dfPowerDiv2 = poOptions->dfPower / 2;

    const double dfSmoothing = poOptions->dfSmoothing;

    const GUInt32 nMaxPoints = poOptions->nMaxPoints;

    double dfNominator = 0.0;

    double dfDenominator = 0.0;

    GUInt32 n = 0;

    for (GUInt32 i = 0; i < nPoints; i++)

    {

        double dfRX = padfX[i] - dfXPoint;

        double dfRY = padfY[i] - dfYPoint;

        const double dfR2 =

            dfRX * dfRX + dfRY * dfRY + dfSmoothing * dfSmoothing;

        if (bRotated)

        {

            const double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;

            const double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;

            dfRX = dfRXRotated;

            dfRY = dfRYRotated;

        }

        // Is this point located inside the search ellipse?

        if (dfRadius2Square * dfRX * dfRX + dfRadius1Square * dfRY * dfRY <=

            dfR12Square)

        {

            // If the test point is close to the grid node, use the point

            // value directly as a node value to avoid singularity.

            if (dfR2 < 0.0000000000001)

            {

                *pdfValue = padfZ[i];

                return CE_None;

            }

            const double dfW = pow(dfR2, dfPowerDiv2);

            const double dfInvW = 1.0 / dfW;

            dfNominator += dfInvW * padfZ[i];

            dfDenominator += dfInvW;

            n++;

            if (nMaxPoints > 0 && n > nMaxPoints)

                break;

        }

    }

    if (n < poOptions->nMinPoints || dfDenominator == 0.0)

    {

        *pdfValue = poOptions->dfNoDataValue;

    }

    else

    {

        *pdfValue = dfNominator / dfDenominator;

    }

    return CE_None;

}

/************************************************************************/

/*           GDALGridInverseDistanceToAPowerNearestNeighbor()           */

/************************************************************************/

/**

 * Inverse distance to a power with nearest neighbor search, ideal when

 * max_points used.

 *

 * The Inverse Distance to a Power gridding method is a weighted average

 * interpolator. You should supply the input arrays with the scattered data

 * values including coordinates of every data point and output grid geometry.

 * The function will compute interpolated value for the given position in

 * output grid.

 *

 * For every grid node the resulting value \f$Z\f$ will be calculated using

 * formula for nearest matches:

 *

 * \f[

 *      Z=\frac{\sum_{i=1}^n{\frac{Z_i}{r_i^p}}}{\sum_{i=1}^n{\frac{1}{r_i^p}}}

 * \f]

 *

 *  where

 *  <ul>

 *      <li> \f$Z_i\f$ is a known value at point \f$i\f$,

 *      <li> \f$r_i\f$ is an Euclidean distance from the grid node

 *           to point \f$i\f$ (with an optional smoothing parameter \f$s\f$),

 *      <li> \f$p\f$ is a weighting power,

 *      <li> \f$n\f$ is a total number of points in search ellipse.

 *  </ul>

 *

 *  In this method the weighting factor \f$w\f$ is

 *

 *  \f[

 *      w=\frac{1}{r^p}

 *  \f]

 *

 * @param poOptionsIn Algorithm parameters. This should point to

 * GDALGridInverseDistanceToAPowerNearestNeighborOptions object.

 * @param nPoints Number of elements in input arrays.

 * @param padfX Input array of X coordinates.

 * @param padfY Input array of Y coordinates.

 * @param padfZ Input array of Z values.

 * @param dfXPoint X coordinate of the point to compute.

 * @param dfYPoint Y coordinate of the point to compute.

 * @param pdfValue Pointer to variable where the computed grid node value

 * will be returned.

 * @param hExtraParamsIn extra parameters.

 *

 * @return CE_None on success or CE_Failure if something goes wrong.

 */

CPLErr GDALGridInverseDistanceToAPowerNearestNeighbor(

    const void *poOptionsIn, GUInt32 nPoints, const double *padfX,

    const double *padfY, const double *padfZ, double dfXPoint, double dfYPoint,

    double *pdfValue, void *hExtraParamsIn)

{

    CPL_IGNORE_RET_VAL(nPoints);

    const GDALGridInverseDistanceToAPowerNearestNeighborOptions *const

        poOptions = static_cast<

            const GDALGridInverseDistanceToAPowerNearestNeighborOptions *>(

            poOptionsIn);

    const double dfRadius = poOptions->dfRadius;

    const double dfSmoothing = poOptions->dfSmoothing;

    const double dfSmoothing2 = dfSmoothing * dfSmoothing;

    const GUInt32 nMaxPoints = poOptions->nMaxPoints;

    GDALGridExtraParameters *psExtraParams =

        static_cast<GDALGridExtraParameters *>(hExtraParamsIn);

    const CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;

    CPLAssert(phQuadTree);

    const double dfRPower2 = psExtraParams->dfRadiusPower2PreComp;

    const double dfPowerDiv2 = psExtraParams->dfPowerDiv2PreComp;

    std::multimap<double, double> oMapDistanceToZValues;

    const double dfSearchRadius = dfRadius;

    CPLRectObj sAoi;

    sAoi.minx = dfXPoint - dfSearchRadius;

    sAoi.miny = dfYPoint - dfSearchRadius;

    sAoi.maxx = dfXPoint + dfSearchRadius;

    sAoi.maxy = dfYPoint + dfSearchRadius;

    int nFeatureCount = 0;

    GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(

        CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));

    if (nFeatureCount != 0)

    {

        for (int k = 0; k < nFeatureCount; k++)

        {

            const int i = papsPoints[k]->i;

            const double dfRX = padfX[i] - dfXPoint;

            const double dfRY = padfY[i] - dfYPoint;

            const double dfR2 = dfRX * dfRX + dfRY * dfRY;

            // real distance + smoothing

            const double dfRsmoothed2 = dfR2 + dfSmoothing2;

            if (dfRsmoothed2 < 0.0000000000001)

            {

                *pdfValue = padfZ[i];

                CPLFree(papsPoints);

                return CE_None;

            }

            // is point within real distance?

            if (dfR2 <= dfRPower2)

            {

                oMapDistanceToZValues.insert(

                    std::make_pair(dfRsmoothed2, padfZ[i]));

            }

        }

    }

    CPLFree(papsPoints);

    double dfNominator = 0.0;

    double dfDenominator = 0.0;

    GUInt32 n = 0;

    // Examine all "neighbors" within the radius (sorted by distance via the

    // multimap), and use the closest n points based on distance until the max

    // is reached.

    for (std::multimap<double, double>::iterator oMapDistanceToZValuesIter =

             oMapDistanceToZValues.begin();

         oMapDistanceToZValuesIter != oMapDistanceToZValues.end();

         ++oMapDistanceToZValuesIter)

    {

        const double dfR2 = oMapDistanceToZValuesIter->first;

        const double dfZ = oMapDistanceToZValuesIter->second;

        const double dfW = pow(dfR2, dfPowerDiv2);

        const double dfInvW = 1.0 / dfW;

        dfNominator += dfInvW * dfZ;

        dfDenominator += dfInvW;

        n++;

        if (nMaxPoints > 0 && n >= nMaxPoints)

        {

            break;

        }

    }

    if (n < poOptions->nMinPoints || dfDenominator == 0.0)

    {

        *pdfValue = poOptions->dfNoDataValue;

    }

    else

    {

        *pdfValue = dfNominator / dfDenominator;

    }

    return CE_None;

}

/************************************************************************/

/*     GDALGridInverseDistanceToAPowerNearestNeighborPerQuadrant()      */

/************************************************************************/

/**

 * Inverse distance to a power with nearest neighbor search, with a per-quadrant

 * search logic.

 */

static CPLErr GDALGridInverseDistanceToAPowerNearestNeighborPerQuadrant(

    const void *poOptionsIn, GUInt32 /*nPoints*/, const double *padfX,

    const double *padfY, const double *padfZ, double dfXPoint, double dfYPoint,

    double *pdfValue, void *hExtraParamsIn)

{

    const GDALGridInverseDistanceToAPowerNearestNeighborOptions *const

        poOptions = static_cast<

            const GDALGridInverseDistanceToAPowerNearestNeighborOptions *>(

            poOptionsIn);

    const double dfRadius = poOptions->dfRadius;

    const double dfSmoothing = poOptions->dfSmoothing;

    const double dfSmoothing2 = dfSmoothing * dfSmoothing;

    const GUInt32 nMaxPoints = poOptions->nMaxPoints;

    const GUInt32 nMinPointsPerQuadrant = poOptions->nMinPointsPerQuadrant;

    const GUInt32 nMaxPointsPerQuadrant = poOptions->nMaxPointsPerQuadrant;

    GDALGridExtraParameters *psExtraParams =

        static_cast<GDALGridExtraParameters *>(hExtraParamsIn);

    const CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;

    CPLAssert(phQuadTree);

    const double dfRPower2 = psExtraParams->dfRadiusPower2PreComp;

    const double dfPowerDiv2 = psExtraParams->dfPowerDiv2PreComp;

    std::multimap<double, double> oMapDistanceToZValuesPerQuadrant[4];

    const double dfSearchRadius = dfRadius;

    CPLRectObj sAoi;

    sAoi.minx = dfXPoint - dfSearchRadius;

    sAoi.miny = dfYPoint - dfSearchRadius;

    sAoi.maxx = dfXPoint + dfSearchRadius;

    sAoi.maxy = dfYPoint + dfSearchRadius;

    int nFeatureCount = 0;

    GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(

        CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));

    if (nFeatureCount != 0)

    {

        for (int k = 0; k < nFeatureCount; k++)

        {

            const int i = papsPoints[k]->i;

            const double dfRX = padfX[i] - dfXPoint;

            const double dfRY = padfY[i] - dfYPoint;

            const double dfR2 = dfRX * dfRX + dfRY * dfRY;

            // real distance + smoothing

            const double dfRsmoothed2 = dfR2 + dfSmoothing2;

            if (dfRsmoothed2 < 0.0000000000001)

            {

                *pdfValue = padfZ[i];

                CPLFree(papsPoints);

                return CE_None;

            }

            // is point within real distance?

            if (dfR2 <= dfRPower2)

            {

                const int iQuadrant =

                    ((dfRX >= 0) ? 1 : 0) | (((dfRY >= 0) ? 1 : 0) << 1);

                oMapDistanceToZValuesPerQuadrant[iQuadrant].insert(

                    std::make_pair(dfRsmoothed2, padfZ[i]));

            }

        }

    }

    CPLFree(papsPoints);

    std::multimap<double, double>::iterator aoIter[] = {

        oMapDistanceToZValuesPerQuadrant[0].begin(),

        oMapDistanceToZValuesPerQuadrant[1].begin(),

        oMapDistanceToZValuesPerQuadrant[2].begin(),

        oMapDistanceToZValuesPerQuadrant[3].begin(),

    };

    constexpr int ALL_QUADRANT_FLAGS = 1 + 2 + 4 + 8;

    // Examine all "neighbors" within the radius (sorted by distance via the

    // multimap), and use the closest n points based on distance until the max

    // is reached.

    // Do that by fetching the nearest point in quadrant 0, then the nearest

    // point in quadrant 1, 2 and 3, and starting again with the next nearest

    // point in quarant 0, etc.

    int nQuadrantIterFinishedFlag = 0;

    GUInt32 anPerQuadrant[4] = {0};

    double dfNominator = 0.0;

    double dfDenominator = 0.0;

    GUInt32 n = 0;

    for (int iQuadrant = 0; /* true */; iQuadrant = (iQuadrant + 1) % 4)

    {

        if (aoIter[iQuadrant] ==

                oMapDistanceToZValuesPerQuadrant[iQuadrant].end() ||

            (nMaxPointsPerQuadrant > 0 &&

             anPerQuadrant[iQuadrant] >= nMaxPointsPerQuadrant))

        {

            nQuadrantIterFinishedFlag |= 1 << iQuadrant;

            if (nQuadrantIterFinishedFlag == ALL_QUADRANT_FLAGS)

                break;

            continue;

        }

        const double dfR2 = aoIter[iQuadrant]->first;

        const double dfZ = aoIter[iQuadrant]->second;

        ++aoIter[iQuadrant];

        const double dfW = pow(dfR2, dfPowerDiv2);

        const double dfInvW = 1.0 / dfW;

        dfNominator += dfInvW * dfZ;

        dfDenominator += dfInvW;

        n++;

        anPerQuadrant[iQuadrant]++;

        if (nMaxPoints > 0 && n >= nMaxPoints)

        {

            break;

        }

    }

    if (nMinPointsPerQuadrant > 0 &&

        (anPerQuadrant[0] < nMinPointsPerQuadrant ||

         anPerQuadrant[1] < nMinPointsPerQuadrant ||

         anPerQuadrant[2] < nMinPointsPerQuadrant ||

         anPerQuadrant[3] < nMinPointsPerQuadrant))

    {

        *pdfValue = poOptions->dfNoDataValue;

    }

    else if (n < poOptions->nMinPoints || dfDenominator == 0.0)

    {

        *pdfValue = poOptions->dfNoDataValue;

    }

    else

    {

        *pdfValue = dfNominator / dfDenominator;

    }

    return CE_None;

}

/************************************************************************/

/*              GDALGridInverseDistanceToAPowerNoSearch()               */

/************************************************************************/

/**

 * Inverse distance to a power for whole data set.

 *

 * This is somewhat optimized version of the Inverse Distance to a Power

 * method. It is used when the search ellips is not set. The algorithm and

 * parameters are the same as in GDALGridInverseDistanceToAPower(), but this

 * implementation works faster, because of no search.

 *

 * @see GDALGridInverseDistanceToAPower()

 */

CPLErr GDALGridInverseDistanceToAPowerNoSearch(

    const void *poOptionsIn, GUInt32 nPoints, const double *padfX,

    const double *padfY, const double *padfZ, double dfXPoint, double dfYPoint,

    double *pdfValue, void * /* hExtraParamsIn */)

{

    const GDALGridInverseDistanceToAPowerOptions *const poOptions =

        static_cast<const GDALGridInverseDistanceToAPowerOptions *>(

            poOptionsIn);

    const double dfPowerDiv2 = poOptions->dfPower / 2.0;

    const double dfSmoothing = poOptions->dfSmoothing;

    const double dfSmoothing2 = dfSmoothing * dfSmoothing;

    double dfNominator = 0.0;

    double dfDenominator = 0.0;

    const bool bPower2 = dfPowerDiv2 == 1.0;

    GUInt32 i = 0;  // Used after if.

    if (bPower2)

    {

        if (dfSmoothing2 > 0)

        {

            for (i = 0; i < nPoints; i++)

            {

                const double dfRX = padfX[i] - dfXPoint;

                const double dfRY = padfY[i] - dfYPoint;

                const double dfR2 = dfRX * dfRX + dfRY * dfRY + dfSmoothing2;

                const double dfInvR2 = 1.0 / dfR2;

                dfNominator += dfInvR2 * padfZ[i];

                dfDenominator += dfInvR2;

            }

        }

        else

        {

            for (i = 0; i < nPoints; i++)

            {

                const double dfRX = padfX[i] - dfXPoint;

                const double dfRY = padfY[i] - dfYPoint;

                const double dfR2 = dfRX * dfRX + dfRY * dfRY;

                // If the test point is close to the grid node, use the point

                // value directly as a node value to avoid singularity.

                if (dfR2 < 0.0000000000001)

                {

                    break;

                }

                const double dfInvR2 = 1.0 / dfR2;

                dfNominator += dfInvR2 * padfZ[i];

                dfDenominator += dfInvR2;

            }

        }

    }

    else

    {

        for (i = 0; i < nPoints; i++)

        {

            const double dfRX = padfX[i] - dfXPoint;

            const double dfRY = padfY[i] - dfYPoint;

            const double dfR2 = dfRX * dfRX + dfRY * dfRY + dfSmoothing2;

            // If the test point is close to the grid node, use the point

            // value directly as a node value to avoid singularity.

            if (dfR2 < 0.0000000000001)

            {

                break;

            }

            const double dfW = pow(dfR2, dfPowerDiv2);

            const double dfInvW = 1.0 / dfW;

            dfNominator += dfInvW * padfZ[i];

            dfDenominator += dfInvW;

        }

    }

    if (i != nPoints)

    {

        *pdfValue = padfZ[i];

    }

    else if (dfDenominator == 0.0)

    {

        *pdfValue = poOptions->dfNoDataValue;

    }

    else

    {

        *pdfValue = dfNominator / dfDenominator;

    }

    return CE_None;

}

/************************************************************************/

/*                       GDALGridMovingAverage()                        */

/************************************************************************/

/**

 * Moving average.

 *

 * The Moving Average is a simple data averaging algorithm. It uses a moving

 * window of elliptic form to search values and averages all data points

 * within the window. Search ellipse can be rotated by specified angle, the

 * center of ellipse located at the grid node. Also the minimum number of data

 * points to average can be set, if there are not enough points in window, the

 * grid node considered empty and will be filled with specified NODATA value.

 *

 * Mathematically it can be expressed with the formula:

 *

 * \f[

 *      Z=\frac{\sum_{i=1}^n{Z_i}}{n}

 * \f]

 *

 *  where

 *  <ul>

 *      <li> \f$Z\f$ is a resulting value at the grid node,

 *      <li> \f$Z_i\f$ is a known value at point \f$i\f$,

 *      <li> \f$n\f$ is a total number of points in search ellipse.

 *  </ul>

 *

 * @param poOptionsIn Algorithm parameters. This should point to

 * GDALGridMovingAverageOptions object.

 * @param nPoints Number of elements in input arrays.

 * @param padfX Input array of X coordinates.

 * @param padfY Input array of Y coordinates.

 * @param padfZ Input array of Z values.

 * @param dfXPoint X coordinate of the point to compute.

 * @param dfYPoint Y coordinate of the point to compute.

 * @param pdfValue Pointer to variable where the computed grid node value

 * will be returned.

 * @param hExtraParamsIn extra parameters (unused)

 *

 * @return CE_None on success or CE_Failure if something goes wrong.

 */

CPLErr GDALGridMovingAverage(const void *poOptionsIn, GUInt32 nPoints,

                             const double *padfX, const double *padfY,

                             const double *padfZ, double dfXPoint,

                             double dfYPoint, double *pdfValue,

                             CPL_UNUSED void *hExtraParamsIn)

{

    // TODO: For optimization purposes pre-computed parameters should be moved

    // out of this routine to the calling function.

    const GDALGridMovingAverageOptions *const poOptions =

        static_cast<const GDALGridMovingAverageOptions *>(poOptionsIn);

    // Pre-compute search ellipse parameters.

    const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;

    const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;

    const double dfSearchRadius =

        std::max(poOptions->dfRadius1, poOptions->dfRadius2);

    const double dfR12Square = dfRadius1Square * dfRadius2Square;

    GDALGridExtraParameters *psExtraParams =

        static_cast<GDALGridExtraParameters *>(hExtraParamsIn);

    const CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;

    // Compute coefficients for coordinate system rotation.

    const double dfAngle = TO_RADIANS * poOptions->dfAngle;

    const bool bRotated = dfAngle != 0.0;

    const double dfCoeff1 = bRotated ? cos(dfAngle) : 0.0;

    const double dfCoeff2 = bRotated ? sin(dfAngle) : 0.0;

    double dfAccumulator = 0.0;

    GUInt32 n = 0;  // Used after for.

    if (phQuadTree != nullptr)

    {

        CPLRectObj sAoi;

        sAoi.minx = dfXPoint - dfSearchRadius;

        sAoi.miny = dfYPoint - dfSearchRadius;

        sAoi.maxx = dfXPoint + dfSearchRadius;

        sAoi.maxy = dfYPoint + dfSearchRadius;

        int nFeatureCount = 0;

        GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(

            CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));

        if (nFeatureCount != 0)

        {

            for (int k = 0; k < nFeatureCount; k++)

            {

                const int i = papsPoints[k]->i;

                double dfRX = padfX[i] - dfXPoint;

                double dfRY = padfY[i] - dfYPoint;

                if (bRotated)

                {

                    const double dfRXRotated =

                        dfRX * dfCoeff1 + dfRY * dfCoeff2;

                    const double dfRYRotated =

                        dfRY * dfCoeff1 - dfRX * dfCoeff2;

                    dfRX = dfRXRotated;

                    dfRY = dfRYRotated;

                }

                if (dfRadius2Square * dfRX * dfRX +

                        dfRadius1Square * dfRY * dfRY <=

                    dfR12Square)

                {

                    dfAccumulator += padfZ[i];

                    n++;

                }

            }

        }

        CPLFree(papsPoints);

    }

    else

    {

        for (GUInt32 i = 0; i < nPoints; i++)

        {

            double dfRX = padfX[i] - dfXPoint;

            double dfRY = padfY[i] - dfYPoint;

            if (bRotated)

            {

                const double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;

                const double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;

                dfRX = dfRXRotated;

                dfRY = dfRYRotated;

            }

            // Is this point located inside the search ellipse?

            if (dfRadius2Square * dfRX * dfRX + dfRadius1Square * dfRY * dfRY <=

                dfR12Square)

            {

                dfAccumulator += padfZ[i];

                n++;

            }

        }

    }

    if (n < poOptions->nMinPoints || n == 0)

    {

        *pdfValue = poOptions->dfNoDataValue;

    }

    else

    {

        *pdfValue = dfAccumulator / n;

    }

    return CE_None;

}

/************************************************************************/

/*                  GDALGridMovingAveragePerQuadrant()                  */

/************************************************************************/

/**

 * Moving average, with a per-quadrant search logic.

 */

static CPLErr GDALGridMovingAveragePerQuadrant(

    const void *poOptionsIn, GUInt32 /*nPoints*/, const double *padfX,

    const double *padfY, const double *padfZ, double dfXPoint, double dfYPoint,

    double *pdfValue, void *hExtraParamsIn)

{

    const GDALGridMovingAverageOptions *const poOptions =

        static_cast<const GDALGridMovingAverageOptions *>(poOptionsIn);

    const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;

    const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;

    const double dfR12Square = dfRadius1Square * dfRadius2Square;

    const GUInt32 nMaxPoints = poOptions->nMaxPoints;

    const GUInt32 nMinPointsPerQuadrant = poOptions->nMinPointsPerQuadrant;

    const GUInt32 nMaxPointsPerQuadrant = poOptions->nMaxPointsPerQuadrant;

    // Compute coefficients for coordinate system rotation.

    const double dfAngle = TO_RADIANS * poOptions->dfAngle;

    const bool bRotated = dfAngle != 0.0;

    const double dfCoeff1 = bRotated ? cos(dfAngle) : 0.0;

    const double dfCoeff2 = bRotated ? sin(dfAngle) : 0.0;

    GDALGridExtraParameters *psExtraParams =

        static_cast<GDALGridExtraParameters *>(hExtraParamsIn);

    const CPLQuadTree *phQuadTree = psExtraParams->hQuadTree;

    CPLAssert(phQuadTree);

    std::multimap<double, double> oMapDistanceToZValuesPerQuadrant[4];

    const double dfSearchRadius =

        std::max(poOptions->dfRadius1, poOptions->dfRadius2);

    CPLRectObj sAoi;

    sAoi.minx = dfXPoint - dfSearchRadius;

    sAoi.miny = dfYPoint - dfSearchRadius;

    sAoi.maxx = dfXPoint + dfSearchRadius;

    sAoi.maxy = dfYPoint + dfSearchRadius;

    int nFeatureCount = 0;

    GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>(

        CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount));

    if (nFeatureCount != 0)

    {

        for (int k = 0; k < nFeatureCount; k++)

        {

            const int i = papsPoints[k]->i;

            double dfRX = padfX[i] - dfXPoint;

            double dfRY = padfY[i] - dfYPoint;

            if (bRotated)

            {

                const double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;

                const double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;

                dfRX = dfRXRotated;

                dfRY = dfRYRotated;

            }

            const double dfRXSquare = dfRX * dfRX;

            const double dfRYSquare = dfRY * dfRY;

            if (dfRadius2Square * dfRXSquare + dfRadius1Square * dfRYSquare <=

                dfR12Square)

            {

                const int iQuadrant =

                    ((dfRX >= 0) ? 1 : 0) | (((dfRY >= 0) ? 1 : 0) << 1);

                oMapDistanceToZValuesPerQuadrant[iQuadrant].insert(

                    std::make_pair(dfRXSquare + dfRYSquare, padfZ[i]));

            }

        }

    }

    CPLFree(papsPoints);

    std::multimap<double, double>::iterator aoIter[] = {

        oMapDistanceToZValuesPerQuadrant[0].begin(),

        oMapDistanceToZValuesPerQuadrant[1].begin(),

        oMapDistanceToZValuesPerQuadrant[2].begin(),

        oMapDistanceToZValuesPerQuadrant[3].begin(),

    };

    constexpr int ALL_QUADRANT_FLAGS = 1 + 2 + 4 + 8;

    // Examine all "neighbors" within the radius (sorted by distance via the

    // multimap), and use the closest n points based on distance until the max

    // is reached.

    // Do that by fetching the nearest point in quadrant 0, then the nearest

    // point in quadrant 1, 2 and 3, and starting again with the next nearest

    // point in quarant 0, etc.

    int nQuadrantIterFinishedFlag = 0;

    GUInt32 anPerQuadrant[4] = {0};

    double dfNominator = 0.0;

    GUInt32 n = 0;

    for (int iQuadrant = 0; /* true */; iQuadrant = (iQuadrant + 1) % 4)

    {

        if (aoIter[iQuadrant] ==

                oMapDistanceToZValuesPerQuadrant[iQuadrant].end() ||

            (nMaxPointsPerQuadrant > 0 &&

             anPerQuadrant[iQuadrant] >= nMaxPointsPerQuadrant))

        {

            nQuadrantIterFinishedFlag |= 1 << iQuadrant;

            if (nQuadrantIterFinishedFlag == ALL_QUADRANT_FLAGS)

                break;

            continue;

        }

        const double dfZ = aoIter[iQuadrant]->second;

        ++aoIter[iQuadrant];

        dfNominator += dfZ;

        n++;

        anPerQuadrant[iQuadrant]++;

        if (nMaxPoints > 0 && n >= nMaxPoints)

        {

            break;

        }

    }

    if (nMinPointsPerQuadrant > 0 &&

        (anPerQuadrant[0] < nMinPointsPerQuadrant ||

         anPerQuadrant[1] < nMinPointsPerQuadrant ||

         anPerQuadrant[2] < nMinPointsPerQuadrant ||

         anPerQuadrant[3] < nMinPointsPerQuadrant))

    {

        *pdfValue = poOptions->dfNoDataValue;

    }

    else if (n < poOptions->nMinPoints || n == 0)

    {

        *pdfValue = poOptions->dfNoDataValue;

    }

    else

    {

        *pdfValue = dfNominator / n;

    }

    return CE_None;

}

/************************************************************************/

/*                      GDALGridNearestNeighbor()                       */

/************************************************************************/

/**

 * Nearest neighbor.

 *

 * The Nearest Neighbor method doesn't perform any interpolation or smoothing,

 * it just takes the value of nearest point found in grid node search ellipse

 * and returns it as a result. If there are no points found, the specified

 * NODATA value will be returned.

 *

 * @param poOptionsIn Algorithm parameters. This should point to

 * GDALGridNearestNeighborOptions object.

 * @param nPoints Number of elements in input arrays.

 * @param padfX Input array of X coordinates.

 * @param padfY Input array of Y coordinates.

 * @param padfZ Input array of Z values.

 * @param dfXPoint X coordinate of the point to compute.

 * @param dfYPoint Y coordinate of the point to compute.

 * @param pdfValue Pointer to variable where the computed grid node value

 * will be returned.

 * @param hExtraParamsIn extra parameters.

 *

 * @return CE_None on success or CE_Failure if something goes wrong.

 */

CPLErr GDALGridNearestNeighbor(const void *poOptionsIn, GUInt32 nPoints,

                               const double *padfX, const double *padfY,

                               const double *padfZ, double dfXPoint,

                               double dfYPoint, double *pdfValue,

                               void *hExtraParamsIn)

{

    // TODO: For optimization purposes pre-computed parameters should be moved

    // out of this routine to the calling function.

    const GDALGridNearestNeighborOptions *const poOptions =

        static_cast<const GDALGridNearestNeighborOptions *>(poOptionsIn);

    // Pre-compute search ellipse parameters.

    const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1;

    const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2;

    const double dfR12Square = dfRadius1Square * dfRadius2Square;

    GDALGridExtraParameters *psExtraParams =

        static_cast<GDALGridExtraParameters *>(hExtraParamsIn);

    CPLQuadTree *hQuadTree = psExtraParams->hQuadTree;

    // Compute coefficients for coordinate system rotation.

    const double dfAngle = TO_RADIANS * poOptions->dfAngle;

    const bool bRotated = dfAngle != 0.0;

    const double dfCoeff1 = bRotated ? cos(dfAngle) : 0.0;

    const double dfCoeff2 = bRotated ? sin(dfAngle) : 0.0;

    // If the nearest point will not be found, its value remains as NODATA.

    double dfNearestValue = poOptions->dfNoDataValue;