/src/gdal/apps/gdalalg_raster_neighbors.cpp

Source
/******************************************************************************
 *
 * Project:  GDAL
 * Purpose:  "neighbors" step of "raster pipeline"
 * Author:   Even Rouault <even dot rouault at spatialys.com>
 *
 ******************************************************************************
 * Copyright (c) 2025, Even Rouault <even dot rouault at spatialys.com>
 *
 * SPDX-License-Identifier: MIT
 ****************************************************************************/

#include "gdalalg_raster_neighbors.h"

#include "gdal_priv.h"
#include "gdal_priv_templates.hpp"
#include "vrtdataset.h"

#include <algorithm>
#include <limits>
#include <map>
#include <optional>
#include <set>
#include <utility>
#include <vector>

//! @cond Doxygen_Suppress

#ifndef _
#define _(x) (x)
#endif

static const std::set<std::string> oSetKernelNames = {
    "u",     "v",       "equal",    "edge1",
    "edge2", "sharpen", "gaussian", "unsharp-masking"};

namespace
{
struct KernelDef
{
    int size = 0;
    std::vector<double> adfCoefficients{};
};
}  // namespace

// clang-format off
// Cf https://en.wikipedia.org/wiki/Kernel_(image_processing)
static const std::map<std::string, std::pair<int, std::vector<int>>> oMapKernelNameToMatrix = {
    { "u", { 3, {  0,  0,  0,
                  -1,  0,  1,
                   0,  0,  0 } } },
    { "v", { 3, {  0, -1,  0,
                   0,  0,  0,
                   0,  1,  0 } } },
    { "edge1", { 3, {  0, -1,  0,
                      -1,  4, -1,
                       0, -1,  0 } } },
    { "edge2", { 3, { -1, -1, -1,
                      -1,  8, -1,
                      -1, -1, -1 } } },
    { "sharpen",  { 3, { 0, -1,  0,
                        -1,  5, -1,
                         0, -1,  0 } } },
    { "gaussian-3x3", { 3, { 1, 2, 1,
                                  2, 4, 2,
                                  1, 2, 1 } } },
    { "gaussian-5x5", { 5, { 1, 4,   6,  4, 1,
                                  4, 16, 24, 16, 4,
                                  6, 24, 36, 24, 6,
                                  4, 16, 24, 16, 4,
                                  1, 4,   6,  4, 1, } } },
    { "unsharp-masking-5x5", { 5, { 1, 4,     6,  4, 1,
                                    4, 16,   24, 16, 4,
                                    6, 24, -476, 24, 6,
                                    4, 16,   24, 16, 4,
                                    1, 4,     6,  4, 1, } } },
};

// clang-format on

/************************************************************************/
/*                        CreateDerivedBandXML()                        */
/************************************************************************/

static bool CreateDerivedBandXML(VRTDataset *poVRTDS, GDALRasterBand *poSrcBand,
                                 GDALDataType eType, const std::string &noData,
                                 const std::string &method,
                                 const KernelDef &kernelDef)
{
    poVRTDS->AddBand(eType, nullptr);

    std::optional<double> dstNoData;
    bool autoSelectNoDataValue = false;
    if (noData.empty())
    {
        autoSelectNoDataValue = true;
    }
    else if (noData != "none")
    {
        // Already validated to be numeric by the validation action
        dstNoData = CPLAtof(noData.c_str());
    }

    auto poVRTBand = cpl::down_cast<VRTSourcedRasterBand *>(
        poVRTDS->GetRasterBand(poVRTDS->GetRasterCount()));

    auto poSource = std::make_unique<VRTKernelFilteredSource>();
    poSrcBand->GetDataset()->Reference();
    poSource->SetSrcBand(poSrcBand);
    poSource->SetKernel(kernelDef.size, /* separable = */ false,
                        kernelDef.adfCoefficients);
    poSource->SetNormalized(method != "sum");
    if (method != "sum" && method != "mean")
        poSource->SetFunction(method.c_str());

    int bSrcHasNoData = false;
    const double dfNoDataValue = poSrcBand->GetNoDataValue(&bSrcHasNoData);
    if (bSrcHasNoData)
    {
        poSource->SetNoDataValue(dfNoDataValue);
        if (autoSelectNoDataValue && !dstNoData.has_value())
        {
            dstNoData = dfNoDataValue;
        }
    }

    if (dstNoData.has_value())
    {
        if (!GDALIsValueExactAs(dstNoData.value(), eType))
        {
            CPLError(CE_Failure, CPLE_AppDefined,
                     "Band output type %s cannot represent NoData value %g",
                     GDALGetDataTypeName(eType), dstNoData.value());
            return false;
        }

        poVRTBand->SetNoDataValue(dstNoData.value());
    }

    poVRTBand->AddSource(std::move(poSource));

    return true;
}

/************************************************************************/
/*                   GDALNeighborsCreateVRTDerived()                    */
/************************************************************************/

static std::unique_ptr<GDALDataset>
GDALNeighborsCreateVRTDerived(GDALDataset *poSrcDS, int nBand,
                              GDALDataType eType, const std::string &noData,
                              const std::vector<std::string> &methods,
                              const std::vector<KernelDef> &aKernelDefs)
{
    CPLAssert(methods.size() == aKernelDefs.size());

    auto ds = std::make_unique<VRTDataset>(poSrcDS->GetRasterXSize(),
                                           poSrcDS->GetRasterYSize());
    GDALGeoTransform gt;
    if (poSrcDS->GetGeoTransform(gt) == CE_None)
        ds->SetGeoTransform(gt);
    if (const OGRSpatialReference *poSRS = poSrcDS->GetSpatialRef())
    {
        ds->SetSpatialRef(poSRS);
    }

    bool ret = true;
    if (nBand != 0)
    {
        for (size_t i = 0; i < aKernelDefs.size() && ret; ++i)
        {
            ret =
                CreateDerivedBandXML(ds.get(), poSrcDS->GetRasterBand(nBand),
                                     eType, noData, methods[i], aKernelDefs[i]);
        }
    }
    else
    {
        for (int iBand = 1; iBand <= poSrcDS->GetRasterCount(); ++iBand)
        {
            for (size_t i = 0; i < aKernelDefs.size() && ret; ++i)
            {
                ret = CreateDerivedBandXML(ds.get(),
                                           poSrcDS->GetRasterBand(iBand), eType,
                                           noData, methods[i], aKernelDefs[i]);
            }
        }
    }
    if (!ret)
        ds.reset();
    return ds;
}

/************************************************************************/
/*     GDALRasterNeighborsAlgorithm::GDALRasterNeighborsAlgorithm()     */
/************************************************************************/

GDALRasterNeighborsAlgorithm::GDALRasterNeighborsAlgorithm(
    bool standaloneStep) noexcept
    : GDALRasterPipelineStepAlgorithm(
          NAME, DESCRIPTION, HELP_URL,
          ConstructorOptions().SetStandaloneStep(standaloneStep))
{
    AddBandArg(&m_band);

    AddArg("method", 0, _("Method to combine weighed source pixels"), &m_method)
        .SetChoices("mean", "sum", "min", "max", "stddev", "median", "mode");

    AddArg("size", 0, _("Neighborhood size"), &m_size)
        .SetMinValueIncluded(3)
        .SetMaxValueIncluded(99)
        .AddValidationAction(
            [this]()
            {
                if ((m_size % 2) != 1)
                {
                    ReportError(CE_Failure, CPLE_IllegalArg,
                                "The value of 'size' must be an odd number.");
                    return false;
                }
                return true;
            });

    AddArg("kernel", 0, _("Convolution kernel(s) to apply"), &m_kernel)
        .SetPackedValuesAllowed(false)
        .SetMinCount(1)
        .SetMinCharCount(1)
        .SetRequired()
        .SetAutoCompleteFunction(
            [](const std::string &v)
            {
                std::vector<std::string> ret;
                if (v.empty() || v.front() != '[')
                {
                    ret.insert(ret.end(), oSetKernelNames.begin(),
                               oSetKernelNames.end());
                    ret.push_back(
                        "[[val00,val10,...,valN0],...,[val0N,val1N,...valNN]]");
                }
                return ret;
            })
        .AddValidationAction(
            [this]()
            {
                for (const std::string &kernel : m_kernel)
                {
                    if (kernel.front() == '[' && kernel.back() == ']')
                    {
                        const CPLStringList aosValues(CSLTokenizeString2(
                            kernel.c_str(), "[] ,",
                            CSLT_STRIPLEADSPACES | CSLT_STRIPENDSPACES));
                        const double dfSize =
                            static_cast<double>(aosValues.size());
                        const int nSqrt = static_cast<int>(sqrt(dfSize) + 0.5);
                        if (!((aosValues.size() % 2) == 1 &&
                              nSqrt * nSqrt == aosValues.size()))
                        {
                            ReportError(
                                CE_Failure, CPLE_IllegalArg,
                                "The number of values in the 'kernel' "
                                "argument must be an odd square number.");
                            return false;
                        }
                        for (int i = 0; i < aosValues.size(); ++i)
                        {
                            if (CPLGetValueType(aosValues[i]) ==
                                CPL_VALUE_STRING)
                            {
                                ReportError(CE_Failure, CPLE_IllegalArg,
                                            "Non-numeric value found in the "
                                            "'kernel' argument");
                                return false;
                            }
                        }
                    }
                    else if (!cpl::contains(oSetKernelNames, kernel))
                    {
                        std::string osMsg =
                            "Valid values for 'kernel' argument are: ";
                        bool bFirst = true;
                        for (const auto &name : oSetKernelNames)
                        {
                            if (!bFirst)
                                osMsg += ", ";
                            bFirst = false;
                            osMsg += '\'';
                            osMsg += name;
                            osMsg += '\'';
                        }
                        osMsg += " or "
                                 "[[val00,val10,...,valN0],...,[val0N,val1N,..."
                                 "valNN]]";
                        ReportError(CE_Failure, CPLE_IllegalArg, "%s",
                                    osMsg.c_str());
                        return false;
                    }
                }
                return true;
            });

    AddOutputDataTypeArg(&m_type).SetDefault("Float64");

    AddNodataArg(&m_nodata, true);

    AddValidationAction(
        [this]()
        {
            if (m_method.size() > 1 && m_method.size() != m_kernel.size())
            {
                ReportError(
                    CE_Failure, CPLE_AppDefined,
                    "The number of values for the 'method' argument should "
                    "be one or exactly the number of values of 'kernel'");
                return false;
            }

            if (m_size > 0)
            {
                for (const std::string &kernel : m_kernel)
                {
                    if (kernel == "gaussian")
                    {
                        if (m_size != 3 && m_size != 5)
                        {
                            ReportError(CE_Failure, CPLE_AppDefined,
                                        "Currently only size = 3 or 5 is "
                                        "supported for kernel '%s'",
                                        kernel.c_str());
                            return false;
                        }
                    }
                    else if (kernel == "unsharp-masking")
                    {
                        if (m_size != 5)
                        {
                            ReportError(CE_Failure, CPLE_AppDefined,
                                        "Currently only size = 5 is supported "
                                        "for kernel '%s'",
                                        kernel.c_str());
                            return false;
                        }
                    }
                    else if (kernel[0] == '[')
                    {
                        const CPLStringList aosValues(CSLTokenizeString2(
                            kernel.c_str(), "[] ,",
                            CSLT_STRIPLEADSPACES | CSLT_STRIPENDSPACES));
                        const double dfSize =
                            static_cast<double>(aosValues.size());
                        const int size =
                            static_cast<int>(std::floor(sqrt(dfSize) + 0.5));
                        if (m_size != size)
                        {
                            ReportError(CE_Failure, CPLE_AppDefined,
                                        "Value of 'size' argument (%d) "
                                        "inconsistent with the one deduced "
                                        "from the kernel matrix (%d)",
                                        m_size, size);
                            return false;
                        }
                    }
                    else if (m_size != 3 && kernel != "equal" &&
                             kernel[0] != '[')
                    {
                        ReportError(CE_Failure, CPLE_AppDefined,
                                    "Currently only size = 3 is supported for "
                                    "kernel '%s'",
                                    kernel.c_str());
                        return false;
                    }
                }
            }

            return true;
        });
}

/************************************************************************/
/*                            GetKernelDef()                            */
/************************************************************************/

static KernelDef GetKernelDef(const std::string &name, bool normalizeCoefs,
                              double weightIfNotNormalized)
{
    auto it = oMapKernelNameToMatrix.find(name);
    CPLAssert(it != oMapKernelNameToMatrix.end());
    KernelDef def;
    def.size = it->second.first;
    int nSum = 0;
    for (const int nVal : it->second.second)
        nSum += nVal;
    const double dfWeight = normalizeCoefs
                                ? 1.0 / (static_cast<double>(nSum) +
                                         std::numeric_limits<double>::min())
                                : weightIfNotNormalized;
    for (const int nVal : it->second.second)
    {
        def.adfCoefficients.push_back(nVal * dfWeight);
    }
    return def;
}

/************************************************************************/
/*               GDALRasterNeighborsAlgorithm::RunStep()                */
/************************************************************************/

bool GDALRasterNeighborsAlgorithm::RunStep(GDALPipelineStepRunContext &)
{
    auto poSrcDS = m_inputDataset[0].GetDatasetRef();
    CPLAssert(!m_outputDataset.GetDatasetRef());

    CPLAssert(m_band <= poSrcDS->GetRasterCount());

    auto eType = GDALGetDataTypeByName(m_type.c_str());
    if (eType == GDT_Unknown)
    {
        eType = GDT_Float64;
    }
    std::vector<KernelDef> aKernelDefs(m_kernel.size());
    std::vector<bool> abNullCoefficientSum(m_kernel.size());
    for (size_t i = 0; i < m_kernel.size(); ++i)
    {
        const std::string &kernel = m_kernel[i];
        if (!kernel.empty() && kernel[0] == '[')
        {
            const CPLStringList aosValues(
                CSLTokenizeString2(kernel.c_str(), "[] ,",
                                   CSLT_STRIPLEADSPACES | CSLT_STRIPENDSPACES));
            const double dfSize = static_cast<double>(aosValues.size());
            KernelDef def;
            def.size = static_cast<int>(std::floor(sqrt(dfSize) + 0.5));
            double dfSum = 0;
            for (const char *pszC : cpl::Iterate(aosValues))
            {
                const double dfV = CPLAtof(pszC);
                dfSum += dfV;
                // Already validated to be numeric by the validation action
                def.adfCoefficients.push_back(dfV);
            }
            abNullCoefficientSum[i] = std::fabs(dfSum) < 1e-10;
            if (abNullCoefficientSum[i] && m_method.size() == m_kernel.size() &&
                m_method[i] == "mean")
            {
                ReportError(
                    CE_Failure, CPLE_AppDefined,
                    "Specifying method = 'mean' for a kernel whose sum of "
                    "coefficients is zero is not allowed. Use 'sum' instead");
                return false;
            }
            aKernelDefs[i] = std::move(def);
        }
    }

    if (m_method.empty())
    {
        for (size_t i = 0; i < m_kernel.size(); ++i)
        {
            const bool bIsZeroSumKernel =
                m_kernel[i] == "u" || m_kernel[i] == "v" ||
                m_kernel[i] == "edge1" || m_kernel[i] == "edge2" ||
                abNullCoefficientSum[i];
            m_method.push_back(bIsZeroSumKernel ? "sum" : "mean");
        }
    }
    else if (m_method.size() == 1)
    {
        const std::string lastValue(m_method.back());
        m_method.resize(m_kernel.size(), lastValue);
    }

    if (m_size == 0 && m_kernel[0][0] != '[')
        m_size = m_kernel[0] == "unsharp-masking" ? 5 : 3;

    for (size_t i = 0; i < m_kernel.size(); ++i)
    {
        const std::string &kernel = m_kernel[i];
        if (aKernelDefs[i].adfCoefficients.empty())
        {
            KernelDef def;
            if (kernel == "edge1" || kernel == "edge2" || kernel == "sharpen")
            {
                CPLAssert(m_size == 3);
                def = GetKernelDef(kernel, false, 1.0);
            }
            else if (kernel == "u" || kernel == "v")
            {
                CPLAssert(m_size == 3);
                def = GetKernelDef(kernel, false, 0.5);
            }
            else if (kernel == "equal")
            {
                def.size = m_size;
                const double dfWeight =
                    m_method[i] == "mean"
                        ? 1.0 / (static_cast<double>(m_size) * m_size +
                                 std::numeric_limits<double>::min())
                        : 1.0;
                def.adfCoefficients.resize(static_cast<size_t>(m_size) * m_size,
                                           dfWeight);
            }
            else if (kernel == "gaussian")
            {
                CPLAssert(m_size == 3 || m_size == 5);
                def = GetKernelDef(
                    m_size == 3 ? "gaussian-3x3" : "gaussian-5x5", true, 0.0);
            }
            else if (kernel == "unsharp-masking")
            {
                CPLAssert(m_size == 5);
                def = GetKernelDef("unsharp-masking-5x5", true, 0.0);
            }
            else
            {
                CPLAssert(false);
            }
            aKernelDefs[i] = std::move(def);
        }
    }

    auto vrt = GDALNeighborsCreateVRTDerived(poSrcDS, m_band, eType, m_nodata,
                                             m_method, aKernelDefs);
    const bool ret = vrt != nullptr;
    if (vrt)
    {
        m_outputDataset.Set(std::move(vrt));
    }
    return ret;
}

GDALRasterNeighborsAlgorithmStandalone::
    ~GDALRasterNeighborsAlgorithmStandalone() = default;

//! @endcond

Coverage Report

Created: 2026-02-14 06:52

Line	Count	Source
1		/******************************************************************************
2		*
3		* Project: GDAL
4		* Purpose: "neighbors" step of "raster pipeline"
5		* Author: Even Rouault <even dot rouault at spatialys.com>
6		*
7		******************************************************************************
8		* Copyright (c) 2025, Even Rouault <even dot rouault at spatialys.com>
9		*
10		* SPDX-License-Identifier: MIT
11		****************************************************************************/
12
13		#include "gdalalg_raster_neighbors.h"
14
15		#include "gdal_priv.h"
16		#include "gdal_priv_templates.hpp"
17		#include "vrtdataset.h"
18
19		#include <algorithm>
20		#include <limits>
21		#include <map>
22		#include <optional>
23		#include <set>
24		#include <utility>
25		#include <vector>
26
27		//! @cond Doxygen_Suppress
28
29		#ifndef _
30	0	#define _(x) (x)
31		#endif
32
33		static const std::set<std::string> oSetKernelNames = {
34		"u", "v", "equal", "edge1",
35		"edge2", "sharpen", "gaussian", "unsharp-masking"};
36
37		namespace
38		{
39		struct KernelDef
40		{
41		int size = 0;
42		std::vector<double> adfCoefficients{};
43		};
44		} // namespace
45
46		// clang-format off
47		// Cf https://en.wikipedia.org/wiki/Kernel_(image_processing)
48		static const std::map<std::string, std::pair<int, std::vector<int>>> oMapKernelNameToMatrix = {
49		{ "u", { 3, { 0, 0, 0,
50		-1, 0, 1,
51		0, 0, 0 } } },
52		{ "v", { 3, { 0, -1, 0,
53		0, 0, 0,
54		0, 1, 0 } } },
55		{ "edge1", { 3, { 0, -1, 0,
56		-1, 4, -1,
57		0, -1, 0 } } },
58		{ "edge2", { 3, { -1, -1, -1,
59		-1, 8, -1,
60		-1, -1, -1 } } },
61		{ "sharpen", { 3, { 0, -1, 0,
62		-1, 5, -1,
63		0, -1, 0 } } },
64		{ "gaussian-3x3", { 3, { 1, 2, 1,
65		2, 4, 2,
66		1, 2, 1 } } },
67		{ "gaussian-5x5", { 5, { 1, 4, 6, 4, 1,
68		4, 16, 24, 16, 4,
69		6, 24, 36, 24, 6,
70		4, 16, 24, 16, 4,
71		1, 4, 6, 4, 1, } } },
72		{ "unsharp-masking-5x5", { 5, { 1, 4, 6, 4, 1,
73		4, 16, 24, 16, 4,
74		6, 24, -476, 24, 6,
75		4, 16, 24, 16, 4,
76		1, 4, 6, 4, 1, } } },
77		};
78
79		// clang-format on
80
81		/************************************************************************/
82		/* CreateDerivedBandXML() */
83		/************************************************************************/
84
85		static bool CreateDerivedBandXML(VRTDataset poVRTDS, GDALRasterBand poSrcBand,
86		GDALDataType eType, const std::string &noData,
87		const std::string &method,
88		const KernelDef &kernelDef)
89	0	{
90	0	poVRTDS->AddBand(eType, nullptr);
91
92	0	std::optional<double> dstNoData;
93	0	bool autoSelectNoDataValue = false;
94	0	if (noData.empty())
95	0	{
96	0	autoSelectNoDataValue = true;
97	0	}
98	0	else if (noData != "none")
99	0	{
100		// Already validated to be numeric by the validation action
101	0	dstNoData = CPLAtof(noData.c_str());
102	0	}
103
104	0	auto poVRTBand = cpl::down_cast<VRTSourcedRasterBand *>(
105	0	poVRTDS->GetRasterBand(poVRTDS->GetRasterCount()));
106
107	0	auto poSource = std::make_unique<VRTKernelFilteredSource>();
108	0	poSrcBand->GetDataset()->Reference();
109	0	poSource->SetSrcBand(poSrcBand);
110	0	poSource->SetKernel(kernelDef.size, /* separable = */ false,
111	0	kernelDef.adfCoefficients);
112	0	poSource->SetNormalized(method != "sum");
113	0	if (method != "sum" && method != "mean")
114	0	poSource->SetFunction(method.c_str());
115
116	0	int bSrcHasNoData = false;
117	0	const double dfNoDataValue = poSrcBand->GetNoDataValue(&bSrcHasNoData);
118	0	if (bSrcHasNoData)
119	0	{
120	0	poSource->SetNoDataValue(dfNoDataValue);
121	0	if (autoSelectNoDataValue && !dstNoData.has_value())
122	0	{
123	0	dstNoData = dfNoDataValue;
124	0	}
125	0	}
126
127	0	if (dstNoData.has_value())
128	0	{
129	0	if (!GDALIsValueExactAs(dstNoData.value(), eType))
130	0	{
131	0	CPLError(CE_Failure, CPLE_AppDefined,
132	0	"Band output type %s cannot represent NoData value %g",
133	0	GDALGetDataTypeName(eType), dstNoData.value());
134	0	return false;
135	0	}
136
137	0	poVRTBand->SetNoDataValue(dstNoData.value());
138	0	}
139
140	0	poVRTBand->AddSource(std::move(poSource));
141
142	0	return true;
143	0	}
144
145		/************************************************************************/
146		/* GDALNeighborsCreateVRTDerived() */
147		/************************************************************************/
148
149		static std::unique_ptr<GDALDataset>
150		GDALNeighborsCreateVRTDerived(GDALDataset *poSrcDS, int nBand,
151		GDALDataType eType, const std::string &noData,
152		const std::vector<std::string> &methods,
153		const std::vector<KernelDef> &aKernelDefs)
154	0	{
155	0	CPLAssert(methods.size() == aKernelDefs.size());
156
157	0	auto ds = std::make_unique<VRTDataset>(poSrcDS->GetRasterXSize(),
158	0	poSrcDS->GetRasterYSize());
159	0	GDALGeoTransform gt;
160	0	if (poSrcDS->GetGeoTransform(gt) == CE_None)
161	0	ds->SetGeoTransform(gt);
162	0	if (const OGRSpatialReference *poSRS = poSrcDS->GetSpatialRef())
163	0	{
164	0	ds->SetSpatialRef(poSRS);
165	0	}
166
167	0	bool ret = true;
168	0	if (nBand != 0)
169	0	{
170	0	for (size_t i = 0; i < aKernelDefs.size() && ret; ++i)
171	0	{
172	0	ret =
173	0	CreateDerivedBandXML(ds.get(), poSrcDS->GetRasterBand(nBand),
174	0	eType, noData, methods[i], aKernelDefs[i]);
175	0	}
176	0	}
177	0	else
178	0	{
179	0	for (int iBand = 1; iBand <= poSrcDS->GetRasterCount(); ++iBand)
180	0	{
181	0	for (size_t i = 0; i < aKernelDefs.size() && ret; ++i)
182	0	{
183	0	ret = CreateDerivedBandXML(ds.get(),
184	0	poSrcDS->GetRasterBand(iBand), eType,
185	0	noData, methods[i], aKernelDefs[i]);
186	0	}
187	0	}
188	0	}
189	0	if (!ret)
190	0	ds.reset();
191	0	return ds;
192	0	}
193
194		/************************************************************************/
195		/* GDALRasterNeighborsAlgorithm::GDALRasterNeighborsAlgorithm() */
196		/************************************************************************/
197
198		GDALRasterNeighborsAlgorithm::GDALRasterNeighborsAlgorithm(
199		bool standaloneStep) noexcept
200	0	: GDALRasterPipelineStepAlgorithm(
201	0	NAME, DESCRIPTION, HELP_URL,
202	0	ConstructorOptions().SetStandaloneStep(standaloneStep))
203	0	{
204	0	AddBandArg(&m_band);
205
206	0	AddArg("method", 0, _("Method to combine weighed source pixels"), &m_method)
207	0	.SetChoices("mean", "sum", "min", "max", "stddev", "median", "mode");
208
209	0	AddArg("size", 0, _("Neighborhood size"), &m_size)
210	0	.SetMinValueIncluded(3)
211	0	.SetMaxValueIncluded(99)
212	0	.AddValidationAction(
213	0	[this]()
214	0	{
215	0	if ((m_size % 2) != 1)
216	0	{
217	0	ReportError(CE_Failure, CPLE_IllegalArg,
218	0	"The value of 'size' must be an odd number.");
219	0	return false;
220	0	}
221	0	return true;
222	0	});
223
224	0	AddArg("kernel", 0, _("Convolution kernel(s) to apply"), &m_kernel)
225	0	.SetPackedValuesAllowed(false)
226	0	.SetMinCount(1)
227	0	.SetMinCharCount(1)
228	0	.SetRequired()
229	0	.SetAutoCompleteFunction(
230	0	[](const std::string &v)
231	0	{
232	0	std::vector<std::string> ret;
233	0	if (v.empty() \|\| v.front() != '[')
234	0	{
235	0	ret.insert(ret.end(), oSetKernelNames.begin(),
236	0	oSetKernelNames.end());
237	0	ret.push_back(
238	0	"[[val00,val10,...,valN0],...,[val0N,val1N,...valNN]]");
239	0	}
240	0	return ret;
241	0	})
242	0	.AddValidationAction(
243	0	[this]()
244	0	{
245	0	for (const std::string &kernel : m_kernel)
246	0	{
247	0	if (kernel.front() == '[' && kernel.back() == ']')
248	0	{
249	0	const CPLStringList aosValues(CSLTokenizeString2(
250	0	kernel.c_str(), "[] ,",
251	0	CSLT_STRIPLEADSPACES \| CSLT_STRIPENDSPACES));
252	0	const double dfSize =
253	0	static_cast<double>(aosValues.size());
254	0	const int nSqrt = static_cast<int>(sqrt(dfSize) + 0.5);
255	0	if (!((aosValues.size() % 2) == 1 &&
256	0	nSqrt * nSqrt == aosValues.size()))
257	0	{
258	0	ReportError(
259	0	CE_Failure, CPLE_IllegalArg,
260	0	"The number of values in the 'kernel' "
261	0	"argument must be an odd square number.");
262	0	return false;
263	0	}
264	0	for (int i = 0; i < aosValues.size(); ++i)
265	0	{
266	0	if (CPLGetValueType(aosValues[i]) ==
267	0	CPL_VALUE_STRING)
268	0	{
269	0	ReportError(CE_Failure, CPLE_IllegalArg,
270	0	"Non-numeric value found in the "
271	0	"'kernel' argument");
272	0	return false;
273	0	}
274	0	}
275	0	}
276	0	else if (!cpl::contains(oSetKernelNames, kernel))
277	0	{
278	0	std::string osMsg =
279	0	"Valid values for 'kernel' argument are: ";
280	0	bool bFirst = true;
281	0	for (const auto &name : oSetKernelNames)
282	0	{
283	0	if (!bFirst)
284	0	osMsg += ", ";
285	0	bFirst = false;
286	0	osMsg += '\'';
287	0	osMsg += name;
288	0	osMsg += '\'';
289	0	}
290	0	osMsg += " or "
291	0	"[[val00,val10,...,valN0],...,[val0N,val1N,..."
292	0	"valNN]]";
293	0	ReportError(CE_Failure, CPLE_IllegalArg, "%s",
294	0	osMsg.c_str());
295	0	return false;
296	0	}
297	0	}
298	0	return true;
299	0	});
300
301	0	AddOutputDataTypeArg(&m_type).SetDefault("Float64");
302
303	0	AddNodataArg(&m_nodata, true);
304
305	0	AddValidationAction(
306	0	[this]()
307	0	{
308	0	if (m_method.size() > 1 && m_method.size() != m_kernel.size())
309	0	{
310	0	ReportError(
311	0	CE_Failure, CPLE_AppDefined,
312	0	"The number of values for the 'method' argument should "
313	0	"be one or exactly the number of values of 'kernel'");
314	0	return false;
315	0	}
316
317	0	if (m_size > 0)
318	0	{
319	0	for (const std::string &kernel : m_kernel)
320	0	{
321	0	if (kernel == "gaussian")
322	0	{
323	0	if (m_size != 3 && m_size != 5)
324	0	{
325	0	ReportError(CE_Failure, CPLE_AppDefined,
326	0	"Currently only size = 3 or 5 is "
327	0	"supported for kernel '%s'",
328	0	kernel.c_str());
329	0	return false;
330	0	}
331	0	}
332	0	else if (kernel == "unsharp-masking")
333	0	{
334	0	if (m_size != 5)
335	0	{
336	0	ReportError(CE_Failure, CPLE_AppDefined,
337	0	"Currently only size = 5 is supported "
338	0	"for kernel '%s'",
339	0	kernel.c_str());
340	0	return false;
341	0	}
342	0	}
343	0	else if (kernel[0] == '[')
344	0	{
345	0	const CPLStringList aosValues(CSLTokenizeString2(
346	0	kernel.c_str(), "[] ,",
347	0	CSLT_STRIPLEADSPACES \| CSLT_STRIPENDSPACES));
348	0	const double dfSize =
349	0	static_cast<double>(aosValues.size());
350	0	const int size =
351	0	static_cast<int>(std::floor(sqrt(dfSize) + 0.5));
352	0	if (m_size != size)
353	0	{
354	0	ReportError(CE_Failure, CPLE_AppDefined,
355	0	"Value of 'size' argument (%d) "
356	0	"inconsistent with the one deduced "
357	0	"from the kernel matrix (%d)",
358	0	m_size, size);
359	0	return false;
360	0	}
361	0	}
362	0	else if (m_size != 3 && kernel != "equal" &&
363	0	kernel[0] != '[')
364	0	{
365	0	ReportError(CE_Failure, CPLE_AppDefined,
366	0	"Currently only size = 3 is supported for "
367	0	"kernel '%s'",
368	0	kernel.c_str());
369	0	return false;
370	0	}
371	0	}
372	0	}
373
374	0	return true;
375	0	});
376	0	}
377
378		/************************************************************************/
379		/* GetKernelDef() */
380		/************************************************************************/
381
382		static KernelDef GetKernelDef(const std::string &name, bool normalizeCoefs,
383		double weightIfNotNormalized)
384	0	{
385	0	auto it = oMapKernelNameToMatrix.find(name);
386	0	CPLAssert(it != oMapKernelNameToMatrix.end());
387	0	KernelDef def;
388	0	def.size = it->second.first;
389	0	int nSum = 0;
390	0	for (const int nVal : it->second.second)
391	0	nSum += nVal;
392	0	const double dfWeight = normalizeCoefs
393	0	? 1.0 / (static_cast<double>(nSum) +
394	0	std::numeric_limits<double>::min())
395	0	: weightIfNotNormalized;
396	0	for (const int nVal : it->second.second)
397	0	{
398	0	def.adfCoefficients.push_back(nVal * dfWeight);
399	0	}
400	0	return def;
401	0	}
402
403		/************************************************************************/
404		/* GDALRasterNeighborsAlgorithm::RunStep() */
405		/************************************************************************/
406
407		bool GDALRasterNeighborsAlgorithm::RunStep(GDALPipelineStepRunContext &)
408	0	{
409	0	auto poSrcDS = m_inputDataset[0].GetDatasetRef();
410	0	CPLAssert(!m_outputDataset.GetDatasetRef());
411
412	0	CPLAssert(m_band <= poSrcDS->GetRasterCount());
413
414	0	auto eType = GDALGetDataTypeByName(m_type.c_str());
415	0	if (eType == GDT_Unknown)
416	0	{
417	0	eType = GDT_Float64;
418	0	}
419	0	std::vector<KernelDef> aKernelDefs(m_kernel.size());
420	0	std::vector<bool> abNullCoefficientSum(m_kernel.size());
421	0	for (size_t i = 0; i < m_kernel.size(); ++i)
422	0	{
423	0	const std::string &kernel = m_kernel[i];
424	0	if (!kernel.empty() && kernel[0] == '[')
425	0	{
426	0	const CPLStringList aosValues(
427	0	CSLTokenizeString2(kernel.c_str(), "[] ,",
428	0	CSLT_STRIPLEADSPACES \| CSLT_STRIPENDSPACES));
429	0	const double dfSize = static_cast<double>(aosValues.size());
430	0	KernelDef def;
431	0	def.size = static_cast<int>(std::floor(sqrt(dfSize) + 0.5));
432	0	double dfSum = 0;
433	0	for (const char *pszC : cpl::Iterate(aosValues))
434	0	{
435	0	const double dfV = CPLAtof(pszC);
436	0	dfSum += dfV;
437		// Already validated to be numeric by the validation action
438	0	def.adfCoefficients.push_back(dfV);
439	0	}
440	0	abNullCoefficientSum[i] = std::fabs(dfSum) < 1e-10;
441	0	if (abNullCoefficientSum[i] && m_method.size() == m_kernel.size() &&
442	0	m_method[i] == "mean")
443	0	{
444	0	ReportError(
445	0	CE_Failure, CPLE_AppDefined,
446	0	"Specifying method = 'mean' for a kernel whose sum of "
447	0	"coefficients is zero is not allowed. Use 'sum' instead");
448	0	return false;
449	0	}
450	0	aKernelDefs[i] = std::move(def);
451	0	}
452	0	}
453
454	0	if (m_method.empty())
455	0	{
456	0	for (size_t i = 0; i < m_kernel.size(); ++i)
457	0	{
458	0	const bool bIsZeroSumKernel =
459	0	m_kernel[i] == "u" \|\| m_kernel[i] == "v" \|\|
460	0	m_kernel[i] == "edge1" \|\| m_kernel[i] == "edge2" \|\|
461	0	abNullCoefficientSum[i];
462	0	m_method.push_back(bIsZeroSumKernel ? "sum" : "mean");
463	0	}
464	0	}
465	0	else if (m_method.size() == 1)
466	0	{
467	0	const std::string lastValue(m_method.back());
468	0	m_method.resize(m_kernel.size(), lastValue);
469	0	}
470
471	0	if (m_size == 0 && m_kernel[0][0] != '[')
472	0	m_size = m_kernel[0] == "unsharp-masking" ? 5 : 3;
473
474	0	for (size_t i = 0; i < m_kernel.size(); ++i)
475	0	{
476	0	const std::string &kernel = m_kernel[i];
477	0	if (aKernelDefs[i].adfCoefficients.empty())
478	0	{
479	0	KernelDef def;
480	0	if (kernel == "edge1" \|\| kernel == "edge2" \|\| kernel == "sharpen")
481	0	{
482	0	CPLAssert(m_size == 3);
483	0	def = GetKernelDef(kernel, false, 1.0);
484	0	}
485	0	else if (kernel == "u" \|\| kernel == "v")
486	0	{
487	0	CPLAssert(m_size == 3);
488	0	def = GetKernelDef(kernel, false, 0.5);
489	0	}
490	0	else if (kernel == "equal")
491	0	{
492	0	def.size = m_size;
493	0	const double dfWeight =
494	0	m_method[i] == "mean"
495	0	? 1.0 / (static_cast<double>(m_size) * m_size +
496	0	std::numeric_limits<double>::min())
497	0	: 1.0;
498	0	def.adfCoefficients.resize(static_cast<size_t>(m_size) * m_size,
499	0	dfWeight);
500	0	}
501	0	else if (kernel == "gaussian")
502	0	{
503	0	CPLAssert(m_size == 3 \|\| m_size == 5);
504	0	def = GetKernelDef(
505	0	m_size == 3 ? "gaussian-3x3" : "gaussian-5x5", true, 0.0);
506	0	}
507	0	else if (kernel == "unsharp-masking")
508	0	{
509	0	CPLAssert(m_size == 5);
510	0	def = GetKernelDef("unsharp-masking-5x5", true, 0.0);
511	0	}
512	0	else
513	0	{
514	0	CPLAssert(false);
515	0	}
516	0	aKernelDefs[i] = std::move(def);
517	0	}
518	0	}
519
520	0	auto vrt = GDALNeighborsCreateVRTDerived(poSrcDS, m_band, eType, m_nodata,
521	0	m_method, aKernelDefs);
522	0	const bool ret = vrt != nullptr;
523	0	if (vrt)
524	0	{
525	0	m_outputDataset.Set(std::move(vrt));
526	0	}
527	0	return ret;
528	0	}
529
530		GDALRasterNeighborsAlgorithmStandalone::
531	0	~GDALRasterNeighborsAlgorithmStandalone() = default;
532
533		//! @endcond