/******************************************************************************
 *
 * Project:  Virtual GDAL Datasets
 * Purpose:  Implementation of some filter types.
 * Author:   Frank Warmerdam <warmerdam@pobox.com>
 *
 ******************************************************************************
 * Copyright (c) 2003, Frank Warmerdam <warmerdam@pobox.com>
 * Copyright (c) 2008-2013, Even Rouault <even dot rouault at spatialys.com>
 *
 * SPDX-License-Identifier: MIT
 ****************************************************************************/

#include "cpl_port.h"
#include "vrtdataset.h"

#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdlib>
#include <cstring>
#include <limits>
#include <map>

#include "cpl_conv.h"
#include "cpl_error.h"
#include "cpl_minixml.h"
#include "cpl_string.h"
#include "cpl_vsi.h"
#include "gdal.h"
#include "gdal_priv.h"

/*! @cond Doxygen_Suppress */

/************************************************************************/
/* ==================================================================== */
/*                          VRTFilteredSource                           */
/* ==================================================================== */
/************************************************************************/

/************************************************************************/
/*                         VRTFilteredSource()                          */
/************************************************************************/

VRTFilteredSource::VRTFilteredSource()
    : m_nSupportedTypesCount(1), m_nExtraEdgePixels(0)
{
    for (size_t i = 0; i < CPL_ARRAYSIZE(m_aeSupportedTypes); ++i)
        m_aeSupportedTypes[i] = GDT_Unknown;

    m_aeSupportedTypes[0] = GDT_Float32;
}

/************************************************************************/
/*                         ~VRTFilteredSource()                         */
/************************************************************************/

VRTFilteredSource::~VRTFilteredSource()
{
}

/************************************************************************/
/*                         SetExtraEdgePixels()                         */
/************************************************************************/

void VRTFilteredSource::SetExtraEdgePixels(int nEdgePixels)

{
    m_nExtraEdgePixels = nEdgePixels;
}

/************************************************************************/
/*                   SetFilteringDataTypesSupported()                   */
/************************************************************************/

void VRTFilteredSource::SetFilteringDataTypesSupported(int nTypeCount,
                                                       GDALDataType *paeTypes)

{
    if (nTypeCount >
        static_cast<int>(sizeof(m_aeSupportedTypes) / sizeof(GDALDataType)))
    {
        CPLAssert(false);
        nTypeCount =
            static_cast<int>(sizeof(m_aeSupportedTypes) / sizeof(GDALDataType));
    }

    m_nSupportedTypesCount = nTypeCount;
    memcpy(m_aeSupportedTypes, paeTypes, sizeof(GDALDataType) * nTypeCount);
}

/************************************************************************/
/*                          IsTypeSupported()                           */
/************************************************************************/

int VRTFilteredSource::IsTypeSupported(GDALDataType eTestType) const

{
    for (int i = 0; i < m_nSupportedTypesCount; i++)
    {
        if (eTestType == m_aeSupportedTypes[i])
            return TRUE;
    }

    return FALSE;
}

/************************************************************************/
/*                              RasterIO()                              */
/************************************************************************/

CPLErr VRTFilteredSource::RasterIO(GDALDataType eVRTBandDataType, int nXOff,
                                   int nYOff, int nXSize, int nYSize,
                                   void *pData, int nBufXSize, int nBufYSize,
                                   GDALDataType eBufType, GSpacing nPixelSpace,
                                   GSpacing nLineSpace,
                                   GDALRasterIOExtraArg *psExtraArg,
                                   WorkingState &oWorkingState)

{
    /* -------------------------------------------------------------------- */
    /*      For now we don't support filtered access to non-full            */
    /*      resolution requests. Just collect the data directly without     */
    /*      any operator.                                                   */
    /* -------------------------------------------------------------------- */
    if (nBufXSize != nXSize || nBufYSize != nYSize)
    {
        return VRTComplexSource::RasterIO(
            eVRTBandDataType, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize,
            nBufYSize, eBufType, nPixelSpace, nLineSpace, psExtraArg,
            oWorkingState);
    }

    double dfXOff = nXOff;
    double dfYOff = nYOff;
    double dfXSize = nXSize;
    double dfYSize = nYSize;
    if (psExtraArg != nullptr && psExtraArg->bFloatingPointWindowValidity)
    {
        dfXOff = psExtraArg->dfXOff;
        dfYOff = psExtraArg->dfYOff;
        dfXSize = psExtraArg->dfXSize;
        dfYSize = psExtraArg->dfYSize;
    }

    // The window we will actually request from the source raster band.
    double dfReqXOff = 0.0;
    double dfReqYOff = 0.0;
    double dfReqXSize = 0.0;
    double dfReqYSize = 0.0;
    int nReqXOff = 0;
    int nReqYOff = 0;
    int nReqXSize = 0;
    int nReqYSize = 0;

    // The window we will actual set _within_ the pData buffer.
    int nOutXOff = 0;
    int nOutYOff = 0;
    int nOutXSize = 0;
    int nOutYSize = 0;

    bool bError = false;
    if (!GetSrcDstWindow(dfXOff, dfYOff, dfXSize, dfYSize, nBufXSize, nBufYSize,
                         psExtraArg ? psExtraArg->eResampleAlg
                                    : GRIORA_NearestNeighbour,
                         &dfReqXOff, &dfReqYOff, &dfReqXSize, &dfReqYSize,
                         &nReqXOff, &nReqYOff, &nReqXSize, &nReqYSize,
                         &nOutXOff, &nOutYOff, &nOutXSize, &nOutYSize, bError))
    {
        return bError ? CE_Failure : CE_None;
    }

    pData = reinterpret_cast<GByte *>(pData) + nPixelSpace * nOutXOff +
            nLineSpace * nOutYOff;

    /* -------------------------------------------------------------------- */
    /*      Determine the data type we want to request.  We try to match    */
    /*      the source or destination request, and if both those fail we    */
    /*      fallback to the first supported type at least as expressive     */
    /*      as the request.                                                 */
    /* -------------------------------------------------------------------- */
    GDALDataType eOperDataType = GDT_Unknown;

    if (IsTypeSupported(eBufType))
        eOperDataType = eBufType;

    auto l_band = GetRasterBand();
    if (!l_band)
        return CE_Failure;

    if (eOperDataType == GDT_Unknown &&
        IsTypeSupported(l_band->GetRasterDataType()))
        eOperDataType = l_band->GetRasterDataType();

    if (eOperDataType == GDT_Unknown)
    {
        for (int i = 0; i < m_nSupportedTypesCount; i++)
        {
            if (GDALDataTypeUnion(m_aeSupportedTypes[i], eBufType) ==
                m_aeSupportedTypes[i])
            {
                eOperDataType = m_aeSupportedTypes[i];
            }
        }
    }

    if (eOperDataType == GDT_Unknown)
    {
        eOperDataType = m_aeSupportedTypes[0];

        for (int i = 1; i < m_nSupportedTypesCount; i++)
        {
            if (GDALGetDataTypeSizeBytes(m_aeSupportedTypes[i]) >
                GDALGetDataTypeSizeBytes(eOperDataType))
            {
                eOperDataType = m_aeSupportedTypes[i];
            }
        }
    }

    /* -------------------------------------------------------------------- */
    /*      Allocate the buffer of data into which our imagery will be      */
    /*      read, with the extra edge pixels as well. This will be the      */
    /*      source data fed into the filter.                                */
    /* -------------------------------------------------------------------- */
    if (nOutXSize > INT_MAX - 2 * m_nExtraEdgePixels ||
        nOutYSize > INT_MAX - 2 * m_nExtraEdgePixels)
    {
        return CE_Failure;
    }
    const int nExtraXSize = nOutXSize + 2 * m_nExtraEdgePixels;
    const int nExtraYSize = nOutYSize + 2 * m_nExtraEdgePixels;

    GByte *pabyWorkData = static_cast<GByte *>(VSI_MALLOC3_VERBOSE(
        nExtraXSize, nExtraYSize, GDALGetDataTypeSizeBytes(eOperDataType)));

    if (pabyWorkData == nullptr)
    {
        return CE_Failure;
    }

    const GPtrDiff_t nPixelOffset = GDALGetDataTypeSizeBytes(eOperDataType);
    const GPtrDiff_t nLineOffset = nPixelOffset * nExtraXSize;

    int bHasNoData = false;
    const double dfSrcNoDataValue = l_band->GetNoDataValue(&bHasNoData);
    if (bHasNoData)
    {
        GDALCopyWords64(&dfSrcNoDataValue, GDT_Float64, 0, pabyWorkData,
                        eOperDataType, static_cast<int>(nPixelOffset),
                        static_cast<size_t>(nExtraXSize) * nExtraYSize);
    }
    else
    {
        memset(pabyWorkData, 0, nLineOffset * nExtraYSize);
    }

    /* -------------------------------------------------------------------- */
    /*      Allocate the output buffer in the same dimensions as the work   */
    /*      buffer. This allows the filter process to write edge pixels     */
    /*      if needed for two-pass (separable) filtering.                   */
    /* -------------------------------------------------------------------- */
    GByte *pabyOutData = static_cast<GByte *>(
        VSI_MALLOC3_VERBOSE(nExtraXSize, nExtraYSize, nPixelOffset));
    if (pabyOutData == nullptr)
    {
        CPLFree(pabyWorkData);
        return CE_Failure;
    }

    /* -------------------------------------------------------------------- */
    /*      Figure out the extended window that we want to load.  Note      */
    /*      that we keep track of the file window as well as the amount     */
    /*      we will need to edge fill past the edge of the source dataset.  */
    /* -------------------------------------------------------------------- */
    int nFileXOff = nReqXOff - m_nExtraEdgePixels;
    int nFileYOff = nReqYOff - m_nExtraEdgePixels;
    int nFileXSize = nExtraXSize;
    int nFileYSize = nExtraYSize;

    int nTopFill = 0;
    int nLeftFill = 0;
    int nRightFill = 0;
    int nBottomFill = 0;

    if (nFileXOff < 0)
    {
        nLeftFill = -nFileXOff;
        nFileXOff = 0;
        nFileXSize -= nLeftFill;
    }

    if (nFileYOff < 0)
    {
        nTopFill = -nFileYOff;
        nFileYOff = 0;
        nFileYSize -= nTopFill;
    }

    if (nFileXOff + nFileXSize > l_band->GetXSize())
    {
        nRightFill = nFileXOff + nFileXSize - l_band->GetXSize();
        nFileXSize -= nRightFill;
    }

    if (nFileYOff + nFileYSize > l_band->GetYSize())
    {
        nBottomFill = nFileYOff + nFileYSize - l_band->GetYSize();
        nFileYSize -= nBottomFill;
    }

    /* -------------------------------------------------------------------- */
    /*      Load the data.                                                  */
    /* -------------------------------------------------------------------- */
    {
        GDALRasterIOExtraArg sExtraArgs;
        INIT_RASTERIO_EXTRA_ARG(sExtraArgs);
        const bool bIsComplex =
            CPL_TO_BOOL(GDALDataTypeIsComplex(eOperDataType));
        const CPLErr eErr = VRTComplexSource::RasterIOInternal<float>(
            l_band, eVRTBandDataType, nFileXOff, nFileYOff, nFileXSize,
            nFileYSize,
            pabyWorkData + nLineOffset * nTopFill + nPixelOffset * nLeftFill,
            nFileXSize, nFileYSize, eOperDataType, nPixelOffset, nLineOffset,
            &sExtraArgs, bIsComplex ? GDT_CFloat32 : GDT_Float32,
            oWorkingState);

        if (eErr != CE_None)
        {
            VSIFree(pabyWorkData);
            VSIFree(pabyOutData);

            return eErr;
        }
    }

    /* -------------------------------------------------------------------- */
    /*      Fill in missing areas.  Note that we replicate the edge         */
    /*      valid values out.  We don't using "mirroring" which might be    */
    /*      more suitable for some times of filters.  We also don't mark    */
    /*      these pixels as "nodata" though perhaps we should.              */
    /* -------------------------------------------------------------------- */
    if (nLeftFill != 0 || nRightFill != 0)
    {
        for (int i = nTopFill; i < nExtraYSize - nBottomFill; i++)
        {
            if (nLeftFill != 0)
                GDALCopyWords(
                    pabyWorkData + nPixelOffset * nLeftFill + i * nLineOffset,
                    eOperDataType, 0, pabyWorkData + i * nLineOffset,
                    eOperDataType, static_cast<int>(nPixelOffset), nLeftFill);

            if (nRightFill != 0)
                GDALCopyWords(pabyWorkData + i * nLineOffset +
                                  nPixelOffset * (nExtraXSize - nRightFill - 1),
                              eOperDataType, 0,
                              pabyWorkData + i * nLineOffset +
                                  nPixelOffset * (nExtraXSize - nRightFill),
                              eOperDataType, static_cast<int>(nPixelOffset),
                              nRightFill);
        }
    }

    for (int i = 0; i < nTopFill; i++)
    {
        memcpy(pabyWorkData + i * nLineOffset,
               pabyWorkData + nTopFill * nLineOffset, nLineOffset);
    }

    for (int i = nExtraYSize - nBottomFill; i < nExtraYSize; i++)
    {
        memcpy(pabyWorkData + i * nLineOffset,
               pabyWorkData + (nExtraYSize - nBottomFill - 1) * nLineOffset,
               nLineOffset);
    }

    /* -------------------------------------------------------------------- */
    /*      Filter the data.                                                */
    /* -------------------------------------------------------------------- */
    const CPLErr eErr = FilterData(nExtraXSize, nExtraYSize, eOperDataType,
                                   pabyWorkData, pabyOutData);

    VSIFree(pabyWorkData);
    if (eErr != CE_None)
    {
        VSIFree(pabyOutData);

        return eErr;
    }

    /* -------------------------------------------------------------------- */
    /*      Copy from work buffer to target buffer.                         */
    /* -------------------------------------------------------------------- */
    GByte *pabySrcRow =
        pabyOutData + (nLineOffset + nPixelOffset) * m_nExtraEdgePixels;
    GByte *pabyDstRow = reinterpret_cast<GByte *>(pData);

    for (int i = 0; i < nOutYSize;
         i++, pabySrcRow += nLineOffset, pabyDstRow += nLineSpace)
    {
        GDALCopyWords(pabySrcRow, eOperDataType, static_cast<int>(nPixelOffset),
                      pabyDstRow, eBufType, static_cast<int>(nPixelSpace),
                      nOutXSize);
    }

    VSIFree(pabyOutData);

    return CE_None;
}

/************************************************************************/
/* ==================================================================== */
/*                       VRTKernelFilteredSource                        */
/* ==================================================================== */
/************************************************************************/

/************************************************************************/
/*                      VRTKernelFilteredSource()                       */
/************************************************************************/

VRTKernelFilteredSource::VRTKernelFilteredSource()
{
    GDALDataType aeSupTypes[] = {GDT_Float32};
    SetFilteringDataTypesSupported(1, aeSupTypes);
}

/************************************************************************/
/*                            GetType()                                 */
/************************************************************************/

const char *VRTKernelFilteredSource::GetType() const
{
    static const char *TYPE = "KernelFilteredSource";
    return TYPE;
}

/************************************************************************/
/*                           SetNormalized()                            */
/************************************************************************/

void VRTKernelFilteredSource::SetNormalized(bool bNormalizedIn)

{
    m_bNormalized = bNormalizedIn;
}

/************************************************************************/
/*                             SetKernel()                              */
/************************************************************************/

CPLErr
VRTKernelFilteredSource::SetKernel(int nNewKernelSize, bool bSeparable,
                                   const std::vector<double> &adfNewCoefs)

{
    if (nNewKernelSize < 1 || (nNewKernelSize % 2) != 1)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Illegal filtering kernel size %d, "
                 "must be odd positive number.",
                 nNewKernelSize);
        return CE_Failure;
    }
    if (adfNewCoefs.size() !=
        static_cast<size_t>(nNewKernelSize) * (bSeparable ? 1 : nNewKernelSize))
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "adfNewCoefs[] is not of expected size");
        return CE_Failure;
    }

    m_nKernelSize = nNewKernelSize;
    m_bSeparable = bSeparable;
    m_adfKernelCoefs = adfNewCoefs;

    SetExtraEdgePixels((nNewKernelSize - 1) / 2);

    return CE_None;
}

/************************************************************************/
/*                             FilterData()                             */
/************************************************************************/

CPLErr VRTKernelFilteredSource::FilterData(int nXSize, int nYSize,
                                           GDALDataType eType,
                                           GByte *pabySrcData,
                                           GByte *pabyDstData)

{
    /* -------------------------------------------------------------------- */
    /*      Validate data type.                                             */
    /* -------------------------------------------------------------------- */
    if (eType != GDT_Float32)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Unsupported data type (%s) in "
                 "VRTKernelFilteredSource::FilterData()",
                 GDALGetDataTypeName(eType));
        return CE_Failure;
    }

    CPLAssert(m_nExtraEdgePixels * 2 + 1 == m_nKernelSize ||
              (m_nKernelSize == 0 && m_nExtraEdgePixels == 0));

    const bool bMin = m_function == "min";
    const bool bMax = m_function == "max";
    const bool bStdDev = m_function == "stddev";
    const bool bMedian = m_function == "median";
    const bool bMode = m_function == "mode";

    /* -------------------------------------------------------------------- */
    /*      Float32 case.                                                   */
    /* -------------------------------------------------------------------- */
    // if( eType == GDT_Float32 )
    {
        float *pafSrcData = reinterpret_cast<float *>(pabySrcData);
        float *pafDstData = reinterpret_cast<float *>(pabyDstData);

        int bHasNoData = FALSE;
        auto l_poBand = GetRasterBand();
        if (!l_poBand)
            return CE_Failure;
        const float fNoData =
            static_cast<float>(l_poBand->GetNoDataValue(&bHasNoData));

        const int nAxisCount = m_bSeparable ? 2 : 1;

        for (int nAxis = 0; nAxis < nAxisCount; ++nAxis)
        {
            const int nISize = nAxis == 0 ? nYSize : nXSize;
            const int nJSize = nAxis == 0 ? nXSize : nYSize;
            const int nIStride = nAxis == 0 ? nXSize : 1;
            const int nJStride = nAxis == 0 ? 1 : nXSize;

            const int nIMin = m_nExtraEdgePixels;
            const int nIMax = nISize - m_nExtraEdgePixels;
            const int nJMin = (m_bSeparable ? 0 : m_nExtraEdgePixels);
            const int nJMax = nJSize - (m_bSeparable ? 0 : m_nExtraEdgePixels);

            for (GPtrDiff_t iJ = nJMin; iJ < nJMax; ++iJ)
            {
                if (nAxis == 1)
                    memcpy(pafSrcData + iJ * nJStride,
                           pafDstData + iJ * nJStride, sizeof(float) * nXSize);

                for (int iI = nIMin; iI < nIMax; ++iI)
                {
                    const GPtrDiff_t iIndex =
                        static_cast<GPtrDiff_t>(iI) * nIStride + iJ * nJStride;

                    if (bHasNoData && pafSrcData[iIndex] == fNoData)
                    {
                        pafDstData[iIndex] = fNoData;
                        continue;
                    }

                    double dfSum = 0.0, dfKernSum = 0.0;
                    size_t nValidCount = 0;
                    double dfRes =
                        bMin   ? std::numeric_limits<double>::infinity()
                        : bMax ? -std::numeric_limits<double>::infinity()
                               : 0.0;
                    double dfMean = 0.0;
                    double dfM2 = 0.0;
                    std::vector<double> adfVals;
                    std::map<double, size_t> mapValToCount;
                    size_t maxCount = 0;

                    for (GPtrDiff_t iII = -m_nExtraEdgePixels, iK = 0;
                         iII <= m_nExtraEdgePixels; ++iII)
                    {
                        for (GPtrDiff_t iJJ =
                                 (m_bSeparable ? 0 : -m_nExtraEdgePixels);
                             iJJ <= (m_bSeparable ? 0 : m_nExtraEdgePixels);
                             ++iJJ, ++iK)
                        {
                            const float *pfData = pafSrcData + iIndex +
                                                  iII * nIStride +
                                                  iJJ * nJStride;
                            if (bHasNoData &&
                                (*pfData == fNoData || std::isnan(*pfData)))
                            {
                                continue;
                            }
                            if (m_adfKernelCoefs[iK] == 0.0)
                            {
                                continue;
                            }
                            const double dfVal = static_cast<double>(*pfData) *
                                                 m_adfKernelCoefs[iK];
                            ++nValidCount;

                            if (bMax)
                            {
                                if (dfVal > dfRes)
                                    dfRes = dfVal;
                            }
                            else if (bMin)
                            {
                                if (dfVal < dfRes)
                                    dfRes = dfVal;
                            }
                            else if (bStdDev)
                            {
                                const double dfDelta = dfVal - dfMean;
                                dfMean += dfDelta / nValidCount;
                                dfM2 += dfDelta * (dfVal - dfMean);
                            }
                            else if (bMedian)
                            {
                                adfVals.push_back(dfVal);
                            }
                            else if (bMode)
                            {
                                const size_t nCountVal = ++mapValToCount[dfVal];
                                if (nCountVal > maxCount)
                                {
                                    maxCount = nCountVal;
                                    dfRes = dfVal;
                                }
                            }
                            else
                            {
                                dfSum += dfVal;
                                dfKernSum += m_adfKernelCoefs[iK];
                            }
                        }
                    }

                    float fResult;
                    if (bMax || bMin || bMode)
                    {
                        fResult = nValidCount  ? static_cast<float>(dfRes)
                                  : bHasNoData ? fNoData
                                               : 0.0f;
                    }
                    else if (bStdDev)
                    {
                        fResult =
                            nValidCount
                                ? static_cast<float>(sqrt(
                                      dfM2 / static_cast<double>(nValidCount)))
                            : bHasNoData ? fNoData
                                         : 0.0f;
                    }
                    else if (bMedian)
                    {
                        if (!adfVals.empty())
                        {
                            if ((adfVals.size() % 2) == 1)
                            {
                                std::nth_element(adfVals.begin(),
                                                 adfVals.begin() +
                                                     adfVals.size() / 2,
                                                 adfVals.end());
                                dfRes = adfVals[adfVals.size() / 2];
                            }
                            else
                            {
                                std::nth_element(adfVals.begin(),
                                                 adfVals.begin() +
                                                     adfVals.size() / 2 - 1,
                                                 adfVals.end());
                                dfRes = adfVals[adfVals.size() / 2 - 1];
                                std::nth_element(adfVals.begin(),
                                                 adfVals.begin() +
                                                     adfVals.size() / 2,
                                                 adfVals.end());
                                dfRes =
                                    (dfRes + adfVals[adfVals.size() / 2]) / 2;
                            }
                            fResult = static_cast<float>(dfRes);
                        }
                        else
                            fResult = bHasNoData ? fNoData : 0.0f;
                    }
                    else if (!m_bNormalized)
                        fResult = static_cast<float>(dfSum);
                    else if (nValidCount == 0 || dfKernSum == 0.0)
                        fResult = bHasNoData ? fNoData : 0.0f;
                    else
                        fResult = static_cast<float>(dfSum / dfKernSum);
                    pafDstData[iIndex] = fResult;
                }
            }
        }
    }

    return CE_None;
}

/************************************************************************/
/*                              XMLInit()                               */
/************************************************************************/

CPLErr
VRTKernelFilteredSource::XMLInit(const CPLXMLNode *psTree,
                                 const char *pszVRTPath,
                                 VRTMapSharedResources &oMapSharedSources)

{
    {
        const CPLErr eErr =
            VRTFilteredSource::XMLInit(psTree, pszVRTPath, oMapSharedSources);
        if (eErr != CE_None)
            return eErr;
    }

    const int nNewKernelSize = atoi(CPLGetXMLValue(psTree, "Kernel.Size", "0"));

    if (nNewKernelSize == 0)
        return CE_None;
    // To prevent a integer overflow when computing nNewKernelSize *
    // nNewKernelSize
    if (nNewKernelSize < 0 ||
        nNewKernelSize > static_cast<int>(std::sqrt(static_cast<double>(
                             std::numeric_limits<int>::max()))))
    {
        CPLError(CE_Failure, CPLE_IllegalArg,
                 "Invalid value for kernel size: %d", nNewKernelSize);
        return CE_Failure;
    }

    const CPLStringList aosCoefItems(
        CSLTokenizeString(CPLGetXMLValue(psTree, "Kernel.Coefs", "")));

    const int nCoefs = aosCoefItems.size();

    const bool bSquare = nCoefs == nNewKernelSize * nNewKernelSize;
    const bool bSeparable = nCoefs == nNewKernelSize && nCoefs != 1;

    if (!bSquare && !bSeparable)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Got wrong number of filter kernel coefficients (%s).  "
                 "Expected %d or %d, got %d.",
                 CPLGetXMLValue(psTree, "Kernel.Coefs", ""),
                 nNewKernelSize * nNewKernelSize, nNewKernelSize, nCoefs);
        return CE_Failure;
    }

    std::vector<double> adfNewCoefs;
    adfNewCoefs.reserve(nCoefs);
    for (int i = 0; i < nCoefs; i++)
        adfNewCoefs.push_back(CPLAtof(aosCoefItems[i]));

    const CPLErr eErr = SetKernel(nNewKernelSize, bSeparable, adfNewCoefs);
    if (eErr == CE_None)
    {
        SetNormalized(atoi(CPLGetXMLValue(psTree, "Kernel.normalized", "0")) !=
                      0);
    }

    const char *pszFunction = CPLGetXMLValue(psTree, "Function", nullptr);
    if (pszFunction)
    {
        if (!EQUAL(pszFunction, "max") && !EQUAL(pszFunction, "min") &&
            !EQUAL(pszFunction, "median") && !EQUAL(pszFunction, "mode") &&
            !EQUAL(pszFunction, "stddev"))
        {
            CPLError(CE_Failure, CPLE_IllegalArg, "Unsupported function: %s",
                     pszFunction);
            return CE_Failure;
        }
        SetFunction(pszFunction);
    }

    return eErr;
}

/************************************************************************/
/*                           SerializeToXML()                           */
/************************************************************************/

CPLXMLNode *VRTKernelFilteredSource::SerializeToXML(const char *pszVRTPath)

{
    CPLXMLNode *psSrc = VRTFilteredSource::SerializeToXML(pszVRTPath);

    if (psSrc == nullptr)
        return nullptr;

    CPLFree(psSrc->pszValue);
    psSrc->pszValue = CPLStrdup("KernelFilteredSource");

    if (m_nKernelSize == 0)
        return psSrc;

    CPLXMLNode *psKernel = CPLCreateXMLNode(psSrc, CXT_Element, "Kernel");

    CPLCreateXMLNode(CPLCreateXMLNode(psKernel, CXT_Attribute, "normalized"),
                     CXT_Text, m_bNormalized ? "1" : "0");

    std::string osCoefs;
    for (auto dfVal : m_adfKernelCoefs)
    {
        if (!osCoefs.empty())
            osCoefs += ' ';
        osCoefs += CPLSPrintf("%.8g", dfVal);
    }

    CPLSetXMLValue(psKernel, "Size", CPLSPrintf("%d", m_nKernelSize));
    CPLSetXMLValue(psKernel, "Coefs", osCoefs.c_str());

    if (!m_function.empty())
        CPLCreateXMLElementAndValue(psSrc, "Function", m_function.c_str());

    return psSrc;
}

/************************************************************************/
/*                       VRTParseFilterSources()                        */
/************************************************************************/

VRTSource *VRTParseFilterSources(const CPLXMLNode *psChild,
                                 const char *pszVRTPath,
                                 VRTMapSharedResources &oMapSharedSources)

{
    if (EQUAL(psChild->pszValue, "KernelFilteredSource"))
    {
        VRTSource *poSrc = new VRTKernelFilteredSource();
        if (poSrc->XMLInit(psChild, pszVRTPath, oMapSharedSources) == CE_None)
            return poSrc;

        delete poSrc;
    }

    return nullptr;
}

/*! @endcond */
