diff --git a/alg/gdalwarper.cpp b/alg/gdalwarper.cpp index 4a9d075b4331..e278979462f8 100644 --- a/alg/gdalwarper.cpp +++ b/alg/gdalwarper.cpp @@ -1249,6 +1249,24 @@ CPLErr GDALWarpDstAlphaMasker(void *pMaskFuncArg, int nBandCount, * an explicit source and target SRS. *
  • MULT_FACTOR_VERTICAL_SHIFT: Multiplication factor for the vertical * shift. Default 1.0
    • + * + *
  • 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. + *
    • + * + *
  • 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 (%)
    • + * * */ diff --git a/alg/gdalwarper.h b/alg/gdalwarper.h index d23e486371fb..90ca632a091b 100644 --- a/alg/gdalwarper.h +++ b/alg/gdalwarper.h @@ -461,6 +461,12 @@ class CPL_DLL GDALWarpKernel bool bApplyVerticalShift = false; double dfMultFactorVerticalShift = 1.0; + + // Tuples of values (e.g. ",," or "(,,),(,,)") that must + // be ignored as contributing source pixels during resampling. Only taken into account by + // Average currently + std::vector> m_aadfExcludedValues{}; + /*! @endcond */ GDALWarpKernel(); diff --git a/alg/gdalwarpkernel.cpp b/alg/gdalwarpkernel.cpp index b8f87b6cda1a..1bfb1fa78191 100644 --- a/alg/gdalwarpkernel.cpp +++ b/alg/gdalwarpkernel.cpp @@ -1407,6 +1407,37 @@ CPLErr GDALWarpKernel::Validate() return CE_Failure; } + // Tuples of values (e.g. ",," or "(,,),(,,)") 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 ,, or " + "(,,),(,,) if there are multiple " + "tuples", + nBands); + return CE_Failure; + } + std::vector 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(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 adfValueReal(poWK->nBands, 0); + std::vector adfValueAveraged(poWK->nBands, 0); + std::vector anCountExcludedValues( + poWK->m_aadfExcludedValues.size(), 0); + + for (int iSrcY = iSrcYMin; iSrcY < iSrcYMax; iSrcY++) + { + const double dfWeightY = COMPUTE_WEIGHT_Y(iSrcY); + iSrcOffset = + iSrcXMin + static_cast(iSrcY) * nSrcXSize; + for (int iSrcX = iSrcXMin; iSrcX < iSrcXMax; + iSrcX++, iSrcOffset++) + { + if (bWrapOverX) + iSrcOffset = + (iSrcX % nSrcXSize) + + static_cast(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) { diff --git a/autotest/alg/warp.py b/autotest/alg/warp.py index 2729fd4f3264..17ba627b8f25 100755 --- a/autotest/alg/warp.py +++ b/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, + ) diff --git a/autotest/pyscripts/test_gdal2tiles.py b/autotest/pyscripts/test_gdal2tiles.py index 00ca7f684070..0b24b8306679 100755 --- a/autotest/pyscripts/test_gdal2tiles.py +++ b/autotest/pyscripts/test_gdal2tiles.py @@ -33,12 +33,13 @@ import os import os.path import shutil +import struct import sys import pytest import test_py_scripts # noqa # pylint: disable=E0401 -from osgeo import gdal # noqa +from osgeo import gdal, osr # noqa from osgeo_utils.gdalcompare import compare_db pytestmark = pytest.mark.skipif( @@ -581,3 +582,44 @@ def test_gdal2tiles_py_webp(script_path, tmp_path, resampling): ) diff_found = compare_db(gdal.Open(webp_filename), gdal.Open(filename)) assert not diff_found, (resampling, filename) + + +@pytest.mark.require_driver("PNG") +def test_gdal2tiles_excluded_values(script_path, tmp_path): + + input_tif = str(tmp_path / "test_gdal2tiles_excluded_values.tif") + output_folder = str(tmp_path / "test_gdal2tiles_excluded_values") + + src_ds = gdal.GetDriverByName("GTiff").Create(input_tif, 256, 256, 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) + ) + srs = osr.SpatialReference() + srs.ImportFromEPSG(3857) + src_ds.SetSpatialRef(srs) + MAX_GM = 20037508.342789244 + RES_Z0 = 2 * MAX_GM / 256 + RES_Z1 = RES_Z0 / 2 + # Spatial extent of tile (0,0) at zoom level 1 + src_ds.SetGeoTransform([-MAX_GM, RES_Z1, 0, MAX_GM, 0, -RES_Z1]) + src_ds = None + + test_py_scripts.run_py_script_as_external_script( + script_path, + "gdal2tiles", + f"-q -z 0-1 --excluded-values=30,31,32 --excluded-values-pct-threshold=50 {input_tif} {output_folder}", + ) + + ds = gdal.Open(f"{output_folder}/0/0/0.png") + assert struct.unpack("B" * 4, ds.ReadRaster(0, 0, 1, 1)) == ( + (10 + 20 + 40) // 3, + (11 + 21 + 41) // 3, + (12 + 22 + 42) // 3, + 255, + ) diff --git a/doc/source/programs/gdal2tiles.rst b/doc/source/programs/gdal2tiles.rst index abc29fc71dd6..79ca7e41e8f6 100644 --- a/doc/source/programs/gdal2tiles.rst +++ b/doc/source/programs/gdal2tiles.rst @@ -22,6 +22,7 @@ Synopsis [-w ] [-t ] [-c <copyright>] [--processes=<NB_PROCESSES>] [--mpi] [--xyz] [--tilesize=<PIXELS>] [--tmscompatible] + [--excluded-values=<EXCLUDED_VALUES>] [--excluded-values-pct-threshold=<EXCLUDED_VALUES_PCT_THRESHOLD>] [-g <googlekey] [-b <bingkey>] <input_file> [<output_dir>] [<COMMON_OPTIONS>] Description @@ -143,6 +144,23 @@ can publish a picture without proper georeferencing too. .. versionadded:: 3.6 +.. option:: --excluded-values=<EXCLUDED_VALUES> + + 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. + +.. option:: --excluded-values-pct-threshold=EXCLUDED_VALUES_PCT_THRESHOLD + + 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(%) + + .. versionadded:: 3.9 .. option:: -h, --help diff --git a/swig/python/gdal-utils/osgeo_utils/gdal2tiles.py b/swig/python/gdal-utils/osgeo_utils/gdal2tiles.py index c33f08f6e537..6d255b1a61e3 100644 --- a/swig/python/gdal-utils/osgeo_utils/gdal2tiles.py +++ b/swig/python/gdal-utils/osgeo_utils/gdal2tiles.py @@ -873,7 +873,27 @@ def scale_query_to_tile(dsquery, dstile, options, tilefilename=""): tile_size = dstile.RasterXSize tilebands = dstile.RasterCount - if options.resampling == "average": + dsquery.SetGeoTransform( + ( + 0.0, + tile_size / float(querysize), + 0.0, + 0.0, + 0.0, + tile_size / float(querysize), + ) + ) + dstile.SetGeoTransform((0.0, 1.0, 0.0, 0.0, 0.0, 1.0)) + + if options.resampling == "average" and options.excluded_values: + + gdal.Warp( + dstile, + dsquery, + options=f"-r average -wo EXCLUDED_VALUES={options.excluded_values} -wo EXCLUDED_VALUES_PCT_THRESHOLD={options.excluded_values_pct_threshold}", + ) + + elif options.resampling == "average": # Function: gdal.RegenerateOverview() for i in range(1, tilebands + 1): @@ -949,18 +969,6 @@ def scale_query_to_tile(dsquery, dstile, options, tilefilename=""): gdal_resampling = gdal.GRA_Q3 # Other algorithms are implemented by gdal.ReprojectImage(). - dsquery.SetGeoTransform( - ( - 0.0, - tile_size / float(querysize), - 0.0, - 0.0, - 0.0, - tile_size / float(querysize), - ) - ) - dstile.SetGeoTransform((0.0, 1.0, 0.0, 0.0, 0.0, 1.0)) - res = gdal.ReprojectImage(dsquery, dstile, None, None, gdal_resampling) if res != 0: exit_with_error( @@ -1346,6 +1354,7 @@ def create_base_tile(tile_job_info: "TileJobInfo", tile_detail: "TileDetail") -> # Tile dataset in memory tilefilename = os.path.join(output, str(tz), str(tx), "%s.%s" % (ty, tileext)) dstile = mem_drv.Create("", tile_size, tile_size, tilebands) + dstile.GetRasterBand(tilebands).SetColorInterpretation(gdal.GCI_AlphaBand) data = alpha = None @@ -1399,6 +1408,8 @@ def create_base_tile(tile_job_info: "TileJobInfo", tile_detail: "TileDetail") -> # Big ReadRaster query in memory scaled to the tile_size - all but 'near' # algo dsquery = mem_drv.Create("", querysize, querysize, tilebands) + dsquery.GetRasterBand(tilebands).SetColorInterpretation(gdal.GCI_AlphaBand) + # TODO: fill the null value in case a tile without alpha is produced (now # only png tiles are supported) dsquery.WriteRaster( @@ -1422,6 +1433,11 @@ def create_base_tile(tile_job_info: "TileJobInfo", tile_detail: "TileDetail") -> tilefilename, dstile, strict=0, options=_get_creation_options(options) ) + # Remove useless side car file + aux_xml = tilefilename + ".aux.xml" + if gdal.VSIStatL(aux_xml) is not None: + gdal.Unlink(aux_xml) + del dstile # Create a KML file for this tile. @@ -1485,10 +1501,12 @@ def create_overview_tile( dsquery = mem_driver.Create( "", 2 * tile_job_info.tile_size, 2 * tile_job_info.tile_size, tilebands ) + dsquery.GetRasterBand(tilebands).SetColorInterpretation(gdal.GCI_AlphaBand) # TODO: fill the null value dstile = mem_driver.Create( "", tile_job_info.tile_size, tile_job_info.tile_size, tilebands ) + dstile.GetRasterBand(tilebands).SetColorInterpretation(gdal.GCI_AlphaBand) usable_base_tiles = [] @@ -1759,6 +1777,19 @@ def optparse_init() -> optparse.OptionParser: type="choice", help="which tile driver to use for the tiles", ) + p.add_option( + "--excluded-values", + dest="excluded_values", + type=str, + help="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 for average resampling", + ) + p.add_option( + "--excluded-values-pct-threshold", + dest="excluded_values_pct_threshold", + type=float, + default=50, + help="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 (%)", + ) # KML options g = optparse.OptionGroup(