From a803e49ed04b4a0abf294be46b651a6247391e18 Mon Sep 17 00:00:00 2001 From: Dan Baston Date: Mon, 13 Jan 2025 20:26:40 -0500 Subject: [PATCH] Warper: Add MODE_TIES warp option --- alg/gdalwarper.cpp | 5 +++ alg/gdalwarper.h | 12 ++++++++ alg/gdalwarpkernel.cpp | 65 +++++++++++++++++++++++++++++++++++---- alg/gdalwarpoperation.cpp | 1 + apps/gdalwarp_lib.cpp | 24 +++++++++++++++ autotest/alg/warp.py | 51 +++++++++++++++++++++++++----- 6 files changed, 144 insertions(+), 14 deletions(-) diff --git a/alg/gdalwarper.cpp b/alg/gdalwarper.cpp index b07664bf1eec..8cddbca37b20 100644 --- a/alg/gdalwarper.cpp +++ b/alg/gdalwarper.cpp @@ -1284,6 +1284,10 @@ CPLErr GDALWarpDstAlphaMasker(void *pMaskFuncArg, int nBandCount, * EXCLUDED_VALUES_PCT_THRESHOLD. * Only taken into account by Average currently. * + *
  • MODE_TIES=FIRST/MIN/MAX: (GDAL >= 3.11) Strategy to use when breaking + * ties with MODE resampling. By default, the first value encountered will be used. + * Alternatively, the minimum or maximum value can be selected.
    • + * * */ @@ -1305,6 +1309,7 @@ GDALWarpOptions *CPL_STDCALL GDALCreateWarpOptions() psOptions->eResampleAlg = GRA_NearestNeighbour; psOptions->pfnProgress = GDALDummyProgress; psOptions->eWorkingDataType = GDT_Unknown; + psOptions->eTieStrategy = GWKTS_First; return psOptions; } diff --git a/alg/gdalwarper.h b/alg/gdalwarper.h index 010fec7afb5f..15eceb561022 100644 --- a/alg/gdalwarper.h +++ b/alg/gdalwarper.h @@ -130,6 +130,14 @@ CPLErr CPL_DLL GDALWarpCutlineMaskerEx(void *pMaskFuncArg, int nBandCount, /*! @endcond */ +/*! GWKMode tie-breaking strategy */ +typedef enum +{ + /* Choose the first value encountered */ GWKTS_First = 1, + /* Choose the minimal value */ GWKTS_Min = 2, + /* Choose the maximum value */ GWKTS_Max = 3, +} GWKTieStrategy; + /************************************************************************/ /* GDALWarpOptions */ /************************************************************************/ @@ -243,6 +251,8 @@ typedef struct * zero. */ double dfCutlineBlendDist; + /** Tie-breaking method */ + GWKTieStrategy eTieStrategy; } GDALWarpOptions; GDALWarpOptions CPL_DLL *CPL_STDCALL GDALCreateWarpOptions(void); @@ -451,6 +461,8 @@ class CPL_DLL GDALWarpKernel // Average currently std::vector> m_aadfExcludedValues{}; + GWKTieStrategy eTieStrategy; + /*! @endcond */ GDALWarpKernel(); diff --git a/alg/gdalwarpkernel.cpp b/alg/gdalwarpkernel.cpp index e5ce23d4b43d..e917cbe24c04 100644 --- a/alg/gdalwarpkernel.cpp +++ b/alg/gdalwarpkernel.cpp @@ -989,7 +989,8 @@ GDALWarpKernel::GDALWarpKernel() nDstXOff(0), nDstYOff(0), pfnTransformer(nullptr), pTransformerArg(nullptr), pfnProgress(GDALDummyProgress), pProgress(nullptr), dfProgressBase(0.0), dfProgressScale(1.0), - padfDstNoDataReal(nullptr), psThreadData(nullptr) + padfDstNoDataReal(nullptr), psThreadData(nullptr), + eTieStrategy(GWKTS_First) { } @@ -6870,6 +6871,9 @@ static void GWKAverageOrModeThread(void *pData) // Only used with nAlgo = 6. float quant = 0.5; + // Only used for GRA_Mode + const GWKTieStrategy eTieStrategy = poWK->eTieStrategy; + // To control array allocation only when data type is complex const bool bIsComplex = GDALDataTypeIsComplex(poWK->eWorkingDataType) != 0; @@ -7599,13 +7603,44 @@ static void GWKAverageOrModeThread(void *pData) // Check array for existing entry. for (i = 0; i < iMaxInd; ++i) - if (pafRealVals[i] == fVal && - ++panRealSums[i] > - panRealSums[iMaxVal]) + { + if (pafRealVals[i] == fVal) { - iMaxVal = i; + bool bValIsMax = + (++panRealSums[i] > + panRealSums[iMaxVal]); + + if (!bValIsMax && + panRealSums[i] == + panRealSums[iMaxVal]) + { + switch (eTieStrategy) + { + case GWKTS_First: + break; + case GWKTS_Min: + bValIsMax = + fVal < + pafRealVals + [iMaxVal]; + break; + case GWKTS_Max: + bValIsMax = + fVal > + pafRealVals + [iMaxVal]; + break; + } + } + + if (bValIsMax) + { + iMaxVal = i; + } + break; } + } // Add to arr if entry not already there. if (i == iMaxInd) @@ -7676,7 +7711,25 @@ static void GWKAverageOrModeThread(void *pData) { const int nVal = static_cast(dfValueRealTmp); - if (++panVals[nVal + nBinsOffset] > nMaxVal) + + bool bValIsMax = + ++panVals[nVal + nBinsOffset] > nMaxVal; + if (!bValIsMax && + panVals[nVal + nBinsOffset] == nMaxVal) + { + switch (eTieStrategy) + { + case GWKTS_First: + break; + case GWKTS_Min: + bValIsMax = nVal < iMaxInd; + break; + case GWKTS_Max: + bValIsMax = nVal > iMaxInd; + break; + } + } + if (bValIsMax) { // Sum the density. // Is it the most common value so far? diff --git a/alg/gdalwarpoperation.cpp b/alg/gdalwarpoperation.cpp index 35a71726faf9..471ecc50040e 100644 --- a/alg/gdalwarpoperation.cpp +++ b/alg/gdalwarpoperation.cpp @@ -1839,6 +1839,7 @@ CPLErr GDALWarpOperation::WarpRegionToBuffer( oWK.eResample = m_bIsTranslationOnPixelBoundaries ? GRA_NearestNeighbour : psOptions->eResampleAlg; + oWK.eTieStrategy = psOptions->eTieStrategy; oWK.nBands = psOptions->nBandCount; oWK.eWorkingDataType = psOptions->eWorkingDataType; diff --git a/apps/gdalwarp_lib.cpp b/apps/gdalwarp_lib.cpp index 802281bc8e67..ae451eadbae7 100644 --- a/apps/gdalwarp_lib.cpp +++ b/apps/gdalwarp_lib.cpp @@ -2946,6 +2946,30 @@ static GDALDatasetH GDALWarpDirect(const char *pszDest, GDALDatasetH hDstDS, psWO->hSrcDS = hWrkSrcDS; psWO->hDstDS = hDstDS; + if (const char *pszTieStrategy = + psOptions->aosWarpOptions.FetchNameValue("MODE_TIES")) + { + if (EQUAL(pszTieStrategy, "FIRST")) + { + psWO->eTieStrategy = GWKTS_First; + } + else if (EQUAL(pszTieStrategy, "MIN")) + { + psWO->eTieStrategy = GWKTS_Min; + } + else if (EQUAL(pszTieStrategy, "MAX")) + { + psWO->eTieStrategy = GWKTS_Max; + } + else + { + CPLError(CE_Failure, CPLE_IllegalArg, + "Unknown value of MODE_TIES: %s", pszTieStrategy); + GDALReleaseDataset(hDstDS); + return nullptr; + } + } + if (!bVRT) { if (psOptions->pfnProgress == GDALDummyProgress) diff --git a/autotest/alg/warp.py b/autotest/alg/warp.py index 1ef52f7cdcac..9eff6a367566 100755 --- a/autotest/alg/warp.py +++ b/autotest/alg/warp.py @@ -1620,21 +1620,56 @@ def test_warp_rms_2(): ) -def test_warp_mode_ties(): - # when performing mode resampling the result in case of a tie will be - # the first value identified as the mode in scanline processing +@pytest.mark.parametrize("tie_strategy", ("FIRST", "MIN", "MAX", "HOPE")) +@pytest.mark.parametrize("dtype", (gdal.GDT_Int16, gdal.GDT_Int32)) +def test_warp_mode_ties(tie_strategy, dtype): numpy = pytest.importorskip("numpy") - src_ds = gdal.GetDriverByName("MEM").Create("", 3, 3, 1, gdal.GDT_Int16) + # 1 and 5 are tied for the mode; 1 encountered first + src_ds = gdal.GetDriverByName("MEM").Create("", 3, 3, 1, dtype) src_ds.SetGeoTransform([1, 1, 0, 1, 0, 1]) src_ds.GetRasterBand(1).WriteArray(numpy.array([[1, 1, 1], [2, 3, 4], [5, 5, 5]])) - out_ds = gdal.Warp("", src_ds, format="MEM", resampleAlg="mode", xRes=3, yRes=3) - assert out_ds.GetRasterBand(1).ReadAsArray()[0, 0] == 1 + with gdaltest.disable_exceptions(): + out_ds = gdal.Warp( + "", + src_ds, + format="MEM", + resampleAlg="mode", + xRes=3, + yRes=3, + warpOptions={"MODE_TIES": tie_strategy}, + ) + + if tie_strategy == "HOPE": + assert out_ds is None + return + result = out_ds.GetRasterBand(1).ReadAsArray()[0, 0] + + if tie_strategy in ("FIRST", "MIN"): + assert result == 1 + else: + assert result == 5 + + # 1 and 5 are tied for the mode; 5 encountered first src_ds.GetRasterBand(1).WriteArray(numpy.array([[1, 5, 1], [2, 5, 4], [5, 1, 0]])) - out_ds = gdal.Warp("", src_ds, format="MEM", resampleAlg="mode", xRes=3, yRes=3) - assert out_ds.GetRasterBand(1).ReadAsArray()[0, 0] == 5 + out_ds = gdal.Warp( + "", + src_ds, + format="MEM", + resampleAlg="mode", + xRes=3, + yRes=3, + warpOptions={"MODE_TIES": tie_strategy}, + ) + + result = out_ds.GetRasterBand(1).ReadAsArray()[0, 0] + + if tie_strategy in ("FIRST", "MAX"): + assert result == 5 + else: + assert result == 1 ###############################################################################