GTiff: endless loop when resampling (bilinear) under some circumstances

[elpaso]

What is the bug?

Discovered downstream in qgis/QGIS#56769

Steps to reproduce the issue

Preliminary test case to reproduce (test.tif file is from downstream QGIS issue mentioned above).

$ gdalinfo test.tif 
Driver: GTiff/GeoTIFF
Files: test.tif
Size is 980, 1012
Coordinate System is:
PROJCRS["ETRS89 / UTM zone 32N",
    BASEGEOGCRS["ETRS89",
        ENSEMBLE["European Terrestrial Reference System 1989 ensemble",
            MEMBER["European Terrestrial Reference Frame 1989"],
            MEMBER["European Terrestrial Reference Frame 1990"],
            MEMBER["European Terrestrial Reference Frame 1991"],
            MEMBER["European Terrestrial Reference Frame 1992"],
            MEMBER["European Terrestrial Reference Frame 1993"],
            MEMBER["European Terrestrial Reference Frame 1994"],
            MEMBER["European Terrestrial Reference Frame 1996"],
            MEMBER["European Terrestrial Reference Frame 1997"],
            MEMBER["European Terrestrial Reference Frame 2000"],
            MEMBER["European Terrestrial Reference Frame 2005"],
            MEMBER["European Terrestrial Reference Frame 2014"],
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[0.1]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4258]],
    CONVERSION["UTM zone 32N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",9,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Europe between 6°E and 12°E: Austria; Belgium; Denmark - onshore and offshore; Germany - onshore and offshore; Norway including - onshore and offshore; Spain - offshore."],
        BBOX[38.76,6,84.33,12]],
    ID["EPSG",25832]]
Data axis to CRS axis mapping: 1,2
Origin = (305000.500000000000000,5681012.500000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  305000.500, 5681012.500) (  6d12'21.32"E, 51d14'49.54"N)
Lower Left  (  305000.500, 5680000.500) (  6d12'23.30"E, 51d14'16.82"N)
Upper Right (  305980.500, 5681012.500) (  6d13'11.80"E, 51d14'50.74"N)
Lower Right (  305980.500, 5680000.500) (  6d13'13.78"E, 51d14'18.02"N)
Center      (  305490.500, 5680506.500) (  6d12'47.55"E, 51d14'33.78"N)
Band 1 Block=980x2 Type=UInt32, ColorInterp=Gray
  NoData Value=-3.4028234663852886e+38
  Metadata:
    RepresentationType=ATHEMATIC

from osgeo import gdal, ogr, osr


def main():

    ds = gdal.Open('test.tif')

    # Get raster band
    band = ds.GetRasterBand(1)

    scanline = band.ReadRaster1(xoff=0, yoff=0,
                        xsize=980, ysize=1012,
                        buf_xsize=956, buf_ysize=988,
                        buf_pixel_space=4, buf_line_space=956*4,
                        buf_type=gdal.GDT_UInt32,
                        resample_alg=gdal.GRIORA_Bilinear)

    print(scanline)

if __name__ == '__main__':
    main()

Versions and provenance

verified in 3.9, will test master soon

Additional context

No response

[elpaso]

incrementing by one does nothing here: