Warper: add a EXCLUDED_VALUES warping option, and use it in gdal2tiles

alg/gdalwarper.cpp

@@ -1249,6 +1249,24 @@ CPLErr GDALWarpDstAlphaMasker(void *pMaskFuncArg, int nBandCount,
  * an explicit source and target SRS.</li>
  * <li>MULT_FACTOR_VERTICAL_SHIFT: Multiplication factor for the vertical
  * shift. Default 1.0</li>
+ *
+ * <li>EXCLUDED_VALUES: (GDAL >= 3.9) Comma-separated tuple of values
+ * (thus typically "R,G,B"), that are ignored as contributing source
+ * pixels during resampling. The number of values in the tuple must be the same
+ * as the number of bands, excluding the alpha band.
+ * Several tuples of excluded values may be specified using the
+ * "(R1,G1,B2),(R2,G2,B2)" syntax.
+ * Only taken into account by Average currently.
+ * This concept is a bit similar to nodata/alpha, but the main difference is
+ * that pixels matching one of the excluded value tuples are still considered
+ * as valid, when determining the target pixel validity/density.
+ * </li>
+ *
+ * <li>EXCLUDED_VALUES_PCT_THRESHOLD=[0-100]: (GDAL >= 3.9) Minimum percentage
+ * of source pixels that must be set at one of the EXCLUDED_VALUES to cause
+ * the excluded value, that is in majority among source pixels, to be used as the
+ * target pixel value. Default value is 50 (%)</li>
+ *
  * </ul>
  */
 

alg/gdalwarper.h

@@ -461,6 +461,12 @@ class CPL_DLL GDALWarpKernel
     bool bApplyVerticalShift = false;
 
     double dfMultFactorVerticalShift = 1.0;
+
+    // Tuples of values (e.g. "<R>,<G>,<B>" or "(<R1>,<G1>,<B1>),(<R2>,<G2>,<B2>)") that must
+    // be ignored as contributing source pixels during resampling. Only taken into account by
+    // Average currently
+    std::vector<std::vector<double>> m_aadfExcludedValues{};
+
     /*! @endcond */
 
     GDALWarpKernel();

alg/gdalwarpkernel.cpp

@@ -1407,6 +1407,37 @@ CPLErr GDALWarpKernel::Validate()
         return CE_Failure;
     }
 
+    // Tuples of values (e.g. "<R>,<G>,<B>" or "(<R1>,<G1>,<B1>),(<R2>,<G2>,<B2>)") that must
+    // be ignored as contributing source pixels during resampling. Only taken into account by
+    // Average currently
+    const char *pszExcludedValues =
+        CSLFetchNameValue(papszWarpOptions, "EXCLUDED_VALUES");
+    if (pszExcludedValues)
+    {
+        const CPLStringList aosTokens(
+            CSLTokenizeString2(pszExcludedValues, "(,)", 0));
+        if ((aosTokens.size() % nBands) != 0)
+        {
+            CPLError(CE_Failure, CPLE_AppDefined,
+                     "EXCLUDED_VALUES should contain one or several tuples of "
+                     "%d values formatted like <R>,<G>,<B> or "
+                     "(<R1>,<G1>,<B1>),(<R2>,<G2>,<B2>) if there are multiple "
+                     "tuples",
+                     nBands);
+            return CE_Failure;
+        }
+        std::vector<double> adfTuple;
+        for (int i = 0; i < aosTokens.size(); ++i)
+        {
+            adfTuple.push_back(CPLAtof(aosTokens[i]));
+            if (((i + 1) % nBands) == 0)
+            {
+                m_aadfExcludedValues.push_back(adfTuple);
+                adfTuple.clear();
+            }
+        }
+    }
+
     return CE_None;
 }
 
@@ -6549,6 +6580,11 @@ static void GWKAverageOrModeThread(void *pData)
     const double dfErrorThreshold = CPLAtof(
         CSLFetchNameValueDef(poWK->papszWarpOptions, "ERROR_THRESHOLD", "0"));
 
+    const double dfExcludedValuesThreshold =
+        CPLAtof(CSLFetchNameValueDef(poWK->papszWarpOptions,
+                                     "EXCLUDED_VALUES_PCT_THRESHOLD", "50")) /
+        100.0;
+
     const int nXMargin =
         2 * std::max(1, static_cast<int>(std::ceil(1. / poWK->dfXScale)));
     const int nYMargin =
@@ -6689,13 +6725,161 @@ static void GWKAverageOrModeThread(void *pData)
             if (iSrcYMin == iSrcYMax && iSrcYMax < nSrcYSize)
                 iSrcYMax++;
 
+#define COMPUTE_WEIGHT_Y(iSrcY)                                                \
+    ((iSrcY == iSrcYMin)                                                       \
+         ? ((iSrcYMin + 1 == iSrcYMax) ? 1.0 : 1 - (dfYMin - iSrcYMin))        \
+     : (iSrcY + 1 == iSrcYMax) ? 1 - (iSrcYMax - dfYMax)                       \
+                               : 1.0)
+
+#define COMPUTE_WEIGHT(iSrcX, dfWeightY)                                       \
+    ((iSrcX == iSrcXMin)       ? ((iSrcXMin + 1 == iSrcXMax)                   \
+                                      ? dfWeightY                              \
+                                      : dfWeightY * (1 - (dfXMin - iSrcXMin))) \
+     : (iSrcX + 1 == iSrcXMax) ? dfWeightY * (1 - (iSrcXMax - dfXMax))         \
+                               : dfWeightY)
+
+            bool bDone = false;
+
+            // Special Average mode where we process all bands together,
+            // to avoid averaging tuples that match an entry of m_aadfExcludedValues
+            if (nAlgo == GWKAOM_Average &&
+                !poWK->m_aadfExcludedValues.empty() &&
+                !poWK->bApplyVerticalShift && !bIsComplex)
+            {
+                double dfTotalWeightInvalid = 0.0;
+                double dfTotalWeightExcluded = 0.0;
+                double dfTotalWeightRegular = 0.0;
+                std::vector<double> adfValueReal(poWK->nBands, 0);
+                std::vector<double> adfValueAveraged(poWK->nBands, 0);
+                std::vector<int> anCountExcludedValues(
+                    poWK->m_aadfExcludedValues.size(), 0);
+
+                for (int iSrcY = iSrcYMin; iSrcY < iSrcYMax; iSrcY++)
+                {
+                    const double dfWeightY = COMPUTE_WEIGHT_Y(iSrcY);
+                    iSrcOffset =
+                        iSrcXMin + static_cast<GPtrDiff_t>(iSrcY) * nSrcXSize;
+                    for (int iSrcX = iSrcXMin; iSrcX < iSrcXMax;
+                         iSrcX++, iSrcOffset++)
+                    {
+                        if (bWrapOverX)
+                            iSrcOffset =
+                                (iSrcX % nSrcXSize) +
+                                static_cast<GPtrDiff_t>(iSrcY) * nSrcXSize;
+
+                        if (poWK->panUnifiedSrcValid != nullptr &&
+                            !CPLMaskGet(poWK->panUnifiedSrcValid, iSrcOffset))
+                        {
+                            continue;
+                        }
+
+                        const double dfWeight =
+                            COMPUTE_WEIGHT(iSrcX, dfWeightY);
+                        if (dfWeight <= 0)
+                            continue;
+
+                        bool bAllValid = true;
+                        for (int iBand = 0; iBand < poWK->nBands; iBand++)
+                        {
+                            double dfBandDensity = 0;
+                            double dfValueImagTmp = 0;
+                            if (!(GWKGetPixelValue(
+                                      poWK, iBand, iSrcOffset, &dfBandDensity,
+                                      &adfValueReal[iBand], &dfValueImagTmp) &&
+                                  dfBandDensity > BAND_DENSITY_THRESHOLD))
+                            {
+                                bAllValid = false;
+                                break;
+                            }
+                        }
+
+                        if (!bAllValid)
+                        {
+                            dfTotalWeightInvalid += dfWeight;
+                            continue;
+                        }
+
+                        bool bExcludedValueFound = false;
+                        for (size_t i = 0;
+                             i < poWK->m_aadfExcludedValues.size(); ++i)
+                        {
+                            if (poWK->m_aadfExcludedValues[i] == adfValueReal)
+                            {
+                                bExcludedValueFound = true;
+                                ++anCountExcludedValues[i];
+                                dfTotalWeightExcluded += dfWeight;
+                                break;
+                            }
+                        }
+                        if (!bExcludedValueFound)
+                        {
+                            // Weighted incremental algorithm mean
+                            // Cf https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Weighted_incremental_algorithm
+                            dfTotalWeightRegular += dfWeight;
+                            for (int iBand = 0; iBand < poWK->nBands; iBand++)
+                            {
+                                adfValueAveraged[iBand] +=
+                                    (dfWeight / dfTotalWeightRegular) *
+                                    (adfValueReal[iBand] -
+                                     adfValueAveraged[iBand]);
+                            }
+                        }
+                    }
+                }
+
+                const double dfTotalWeight = dfTotalWeightInvalid +
+                                             dfTotalWeightExcluded +
+                                             dfTotalWeightRegular;
+                if (dfTotalWeightExcluded > 0 &&
+                    dfTotalWeightExcluded >=
+                        dfExcludedValuesThreshold * dfTotalWeight)
+                {
+                    // Find the most represented excluded value tuple
+                    size_t iExcludedValue = 0;
+                    int nExcludedValueCount = 0;
+                    for (size_t i = 0; i < poWK->m_aadfExcludedValues.size();
+                         ++i)
+                    {
+                        if (anCountExcludedValues[i] > nExcludedValueCount)
+                        {
+                            iExcludedValue = i;
+                            nExcludedValueCount = anCountExcludedValues[i];
+                        }
+                    }
+
+                    bHasFoundDensity = true;
+
+                    for (int iBand = 0; iBand < poWK->nBands; iBand++)
+                    {
+                        GWKSetPixelValue(
+                            poWK, iBand, iDstOffset, /* dfBandDensity = */ 1.0,
+                            poWK->m_aadfExcludedValues[iExcludedValue][iBand],
+                            0);
+                    }
+                }
+                else if (dfTotalWeightRegular > 0)
+                {
+                    bHasFoundDensity = true;
+
+                    for (int iBand = 0; iBand < poWK->nBands; iBand++)
+                    {
+                        GWKSetPixelValue(poWK, iBand, iDstOffset,
+                                         /* dfBandDensity = */ 1.0,
+                                         adfValueAveraged[iBand], 0);
+                    }
+                }
+
+                // Skip below loop on bands
+                bDone = true;
+            }
+
             /* ====================================================================
              */
             /*      Loop processing each band. */
             /* ====================================================================
              */
 
-            for (int iBand = 0; iBand < poWK->nBands; iBand++)
+            for (int iBand = 0; !bDone && iBand < poWK->nBands; iBand++)
             {
                 double dfBandDensity = 0.0;
                 double dfValueReal = 0.0;
@@ -6711,19 +6895,6 @@ static void GWKAverageOrModeThread(void *pData)
 
                 // Loop over source lines and pixels - 3 possible algorithms.
 
-#define COMPUTE_WEIGHT_Y(iSrcY)                                                \
-    ((iSrcY == iSrcYMin)                                                       \
-         ? ((iSrcYMin + 1 == iSrcYMax) ? 1.0 : 1 - (dfYMin - iSrcYMin))        \
-     : (iSrcY + 1 == iSrcYMax) ? 1 - (iSrcYMax - dfYMax)                       \
-                               : 1.0)
-
-#define COMPUTE_WEIGHT(iSrcX, dfWeightY)                                       \
-    ((iSrcX == iSrcXMin)       ? ((iSrcXMin + 1 == iSrcXMax)                   \
-                                      ? dfWeightY                              \
-                                      : dfWeightY * (1 - (dfXMin - iSrcXMin))) \
-     : (iSrcX + 1 == iSrcXMax) ? dfWeightY * (1 - (iSrcXMax - dfXMax))         \
-                               : dfWeightY)
-
                 // poWK->eResample == GRA_Average.
                 if (nAlgo == GWKAOM_Average)
                 {

autotest/alg/warp.py

@@ -1740,3 +1740,78 @@ def test_non_square():
         assert res.ymax == pytest.approx(30.25)
         assert res.xmin == pytest.approx(9.9)
         assert res.xmax == pytest.approx(10.5)
+
+
+###############################################################################
+# Test EXCLUDED_VALUES warping option with average resampling
+
+
+def test_warp_average_excluded_values():
+
+    src_ds = gdal.GetDriverByName("MEM").Create("", 2, 2, 3, gdal.GDT_Byte)
+    src_ds.GetRasterBand(1).WriteRaster(
+        0, 0, 2, 2, struct.pack("B" * 4, 10, 20, 30, 40)
+    )
+    src_ds.GetRasterBand(2).WriteRaster(
+        0, 0, 2, 2, struct.pack("B" * 4, 11, 21, 31, 41)
+    )
+    src_ds.GetRasterBand(3).WriteRaster(
+        0, 0, 2, 2, struct.pack("B" * 4, 12, 22, 32, 42)
+    )
+    src_ds.SetGeoTransform([1, 1, 0, 1, 0, 1])
+
+    with pytest.raises(
+        Exception,
+        match="EXCLUDED_VALUES should contain one or several tuples of 3 values",
+    ):
+        out_ds = gdal.Warp(
+            "", src_ds, options="-of MEM -ts 1 1 -r average -wo EXCLUDED_VALUES=30,31"
+        )
+
+    # The excluded value is just ignored in contributing source pixels that are average, as it represents only 25% of contributing pixels
+    out_ds = gdal.Warp(
+        "", src_ds, options="-of MEM -ts 1 1 -r average -wo EXCLUDED_VALUES=(30,31,32)"
+    )
+    assert struct.unpack("B" * 3, out_ds.ReadRaster()) == (
+        (10 + 20 + 40) // 3,
+        (11 + 21 + 41) // 3,
+        (12 + 22 + 42) // 3,
+    )
+
+    # The excluded value is selected because its contributing 25% is >= 0%
+    out_ds = gdal.Warp(
+        "",
+        src_ds,
+        options="-of MEM -ts 1 1 -r average -wo EXCLUDED_VALUES=(30,31,32) -wo EXCLUDED_VALUES_PCT_THRESHOLD=0",
+    )
+    assert struct.unpack("B" * 3, out_ds.ReadRaster()) == (30, 31, 32)
+
+    # The excluded value is selected because its contributing 25% is >= 24%
+    out_ds = gdal.Warp(
+        "",
+        src_ds,
+        options="-of MEM -ts 1 1 -r average -wo EXCLUDED_VALUES=(30,31,32) -wo EXCLUDED_VALUES_PCT_THRESHOLD=24",
+    )
+    assert struct.unpack("B" * 3, out_ds.ReadRaster()) == (30, 31, 32)
+
+    # The excluded value is selected because its contributing 25% is < 26%
+    out_ds = gdal.Warp(
+        "",
+        src_ds,
+        options="-of MEM -ts 1 1 -r average -wo EXCLUDED_VALUES=(30,31,32) -wo EXCLUDED_VALUES_PCT_THRESHOLD=26",
+    )
+    assert struct.unpack("B" * 3, out_ds.ReadRaster()) == (
+        (10 + 20 + 40) // 3,
+        (11 + 21 + 41) // 3,
+        (12 + 22 + 42) // 3,
+    )
+
+    # No match of excluded value
+    out_ds = gdal.Warp(
+        "", src_ds, options="-of MEM -ts 1 1 -r average -wo EXCLUDED_VALUES=(30,31,0)"
+    )
+    assert struct.unpack("B" * 3, out_ds.ReadRaster()) == (
+        (10 + 20 + 30 + 40) // 4,
+        (11 + 21 + 31 + 41) // 4,
+        (12 + 22 + 32 + 42) // 4,
+    )