From 9a0494a93461e852d1bcf77a7678374f8715b701 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Thu, 18 Jan 2024 15:02:49 +0100
Subject: [PATCH] GDALBandGetBestOverviewLevel2(): avoid numerical instability
 with GDAL_OVERVIEW_OVERSAMPLING_THRESHOLD = 1.0

Fixes https://github.com/OSGeo/gdal/pull/9040#issuecomment-1898524693
---
 gcore/rasterio.cpp | 20 +++++++++++++-------
 1 file changed, 13 insertions(+), 7 deletions(-)

diff --git a/gcore/rasterio.cpp b/gcore/rasterio.cpp
index fe1500da201c..b3c7f22d29fb 100644
--- a/gcore/rasterio.cpp
+++ b/gcore/rasterio.cpp
@@ -3602,11 +3602,20 @@ int GDALBandGetBestOverviewLevel2(GDALRasterBand *poBand, int &nXOff,
 
     const char *pszOversampligThreshold =
         CPLGetConfigOption("GDAL_OVERVIEW_OVERSAMPLING_THRESHOLD", nullptr);
-    const double dfOversamplingThreshold =
+
+    // 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(
         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);
@@ -3624,11 +3633,8 @@ int GDALBandGetBestOverviewLevel2(GDALRasterBand *poBand, int &nXOff,
 
         // Is it nearly the requested factor and better (lower) than
         // the current best factor?
-        if ((dfOversamplingThreshold == 1.0 &&
-             dfDownsamplingFactor > dfDesiredDownsamplingFactor) ||
-            (dfOversamplingThreshold > 1.0 &&
-             dfDownsamplingFactor >=
-                 dfDesiredDownsamplingFactor * dfOversamplingThreshold) ||
+        if (dfDownsamplingFactor >=
+                dfDesiredDownsamplingFactor * dfOversamplingThreshold ||
             dfDownsamplingFactor <= dfBestDownsamplingFactor)
         {
             continue;