From 94e69abd2e96e8f0dc333be9fe42509da4cc98f2 Mon Sep 17 00:00:00 2001 From: Dan Baston Date: Mon, 13 Jan 2025 20:26:40 -0500 Subject: [PATCH 1/2] Warper: Add MODE_TIES warp option --- alg/gdalwarper.cpp | 5 +++ alg/gdalwarper.h | 12 ++++++++ alg/gdalwarpkernel.cpp | 65 +++++++++++++++++++++++++++++++++++---- alg/gdalwarpoperation.cpp | 37 ++++++++++++++++++++-- autotest/alg/warp.py | 51 +++++++++++++++++++++++++----- 5 files changed, 154 insertions(+), 16 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..f4997c01b16a 100644 --- a/alg/gdalwarpoperation.cpp +++ b/alg/gdalwarpoperation.cpp @@ -473,6 +473,36 @@ static void SetAlphaMax(GDALWarpOptions *psOptions, GDALRasterBandH hBand, CPLDebug("WARP", "SetAlphaMax: AlphaMax not set."); } +/************************************************************************/ +/* SetTieStrategy() */ +/************************************************************************/ + +static void SetTieStrategy(GDALWarpOptions *psOptions, CPLErr *peErr) +{ + if (const char *pszTieStrategy = + CSLFetchNameValue(psOptions->papszWarpOptions, "MODE_TIES")) + { + if (EQUAL(pszTieStrategy, "FIRST")) + { + psOptions->eTieStrategy = GWKTS_First; + } + else if (EQUAL(pszTieStrategy, "MIN")) + { + psOptions->eTieStrategy = GWKTS_Min; + } + else if (EQUAL(pszTieStrategy, "MAX")) + { + psOptions->eTieStrategy = GWKTS_Max; + } + else + { + CPLError(CE_Failure, CPLE_IllegalArg, + "Unknown value of MODE_TIES: %s", pszTieStrategy); + *peErr = CE_Failure; + } + } +} + /************************************************************************/ /* Initialize() */ /************************************************************************/ @@ -483,7 +513,7 @@ static void SetAlphaMax(GDALWarpOptions *psOptions, GDALRasterBandH hBand, * This method initializes the GDALWarpOperation's concept of the warp * options in effect. It creates an internal copy of the GDALWarpOptions * structure and defaults a variety of additional fields in the internal - * copy if not set in the provides warp options. + * copy if not set in the provided warp options. * * Defaulting operations include: * - If the nBandCount is 0, it will be set to the number of bands in the @@ -505,6 +535,8 @@ CPLErr GDALWarpOperation::Initialize(const GDALWarpOptions *psNewOptions) if (psOptions != nullptr) WipeOptions(); + CPLErr eErr = CE_None; + psOptions = GDALCloneWarpOptions(psNewOptions); psOptions->papszWarpOptions = CSLSetNameValue(psOptions->papszWarpOptions, "EXTRA_ELTS", @@ -523,6 +555,7 @@ CPLErr GDALWarpOperation::Initialize(const GDALWarpOptions *psNewOptions) } GDALWarpResolveWorkingDataType(psOptions); + SetTieStrategy(psOptions, &eErr); /* -------------------------------------------------------------------- */ /* Default memory available. */ @@ -548,7 +581,6 @@ CPLErr GDALWarpOperation::Initialize(const GDALWarpOptions *psNewOptions) const char *pszCutlineWKT = CSLFetchNameValue(psOptions->papszWarpOptions, "CUTLINE"); - CPLErr eErr = CE_None; if (pszCutlineWKT && psOptions->hCutline == nullptr) { char *pszWKTTmp = const_cast(pszCutlineWKT); @@ -1839,6 +1871,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/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 ############################################################################### From b075d4914360c21fc674d3afa97db3dae0236220 Mon Sep 17 00:00:00 2001 From: Daniel Baston Date: Wed, 15 Jan 2025 10:22:16 -0500 Subject: [PATCH 2/2] Warper: Mode resampling: Use src pixel coverage and allow mode of -1 --- alg/gdalwarpkernel.cpp | 134 +++++++++++++++++++++-------------------- autotest/alg/warp.py | 30 +++++++++ 2 files changed, 98 insertions(+), 66 deletions(-) diff --git a/alg/gdalwarpkernel.cpp b/alg/gdalwarpkernel.cpp index e917cbe24c04..c233d84ae107 100644 --- a/alg/gdalwarpkernel.cpp +++ b/alg/gdalwarpkernel.cpp @@ -6857,23 +6857,16 @@ static void GWKAverageOrModeThread(void *pData) /* -------------------------------------------------------------------- */ int nAlgo = 0; - // These vars only used with nAlgo == 3. - int *panVals = nullptr; + // Only used for GRA_Mode + float *pafRealVals = nullptr; + float *pafCounts = nullptr; int nBins = 0; int nBinsOffset = 0; - - // Only used with nAlgo = 2. - float *pafRealVals = nullptr; - float *pafImagVals = nullptr; - int *panRealSums = nullptr; - int *panImagSums = nullptr; + const GWKTieStrategy eTieStrategy = poWK->eTieStrategy; // 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; @@ -6917,9 +6910,9 @@ static void GWKAverageOrModeThread(void *pData) { nBins = 65536; } - panVals = - static_cast(VSI_MALLOC_VERBOSE(nBins * sizeof(int))); - if (panVals == nullptr) + pafCounts = + static_cast(VSI_MALLOC_VERBOSE(nBins * sizeof(float))); + if (pafCounts == nullptr) return; } else @@ -6930,12 +6923,12 @@ static void GWKAverageOrModeThread(void *pData) { pafRealVals = static_cast( VSI_MALLOC3_VERBOSE(nSrcXSize, nSrcYSize, sizeof(float))); - panRealSums = static_cast( - VSI_MALLOC3_VERBOSE(nSrcXSize, nSrcYSize, sizeof(int))); - if (pafRealVals == nullptr || panRealSums == nullptr) + pafCounts = static_cast( + VSI_MALLOC3_VERBOSE(nSrcXSize, nSrcYSize, sizeof(float))); + if (pafRealVals == nullptr || pafCounts == nullptr) { VSIFree(pafRealVals); - VSIFree(panRealSums); + VSIFree(pafCounts); return; } } @@ -7570,12 +7563,12 @@ static void GWKAverageOrModeThread(void *pData) // majority filter on floating point data? But, here it // is for the sake of compatibility. It won't look // right on RGB images by the nature of the filter. - int iMaxInd = 0; - int iMaxVal = -1; - int i = 0; + nBins = 0; + int iModeIndex = -1; for (int iSrcY = iSrcYMin; iSrcY < iSrcYMax; iSrcY++) { + const double dfWeightY = COMPUTE_WEIGHT_Y(iSrcY); iSrcOffset = iSrcXMin + static_cast(iSrcY) * nSrcXSize; @@ -7600,42 +7593,48 @@ static void GWKAverageOrModeThread(void *pData) { const float fVal = static_cast(dfValueRealTmp); + const double dfWeight = + COMPUTE_WEIGHT(iSrcX, dfWeightY); // Check array for existing entry. - for (i = 0; i < iMaxInd; ++i) + int i = 0; + for (i = 0; i < nBins; ++i) { if (pafRealVals[i] == fVal) { - bool bValIsMax = - (++panRealSums[i] > - panRealSums[iMaxVal]); - if (!bValIsMax && - panRealSums[i] == - panRealSums[iMaxVal]) + pafCounts[i] += + static_cast(dfWeight); + bool bValIsMaxCount = + (pafCounts[i] > + pafCounts[iModeIndex]); + + if (!bValIsMaxCount && + pafCounts[i] == + pafCounts[iModeIndex]) { switch (eTieStrategy) { case GWKTS_First: break; case GWKTS_Min: - bValIsMax = + bValIsMaxCount = fVal < pafRealVals - [iMaxVal]; + [iModeIndex]; break; case GWKTS_Max: - bValIsMax = + bValIsMaxCount = fVal > pafRealVals - [iMaxVal]; + [iModeIndex]; break; } } - if (bValIsMax) + if (bValIsMaxCount) { - iMaxVal = i; + iModeIndex = i; } break; @@ -7643,23 +7642,24 @@ static void GWKAverageOrModeThread(void *pData) } // Add to arr if entry not already there. - if (i == iMaxInd) + if (i == nBins) { - pafRealVals[iMaxInd] = fVal; - panRealSums[iMaxInd] = 1; + pafRealVals[i] = fVal; + pafCounts[i] = + static_cast(dfWeight); - if (iMaxVal < 0) - iMaxVal = iMaxInd; + if (iModeIndex < 0) + iModeIndex = i; - ++iMaxInd; + ++nBins; } } } } - if (iMaxVal != -1) + if (iModeIndex != -1) { - dfValueReal = pafRealVals[iMaxVal]; + dfValueReal = pafRealVals[iModeIndex]; if (poWK->bApplyVerticalShift) { @@ -7680,13 +7680,15 @@ static void GWKAverageOrModeThread(void *pData) } else // byte or int16. { - int nMaxVal = 0; - int iMaxInd = -1; + float fMaxCount = 0.0f; + int nMode = -1; + bool bHasSourceValues = false; - memset(panVals, 0, nBins * sizeof(int)); + memset(pafCounts, 0, nBins * sizeof(float)); for (int iSrcY = iSrcYMin; iSrcY < iSrcYMax; iSrcY++) { + const double dfWeightY = COMPUTE_WEIGHT_Y(iSrcY); iSrcOffset = iSrcXMin + static_cast(iSrcY) * nSrcXSize; @@ -7709,40 +7711,46 @@ static void GWKAverageOrModeThread(void *pData) &dfValueRealTmp, &dfValueImagTmp) && dfBandDensity > BAND_DENSITY_THRESHOLD) { + bHasSourceValues = true; const int nVal = static_cast(dfValueRealTmp); - - bool bValIsMax = - ++panVals[nVal + nBinsOffset] > nMaxVal; - if (!bValIsMax && - panVals[nVal + nBinsOffset] == nMaxVal) + const int iBin = nVal + nBinsOffset; + const double dfWeight = + COMPUTE_WEIGHT(iSrcX, dfWeightY); + + // Sum the density. + pafCounts[iBin] += + static_cast(dfWeight); + // Is it the most common value so far? + bool bUpdateMode = + pafCounts[iBin] > fMaxCount; + if (!bUpdateMode && + pafCounts[iBin] == fMaxCount) { switch (eTieStrategy) { case GWKTS_First: break; case GWKTS_Min: - bValIsMax = nVal < iMaxInd; + bUpdateMode = nVal < nMode; break; case GWKTS_Max: - bValIsMax = nVal > iMaxInd; + bUpdateMode = nVal > nMode; break; } } - if (bValIsMax) + if (bUpdateMode) { - // Sum the density. - // Is it the most common value so far? - iMaxInd = nVal; - nMaxVal = panVals[nVal + nBinsOffset]; + nMode = nVal; + fMaxCount = pafCounts[iBin]; } } } } - if (iMaxInd != -1) + if (bHasSourceValues) { - dfValueReal = iMaxInd; + dfValueReal = nMode; if (poWK->bApplyVerticalShift) { @@ -8003,14 +8011,8 @@ static void GWKAverageOrModeThread(void *pData) CPLFree(padfZ2); CPLFree(pabSuccess); CPLFree(pabSuccess2); - VSIFree(panVals); + VSIFree(pafCounts); VSIFree(pafRealVals); - VSIFree(panRealSums); - if (bIsComplex) - { - VSIFree(pafImagVals); - VSIFree(panImagSums); - } } /************************************************************************/ diff --git a/autotest/alg/warp.py b/autotest/alg/warp.py index 9eff6a367566..24158ceaec6a 100755 --- a/autotest/alg/warp.py +++ b/autotest/alg/warp.py @@ -1160,6 +1160,36 @@ def test_warp_weighted_average(): assert maxdiff <= 1, "Image too different from reference" +############################################################################### +# test weighted mode + + +@pytest.mark.parametrize( + "dtype", (gdal.GDT_Int16, gdal.GDT_Int32), ids=gdal.GetDataTypeName +) +def test_warp_weighted_mode(dtype): + + np = pytest.importorskip("numpy") + + with gdal.GetDriverByName("MEM").Create("", 3, 3, eType=dtype) as src_ds: + src_ds.SetGeoTransform([0, 1, 0, 3, 0, -1]) + src_ds.WriteArray(np.array([[1, 1, 1], [-1, -1, -1], [3, 3, 3]])) + + dst_ds = gdal.Warp( + "", + src_ds, + format="MEM", + resampleAlg="mode", + width=1, + height=1, + outputBounds=(0.5, 0.5, 2.5, 2.5), + ) + + result = dst_ds.ReadAsArray()[0, 0] + + assert result == -1 + + ############################################################################### # test weighted average, with src offset (fix for #2665)