Fix round-up error in 16-bit filtering

Author: akx

Tests/test_image_filter.py

@@ -214,10 +214,11 @@ def test_consistency_5x5(mode: str) -> None:
 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)
-    pixel = im.filter(ImageFilter.SMOOTH).getpixel((4, 4))
-    assert isinstance(pixel, (int, float))
-    # Allow for kernel rounding errors.
-    assert abs(pixel - 1000) <= 1
+    # 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(

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).
@@ -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);