Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Warper: Add TIE_STRATEGY warp option #11649

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
+142 −12
Open

Warper: Add TIE_STRATEGY warp option #11649

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions alg/gdalwarper.cpp
View file
Open in desktop
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>TIE_STRATEGY=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
Open in desktop
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
Open in desktop
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;
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@rouault probably not worth backporting, but this is doing extra work by not hitting break on pafRealVals[i] == fVal regardless of whether ++panRealSums[i] > panRealSums[iMaxVal]

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 =
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I know this isn't due to your changes, but the current variable naming is terrible and makes it hard to follow the algorithm. Would you mind renaming:

  • panVals to something like panCount
  • nMaxVal to nMaxCount
  • iMaxInd to nModeValCandidate ?
  • bValIsMax to bValIsMaxCount

I believe panRealSums could be removed, and just replaced by panCount shared between int and float implementations

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, the variable names caused me some confusion when working on this. I'll try to improve them in a separate commit.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually, I was surprised to see that "mode" resampling isn't accounting for pixel coverage fractions in the way that "sum" and "average" are. I'm planning to address this as I rename the variables.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm planning to address this as I rename the variables.

yes, that would definitely make sense for "mode" as well

++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
Open in desktop
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
Open in desktop
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("TIE_STRATEGY"))
{
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 TIE_STRATEGY: %s", pszTieStrategy);
GDALReleaseDataset(hDstDS);
return nullptr;
}
}

if (!bVRT)
{
if (psOptions->pfnProgress == GDALDummyProgress)
Expand Down
47 changes: 41 additions & 6 deletions autotest/alg/warp.py
View file
Open in desktop
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():
@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):
# when performing mode resampling the result in case of a tie will be
# the first value identified as the mode in scanline processing
numpy = pytest.importorskip("numpy")

src_ds = gdal.GetDriverByName("MEM").Create("", 3, 3, 1, gdal.GDT_Int16)
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={"TIE_STRATEGY": 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

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={"TIE_STRATEGY": 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
Loading