Fix 16-bit ImageFilter accuracy

Tests/test_image_filter.py

@@ -210,6 +210,17 @@ def test_consistency_5x5(mode: str) -> None:
             assert_image_equal(source.filter(kernel), reference)
 
 
+@pytest.mark.parametrize("mode", ("I;16", "I;16L", "I;16B", "I;16N"))
+def test_consistency_i16_high_byte(mode: str) -> None:
+    # Exercise filters with a 16bpc image that has content in the high byte, too.
+    im = Image.new(mode, (8, 8), 1000)
+    # Ensure repeated smoothing retains the exact same color value,
+    # rather than drifting due to rounding errors.
+    for _ in range(5):
+        im = im.filter(ImageFilter.SMOOTH)
+        assert im.getpixel((4, 4)) == 1000
+
+
 @pytest.mark.parametrize(
     "radius",
     (

src/libImaging/Filter.c

@@ -27,7 +27,6 @@
 #include <stdint.h>
 #include "Imaging.h"
 
-#define ROUND_UP(f) ((int)((f) >= 0.0 ? (f) + 0.5F : (f) - 0.5F))
 #define INT32_MAX_F 2147483647.0F
 
 static inline UINT8
@@ -38,6 +37,14 @@ clip8(float in) {
     return (UINT8)in;
 }
 
+static inline int
+clip16(float in) {
+    // Branchless clamp to [0, 65535].
+    in = in < 0.0f ? 0.0f : in;
+    in = in > 65535.0f ? 65535.0f : in;
+    return (int)in;
+}
+
 static inline INT32
 clip32(float in) {
     // Clamp the low bound branchlessly (maxss).
@@ -114,7 +121,7 @@ kernel_i16(int size, UINT8 *in0, int x, const float *kernel, int bigendian) {
     for (i = 0; i < size; i++) {
         int x1 = x + i - half_size;
         result += (float)(in0[x1 * 2 + (bigendian ? 1 : 0)] +
-                          (in0[x1 * 2 + (bigendian ? 0 : 1)] >> 8)) *
+                          (in0[x1 * 2 + (bigendian ? 0 : 1)] << 8)) *
                   kernel[i];
     }
     return result;
@@ -182,9 +189,10 @@ ImagingFilter3x3(Imaging imOut, Imaging im, const float *kernel, float offset) {
                         ss += kernel_i16(3, in1, x, &kernel[0], bigendian);
                         ss += kernel_i16(3, in0, x, &kernel[3], bigendian);
                         ss += kernel_i16(3, in_1, x, &kernel[6], bigendian);
-                        int ss_int = ROUND_UP(ss);
-                        out[x * 2 + (bigendian ? 1 : 0)] = clip8(ss_int % 256);
-                        out[x * 2 + (bigendian ? 0 : 1)] = clip8(ss_int >> 8);
+                        // NOT rounding here because `offset` already has a +0.5 bias.
+                        int ss_int = clip16(ss);
+                        out[x * 2 + (bigendian ? 1 : 0)] = (UINT8)(ss_int & 0xff);
+                        out[x * 2 + (bigendian ? 0 : 1)] = (UINT8)(ss_int >> 8);
                     } else {
                         ss += KERNEL1x3(in1, x, &kernel[0], 1);
                         ss += KERNEL1x3(in0, x, &kernel[3], 1);
@@ -348,9 +356,10 @@ ImagingFilter5x5(Imaging imOut, Imaging im, const float *kernel, float offset) {
                         ss += kernel_i16(5, in0, x, &kernel[10], bigendian);
                         ss += kernel_i16(5, in_1, x, &kernel[15], bigendian);
                         ss += kernel_i16(5, in_2, x, &kernel[20], bigendian);
-                        int ss_int = ROUND_UP(ss);
-                        out[x * 2 + (bigendian ? 1 : 0)] = clip8(ss_int % 256);
-                        out[x * 2 + (bigendian ? 0 : 1)] = clip8(ss_int >> 8);
+                        // NOT rounding here because `offset` already has a +0.5 bias.
+                        int ss_int = clip16(ss);
+                        out[x * 2 + (bigendian ? 1 : 0)] = (UINT8)(ss_int & 0xff);
+                        out[x * 2 + (bigendian ? 0 : 1)] = (UINT8)(ss_int >> 8);
                     } else {
                         ss += KERNEL1x5(in2, x, &kernel[0], 1);
                         ss += KERNEL1x5(in1, x, &kernel[5], 1);