Fix integer overflow in 16-bit resampling

This commit fixes a bug in Resample.c where downsampling 16-bit
images (I;16) using filters with negative lobes (such as
Image.Resampling.LANCZOS) could result in byte corruption.

Because Lanczos weighting can create overshoots (ringing artifacts)
near sharp edges, the accumulated floating-point sum can sometimes
exceed the 16-bit maximum (65535) or fall below zero. Previously,
these out-of-bounds values were not correctly clamped before being
cast or packed into the 16-bit output buffer, leading to integer
overflow/underflow and corrupted pixels.

This update correctly clamps the accumulated float values to the
[0, 65535] range for I;16 images during resampling.

Author: hayatoikoma

Tests/test_image_resample.py

@@ -627,3 +627,37 @@ def test_skip_vertical(self, flt: Image.Resampling) -> None:
                 0.4,
                 f">>> {size} {box} {flt}",
             )
+
+
+class TestCoreResample16bpc:
+    def test_resampling_clamp(self) -> None:
+        # Lanczos weighting during downsampling can push accumulated float sums
+        # above 65535. These must be clamped to 65535, not corrupted byte-by-byte.
+        width, height = 100, 10
+        # I;16 image: left half = 0, right half = 65535
+        im_16 = Image.new("I;16", (width, height))
+        for y in range(height):
+            for x in range(width // 2, width):
+                im_16.putpixel((x, y), 65535)
+        # F image: same values as float reference
+        im_f = Image.new("F", (width, height))
+        for y in range(height):
+            for x in range(width // 2, width):
+                im_f.putpixel((x, y), 65535.0)
+
+        # 5x downsampling with Lanczos creates ~8.7% overshoot at the step edge
+        result_16 = im_16.resize((20, height), Image.Resampling.LANCZOS)
+        result_f = im_f.resize((20, height), Image.Resampling.LANCZOS)
+
+        px_16 = result_16.load()
+        px_f = result_f.load()
+        assert px_16 is not None
+        assert px_f is not None
+        for y in range(height):
+            for x in range(20):
+                v = px_f[x, y]
+                assert isinstance(v, float)
+                expected = max(0, min(65535, round(v)))
+                assert (
+                    px_16[x, y] == expected
+                ), f"Pixel ({x}, {y}): expected {expected}, got {px_16[x, y]}"

src/libImaging/Resample.c

@@ -493,8 +493,13 @@ ImagingResampleHorizontal_16bpc(
                       k[x];
             }
             ss_int = ROUND_UP(ss);
-            imOut->image8[yy][xx * 2 + (bigendian ? 1 : 0)] = CLIP8(ss_int % 256);
-            imOut->image8[yy][xx * 2 + (bigendian ? 0 : 1)] = CLIP8(ss_int >> 8);
+            if (ss_int < 0) {
+                ss_int = 0;
+            } else if (ss_int > 65535) {
+                ss_int = 65535;
+            }
+            imOut->image8[yy][xx * 2 + (bigendian ? 1 : 0)] = ss_int & 0xFF;
+            imOut->image8[yy][xx * 2 + (bigendian ? 0 : 1)] = ss_int >> 8;
         }
     }
     ImagingSectionLeave(&cookie);
@@ -532,8 +537,13 @@ ImagingResampleVertical_16bpc(
                       k[y];
             }
             ss_int = ROUND_UP(ss);
-            imOut->image8[yy][xx * 2 + (bigendian ? 1 : 0)] = CLIP8(ss_int % 256);
-            imOut->image8[yy][xx * 2 + (bigendian ? 0 : 1)] = CLIP8(ss_int >> 8);
+            if (ss_int < 0) {
+                ss_int = 0;
+            } else if (ss_int > 65535) {
+                ss_int = 65535;
+            }
+            imOut->image8[yy][xx * 2 + (bigendian ? 1 : 0)] = ss_int & 0xFF;
+            imOut->image8[yy][xx * 2 + (bigendian ? 0 : 1)] = ss_int >> 8;
         }
     }
     ImagingSectionLeave(&cookie);