Speed up Image.filter() for I;16* (#9760)

Co-authored-by: Andrew Murray <radarhere@users.noreply.github.com>

Author: akx

Tests/images/hopper_emboss.bmp Tests/images/hopper_emboss_3x3.bmpTests/images/hopper_emboss.bmp renamed to Tests/images/hopper_emboss_3x3.bmp

@@ -8,6 +8,22 @@
 
 from .helper import assert_image_equal, hopper
 
+MODES = (
+    "L",
+    "LA",
+    "La",
+    "I",
+    "I;16",
+    "I;16B",
+    "I;16L",
+    "I;16N",
+    "RGB",
+    "RGBA",
+    "RGBa",
+    "RGBX",
+    "CMYK",
+)
+
 
 @pytest.mark.parametrize(
     "filter_to_apply",
@@ -37,9 +53,7 @@
         ImageFilter.UnsharpMask(10),
     ),
 )
-@pytest.mark.parametrize(
-    "mode", ("L", "I", "I;16", "I;16L", "I;16B", "I;16N", "RGB", "CMYK")
-)
+@pytest.mark.parametrize("mode", MODES)
 def test_sanity(
     filter_to_apply: ImageFilter.Filter | type[ImageFilter.Filter], mode: str
 ) -> None:
@@ -53,20 +67,24 @@ def test_sanity(
         assert out.size == im.size
 
 
-@pytest.mark.parametrize(
-    "mode", ("L", "I", "I;16", "I;16L", "I;16B", "I;16N", "RGB", "CMYK")
-)
+@pytest.mark.parametrize("mode", MODES)
 def test_sanity_error(mode: str) -> None:
     im = hopper(mode)
     with pytest.raises(TypeError):
         im.filter("hello")  # type: ignore[arg-type]
 
 
-# crashes on small images
-@pytest.mark.parametrize("size", ((1, 1), (2, 2), (3, 3)))
-def test_crash(size: tuple[int, int]) -> None:
-    im = Image.new("RGB", size)
-    im.filter(ImageFilter.SMOOTH)
+def test_noop_on_small_images() -> None:
+    # If image is smaller than the kernel size, return it as-is.
+    kernel_size: tuple[int, int] = ImageFilter.SMOOTH_MORE.filterargs[0]
+    kernel_w, kernel_h = kernel_size
+    for w in range(1, kernel_w):
+        for h in range(1, kernel_h):
+            im = hopper("RGB").resize((w, h))
+            # Precondition for the below equality test:
+            # filter is larger or equal to image.
+            assert im.size < kernel_size
+            assert_image_equal(im.filter(ImageFilter.SMOOTH_MORE), im)
 
 
 @pytest.mark.parametrize(
@@ -174,38 +192,39 @@ def test_kernel_not_enough_coefficients() -> None:
         ImageFilter.Kernel((3, 3), (0, 0))
 
 
-@pytest.mark.parametrize(
-    "mode", ("L", "LA", "I", "I;16", "I;16L", "I;16B", "I;16N", "RGB", "CMYK")
-)
-def test_consistency_3x3(mode: str) -> None:
-    matrix = (
+EMBOSS_MATRIX = {
+    3: (
         -1, -1,  0,
         -1,  0,  1,
          0,  1,  1,
-    )  # fmt: skip
-    with Image.open("Tests/images/hopper.bmp") as source:
-        with Image.open("Tests/images/hopper_emboss.bmp") as reference:
-            kernel = ImageFilter.Kernel((3, 3), matrix, 0.3)
-            assert_image_equal(source.filter(kernel), reference)
-
-
-@pytest.mark.parametrize(
-    "mode", ("L", "LA", "I", "I;16", "I;16L", "I;16B", "I;16N", "RGB", "CMYK")
-)
-def test_consistency_5x5(mode: str) -> None:
-    matrix = (
+    ),
+    5: (
         -1, -1, -1, -1,  0,
         -1, -1, -1,  0,  1,
         -1, -1,  0,  1,  1,
         -1,  0,  1,  1,  1,
          0,  1,  1,  1,  1,
-    )  # fmt: skip
+    ),
+}  # fmt: skip
+
+
+@pytest.mark.parametrize("size", (3, 5))
+def test_consistency(size: int) -> None:
+    kernel = ImageFilter.Kernel((size, size), EMBOSS_MATRIX[size], 0.3)
     with Image.open("Tests/images/hopper.bmp") as source:
-        with Image.open("Tests/images/hopper_emboss_more.bmp") as reference:
-            kernel = ImageFilter.Kernel((5, 5), matrix, 0.3)
+        with Image.open(f"Tests/images/hopper_emboss_{size}x{size}.bmp") as reference:
             assert_image_equal(source.filter(kernel), reference)
 
 
+@pytest.mark.parametrize("size", (3, 5))
+@pytest.mark.parametrize("mode", ("I;16", "I;16L", "I;16B", "I;16N"))
+def test_consistency_i16(size: int, mode: str) -> None:
+    kernel = ImageFilter.Kernel((size, size), EMBOSS_MATRIX[size], 0.3)
+    reference = hopper("I").filter(kernel)
+    result = hopper(mode).filter(kernel)
+    assert_image_equal(result.convert("I"), 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.

Tests/images/hopper_emboss_more.bmp Tests/images/hopper_emboss_5x5.bmpTests/images/hopper_emboss_more.bmp renamed to Tests/images/hopper_emboss_5x5.bmp

@@ -53,10 +53,8 @@ def filter(self, image: _imaging.ImagingCore) -> _imaging.ImagingCore:
 
 class Kernel(BuiltinFilter):
     """
-    Create a convolution kernel. This only supports 3x3 and 5x5 integer and floating
-    point kernels.
-
-    Kernels can only be applied to "L" and "RGB" images.
+    Create a convolution kernel.
+    This only supports 3x3 and 5x5 integer and floating point kernels.
 
     :param size: Kernel size, given as (width, height). This must be (3,3) or (5,5).
     :param kernel: A sequence containing kernel weights. The kernel will be flipped

Tests/test_image_filter.py

@@ -161,17 +161,15 @@ ImagingFilter3x3(Imaging imOut, Imaging im, const float *kernel, float offset) {
                 }
                 out[x] = in0[x];
             }
-        } else {
+        } else if (im->type == IMAGING_TYPE_I16) {
             int bigendian = 0;
-            if (im->type == IMAGING_TYPE_I16) {
-                if (
-                    im->mode == IMAGING_MODE_I_16B
+            if (
+                im->mode == IMAGING_MODE_I_16B
 #ifdef WORDS_BIGENDIAN
-                    || im->mode == IMAGING_MODE_I_16N
+                || im->mode == IMAGING_MODE_I_16N
 #endif
-                ) {
-                    bigendian = 1;
-                }
+            ) {
+                bigendian = 1;
             }
             for (y = 1; y < ysize - 1; y++) {
                 UINT8 *restrict in_1 = (UINT8 *)im->image[y - 1];
@@ -180,32 +178,36 @@ ImagingFilter3x3(Imaging imOut, Imaging im, const float *kernel, float offset) {
                 UINT8 *restrict out = (UINT8 *)imOut->image[y];
 
                 out[0] = in0[0];
-                if (im->type == IMAGING_TYPE_I16) {
-                    out[1] = in0[1];
-                }
+                out[1] = in0[1];
                 for (x = 1; x < xsize - 1; x++) {
                     float ss = offset;
-                    if (im->type == IMAGING_TYPE_I16) {
-                        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);
-                        // 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);
-                        ss += KERNEL1x3(in_1, x, &kernel[6], 1);
-                        out[x] = clip8(ss);
-                    }
+                    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);
+                    // 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);
                 }
-                if (im->type == IMAGING_TYPE_I16) {
-                    out[x * 2] = in0[x * 2];
-                    out[x * 2 + 1] = in0[x * 2 + 1];
-                } else {
-                    out[x] = in0[x];
+                out[x * 2] = in0[x * 2];
+                out[x * 2 + 1] = in0[x * 2 + 1];
+            }
+        } else {
+            for (y = 1; y < ysize - 1; y++) {
+                UINT8 *restrict in_1 = (UINT8 *)im->image[y - 1];
+                UINT8 *restrict in0 = (UINT8 *)im->image[y];
+                UINT8 *restrict in1 = (UINT8 *)im->image[y + 1];
+                UINT8 *restrict out = (UINT8 *)imOut->image[y];
+
+                out[0] = in0[0];
+                for (x = 1; x < xsize - 1; x++) {
+                    float ss = offset;
+                    ss += KERNEL1x3(in1, x, &kernel[0], 1);
+                    ss += KERNEL1x3(in0, x, &kernel[3], 1);
+                    ss += KERNEL1x3(in_1, x, &kernel[6], 1);
+                    out[x] = clip8(ss);
                 }
+                out[x] = in0[x];
             }
         }
     } else {
@@ -322,17 +324,15 @@ ImagingFilter5x5(Imaging imOut, Imaging im, const float *kernel, float offset) {
                 out[x + 0] = in0[x + 0];
                 out[x + 1] = in0[x + 1];
             }
-        } else {
+        } else if (im->type == IMAGING_TYPE_I16) {
             int bigendian = 0;
-            if (im->type == IMAGING_TYPE_I16) {
-                if (
-                    im->mode == IMAGING_MODE_I_16B
+            if (
+                im->mode == IMAGING_MODE_I_16B
 #ifdef WORDS_BIGENDIAN
-                    || im->mode == IMAGING_MODE_I_16N
+                || im->mode == IMAGING_MODE_I_16N
 #endif
-                ) {
-                    bigendian = 1;
-                }
+            ) {
+                bigendian = 1;
             }
             for (y = 2; y < ysize - 2; y++) {
                 UINT8 *restrict in_2 = (UINT8 *)im->image[y - 2];
@@ -344,40 +344,47 @@ ImagingFilter5x5(Imaging imOut, Imaging im, const float *kernel, float offset) {
 
                 out[0] = in0[0];
                 out[1] = in0[1];
-                if (im->type == IMAGING_TYPE_I16) {
-                    out[2] = in0[2];
-                    out[3] = in0[3];
-                }
+                out[2] = in0[2];
+                out[3] = in0[3];
                 for (x = 2; x < xsize - 2; x++) {
                     float ss = offset;
-                    if (im->type == IMAGING_TYPE_I16) {
-                        ss += kernel_i16(5, in2, x, &kernel[0], bigendian);
-                        ss += kernel_i16(5, in1, x, &kernel[5], bigendian);
-                        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);
-                        // 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);
-                        ss += KERNEL1x5(in0, x, &kernel[10], 1);
-                        ss += KERNEL1x5(in_1, x, &kernel[15], 1);
-                        ss += KERNEL1x5(in_2, x, &kernel[20], 1);
-                        out[x] = clip8(ss);
-                    }
+                    ss += kernel_i16(5, in2, x, &kernel[0], bigendian);
+                    ss += kernel_i16(5, in1, x, &kernel[5], bigendian);
+                    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);
+                    // 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);
                 }
-                if (im->type == IMAGING_TYPE_I16) {
-                    out[x * 2 + 0] = in0[x * 2 + 0];
-                    out[x * 2 + 1] = in0[x * 2 + 1];
-                    out[x * 2 + 2] = in0[x * 2 + 2];
-                    out[x * 2 + 3] = in0[x * 2 + 3];
-                } else {
-                    out[x + 0] = in0[x + 0];
-                    out[x + 1] = in0[x + 1];
+                out[x * 2 + 0] = in0[x * 2 + 0];
+                out[x * 2 + 1] = in0[x * 2 + 1];
+                out[x * 2 + 2] = in0[x * 2 + 2];
+                out[x * 2 + 3] = in0[x * 2 + 3];
+            }
+        } else {
+            for (y = 2; y < ysize - 2; y++) {
+                UINT8 *restrict in_2 = (UINT8 *)im->image[y - 2];
+                UINT8 *restrict in_1 = (UINT8 *)im->image[y - 1];
+                UINT8 *restrict in0 = (UINT8 *)im->image[y];
+                UINT8 *restrict in1 = (UINT8 *)im->image[y + 1];
+                UINT8 *restrict in2 = (UINT8 *)im->image[y + 2];
+                UINT8 *restrict out = (UINT8 *)imOut->image[y];
+
+                out[0] = in0[0];
+                out[1] = in0[1];
+                for (x = 2; x < xsize - 2; x++) {
+                    float ss = offset;
+                    ss += KERNEL1x5(in2, x, &kernel[0], 1);
+                    ss += KERNEL1x5(in1, x, &kernel[5], 1);
+                    ss += KERNEL1x5(in0, x, &kernel[10], 1);
+                    ss += KERNEL1x5(in_1, x, &kernel[15], 1);
+                    ss += KERNEL1x5(in_2, x, &kernel[20], 1);
+                    out[x] = clip8(ss);
                 }
+                out[x + 0] = in0[x + 0];
+                out[x + 1] = in0[x + 1];
             }
         }
     } else {