From 7ec0dca84982048203a20123d60a971ca4579edc Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 26 Sep 2024 18:06:30 +0200 Subject: [PATCH] gdalwarp: be more robust to numerical instability when selecting overviews Fixes #10873 Uniformize logic with the one of GDALBandGetBestOverviewLevel2() --- apps/gdalwarp_lib.cpp | 42 ++++++++++++++++++++--------- autotest/utilities/test_gdalwarp.py | 19 +++++++++++++ gcore/rasterio.cpp | 24 ++++++++++------- 3 files changed, 63 insertions(+), 22 deletions(-) diff --git a/apps/gdalwarp_lib.cpp b/apps/gdalwarp_lib.cpp index 53fa14271b03..ed0345b01dee 100644 --- a/apps/gdalwarp_lib.cpp +++ b/apps/gdalwarp_lib.cpp @@ -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 @@ -2827,18 +2836,27 @@ static GDALDatasetH GDALWarpDirect(const char *pszDest, GDALDatasetH hDstDS, poSrcDS->GetRasterBand(1) ->GetOverview(iOvr) ->GetXSize(); - const double dfNextOvrRatio = - static_cast(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, diff --git a/autotest/utilities/test_gdalwarp.py b/autotest/utilities/test_gdalwarp.py index ee749fed657d..5b09b8fc0fdf 100755 --- a/autotest/utilities/test_gdalwarp.py +++ b/autotest/utilities/test_gdalwarp.py @@ -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") diff --git a/gcore/rasterio.cpp b/gcore/rasterio.cpp index 8f6a99ee48eb..995a4318c99a 100644 --- a/gcore/rasterio.cpp +++ b/gcore/rasterio.cpp @@ -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); @@ -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; @@ -3704,6 +3702,12 @@ int GDALBandGetBestOverviewLevel2(GDALRasterBand *poBand, int &nXOff, poBestOverview = poOverview; nBestOverviewLevel = iOverview; dfBestDownsamplingFactor = dfDownsamplingFactor; + + if (std::abs(dfDesiredDownsamplingFactor - dfDownsamplingFactor) < + EPSILON) + { + break; + } } /* -------------------------------------------------------------------- */