Skip to content

Commit

Permalink
Warper: Add MODE_TIES warp option
Browse files Browse the repository at this point in the history
  • Loading branch information
@dbaston
dbaston committed Jan 15, 2025
1 parent b4ddb34 commit a803e49
Show file tree
Hide file tree
Showing 6 changed files with 144 additions and 14 deletions.
5 changes: 5 additions & 0 deletions alg/gdalwarper.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1284,6 +1284,10 @@ CPLErr GDALWarpDstAlphaMasker(void *pMaskFuncArg, int nBandCount,
* EXCLUDED_VALUES_PCT_THRESHOLD.
* Only taken into account by Average currently.</li>
*
* <li>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.</li>
*
* </ul>
*/

Expand All @@ -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;
}
Expand Down
12 changes: 12 additions & 0 deletions alg/gdalwarper.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 */
/************************************************************************/
Expand Down Expand Up @@ -243,6 +251,8 @@ typedef struct
* zero. */
double dfCutlineBlendDist;

/** Tie-breaking method */
GWKTieStrategy eTieStrategy;
} GDALWarpOptions;

GDALWarpOptions CPL_DLL *CPL_STDCALL GDALCreateWarpOptions(void);
Expand Down Expand Up @@ -451,6 +461,8 @@ class CPL_DLL GDALWarpKernel
// Average currently
std::vector<std::vector<double>> m_aadfExcludedValues{};

GWKTieStrategy eTieStrategy;

/*! @endcond */

GDALWarpKernel();
Expand Down
65 changes: 59 additions & 6 deletions alg/gdalwarpkernel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
{
}

Expand Down Expand Up @@ -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;

Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -7676,7 +7711,25 @@ static void GWKAverageOrModeThread(void *pData)
{
const int nVal =
static_cast<int>(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?
Expand Down
1 change: 1 addition & 0 deletions alg/gdalwarpoperation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down
24 changes: 24 additions & 0 deletions apps/gdalwarp_lib.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
51 changes: 43 additions & 8 deletions autotest/alg/warp.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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


###############################################################################
Expand Down

0 comments on commit a803e49

Please sign in to comment.