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

 *

 * Project:  GDAL

 * Purpose:  Interpolate in nodata areas.

 * Author:   Frank Warmerdam, warmerdam@pobox.com

 *

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

 * Copyright (c) 2008, Frank Warmerdam

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

 * Copyright (c) 2015, Sean Gillies <sean@mapbox.com>

 *

 * SPDX-License-Identifier: MIT

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

#include "cpl_port.h"

#include "gdal_alg.h"

#include <cmath>

#include <cstring>

#include <algorithm>

#include <string>

#include <utility>

#include "cpl_conv.h"

#include "cpl_error.h"

#include "cpl_progress.h"

#include "cpl_string.h"

#include "cpl_vsi.h"

#include "gdal.h"

#include "gdal_priv.h"

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

/*                           GDALFilterLine()                           */

/*                                                                      */

/*      Apply 3x3 filtering one one scanline with masking for which     */

/*      pixels are to be interpolated (ThisFMask) and which window      */

/*      pixels are valid to include in the interpolation (TMask).       */

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

static void GDALFilterLine(const float *pafLastLine, const float *pafThisLine,

                           const float *pafNextLine, float *pafOutLine,

                           const GByte *pabyLastTMask,

                           const GByte *pabyThisTMask,

                           const GByte *pabyNextTMask,

                           const GByte *pabyThisFMask, int nXSize)

{

    for (int iX = 0; iX < nXSize; iX++)

    {

        if (!pabyThisFMask[iX])

        {

            pafOutLine[iX] = pafThisLine[iX];

            continue;

        }

        CPLAssert(pabyThisTMask[iX]);

        double dfValSum = 0.0;

        double dfWeightSum = 0.0;

        // Previous line.

        if (pafLastLine != nullptr)

        {

            if (iX > 0 && pabyLastTMask[iX - 1])

            {

                dfValSum += pafLastLine[iX - 1];

                dfWeightSum += 1.0;

            }

            if (pabyLastTMask[iX])

            {

                dfValSum += pafLastLine[iX];

                dfWeightSum += 1.0;

            }

            if (iX < nXSize - 1 && pabyLastTMask[iX + 1])

            {

                dfValSum += pafLastLine[iX + 1];

                dfWeightSum += 1.0;

            }

        }

        // Current Line.

        if (iX > 0 && pabyThisTMask[iX - 1])

        {

            dfValSum += pafThisLine[iX - 1];

            dfWeightSum += 1.0;

        }

        if (pabyThisTMask[iX])

        {

            dfValSum += pafThisLine[iX];

            dfWeightSum += 1.0;

        }

        if (iX < nXSize - 1 && pabyThisTMask[iX + 1])

        {

            dfValSum += pafThisLine[iX + 1];

            dfWeightSum += 1.0;

        }

        // Next line.

        if (pafNextLine != nullptr)

        {

            if (iX > 0 && pabyNextTMask[iX - 1])

            {

                dfValSum += pafNextLine[iX - 1];

                dfWeightSum += 1.0;

            }

            if (pabyNextTMask[iX])

            {

                dfValSum += pafNextLine[iX];

                dfWeightSum += 1.0;

            }

            if (iX < nXSize - 1 && pabyNextTMask[iX + 1])

            {

                dfValSum += pafNextLine[iX + 1];

                dfWeightSum += 1.0;

            }

        }

        pafOutLine[iX] = static_cast<float>(dfValSum / dfWeightSum);

    }

}

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

/*                          GDALMultiFilter()                           */

/*                                                                      */

/*      Apply multiple iterations of a 3x3 smoothing filter over a      */

/*      band with masking controlling what pixels should be             */

/*      filtered (FiltMaskBand non zero) and which pixels can be        */

/*      considered valid contributors to the filter                     */

/*      (TargetMaskBand non zero).                                      */

/*                                                                      */

/*      This implementation attempts to apply many iterations in        */

/*      one IO pass by managing the filtering over a rolling buffer     */

/*      of nIterations+2 scanlines.  While possibly clever this        */

/*      makes the algorithm implementation largely                      */

/*      incomprehensible.                                               */

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

static CPLErr GDALMultiFilter(GDALRasterBandH hTargetBand,

                              GDALRasterBandH hTargetMaskBand,

                              GDALRasterBandH hFiltMaskBand, int nIterations,

                              GDALProgressFunc pfnProgress, void *pProgressArg)

{

    const int nXSize = GDALGetRasterBandXSize(hTargetBand);

    const int nYSize = GDALGetRasterBandYSize(hTargetBand);

    /* -------------------------------------------------------------------- */

    /*      Report starting progress value.                                 */

    /* -------------------------------------------------------------------- */

    if (!pfnProgress(0.0, "Smoothing Filter...", pProgressArg))

    {

        CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");

        return CE_Failure;

    }

    /* -------------------------------------------------------------------- */

    /*      Allocate rotating buffers.                                      */

    /* -------------------------------------------------------------------- */

    const int nBufLines = nIterations + 2;

    GByte *pabyTMaskBuf =

        static_cast<GByte *>(VSI_MALLOC2_VERBOSE(nXSize, nBufLines));

    GByte *pabyFMaskBuf =

        static_cast<GByte *>(VSI_MALLOC2_VERBOSE(nXSize, nBufLines));

    float *paf3PassLineBuf = static_cast<float *>(

        VSI_MALLOC3_VERBOSE(nXSize, nBufLines, 3 * sizeof(float)));

    if (pabyTMaskBuf == nullptr || pabyFMaskBuf == nullptr ||

        paf3PassLineBuf == nullptr)

    {

        CPLFree(pabyTMaskBuf);

        CPLFree(pabyFMaskBuf);

        CPLFree(paf3PassLineBuf);

        return CE_Failure;

    }

    /* -------------------------------------------------------------------- */

    /*      Process rotating buffers.                                       */

    /* -------------------------------------------------------------------- */

    CPLErr eErr = CE_None;

    int iPassCounter = 0;

    for (int nNewLine = 0;  // Line being loaded (zero based scanline).

         eErr == CE_None && nNewLine < nYSize + nIterations; nNewLine++)

    {

        /* --------------------------------------------------------------------

         */

        /*      Rotate pass buffers. */

        /* --------------------------------------------------------------------

         */

        iPassCounter = (iPassCounter + 1) % 3;

        float *const pafSLastPass =

            paf3PassLineBuf + ((iPassCounter + 0) % 3) * nXSize * nBufLines;

        float *const pafLastPass =

            paf3PassLineBuf + ((iPassCounter + 1) % 3) * nXSize * nBufLines;

        float *const pafThisPass =

            paf3PassLineBuf + ((iPassCounter + 2) % 3) * nXSize * nBufLines;

        /* --------------------------------------------------------------------

         */

        /*      Where does the new line go in the rotating buffer? */

        /* --------------------------------------------------------------------

         */

        const int iBufOffset = nNewLine % nBufLines;

        /* --------------------------------------------------------------------

         */

        /*      Read the new data line if it is't off the bottom of the */

        /*      image. */

        /* --------------------------------------------------------------------

         */

        if (nNewLine < nYSize)

        {

            eErr = GDALRasterIO(hTargetMaskBand, GF_Read, 0, nNewLine, nXSize,

                                1, pabyTMaskBuf + nXSize * iBufOffset, nXSize,

                                1, GDT_Byte, 0, 0);

            if (eErr != CE_None)

                break;

            eErr = GDALRasterIO(hFiltMaskBand, GF_Read, 0, nNewLine, nXSize, 1,

                                pabyFMaskBuf + nXSize * iBufOffset, nXSize, 1,

                                GDT_Byte, 0, 0);

            if (eErr != CE_None)

                break;

            eErr = GDALRasterIO(hTargetBand, GF_Read, 0, nNewLine, nXSize, 1,

                                pafThisPass + nXSize * iBufOffset, nXSize, 1,

                                GDT_Float32, 0, 0);

            if (eErr != CE_None)

                break;

        }

        /* --------------------------------------------------------------------

         */

        /*      Loop over the loaded data, applying the filter to all loaded */

        /*      lines with neighbours. */

        /* --------------------------------------------------------------------

         */

        for (int iFLine = nNewLine - 1;

             eErr == CE_None && iFLine >= nNewLine - nIterations; iFLine--)

        {

            const int iLastOffset = (iFLine - 1) % nBufLines;

            const int iThisOffset = (iFLine) % nBufLines;

            const int iNextOffset = (iFLine + 1) % nBufLines;

            // Default to preserving the old value.

            if (iFLine >= 0)

                memcpy(pafThisPass + iThisOffset * nXSize,

                       pafLastPass + iThisOffset * nXSize,

                       sizeof(float) * nXSize);

            // TODO: Enable first and last line.

            // Skip the first and last line.

            if (iFLine < 1 || iFLine >= nYSize - 1)

            {

                continue;

            }

            GDALFilterLine(pafSLastPass + iLastOffset * nXSize,

                           pafLastPass + iThisOffset * nXSize,

                           pafThisPass + iNextOffset * nXSize,

                           pafThisPass + iThisOffset * nXSize,

                           pabyTMaskBuf + iLastOffset * nXSize,

                           pabyTMaskBuf + iThisOffset * nXSize,

                           pabyTMaskBuf + iNextOffset * nXSize,

                           pabyFMaskBuf + iThisOffset * nXSize, nXSize);

        }

        /* --------------------------------------------------------------------

         */

        /*      Write out the top data line that will be rolling out of our */

        /*      buffer. */

        /* --------------------------------------------------------------------

         */

        const int iLineToSave = nNewLine - nIterations;

        if (iLineToSave >= 0 && eErr == CE_None)

        {

            const int iBufOffset2 = iLineToSave % nBufLines;

            eErr = GDALRasterIO(hTargetBand, GF_Write, 0, iLineToSave, nXSize,

                                1, pafThisPass + nXSize * iBufOffset2, nXSize,

                                1, GDT_Float32, 0, 0);

        }

        /* --------------------------------------------------------------------

         */

        /*      Report progress. */

        /* --------------------------------------------------------------------

         */

        if (eErr == CE_None &&

            !pfnProgress((nNewLine + 1) /

                             static_cast<double>(nYSize + nIterations),

                         "Smoothing Filter...", pProgressArg))

        {

            CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");

            eErr = CE_Failure;

        }

    }

    /* -------------------------------------------------------------------- */

    /*      Cleanup                                                         */

    /* -------------------------------------------------------------------- */

    CPLFree(pabyTMaskBuf);

    CPLFree(pabyFMaskBuf);

    CPLFree(paf3PassLineBuf);

    return eErr;

}

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

/*                             QUAD_CHECK()                             */

/*                                                                      */

/*      macro for checking whether a point is nearer than the           */

/*      existing closest point.                                         */

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

inline void QUAD_CHECK(double &dfQuadDist, float &fQuadValue, int target_x,

                       GUInt32 target_y, int origin_x, int origin_y,

                       float fTargetValue, GUInt32 nNoDataVal)

{

    if (target_y != nNoDataVal)

    {

        const double dfDx =

            static_cast<double>(target_x) - static_cast<double>(origin_x);

        const double dfDy =

            static_cast<double>(target_y) - static_cast<double>(origin_y);

        double dfDistSq = dfDx * dfDx + dfDy * dfDy;

        if (dfDistSq < dfQuadDist * dfQuadDist)

        {

            CPLAssert(dfDistSq > 0.0);

            dfQuadDist = sqrt(dfDistSq);

            fQuadValue = fTargetValue;

        }

    }

}

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

/*                           GDALFillNodata()                           */

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

/**

 * Fill selected raster regions by interpolation from the edges.

 *

 * This algorithm will interpolate values for all designated

 * nodata pixels (marked by zeros in hMaskBand). For each pixel

 * a four direction conic search is done to find values to interpolate

 * from (using inverse distance weighting by default). Once all values are

 * interpolated, zero or more smoothing iterations (3x3 average

 * filters on interpolated pixels) are applied to smooth out

 * artifacts.

 *

 * This algorithm is generally suitable for interpolating missing

 * regions of fairly continuously varying rasters (such as elevation

 * models for instance). It is also suitable for filling small holes

 * and cracks in more irregularly varying images (like airphotos). It

 * is generally not so great for interpolating a raster from sparse

 * point data - see the algorithms defined in gdal_grid.h for that case.

 *

 * @param hTargetBand the raster band to be modified in place.

 * @param hMaskBand a mask band indicating pixels to be interpolated

 * (zero valued). If hMaskBand is set to NULL, this method will internally use

 * the mask band returned by GDALGetMaskBand(hTargetBand).

 * @param dfMaxSearchDist the maximum number of pixels to search in all

 * directions to find values to interpolate from.

 * @param bDeprecatedOption unused argument, should be zero.

 * @param nSmoothingIterations the number of 3x3 smoothing filter passes to

 * run (0 or more).

 * @param papszOptions additional name=value options in a string list.

 * <ul>

 * <li>TEMP_FILE_DRIVER=gdal_driver_name. For example MEM.</li>

 * <li>NODATA=value (starting with GDAL 2.4).

 * Source pixels at that value will be ignored by the interpolator. Warning:

 * currently this will not be honored by smoothing passes.</li>

 * <li>INTERPOLATION=INV_DIST/NEAREST (GDAL >= 3.9). By default, pixels are

 * interpolated using an inverse distance weighting (INV_DIST). It is also

 * possible to choose a nearest neighbour (NEAREST) strategy.</li>

 * </ul>

 * @param pfnProgress the progress function to report completion.

 * @param pProgressArg callback data for progress function.

 *

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

 */

CPLErr CPL_STDCALL GDALFillNodata(GDALRasterBandH hTargetBand,

                                  GDALRasterBandH hMaskBand,

                                  double dfMaxSearchDist,

                                  CPL_UNUSED int bDeprecatedOption,

                                  int nSmoothingIterations, char **papszOptions,

                                  GDALProgressFunc pfnProgress,

                                  void *pProgressArg)

{

    VALIDATE_POINTER1(hTargetBand, "GDALFillNodata", CE_Failure);

    const int nXSize = GDALGetRasterBandXSize(hTargetBand);

    const int nYSize = GDALGetRasterBandYSize(hTargetBand);

    if (dfMaxSearchDist == 0.0)

        dfMaxSearchDist = std::max(nXSize, nYSize) + 1;

    const int nMaxSearchDist = static_cast<int>(floor(dfMaxSearchDist));

    const char *pszInterpolation =

        CSLFetchNameValueDef(papszOptions, "INTERPOLATION", "INV_DIST");

    const bool bNearest = EQUAL(pszInterpolation, "NEAREST");

    if (!EQUAL(pszInterpolation, "INV_DIST") &&

        !EQUAL(pszInterpolation, "NEAREST"))

    {

        CPLError(CE_Failure, CPLE_NotSupported,

                 "Unsupported interpolation method: %s", pszInterpolation);

        return CE_Failure;

    }

    // Special "x" pixel values identifying pixels as special.

    GDALDataType eType = GDT_UInt16;

    GUInt32 nNoDataVal = 65535;

    if (nXSize > 65533 || nYSize > 65533)

    {

        eType = GDT_UInt32;

        nNoDataVal = 4000002;

    }

    /* -------------------------------------------------------------------- */

    /*      Determine format driver for temp work files.                    */

    /* -------------------------------------------------------------------- */

    CPLString osTmpFileDriver =

        CSLFetchNameValueDef(papszOptions, "TEMP_FILE_DRIVER", "GTiff");

    GDALDriverH hDriver = GDALGetDriverByName(osTmpFileDriver.c_str());

    if (hDriver == nullptr)

    {

        CPLError(CE_Failure, CPLE_AppDefined,

                 "TEMP_FILE_DRIVER=%s driver is not registered",

                 osTmpFileDriver.c_str());

        return CE_Failure;

    }

    if (GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATE, nullptr) == nullptr)

    {

        CPLError(CE_Failure, CPLE_AppDefined,

                 "TEMP_FILE_DRIVER=%s driver is incapable of creating "

                 "temp work files",

                 osTmpFileDriver.c_str());

        return CE_Failure;

    }

    CPLStringList aosWorkFileOptions;

    if (osTmpFileDriver == "GTiff")

    {

        aosWorkFileOptions.SetNameValue("COMPRESS", "LZW");

        aosWorkFileOptions.SetNameValue("BIGTIFF", "IF_SAFER");

    }

    const CPLString osTmpFile = CPLGenerateTempFilenameSafe("");

    std::unique_ptr<GDALDataset> poTmpMaskDS;

    if (hMaskBand == nullptr)

    {

        hMaskBand = GDALGetMaskBand(hTargetBand);

    }

    else if (nSmoothingIterations > 0 &&

             hMaskBand != GDALGetMaskBand(hTargetBand))

    {

        // If doing smoothing operations and the user provided its own

        // mask band, we must make a copy of it to be able to update it

        // when we fill pixels during the initial pass.

        const CPLString osMaskTmpFile = osTmpFile + "fill_mask_work.tif";

        poTmpMaskDS.reset(GDALDataset::FromHandle(

            GDALCreate(hDriver, osMaskTmpFile, nXSize, nYSize, 1, GDT_Byte,

                       aosWorkFileOptions.List())));

        if (poTmpMaskDS == nullptr)

        {

            CPLError(CE_Failure, CPLE_AppDefined,

                     "Could not create poTmpMaskDS work file. Check driver "

                     "capabilities.");

            return CE_Failure;

        }

        poTmpMaskDS->MarkSuppressOnClose();

        auto hTmpMaskBand =

            GDALRasterBand::ToHandle(poTmpMaskDS->GetRasterBand(1));

        if (GDALRasterBandCopyWholeRaster(hMaskBand, hTmpMaskBand, nullptr,

                                          nullptr, nullptr) != CE_None)

        {

            return CE_Failure;

        }

        hMaskBand = hTmpMaskBand;

    }

    // If there are smoothing iterations, reserve 10% of the progress for them.

    const double dfProgressRatio = nSmoothingIterations > 0 ? 0.9 : 1.0;

    const char *pszNoData = CSLFetchNameValue(papszOptions, "NODATA");

    bool bHasNoData = false;

    float fNoData = 0.0f;

    if (pszNoData)

    {

        bHasNoData = true;

        fNoData = static_cast<float>(CPLAtof(pszNoData));

    }

    /* -------------------------------------------------------------------- */

    /*      Initialize progress counter.                                    */

    /* -------------------------------------------------------------------- */

    if (pfnProgress == nullptr)

        pfnProgress = GDALDummyProgress;

    if (!pfnProgress(0.0, "Filling...", pProgressArg))

    {

        CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");

        return CE_Failure;

    }

    /* -------------------------------------------------------------------- */

    /*      Create a work file to hold the Y "last value" indices.          */

    /* -------------------------------------------------------------------- */

    const CPLString osYTmpFile = osTmpFile + "fill_y_work.tif";

    auto poYDS = std::unique_ptr<GDALDataset>(GDALDataset::FromHandle(

        GDALCreate(hDriver, osYTmpFile, nXSize, nYSize, 1, eType,

                   aosWorkFileOptions.List())));

    if (poYDS == nullptr)

    {

        CPLError(

            CE_Failure, CPLE_AppDefined,

            "Could not create Y index work file. Check driver capabilities.");

        return CE_Failure;

    }

    poYDS->MarkSuppressOnClose();

    GDALRasterBandH hYBand =

        GDALRasterBand::FromHandle(poYDS->GetRasterBand(1));

    /* -------------------------------------------------------------------- */

    /*      Create a work file to hold the pixel value associated with      */

    /*      the "last xy value" pixel.                                      */

    /* -------------------------------------------------------------------- */

    const CPLString osValTmpFile = osTmpFile + "fill_val_work.tif";

    auto poValDS =

        std::unique_ptr<GDALDataset>(GDALDataset::FromHandle(GDALCreate(

            hDriver, osValTmpFile, nXSize, nYSize, 1,

            GDALGetRasterDataType(hTargetBand), aosWorkFileOptions.List())));

    if (poValDS == nullptr)

    {

        CPLError(

            CE_Failure, CPLE_AppDefined,

            "Could not create XY value work file. Check driver capabilities.");

        return CE_Failure;

    }

    poValDS->MarkSuppressOnClose();

    GDALRasterBandH hValBand =

        GDALRasterBand::FromHandle(poValDS->GetRasterBand(1));

    /* -------------------------------------------------------------------- */

    /*      Create a mask file to make it clear what pixels can be filtered */

    /*      on the filtering pass.                                          */

    /* -------------------------------------------------------------------- */

    const CPLString osFiltMaskTmpFile = osTmpFile + "fill_filtmask_work.tif";

    auto poFiltMaskDS = std::unique_ptr<GDALDataset>(GDALDataset::FromHandle(

        GDALCreate(hDriver, osFiltMaskTmpFile, nXSize, nYSize, 1, GDT_Byte,

                   aosWorkFileOptions.List())));

    if (poFiltMaskDS == nullptr)

    {

        CPLError(CE_Failure, CPLE_AppDefined,

                 "Could not create mask work file. Check driver capabilities.");

        return CE_Failure;

    }

    poFiltMaskDS->MarkSuppressOnClose();

    GDALRasterBandH hFiltMaskBand =

        GDALRasterBand::FromHandle(poFiltMaskDS->GetRasterBand(1));

    /* -------------------------------------------------------------------- */

    /*      Allocate buffers for last scanline and this scanline.           */

    /* -------------------------------------------------------------------- */

    GUInt32 *panLastY =

        static_cast<GUInt32 *>(VSI_CALLOC_VERBOSE(nXSize, sizeof(GUInt32)));

    GUInt32 *panThisY =

        static_cast<GUInt32 *>(VSI_CALLOC_VERBOSE(nXSize, sizeof(GUInt32)));

    GUInt32 *panTopDownY =

        static_cast<GUInt32 *>(VSI_CALLOC_VERBOSE(nXSize, sizeof(GUInt32)));

    float *pafLastValue =

        static_cast<float *>(VSI_CALLOC_VERBOSE(nXSize, sizeof(float)));

    float *pafThisValue =

        static_cast<float *>(VSI_CALLOC_VERBOSE(nXSize, sizeof(float)));

    float *pafTopDownValue =

        static_cast<float *>(VSI_CALLOC_VERBOSE(nXSize, sizeof(float)));

    float *pafScanline =

        static_cast<float *>(VSI_CALLOC_VERBOSE(nXSize, sizeof(float)));

    GByte *pabyMask = static_cast<GByte *>(VSI_CALLOC_VERBOSE(nXSize, 1));

    GByte *pabyFiltMask = static_cast<GByte *>(VSI_CALLOC_VERBOSE(nXSize, 1));

    CPLErr eErr = CE_None;

    if (panLastY == nullptr || panThisY == nullptr || panTopDownY == nullptr ||

        pafLastValue == nullptr || pafThisValue == nullptr ||

        pafTopDownValue == nullptr || pafScanline == nullptr ||

        pabyMask == nullptr || pabyFiltMask == nullptr)

    {

        eErr = CE_Failure;

        goto end;

    }

    for (int iX = 0; iX < nXSize; iX++)

    {

        panLastY[iX] = nNoDataVal;

    }

    /* ==================================================================== */

    /*      Make first pass from top to bottom collecting the "last         */

    /*      known value" for each column and writing it out to the work     */

    /*      files.                                                          */

    /* ==================================================================== */

    for (int iY = 0; iY < nYSize && eErr == CE_None; iY++)

    {

        /* --------------------------------------------------------------------

         */

        /*      Read data and mask for this line. */

        /* --------------------------------------------------------------------

         */

        eErr = GDALRasterIO(hMaskBand, GF_Read, 0, iY, nXSize, 1, pabyMask,

                            nXSize, 1, GDT_Byte, 0, 0);

        if (eErr != CE_None)

            break;

        eErr = GDALRasterIO(hTargetBand, GF_Read, 0, iY, nXSize, 1, pafScanline,

                            nXSize, 1, GDT_Float32, 0, 0);

        if (eErr != CE_None)

            break;

        /* --------------------------------------------------------------------

         */

        /*      Figure out the most recent pixel for each column. */

        /* --------------------------------------------------------------------

         */

        for (int iX = 0; iX < nXSize; iX++)

        {

            if (pabyMask[iX])

            {

                pafThisValue[iX] = pafScanline[iX];

                panThisY[iX] = iY;

            }

            else if (iY <= dfMaxSearchDist + panLastY[iX])

            {

                pafThisValue[iX] = pafLastValue[iX];

                panThisY[iX] = panLastY[iX];

            }

            else

            {

                panThisY[iX] = nNoDataVal;

            }

        }

        /* --------------------------------------------------------------------

         */

        /*      Write out best index/value to working files. */

        /* --------------------------------------------------------------------

         */

        eErr = GDALRasterIO(hYBand, GF_Write, 0, iY, nXSize, 1, panThisY,

                            nXSize, 1, GDT_UInt32, 0, 0);

        if (eErr != CE_None)

            break;

        eErr = GDALRasterIO(hValBand, GF_Write, 0, iY, nXSize, 1, pafThisValue,

                            nXSize, 1, GDT_Float32, 0, 0);

        if (eErr != CE_None)

            break;

        /* --------------------------------------------------------------------

         */

        /*      Flip this/last buffers. */

        /* --------------------------------------------------------------------

         */

        std::swap(pafThisValue, pafLastValue);

        std::swap(panThisY, panLastY);

        /* --------------------------------------------------------------------

         */

        /*      report progress. */

        /* --------------------------------------------------------------------

         */

        if (!pfnProgress(dfProgressRatio *

                             (0.5 * (iY + 1) / static_cast<double>(nYSize)),

                         "Filling...", pProgressArg))

        {

            CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");

            eErr = CE_Failure;

        }

    }

    for (int iX = 0; iX < nXSize; iX++)

    {

        panLastY[iX] = nNoDataVal;

    }

    /* ==================================================================== */

    /*      Now we will do collect similar this/last information from       */

    /*      bottom to top and use it in combination with the top to         */

    /*      bottom search info to interpolate.                              */

    /* ==================================================================== */

    for (int iY = nYSize - 1; iY >= 0 && eErr == CE_None; iY--)

    {

        eErr = GDALRasterIO(hMaskBand, GF_Read, 0, iY, nXSize, 1, pabyMask,

                            nXSize, 1, GDT_Byte, 0, 0);

        if (eErr != CE_None)

            break;

        eErr = GDALRasterIO(hTargetBand, GF_Read, 0, iY, nXSize, 1, pafScanline,

                            nXSize, 1, GDT_Float32, 0, 0);

        if (eErr != CE_None)

            break;

        /* --------------------------------------------------------------------

         */

        /*      Figure out the most recent pixel for each column. */

        /* --------------------------------------------------------------------

         */

        for (int iX = 0; iX < nXSize; iX++)

        {

            if (pabyMask[iX])

            {

                pafThisValue[iX] = pafScanline[iX];

                panThisY[iX] = iY;

            }

            else if (panLastY[iX] - iY <= dfMaxSearchDist)

            {

                pafThisValue[iX] = pafLastValue[iX];

                panThisY[iX] = panLastY[iX];

            }

            else

            {

                panThisY[iX] = nNoDataVal;

            }

        }

        /* --------------------------------------------------------------------

         */

        /*      Load the last y and corresponding value from the top down pass.

         */

        /* --------------------------------------------------------------------

         */

        eErr = GDALRasterIO(hYBand, GF_Read, 0, iY, nXSize, 1, panTopDownY,

                            nXSize, 1, GDT_UInt32, 0, 0);

        if (eErr != CE_None)

            break;

        eErr = GDALRasterIO(hValBand, GF_Read, 0, iY, nXSize, 1,

                            pafTopDownValue, nXSize, 1, GDT_Float32, 0, 0);

        if (eErr != CE_None)

            break;

        /* --------------------------------------------------------------------

         */

        /*      Attempt to interpolate any pixels that are nodata. */

        /* --------------------------------------------------------------------

         */

        memset(pabyFiltMask, 0, nXSize);

        for (int iX = 0; iX < nXSize; iX++)

        {

            int nThisMaxSearchDist = nMaxSearchDist;

            // If this was a valid target - no change.

            if (pabyMask[iX])

                continue;

            enum Quadrants

            {

                QUAD_TOP_LEFT = 0,

                QUAD_BOTTOM_LEFT = 1,

                QUAD_TOP_RIGHT = 2,

                QUAD_BOTTOM_RIGHT = 3,

            };

            constexpr int QUAD_COUNT = 4;

            double adfQuadDist[QUAD_COUNT] = {};

            float afQuadValue[QUAD_COUNT] = {};

            for (int iQuad = 0; iQuad < QUAD_COUNT; iQuad++)

            {

                adfQuadDist[iQuad] = dfMaxSearchDist + 1.0;

                afQuadValue[iQuad] = 0.0;

            }

            // Step left and right by one pixel searching for the closest

            // target value for each quadrant.

            for (int iStep = 0; iStep <= nThisMaxSearchDist; iStep++)

            {

                const int iLeftX = std::max(0, iX - iStep);

                const int iRightX = std::min(nXSize - 1, iX + iStep);

                // Top left includes current line.

                QUAD_CHECK(adfQuadDist[QUAD_TOP_LEFT],

                           afQuadValue[QUAD_TOP_LEFT], iLeftX,

                           panTopDownY[iLeftX], iX, iY, pafTopDownValue[iLeftX],

                           nNoDataVal);

                // Bottom left.

                QUAD_CHECK(adfQuadDist[QUAD_BOTTOM_LEFT],

                           afQuadValue[QUAD_BOTTOM_LEFT], iLeftX,

                           panLastY[iLeftX], iX, iY, pafLastValue[iLeftX],

                           nNoDataVal);

                // Top right and bottom right do no include center pixel.

                if (iStep == 0)

                    continue;

                // Top right includes current line.

                QUAD_CHECK(adfQuadDist[QUAD_TOP_RIGHT],

                           afQuadValue[QUAD_TOP_RIGHT], iRightX,

                           panTopDownY[iRightX], iX, iY,

                           pafTopDownValue[iRightX], nNoDataVal);

                // Bottom right.

                QUAD_CHECK(adfQuadDist[QUAD_BOTTOM_RIGHT],

                           afQuadValue[QUAD_BOTTOM_RIGHT], iRightX,

                           panLastY[iRightX], iX, iY, pafLastValue[iRightX],

                           nNoDataVal);

                // Every four steps, recompute maximum distance.

                if ((iStep & 0x3) == 0)

                    nThisMaxSearchDist = static_cast<int>(floor(

                        std::max(std::max(adfQuadDist[0], adfQuadDist[1]),

                                 std::max(adfQuadDist[2], adfQuadDist[3]))));

            }

            bool bHasSrcValues = false;

            if (bNearest)

            {

                double dfNearestDist = dfMaxSearchDist + 1;

                float fNearestValue = 0.0f;

                for (int iQuad = 0; iQuad < QUAD_COUNT; iQuad++)

                {

                    if (adfQuadDist[iQuad] < dfNearestDist)

                    {

                        bHasSrcValues = true;

                        if (!bHasNoData || afQuadValue[iQuad] != fNoData)

                        {

                            fNearestValue = afQuadValue[iQuad];

                            dfNearestDist = adfQuadDist[iQuad];

                        }

                    }

                }

                if (bHasSrcValues)

                {

                    pabyFiltMask[iX] = 255;

                    if (dfNearestDist <= dfMaxSearchDist)

                    {

                        pabyMask[iX] = 255;

                        pafScanline[iX] = fNearestValue;

                    }

                    else

                        pafScanline[iX] = fNoData;

                }

            }

            else

            {

                double dfWeightSum = 0.0;

                double dfValueSum = 0.0;

                for (int iQuad = 0; iQuad < QUAD_COUNT; iQuad++)

                {

                    if (adfQuadDist[iQuad] <= dfMaxSearchDist)

                    {

                        bHasSrcValues = true;

                        if (!bHasNoData || afQuadValue[iQuad] != fNoData)

                        {

                            const double dfWeight = 1.0 / adfQuadDist[iQuad];

                            dfWeightSum += dfWeight;

                            dfValueSum += afQuadValue[iQuad] * dfWeight;

                        }

                    }

                }

                if (bHasSrcValues)

                {

                    pabyFiltMask[iX] = 255;

                    if (dfWeightSum > 0.0)

                    {

                        pabyMask[iX] = 255;

                        pafScanline[iX] =

                            static_cast<float>(dfValueSum / dfWeightSum);

                    }

                    else

                        pafScanline[iX] = fNoData;