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

 *

 * Project:  GDAL

 * Purpose:  Compute each pixel's proximity to a set of target pixels.

 * Author:   Frank Warmerdam, warmerdam@pobox.com

 *

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

 * Copyright (c) 2008, Frank Warmerdam

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

 *

 * SPDX-License-Identifier: MIT

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

#include "cpl_port.h"

#include "gdal_alg.h"

#include <cmath>

#include <cstdlib>

#include <algorithm>

#include "cpl_conv.h"

#include "cpl_error.h"

#include "cpl_progress.h"

#include "cpl_string.h"

#include "cpl_vsi.h"

#include "gdal.h"

static CPLErr ProcessProximityLine(GInt32 *panSrcScanline, int *panNearX,

                                   int *panNearY, int bForward, int iLine,

                                   int nXSize, double nMaxDist,

                                   float *pafProximity,

                                   double *pdfSrcNoDataValue, int nTargetValues,

                                   int *panTargetValues);

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

/*                        GDALComputeProximity()                        */

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

/**

Compute the proximity of all pixels in the image to a set of pixels in

the source image.

This function attempts to compute the proximity of all pixels in

the image to a set of pixels in the source image.  The following

options are used to define the behavior of the function.  By

default all non-zero pixels in hSrcBand will be considered the

"target", and all proximities will be computed in pixels.  Note

that target pixels are set to the value corresponding to a distance

of zero.

The progress function args may be NULL or a valid progress reporting function

such as GDALTermProgress/NULL.

Options:

  VALUES=n[,n]*

A list of target pixel values to measure the distance from.  If this

option is not provided proximity will be computed from non-zero

pixel values.  Currently pixel values are internally processed as

integers.

  DISTUNITS=[PIXEL]/GEO

Indicates whether distances will be computed in pixel units or

in georeferenced units.  The default is pixel units.  This also

determines the interpretation of MAXDIST.

  MAXDIST=n

The maximum distance to search.  Proximity distances greater than

this value will not be computed.  Instead output pixels will be

set to a nodata value.

  NODATA=n

The NODATA value to use on the output band for pixels that are

beyond MAXDIST.  If not provided, the hProximityBand will be

queried for a nodata value.  If one is not found, 65535 will be used.

  USE_INPUT_NODATA=YES/NO

If this option is set, the input data set no-data value will be

respected. Leaving no data pixels in the input as no data pixels in

the proximity output.

  FIXED_BUF_VAL=n

If this option is set, all pixels within the MAXDIST threshold are

set to this fixed value instead of to a proximity distance.

*/

CPLErr CPL_STDCALL GDALComputeProximity(GDALRasterBandH hSrcBand,

                                        GDALRasterBandH hProximityBand,

                                        char **papszOptions,

                                        GDALProgressFunc pfnProgress,

                                        void *pProgressArg)

{

    VALIDATE_POINTER1(hSrcBand, "GDALComputeProximity", CE_Failure);

    VALIDATE_POINTER1(hProximityBand, "GDALComputeProximity", CE_Failure);

    if (pfnProgress == nullptr)

        pfnProgress = GDALDummyProgress;

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

    /*      Are we using pixels or georeferenced coordinates for distances? */

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

    double dfDistMult = 1.0;

    const char *pszOpt = CSLFetchNameValue(papszOptions, "DISTUNITS");

    if (pszOpt)

    {

        if (EQUAL(pszOpt, "GEO"))

        {

            GDALDatasetH hSrcDS = GDALGetBandDataset(hSrcBand);

            if (hSrcDS)

            {

                double adfGeoTransform[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};

                GDALGetGeoTransform(hSrcDS, adfGeoTransform);

                if (std::abs(adfGeoTransform[1]) !=

                    std::abs(adfGeoTransform[5]))

                    CPLError(

                        CE_Warning, CPLE_AppDefined,

                        "Pixels not square, distances will be inaccurate.");

                dfDistMult = std::abs(adfGeoTransform[1]);

            }

        }

        else if (!EQUAL(pszOpt, "PIXEL"))

        {

            CPLError(

                CE_Failure, CPLE_AppDefined,

                "Unrecognized DISTUNITS value '%s', should be GEO or PIXEL.",

                pszOpt);

            return CE_Failure;

        }

    }

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

    /*      What is our maxdist value?                                      */

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

    pszOpt = CSLFetchNameValue(papszOptions, "MAXDIST");

    const double dfMaxDist = pszOpt ? CPLAtof(pszOpt) / dfDistMult

                                    : GDALGetRasterBandXSize(hSrcBand) +

                                          GDALGetRasterBandYSize(hSrcBand);

    CPLDebug("GDAL", "MAXDIST=%g, DISTMULT=%g", dfMaxDist, dfDistMult);

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

    /*      Verify the source and destination are compatible.               */

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

    const int nXSize = GDALGetRasterBandXSize(hSrcBand);

    const int nYSize = GDALGetRasterBandYSize(hSrcBand);

    if (nXSize != GDALGetRasterBandXSize(hProximityBand) ||

        nYSize != GDALGetRasterBandYSize(hProximityBand))

    {

        CPLError(CE_Failure, CPLE_AppDefined,

                 "Source and proximity bands are not the same size.");

        return CE_Failure;

    }

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

    /*      Get input NODATA value.                                         */

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

    double dfSrcNoDataValue = 0.0;

    double *pdfSrcNoData = nullptr;

    if (CPLFetchBool(papszOptions, "USE_INPUT_NODATA", false))

    {

        int bSrcHasNoData = 0;

        dfSrcNoDataValue = GDALGetRasterNoDataValue(hSrcBand, &bSrcHasNoData);

        if (bSrcHasNoData)

            pdfSrcNoData = &dfSrcNoDataValue;

    }

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

    /*      Get output NODATA value.                                        */

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

    float fNoDataValue = 0.0f;

    pszOpt = CSLFetchNameValue(papszOptions, "NODATA");

    if (pszOpt != nullptr)

    {

        fNoDataValue = static_cast<float>(CPLAtof(pszOpt));

    }

    else

    {

        int bSuccess = FALSE;

        fNoDataValue = static_cast<float>(

            GDALGetRasterNoDataValue(hProximityBand, &bSuccess));

        if (!bSuccess)

            fNoDataValue = 65535.0;

    }

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

    /*      Is there a fixed value we wish to force the buffer area to?     */

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

    double dfFixedBufVal = 0.0;

    bool bFixedBufVal = false;

    pszOpt = CSLFetchNameValue(papszOptions, "FIXED_BUF_VAL");

    if (pszOpt)

    {

        dfFixedBufVal = CPLAtof(pszOpt);

        bFixedBufVal = true;

    }

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

    /*      Get the target value(s).                                        */

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

    int *panTargetValues = nullptr;

    int nTargetValues = 0;

    pszOpt = CSLFetchNameValue(papszOptions, "VALUES");

    if (pszOpt != nullptr)

    {

        char **papszValuesTokens =

            CSLTokenizeStringComplex(pszOpt, ",", FALSE, FALSE);

        nTargetValues = CSLCount(papszValuesTokens);

        panTargetValues =

            static_cast<int *>(CPLCalloc(sizeof(int), nTargetValues));

        for (int i = 0; i < nTargetValues; i++)

            panTargetValues[i] = atoi(papszValuesTokens[i]);

        CSLDestroy(papszValuesTokens);

    }

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

    /*      Initialize progress counter.                                    */

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

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

    {

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

        CPLFree(panTargetValues);

        return CE_Failure;

    }

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

    /*      We need a signed type for the working proximity values kept     */

    /*      on disk.  If our proximity band is not signed, then create a    */

    /*      temporary file for this purpose.                                */

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

    GDALRasterBandH hWorkProximityBand = hProximityBand;

    GDALDatasetH hWorkProximityDS = nullptr;

    const GDALDataType eProxType = GDALGetRasterDataType(hProximityBand);

    CPLErr eErr = CE_None;

    // TODO(schwehr): Localize after removing gotos.

    float *pafProximity = nullptr;

    int *panNearX = nullptr;

    int *panNearY = nullptr;

    GInt32 *panSrcScanline = nullptr;

    bool bTempFileAlreadyDeleted = false;

    if (eProxType == GDT_UInt8 || eProxType == GDT_UInt16 ||

        eProxType == GDT_UInt32)

    {

        GDALDriverH hDriver = GDALGetDriverByName("GTiff");

        if (hDriver == nullptr)

        {

            CPLError(CE_Failure, CPLE_AppDefined,

                     "GDALComputeProximity needs GTiff driver");

            eErr = CE_Failure;

            goto end;

        }

        CPLString osTmpFile = CPLGenerateTempFilenameSafe("proximity");

        hWorkProximityDS = GDALCreate(hDriver, osTmpFile, nXSize, nYSize, 1,

                                      GDT_Float32, nullptr);

        if (hWorkProximityDS == nullptr)

        {

            eErr = CE_Failure;

            goto end;

        }

        // On Unix, attempt at deleting the temporary file now, so that

        // if the process gets interrupted, it is automatically destroyed

        // by the operating system.

        bTempFileAlreadyDeleted = VSIUnlink(osTmpFile) == 0;

        hWorkProximityBand = GDALGetRasterBand(hWorkProximityDS, 1);

    }

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

    /*      Allocate buffer for two scanlines of distances as floats        */

    /*      (the current and last line).                                    */

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

    pafProximity =

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

    panNearX = static_cast<int *>(VSI_MALLOC2_VERBOSE(sizeof(int), nXSize));

    panNearY = static_cast<int *>(VSI_MALLOC2_VERBOSE(sizeof(int), nXSize));

    panSrcScanline =

        static_cast<GInt32 *>(VSI_MALLOC2_VERBOSE(sizeof(GInt32), nXSize));

    if (pafProximity == nullptr || panNearX == nullptr || panNearY == nullptr ||

        panSrcScanline == nullptr)

    {

        eErr = CE_Failure;

        goto end;

    }

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

    /*      Loop from top to bottom of the image.                           */

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

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

    {

        panNearX[i] = -1;

        panNearY[i] = -1;

    }

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

    {

        // Read for target values.

        eErr = GDALRasterIO(hSrcBand, GF_Read, 0, iLine, nXSize, 1,

                            panSrcScanline, nXSize, 1, GDT_Int32, 0, 0);

        if (eErr != CE_None)

            break;

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

            pafProximity[i] = -1.0;

        // Left to right.

        ProcessProximityLine(panSrcScanline, panNearX, panNearY, TRUE, iLine,

                             nXSize, dfMaxDist, pafProximity, pdfSrcNoData,

                             nTargetValues, panTargetValues);

        // Right to Left.

        ProcessProximityLine(panSrcScanline, panNearX, panNearY, FALSE, iLine,

                             nXSize, dfMaxDist, pafProximity, pdfSrcNoData,

                             nTargetValues, panTargetValues);

        // Write out results.

        eErr = GDALRasterIO(hWorkProximityBand, GF_Write, 0, iLine, nXSize, 1,

                            pafProximity, nXSize, 1, GDT_Float32, 0, 0);

        if (eErr != CE_None)

            break;

        if (!pfnProgress(0.5 * (iLine + 1) / static_cast<double>(nYSize), "",

                         pProgressArg))

        {

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

            eErr = CE_Failure;

        }

    }

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

    /*      Loop from bottom to top of the image.                           */

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

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

    {

        panNearX[i] = -1;

        panNearY[i] = -1;

    }

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

    {

        // Read first pass proximity.

        eErr = GDALRasterIO(hWorkProximityBand, GF_Read, 0, iLine, nXSize, 1,

                            pafProximity, nXSize, 1, GDT_Float32, 0, 0);

        if (eErr != CE_None)

            break;

        // Read pixel values.

        eErr = GDALRasterIO(hSrcBand, GF_Read, 0, iLine, nXSize, 1,

                            panSrcScanline, nXSize, 1, GDT_Int32, 0, 0);

        if (eErr != CE_None)

            break;

        // Right to left.

        ProcessProximityLine(panSrcScanline, panNearX, panNearY, FALSE, iLine,

                             nXSize, dfMaxDist, pafProximity, pdfSrcNoData,

                             nTargetValues, panTargetValues);

        // Left to right.

        ProcessProximityLine(panSrcScanline, panNearX, panNearY, TRUE, iLine,

                             nXSize, dfMaxDist, pafProximity, pdfSrcNoData,

                             nTargetValues, panTargetValues);

        // Final post processing of distances.

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

        {

            if (pafProximity[i] < 0.0f)

                pafProximity[i] = fNoDataValue;

            else if (pafProximity[i] > 0.0f)

            {

                if (bFixedBufVal)

                    pafProximity[i] = static_cast<float>(dfFixedBufVal);

                else

                    pafProximity[i] =

                        pafProximity[i] * static_cast<float>(dfDistMult);

            }

        }

        // Write out results.

        eErr = GDALRasterIO(hProximityBand, GF_Write, 0, iLine, nXSize, 1,

                            pafProximity, nXSize, 1, GDT_Float32, 0, 0);

        if (eErr != CE_None)

            break;

        if (!pfnProgress(0.5 + 0.5 * (nYSize - iLine) /

                                   static_cast<double>(nYSize),

                         "", pProgressArg))

        {

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

            eErr = CE_Failure;

        }

    }

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

/*      Cleanup                                                         */

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

end:

    CPLFree(panNearX);

    CPLFree(panNearY);

    CPLFree(panSrcScanline);

    CPLFree(pafProximity);

    CPLFree(panTargetValues);

    if (hWorkProximityDS != nullptr)

    {

        CPLString osProxFile = GDALGetDescription(hWorkProximityDS);

        GDALClose(hWorkProximityDS);

        if (!bTempFileAlreadyDeleted)

        {

            GDALDeleteDataset(GDALGetDriverByName("GTiff"), osProxFile);

        }

    }

    return eErr;

}

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

/*                           SquareDistance()                           */

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

static double SquareDistance(double dfX1, double dfX2, double dfY1, double dfY2)

{

    const double dfDX = dfX1 - dfX2;

    const double dfDY = dfY1 - dfY2;

    return dfDX * dfDX + dfDY * dfDY;

}

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

/*                        ProcessProximityLine()                        */

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

static CPLErr ProcessProximityLine(GInt32 *panSrcScanline, int *panNearX,

                                   int *panNearY, int bForward, int iLine,

                                   int nXSize, double dfMaxDist,

                                   float *pafProximity,

                                   double *pdfSrcNoDataValue, int nTargetValues,

                                   int *panTargetValues)

{

    const int iStart = bForward ? 0 : nXSize - 1;

    const int iEnd = bForward ? nXSize : -1;

    const int iStep = bForward ? 1 : -1;

    for (int iPixel = iStart; iPixel != iEnd; iPixel += iStep)

    {

        bool bIsTarget = false;

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

         */

        /*      Is the current pixel a target pixel? */

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

         */

        if (nTargetValues == 0)

        {

            bIsTarget = panSrcScanline[iPixel] != 0;

        }

        else

        {

            for (int i = 0; i < nTargetValues; i++)

            {

                if (panSrcScanline[iPixel] == panTargetValues[i])

                    bIsTarget = TRUE;

            }

        }

        if (bIsTarget)

        {

            pafProximity[iPixel] = 0.0;

            panNearX[iPixel] = iPixel;

            panNearY[iPixel] = iLine;

            continue;

        }

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

         */

        /*      Are we near(er) to the closest target to the above (below) */

        /*      pixel? */

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

         */

        double dfNearDistSq = std::max(dfMaxDist, static_cast<double>(nXSize)) *

                              std::max(dfMaxDist, static_cast<double>(nXSize)) *

                              2.0;

        if (panNearX[iPixel] != -1)

        {

            const double dfDistSq = SquareDistance(panNearX[iPixel], iPixel,

                                                   panNearY[iPixel], iLine);

            if (dfDistSq < dfNearDistSq)

            {

                dfNearDistSq = dfDistSq;

            }

            else

            {

                panNearX[iPixel] = -1;

                panNearY[iPixel] = -1;

            }

        }

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

         */

        /*      Are we near(er) to the closest target to the left (right) */

        /*      pixel? */

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

         */

        const int iLast = iPixel - iStep;

        if (iPixel != iStart && panNearX[iLast] != -1)

        {

            const double dfDistSq =

                SquareDistance(panNearX[iLast], iPixel, panNearY[iLast], iLine);

            if (dfDistSq < dfNearDistSq)

            {

                dfNearDistSq = dfDistSq;

                panNearX[iPixel] = panNearX[iLast];

                panNearY[iPixel] = panNearY[iLast];

            }

        }

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

         */

        /*      Are we near(er) to the closest target to the topright */

        /*      (bottom left) pixel? */

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

         */

        const int iTR = iPixel + iStep;

        if (iTR != iEnd && panNearX[iTR] != -1)

        {

            const double dfDistSq =

                SquareDistance(panNearX[iTR], iPixel, panNearY[iTR], iLine);

            if (dfDistSq < dfNearDistSq)

            {

                dfNearDistSq = dfDistSq;

                panNearX[iPixel] = panNearX[iTR];

                panNearY[iPixel] = panNearY[iTR];

            }

        }

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

         */

        /*      Update our proximity value. */

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

         */

        if (panNearX[iPixel] != -1 &&

            (pdfSrcNoDataValue == nullptr ||

             panSrcScanline[iPixel] != *pdfSrcNoDataValue) &&

            dfNearDistSq <= dfMaxDist * dfMaxDist &&

            (pafProximity[iPixel] < 0 ||

             dfNearDistSq < static_cast<double>(pafProximity[iPixel]) *

                                static_cast<double>(pafProximity[iPixel])))

        {

            pafProximity[iPixel] = static_cast<float>(sqrt(dfNearDistSq));

        }

    }

    return CE_None;

}