From 57d953eae65648484031ea1e8a274a2c8ecc2bc1 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Mon, 30 Sep 2024 14:09:58 +0200
Subject: [PATCH] VRTProcessedDataset: fix issue when computing RasterIO window
 on auxiliary datasets on right-most/bottom-most tiles

---
 autotest/gdrivers/vrtprocesseddataset.py   | 52 ++++++++++++++++++++++
 frmts/vrt/vrtprocesseddatasetfunctions.cpp |  6 ++-
 2 files changed, 56 insertions(+), 2 deletions(-)

diff --git a/autotest/gdrivers/vrtprocesseddataset.py b/autotest/gdrivers/vrtprocesseddataset.py
index b11483b16df0..fd64c14e180e 100755
--- a/autotest/gdrivers/vrtprocesseddataset.py
+++ b/autotest/gdrivers/vrtprocesseddataset.py
@@ -719,6 +719,58 @@ def test_vrtprocesseddataset_dehazing_different_resolution(tmp_vsimem):
     )
 
 
+###############################################################################
+# Test we properly request auxiliary datasets on the right-most/bottom-most
+# truncated tile
+
+
+def test_vrtprocesseddataset_dehazing_edge_effects(tmp_vsimem):
+
+    src_filename = str(tmp_vsimem / "src.tif")
+    src_ds = gdal.GetDriverByName("GTiff").Create(
+        src_filename,
+        257,
+        257,
+        1,
+        gdal.GDT_Byte,
+        ["TILED=YES", "BLOCKXSIZE=256", "BLOCKYSIZE=256"],
+    )
+    src_ds.GetRasterBand(1).Fill(10)
+    src_ds.SetGeoTransform([0, 1, 0, 0, 0, -1])
+    src_ds.Close()
+
+    gain_filename = str(tmp_vsimem / "gain.tif")
+    gain_ds = gdal.GetDriverByName("GTiff").Create(gain_filename, 1, 1)
+    gain_ds.GetRasterBand(1).Fill(2)
+    gain_ds.SetGeoTransform([0, 257, 0, 0, 0, -257])
+    gain_ds.Close()
+
+    offset_filename = str(tmp_vsimem / "offset.tif")
+    offset_ds = gdal.GetDriverByName("GTiff").Create(offset_filename, 1, 1)
+    offset_ds.GetRasterBand(1).Fill(3)
+    offset_ds.SetGeoTransform([0, 257, 0, 0, 0, -257])
+    offset_ds.Close()
+
+    ds = gdal.Open(
+        f"""<VRTDataset subclass='VRTProcessedDataset'>
+    <Input>
+        <SourceFilename>{src_filename}</SourceFilename>
+    </Input>
+    <ProcessingSteps>
+        <Step>
+            <Algorithm>LocalScaleOffset</Algorithm>
+            <Argument name="gain_dataset_filename_1">{gain_filename}</Argument>
+            <Argument name="gain_dataset_band_1">1</Argument>
+            <Argument name="offset_dataset_filename_1">{offset_filename}</Argument>
+            <Argument name="offset_dataset_band_1">1</Argument>
+        </Step>
+    </ProcessingSteps>
+    </VRTDataset>
+        """
+    )
+    assert ds.GetRasterBand(1).ComputeRasterMinMax() == (17, 17)
+
+
 ###############################################################################
 # Test error cases of LocalScaleOffset algorithm
 
diff --git a/frmts/vrt/vrtprocesseddatasetfunctions.cpp b/frmts/vrt/vrtprocesseddatasetfunctions.cpp
index a34aca2919b9..e5e5d86c25cc 100644
--- a/frmts/vrt/vrtprocesseddatasetfunctions.cpp
+++ b/frmts/vrt/vrtprocesseddatasetfunctions.cpp
@@ -963,8 +963,10 @@ static bool LoadAuxData(double dfULX, double dfULY, double dfLRX, double dfLRY,
         return false;
     }
 
-    const int nAuxXOff = std::max(0, static_cast<int>(std::round(dfULPixel)));
-    const int nAuxYOff = std::max(0, static_cast<int>(std::round(dfULLine)));
+    const int nAuxXOff = std::clamp(static_cast<int>(std::round(dfULPixel)), 0,
+                                    poAuxBand->GetXSize() - 1);
+    const int nAuxYOff = std::clamp(static_cast<int>(std::round(dfULLine)), 0,
+                                    poAuxBand->GetYSize() - 1);
     const int nAuxX2Off = std::min(poAuxBand->GetXSize(),
                                    static_cast<int>(std::round(dfLRPixel)));
     const int nAuxY2Off =