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

gdalwarp: be more robust to numerical instability when selecting overviews #10879

Merged
merged 1 commit into from
Oct 4, 2024
Merged

gdalwarp: be more robust to numerical instability when selecting overviews #10879

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
42 changes: 30 additions & 12 deletions apps/gdalwarp_lib.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2817,8 +2817,17 @@ static GDALDatasetH GDALWarpDirect(const char *pszDest, GDALDatasetH hDstDS,

if (dfTargetRatio > 1.0)
{
int iOvr = -1;
for (; iOvr < nOvCount - 1; iOvr++)
// Note: keep this logic for overview selection in sync between
// gdalwarp_lib.cpp and rasterio.cpp
const char *pszOversampligThreshold = CPLGetConfigOption(
"GDALWARP_OVERSAMPLING_THRESHOLD", nullptr);
const double dfOversamplingThreshold =
pszOversampligThreshold ? CPLAtof(pszOversampligThreshold)
: 1.0;

int iBestOvr = -1;
double dfBestRatio = 0;
for (int iOvr = -1; iOvr < nOvCount; iOvr++)
{
const double dfOvrRatio =
iOvr < 0
Expand All @@ -2827,18 +2836,27 @@ static GDALDatasetH GDALWarpDirect(const char *pszDest, GDALDatasetH hDstDS,
poSrcDS->GetRasterBand(1)
->GetOverview(iOvr)
->GetXSize();
const double dfNextOvrRatio =
static_cast<double>(poSrcDS->GetRasterXSize()) /
poSrcDS->GetRasterBand(1)
->GetOverview(iOvr + 1)
->GetXSize();
if (dfOvrRatio < dfTargetRatio &&
dfNextOvrRatio > dfTargetRatio)
break;
if (fabs(dfOvrRatio - dfTargetRatio) < 1e-1)

// Is it nearly the requested factor and better (lower) than
// the current best factor?
// Use an epsilon because of numerical instability.
constexpr double EPSILON = 1e-1;
if (dfOvrRatio >=
dfTargetRatio * dfOversamplingThreshold + EPSILON ||
dfOvrRatio <= dfBestRatio)
{
continue;
}

iBestOvr = iOvr;
dfBestRatio = dfOvrRatio;
if (std::abs(dfTargetRatio - dfOvrRatio) < EPSILON)
{
break;
}
}
iOvr += (psOptions->nOvLevel - OVR_LEVEL_AUTO);
const int iOvr =
iBestOvr + (psOptions->nOvLevel - OVR_LEVEL_AUTO);
if (iOvr >= 0)
{
CPLDebug("WARP", "Selecting overview level %d for %s", iOvr,
Expand Down
4 changes: 2 additions & 2 deletions autotest/gcore/rasterio.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3318,11 +3318,11 @@ def test_rasterio_overview_selection():
assert ds.GetRasterBand(1).ReadRaster(0, 0, 100, 100, 101, 101)[0] == 1
assert ds.GetRasterBand(1).ReadRaster(0, 0, 100, 100, 100, 100)[0] == 1
assert ds.GetRasterBand(1).ReadRaster(0, 0, 100, 100, 99, 99)[0] == 1
assert ds.GetRasterBand(1).ReadRaster(0, 0, 100, 100, 60, 60)[0] == 1
assert ds.GetRasterBand(1).ReadRaster(0, 0, 100, 100, 60, 60)[0] == 2
assert ds.GetRasterBand(1).ReadRaster(0, 0, 100, 100, 59, 59)[0] == 2
assert ds.GetRasterBand(1).ReadRaster(0, 0, 100, 100, 50, 50)[0] == 2
assert ds.GetRasterBand(1).ReadRaster(0, 0, 100, 100, 49, 49)[0] == 2
assert ds.GetRasterBand(1).ReadRaster(0, 0, 100, 100, 30, 30)[0] == 2
assert ds.GetRasterBand(1).ReadRaster(0, 0, 100, 100, 30, 30)[0] == 3
assert ds.GetRasterBand(1).ReadRaster(0, 0, 100, 100, 29, 29)[0] == 3
assert ds.GetRasterBand(1).ReadRaster(0, 0, 100, 100, 25, 25)[0] == 3
assert ds.GetRasterBand(1).ReadRaster(0, 0, 100, 100, 24, 24)[0] == 3
Expand Down
19 changes: 19 additions & 0 deletions autotest/utilities/test_gdalwarp.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1164,6 +1164,25 @@ def test_gdalwarp_40(gdalwarp_path, tmp_path):
expected_cs = ds.GetRasterBand(1).Checksum()
ds = None

# Test that tiny variations in -te that result in a target resampling factor
# very close to the one of overview 0 lead to overview 0 been selected

gdaltest.runexternal(
f"{gdalwarp_path} {src_tif} {dst_vrt} -overwrite -ts 10 10 -te 440721 3750120 441920 3751320 -of VRT"
)

ds = gdal.Open(dst_vrt)
assert ds.GetRasterBand(1).Checksum() == cs_ov0
ds = None

gdaltest.runexternal(
f"{gdalwarp_path} {src_tif} {dst_vrt} -overwrite -ts 10 10 -te 440719 3750120 441920 3751320 -of VRT"
)

ds = gdal.Open(dst_vrt)
assert ds.GetRasterBand(1).Checksum() == cs_ov0
ds = None

# Should select overview 0 too
gdaltest.runexternal(f"{gdalwarp_path} {src_tif} {dst_tif} -overwrite -ts 7 7")

Expand Down
24 changes: 14 additions & 10 deletions gcore/rasterio.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3656,19 +3656,14 @@ int GDALBandGetBestOverviewLevel2(GDALRasterBand *poBand, int &nXOff,
const char *pszOversampligThreshold =
CPLGetConfigOption("GDAL_OVERVIEW_OVERSAMPLING_THRESHOLD", nullptr);

// Note: keep this logic for overview selection in sync between
// gdalwarp_lib.cpp and rasterio.cpp
// Cf https://github.com/OSGeo/gdal/pull/9040#issuecomment-1898524693
// Do not exactly use a oversampling threshold of 1.0 because of numerical
// instability.
const auto AdjustThreshold = [](double x)
{
constexpr double EPS = 1e-2;
return x == 1.0 ? x + EPS : x;
};
const double dfOversamplingThreshold = AdjustThreshold(
const double dfOversamplingThreshold =
pszOversampligThreshold ? CPLAtof(pszOversampligThreshold)
: psExtraArg && psExtraArg->eResampleAlg != GRIORA_NearestNeighbour
? 1.0
: 1.2);
: 1.2;
for (int iOverview = 0; iOverview < nOverviewCount; iOverview++)
{
GDALRasterBand *poOverview = poBand->GetOverview(iOverview);
Expand All @@ -3686,8 +3681,11 @@ int GDALBandGetBestOverviewLevel2(GDALRasterBand *poBand, int &nXOff,

// Is it nearly the requested factor and better (lower) than
// the current best factor?
// Use an epsilon because of numerical instability.
constexpr double EPSILON = 1e-1;
if (dfDownsamplingFactor >=
dfDesiredDownsamplingFactor * dfOversamplingThreshold ||
dfDesiredDownsamplingFactor * dfOversamplingThreshold +
EPSILON ||
dfDownsamplingFactor <= dfBestDownsamplingFactor)
{
continue;
Expand All @@ -3704,6 +3702,12 @@ int GDALBandGetBestOverviewLevel2(GDALRasterBand *poBand, int &nXOff,
poBestOverview = poOverview;
nBestOverviewLevel = iOverview;
dfBestDownsamplingFactor = dfDownsamplingFactor;

if (std::abs(dfDesiredDownsamplingFactor - dfDownsamplingFactor) <
EPSILON)
{
break;
}
}

/* -------------------------------------------------------------------- */
Expand Down
Loading