diff --git a/Tests/test_image_resample.py b/Tests/test_image_resample.py
index 73b25ed51b2..f188be81ecb 100644
--- a/Tests/test_image_resample.py
+++ b/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:
+    # Lanczos weighting during downsampling can push accumulated float sums
+    @pytest.mark.parametrize(
+        "offset",
+        (
+            # below 0. These must be clamped to 0, not corrupted byte-by-byte.
+            0,  # Left half = 65535, right half = 0
+            # above 65535. These must be clamped to 65535, not corrupted byte-by-byte.
+            50,  # # Left half = 0, right half = 65535
+        ),
+    )
+    def test_resampling_clamp_overflow(self, offset: int) -> None:
+        ims = {}
+        width, height = 100, 10
+        for mode in ("I;16", "F"):
+            im = Image.new(mode, (width, height))
+            im.paste(65535, (offset, 0, offset + width // 2, height))
+
+            # 5x downsampling with Lanczos
+            # creates ~8.7% overshoot or undershoot at the step edge
+            ims[mode] = im.resize((20, height), Image.Resampling.LANCZOS)
+
+        for y in range(height):
+            for x in range(20):
+                v = ims["F"].getpixel((x, y))
+                assert isinstance(v, float)
+                expected = max(0, min(65535, round(v)))
+
+                value = ims["I;16"].getpixel((x, y))
+                assert (
+                    value == expected
+                ), f"Pixel ({x}, {y}): expected {expected}, got {value}"
diff --git a/src/libImaging/ImagingUtils.h b/src/libImaging/ImagingUtils.h
index 714458ad02a..a362780d0ee 100644
--- a/src/libImaging/ImagingUtils.h
+++ b/src/libImaging/ImagingUtils.h
@@ -27,6 +27,8 @@
 
 #define CLIP8(v) ((v) <= 0 ? 0 : (v) < 256 ? (v) : 255)
 
+#define CLIP16(v) ((v) <= 0 ? 0 : (v) < 65536 ? (v) : 65535)
+
 /* This is to work around a bug in GCC prior 4.9 in 64 bit mode.
    GCC generates code with partial dependency which is 3 times slower.
    See: https://stackoverflow.com/a/26588074/253146 */
diff --git a/src/libImaging/Resample.c b/src/libImaging/Resample.c
index cbd18d0c116..a3625701455 100644
--- a/src/libImaging/Resample.c
+++ b/src/libImaging/Resample.c
@@ -492,9 +492,9 @@ ImagingResampleHorizontal_16bpc(
                         << 8)) *
                       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);
+            ss_int = CLIP16(ROUND_UP(ss));
+            imOut->image8[yy][xx * 2 + (bigendian ? 1 : 0)] = ss_int & 0xFF;
+            imOut->image8[yy][xx * 2 + (bigendian ? 0 : 1)] = ss_int >> 8;
         }
     }
     ImagingSectionLeave(&cookie);
@@ -531,9 +531,9 @@ ImagingResampleVertical_16bpc(
                        (imIn->image8[y + ymin][xx * 2 + (bigendian ? 0 : 1)] << 8)) *
                       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);
+            ss_int = CLIP16(ROUND_UP(ss));
+            imOut->image8[yy][xx * 2 + (bigendian ? 1 : 0)] = ss_int & 0xFF;
+            imOut->image8[yy][xx * 2 + (bigendian ? 0 : 1)] = ss_int >> 8;
         }
     }
     ImagingSectionLeave(&cookie);