Skip to content

Commit

Permalink
Merge pull request #11661 from dbaston/warp-mode-fixes
Browse files Browse the repository at this point in the history 
Warper: Fixes related to mode resampling
  • Loading branch information
@rouault
rouault authored Jan 16, 2025
2 parents a16d29c + b075d49 commit b68bf74
Show file tree
Hide file tree
Showing 5 changed files with 229 additions and 59 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
153 changes: 104 additions & 49 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 @@ -6856,16 +6857,12 @@ 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;
Expand Down Expand Up @@ -6913,9 +6910,9 @@ static void GWKAverageOrModeThread(void *pData)
{
nBins = 65536;
}
panVals =
static_cast<int *>(VSI_MALLOC_VERBOSE(nBins * sizeof(int)));
if (panVals == nullptr)
pafCounts =
static_cast<float *>(VSI_MALLOC_VERBOSE(nBins * sizeof(float)));
if (pafCounts == nullptr)
return;
}
else
Expand All @@ -6926,12 +6923,12 @@ static void GWKAverageOrModeThread(void *pData)
{
pafRealVals = static_cast<float *>(
VSI_MALLOC3_VERBOSE(nSrcXSize, nSrcYSize, sizeof(float)));
panRealSums = static_cast<int *>(
VSI_MALLOC3_VERBOSE(nSrcXSize, nSrcYSize, sizeof(int)));
if (pafRealVals == nullptr || panRealSums == nullptr)
pafCounts = static_cast<float *>(
VSI_MALLOC3_VERBOSE(nSrcXSize, nSrcYSize, sizeof(float)));
if (pafRealVals == nullptr || pafCounts == nullptr)
{
VSIFree(pafRealVals);
VSIFree(panRealSums);
VSIFree(pafCounts);
return;
}
}
Expand Down Expand Up @@ -7566,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<GPtrDiff_t>(iSrcY) * nSrcXSize;
Expand All @@ -7596,35 +7593,73 @@ static void GWKAverageOrModeThread(void *pData)
{
const float fVal =
static_cast<float>(dfValueRealTmp);
const double dfWeight =
COMPUTE_WEIGHT(iSrcX, dfWeightY);

// Check array for existing entry.
for (i = 0; i < iMaxInd; ++i)
if (pafRealVals[i] == fVal &&
++panRealSums[i] >
panRealSums[iMaxVal])
int i = 0;
for (i = 0; i < nBins; ++i)
{
if (pafRealVals[i] == fVal)
{
iMaxVal = i;

pafCounts[i] +=
static_cast<float>(dfWeight);
bool bValIsMaxCount =
(pafCounts[i] >
pafCounts[iModeIndex]);

if (!bValIsMaxCount &&
pafCounts[i] ==
pafCounts[iModeIndex])
{
switch (eTieStrategy)
{
case GWKTS_First:
break;
case GWKTS_Min:
bValIsMaxCount =
fVal <
pafRealVals
[iModeIndex];
break;
case GWKTS_Max:
bValIsMaxCount =
fVal >
pafRealVals
[iModeIndex];
break;
}
}

if (bValIsMaxCount)
{
iModeIndex = i;
}

break;
}
}

// 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<float>(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)
{
Expand All @@ -7645,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<GPtrDiff_t>(iSrcY) * nSrcXSize;
Expand All @@ -7674,22 +7711,46 @@ static void GWKAverageOrModeThread(void *pData)
&dfValueRealTmp, &dfValueImagTmp) &&
dfBandDensity > BAND_DENSITY_THRESHOLD)
{
bHasSourceValues = true;
const int nVal =
static_cast<int>(dfValueRealTmp);
if (++panVals[nVal + nBinsOffset] > nMaxVal)
const int iBin = nVal + nBinsOffset;
const double dfWeight =
COMPUTE_WEIGHT(iSrcX, dfWeightY);

// Sum the density.
pafCounts[iBin] +=
static_cast<float>(dfWeight);
// Is it the most common value so far?
bool bUpdateMode =
pafCounts[iBin] > fMaxCount;
if (!bUpdateMode &&
pafCounts[iBin] == fMaxCount)
{
// Sum the density.
// Is it the most common value so far?
iMaxInd = nVal;
nMaxVal = panVals[nVal + nBinsOffset];
switch (eTieStrategy)
{
case GWKTS_First:
break;
case GWKTS_Min:
bUpdateMode = nVal < nMode;
break;
case GWKTS_Max:
bUpdateMode = nVal > nMode;
break;
}
}
if (bUpdateMode)
{
nMode = nVal;
fMaxCount = pafCounts[iBin];
}
}
}
}

if (iMaxInd != -1)
if (bHasSourceValues)
{
dfValueReal = iMaxInd;
dfValueReal = nMode;

if (poWK->bApplyVerticalShift)
{
Expand Down Expand Up @@ -7950,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);
}
}

/************************************************************************/
Expand Down
37 changes: 35 additions & 2 deletions alg/gdalwarpoperation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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() */
/************************************************************************/
Expand All @@ -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
Expand All @@ -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",
Expand All @@ -523,6 +555,7 @@ CPLErr GDALWarpOperation::Initialize(const GDALWarpOptions *psNewOptions)
}

GDALWarpResolveWorkingDataType(psOptions);
SetTieStrategy(psOptions, &eErr);

/* -------------------------------------------------------------------- */
/* Default memory available. */
Expand All @@ -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<char *>(pszCutlineWKT);
Expand Down Expand Up @@ -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;

Expand Down
Loading

1 comment on commit b68bf74

@schwehr
Copy link
Member

Choose a reason for hiding this comment

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

rasterio/rasterio#3293 to handle a new value for mode in test_reproject_resampling

Please sign in to comment.