Merge pull request #961 from homm/fast-box-blur

Merge Fast Gaussian Blur

PIL/ImageOps.py

@@ -20,6 +20,7 @@
 from PIL import Image
 from PIL._util import isStringType
 import operator
+import math
 from functools import reduce
 
 
@@ -441,3 +442,22 @@ def unsharp_mask(im, radius=None, percent=None, threshold=None):
     return im.im.unsharp_mask(radius, percent, threshold)
 
 usm = unsharp_mask
+
+
+def box_blur(image, radius):
+    """
+    Blur the image by setting each pixel to the average value of the pixels
+    in a square box extending radius pixels in each direction.
+    Supports float radius of arbitrary size. Uses an optimized implementation
+    which runs in linear time relative to the size of the image
+    for any radius value.
+
+    :param image: The image to blur.
+    :param radius: Size of the box in one direction. Radius 0 does not blur,
+                   returns an identical image. Radius 1 takes 1 pixel
+                   in each direction, i.e. 9 pixels in total.
+    :return: An image.
+    """
+    image.load()
+
+    return image._new(image.im.box_blur(radius))

Tests/test_box_blur.py

@@ -0,0 +1,229 @@
+from helper import unittest, PillowTestCase
+
+from PIL import Image, ImageOps
+
+
+sample = Image.new("L", (7, 5))
+sample.putdata(sum([
+    [210, 50,  20,  10,  220, 230, 80],
+    [190, 210, 20,  180, 170, 40,  110],
+    [120, 210, 250, 60,  220, 0,   220],
+    [220, 40,  230, 80,  130, 250, 40],
+    [250, 0,   80,  30,  60,  20,  110],
+], []))
+
+
+class ImageMock(object):
+    def __init__(self):
+        self.im = self
+
+    def load(self):
+        pass
+
+    def _new(self, im):
+        return im
+
+    def box_blur(self, radius, n):
+        return radius, n
+
+
+class TestBoxBlurApi(PillowTestCase):
+
+    def test_imageops_box_blur(self):
+        i = ImageOps.box_blur(sample, 1)
+        self.assertEqual(i.mode, sample.mode)
+        self.assertEqual(i.size, sample.size)
+        self.assertIsInstance(i, Image.Image)
+
+
+class TestBoxBlur(PillowTestCase):
+
+    def box_blur(self, image, radius=1, n=1):
+        return image._new(image.im.box_blur(radius, n))
+
+    def assertImage(self, im, data, delta=0):
+        it = iter(im.getdata())
+        for data_row in data:
+            im_row = [next(it) for _ in range(im.size[0])]
+            if any(
+                abs(data_v - im_v) > delta
+                for data_v, im_v in zip(data_row, im_row)
+            ):
+                self.assertEqual(im_row, data_row)
+        self.assertRaises(StopIteration, next, it)
+
+    def assertBlur(self, im, radius, data, passes=1, delta=0):
+        # check grayscale image
+        self.assertImage(self.box_blur(im, radius, passes), data, delta)
+        rgba = Image.merge('RGBA', (im, im, im, im))
+        for band in self.box_blur(rgba, radius, passes).split():
+            self.assertImage(band, data, delta)
+
+    def test_color_modes(self):
+        self.assertRaises(ValueError, self.box_blur, sample.convert("1"))
+        self.assertRaises(ValueError, self.box_blur, sample.convert("P"))
+        self.box_blur(sample.convert("L"))
+        self.box_blur(sample.convert("LA"))
+        self.assertRaises(ValueError, self.box_blur, sample.convert("I"))
+        self.assertRaises(ValueError, self.box_blur, sample.convert("F"))
+        self.box_blur(sample.convert("RGB"))
+        self.box_blur(sample.convert("RGBA"))
+        self.box_blur(sample.convert("CMYK"))
+        self.assertRaises(ValueError, self.box_blur, sample.convert("YCbCr"))
+
+    def test_radius_0(self):
+        self.assertBlur(
+            sample, 0,
+            [
+                [210, 50,  20,  10,  220, 230, 80],
+                [190, 210, 20,  180, 170, 40,  110],
+                [120, 210, 250, 60,  220, 0,   220],
+                [220, 40,  230, 80,  130, 250, 40],
+                [250, 0,   80,  30,  60,  20,  110],
+            ]
+        )
+
+    def test_radius_0_02(self):
+        self.assertBlur(
+            sample, 0.02,
+            [
+                [206, 55,  20,  17,  215, 223, 83],
+                [189, 203, 31,  171, 169, 46,  110],
+                [125, 206, 241, 69,  210, 13,  210],
+                [215, 49,  221, 82,  131, 235, 48],
+                [244, 7,   80,  32,  60,  27,  107],
+            ],
+            delta=2,
+        )
+
+    def test_radius_0_05(self):
+        self.assertBlur(
+            sample, 0.05,
+            [
+                [202, 62,  22,  27,  209, 215, 88],
+                [188, 194, 44,  161, 168, 56,  111],
+                [131, 201, 229, 81,  198, 31,  198],
+                [209, 62,  209, 86,  133, 216, 59],
+                [237, 17,  80,  36,  60,  35,  103],
+            ],
+            delta=2,
+        )
+
+    def test_radius_0_1(self):
+        self.assertBlur(
+            sample, 0.1,
+            [
+                [196, 72,  24,  40,  200, 203, 93],
+                [187, 183, 62,  148, 166, 68,  111],
+                [139, 193, 213, 96,  182, 54,  182],
+                [201, 78,  193, 91,  133, 191, 73],
+                [227, 31,  80,  42,  61,  47,  99],
+            ],
+            delta=1,
+        )
+
+    def test_radius_0_5(self):
+        self.assertBlur(
+            sample, 0.5,
+            [
+                [176, 101, 46,  83,  163, 165, 111],
+                [176, 149, 108, 122, 144, 120, 117],
+                [164, 171, 159, 141, 134, 119, 129],
+                [170, 136, 133, 114, 116, 124, 109],
+                [184, 95,  72,  70,  69,  81,  89],
+            ],
+            delta=1,
+        )
+
+    def test_radius_1(self):
+        self.assertBlur(
+            sample, 1,
+            [
+                [170, 109, 63,  97,  146, 153, 116],
+                [168, 142, 112, 128, 126, 143, 121],
+                [169, 166, 142, 149, 126, 131, 114],
+                [159, 156, 109, 127, 94,  117, 112],
+                [164, 128, 63,  87,  76,  89,  90],
+            ],
+            delta=1,
+        )
+
+    def test_radius_1_5(self):
+        self.assertBlur(
+            sample, 1.5,
+            [
+                [155, 120, 105, 112, 124, 137, 130],
+                [160, 136, 124, 125, 127, 134, 130],
+                [166, 147, 130, 125, 120, 121, 119],
+                [168, 145, 119, 109, 103, 105, 110],
+                [168, 134, 96,  85,  85,  89,  97],
+            ],
+            delta=1,
+        )
+
+    def test_radius_bigger_then_half(self):
+        self.assertBlur(
+            sample, 3,
+            [
+                [144, 145, 142, 128, 114, 115, 117],
+                [148, 145, 137, 122, 109, 111, 112],
+                [152, 145, 131, 117, 103, 107, 108],
+                [156, 144, 126, 111, 97,  102, 103],
+                [160, 144, 121, 106, 92,  98,  99],
+            ],
+            delta=1,
+        )
+
+    def test_radius_bigger_then_width(self):
+        self.assertBlur(
+            sample, 10,
+            [
+                [158, 153, 147, 141, 135, 129, 123],
+                [159, 153, 147, 141, 136, 130, 124],
+                [159, 154, 148, 142, 136, 130, 124],
+                [160, 154, 148, 142, 137, 131, 125],
+                [160, 155, 149, 143, 137, 131, 125],
+            ],
+            delta=0,
+        )
+
+    def test_exteme_large_radius(self):
+        self.assertBlur(
+            sample, 600,
+            [
+                [162, 162, 162, 162, 162, 162, 162],
+                [162, 162, 162, 162, 162, 162, 162],
+                [162, 162, 162, 162, 162, 162, 162],
+                [162, 162, 162, 162, 162, 162, 162],
+                [162, 162, 162, 162, 162, 162, 162],
+            ],
+            delta=1,
+        )
+
+    def test_two_passes(self):
+        self.assertBlur(
+            sample, 1,
+            [
+                [153, 123, 102, 109, 132, 135, 129],
+                [159, 138, 123, 121, 133, 131, 126],
+                [162, 147, 136, 124, 127, 121, 121],
+                [159, 140, 125, 108, 111, 106, 108],
+                [154, 126, 105, 87,  94,  93,  97],
+            ],
+            passes=2,
+            delta=1,
+        )
+
+    def test_three_passes(self):
+        self.assertBlur(
+            sample, 1,
+            [
+                [146, 131, 116, 118, 126, 131, 130],
+                [151, 138, 125, 123, 126, 128, 127],
+                [154, 143, 129, 123, 120, 120, 119],
+                [152, 139, 122, 113, 108, 108, 108],
+                [148, 132, 112, 102, 97,  99,  100],
+            ],
+            passes=3,
+            delta=1,
+        )

Tests/test_imageops_usm.py

@@ -5,6 +5,7 @@ from PIL import ImageOps
 from PIL import ImageFilter
 
 im = Image.open("Tests/images/hopper.ppm")
+snakes = Image.open("Tests/images/color_snakes.png")
 
 
 class TestImageOpsUsm(PillowTestCase):
@@ -16,7 +17,7 @@ class TestImageOpsUsm(PillowTestCase):
         self.assertEqual(i.size, (128, 128))
         # i.save("blur.bmp")
 
-        i = ImageOps.usm(im, 2.0, 125, 8)
+        i = ImageOps.unsharp_mask(im, 2.0, 125, 8)
         self.assertEqual(i.mode, "RGB")
         self.assertEqual(i.size, (128, 128))
         # i.save("usm.bmp")
@@ -33,7 +34,7 @@ class TestImageOpsUsm(PillowTestCase):
         self.assertEqual(i.mode, "RGB")
         self.assertEqual(i.size, (128, 128))
 
-    def test_usm(self):
+    def test_usm_formats(self):
 
         usm = ImageOps.unsharp_mask
         self.assertRaises(ValueError, lambda: usm(im.convert("1")))
@@ -45,7 +46,7 @@ class TestImageOpsUsm(PillowTestCase):
         usm(im.convert("CMYK"))
         self.assertRaises(ValueError, lambda: usm(im.convert("YCbCr")))
 
-    def test_blur(self):
+    def test_blur_formats(self):
 
         blur = ImageOps.gaussian_blur
         self.assertRaises(ValueError, lambda: blur(im.convert("1")))
@@ -57,6 +58,32 @@ class TestImageOpsUsm(PillowTestCase):
         blur(im.convert("CMYK"))
         self.assertRaises(ValueError, lambda: blur(im.convert("YCbCr")))
 
+    def test_usm_accuracy(self):
+
+        src = snakes.convert('RGB')
+        i = src._new(ImageOps.unsharp_mask(src, 5, 1024, 0))
+        # Image should not be changed because it have only 0 and 255 levels.
+        self.assertEqual(i.tobytes(), src.tobytes())
+
+    def test_blur_accuracy(self):
+
+        i = snakes._new(ImageOps.gaussian_blur(snakes, .4))
+        # These pixels surrounded with pixels with 255 intensity.
+        # They must be very close to 255.
+        for x, y, c in [(1, 0, 1), (2, 0, 1), (7, 8, 1), (8, 8, 1), (2, 9, 1),
+                        (7, 3, 0), (8, 3, 0), (5, 8, 0), (5, 9, 0), (1, 3, 0),
+                        (4, 3, 2), (4, 2, 2)]:
+            self.assertGreaterEqual(i.im.getpixel((x, y))[c], 250)
+        # Fuzzy match.
+        gp = lambda x, y: i.im.getpixel((x, y))
+        self.assertTrue(236 <= gp(7, 4)[0] <= 239)
+        self.assertTrue(236 <= gp(7, 5)[2] <= 239)
+        self.assertTrue(236 <= gp(7, 6)[2] <= 239)
+        self.assertTrue(236 <= gp(7, 7)[1] <= 239)
+        self.assertTrue(236 <= gp(8, 4)[0] <= 239)
+        self.assertTrue(236 <= gp(8, 5)[2] <= 239)
+        self.assertTrue(236 <= gp(8, 6)[2] <= 239)
+        self.assertTrue(236 <= gp(8, 7)[1] <= 239)
 
 if __name__ == '__main__':
     unittest.main()

_imaging.c

@@ -863,7 +863,8 @@ _gaussian_blur(ImagingObject* self, PyObject* args)
     Imaging imOut;
 
     float radius = 0;
-    if (!PyArg_ParseTuple(args, "f", &radius))
+    int passes = 3;
+    if (!PyArg_ParseTuple(args, "f|i", &radius, &passes))
         return NULL;
 
     imIn = self->image;
@@ -871,7 +872,7 @@ _gaussian_blur(ImagingObject* self, PyObject* args)
     if (!imOut)
         return NULL;
 
-    if (!ImagingGaussianBlur(imIn, imOut, radius))
+    if (!ImagingGaussianBlur(imOut, imIn, radius, passes))
         return NULL;
 
     return PyImagingNew(imOut);
@@ -1769,18 +1770,39 @@ _unsharp_mask(ImagingObject* self, PyObject* args)
     if (!PyArg_ParseTuple(args, "fii", &radius, &percent, &threshold))
         return NULL;
 
+    imIn = self->image;
+    imOut = ImagingNew(imIn->mode, imIn->xsize, imIn->ysize);
+    if (!imOut)
+        return NULL;
+
+    if (!ImagingUnsharpMask(imOut, imIn, radius, percent, threshold))
+        return NULL;
+
+    return PyImagingNew(imOut);
+}
+#endif
+
+static PyObject*
+_box_blur(ImagingObject* self, PyObject* args)
+{
+    Imaging imIn;
+    Imaging imOut;
+
+    float radius;
+    int n = 1;
+    if (!PyArg_ParseTuple(args, "f|i", &radius, &n))
+        return NULL;
 
     imIn = self->image;
     imOut = ImagingNew(imIn->mode, imIn->xsize, imIn->ysize);
     if (!imOut)
         return NULL;
 
-    if (!ImagingUnsharpMask(imIn, imOut, radius, percent, threshold))
+    if (!ImagingBoxBlur(imOut, imIn, radius, n))
         return NULL;
 
     return PyImagingNew(imOut);
 }
-#endif
 
 /* -------------------------------------------------------------------- */
 
@@ -3056,6 +3078,8 @@ static struct PyMethodDef methods[] = {
     {"unsharp_mask", (PyCFunction)_unsharp_mask, 1},
 #endif
 
+    {"box_blur", (PyCFunction)_box_blur, 1},    
+
 #ifdef WITH_EFFECTS
     /* Special effects */
     {"effect_spread", (PyCFunction)_effect_spread, 1},

libImaging/BoxBlur.c

@@ -0,0 +1,308 @@
+#include "Python.h"
+#include "Imaging.h"
+
+
+#define MAX(x, y) (((x) > (y)) ? (x) : (y))
+#define MIN(x, y) (((x) < (y)) ? (x) : (y))
+
+
+typedef UINT8 pixel[4];
+
+void inline
+ImagingLineBoxBlur32(pixel *lineOut, pixel *lineIn, int lastx, int radius, int edgeA,
+    int edgeB, UINT32 ww, UINT32 fw)
+{
+    int x;
+    UINT32 acc[4];
+    UINT32 bulk[4];
+
+    #define MOVE_ACC(acc, subtract, add) \
+        acc[0] += lineIn[add][0] - lineIn[subtract][0]; \
+        acc[1] += lineIn[add][1] - lineIn[subtract][1]; \
+        acc[2] += lineIn[add][2] - lineIn[subtract][2]; \
+        acc[3] += lineIn[add][3] - lineIn[subtract][3];
+
+    #define ADD_FAR(bulk, acc, left, right) \
+        bulk[0] = (acc[0] * ww) + (lineIn[left][0] + lineIn[right][0]) * fw; \
+        bulk[1] = (acc[1] * ww) + (lineIn[left][1] + lineIn[right][1]) * fw; \
+        bulk[2] = (acc[2] * ww) + (lineIn[left][2] + lineIn[right][2]) * fw; \
+        bulk[3] = (acc[3] * ww) + (lineIn[left][3] + lineIn[right][3]) * fw;
+
+    #define SAVE(x, bulk) \
+        lineOut[x][0] = (UINT8)((bulk[0] + (1 << 23)) >> 24); \
+        lineOut[x][1] = (UINT8)((bulk[1] + (1 << 23)) >> 24); \
+        lineOut[x][2] = (UINT8)((bulk[2] + (1 << 23)) >> 24); \
+        lineOut[x][3] = (UINT8)((bulk[3] + (1 << 23)) >> 24);
+
+    /* Compute acc for -1 pixel (outside of image):
+       From "-radius-1" to "-1" get first pixel,
+       then from "0" to "radius-1". */
+    acc[0] = lineIn[0][0] * (radius + 1);
+    acc[1] = lineIn[0][1] * (radius + 1);
+    acc[2] = lineIn[0][2] * (radius + 1);
+    acc[3] = lineIn[0][3] * (radius + 1);
+    /* As radius can be bigger than xsize, iterate to edgeA -1. */
+    for (x = 0; x < edgeA - 1; x++) {
+        acc[0] += lineIn[x][0];
+        acc[1] += lineIn[x][1];
+        acc[2] += lineIn[x][2];
+        acc[3] += lineIn[x][3];
+    }
+    /* Then multiply remainder to last x. */
+    acc[0] += lineIn[lastx][0] * (radius - edgeA + 1);
+    acc[1] += lineIn[lastx][1] * (radius - edgeA + 1);
+    acc[2] += lineIn[lastx][2] * (radius - edgeA + 1);
+    acc[3] += lineIn[lastx][3] * (radius - edgeA + 1);
+
+    if (edgeA <= edgeB)
+    {
+        /* Subtract pixel from left ("0").
+           Add pixels from radius. */
+        for (x = 0; x < edgeA; x++) {
+            MOVE_ACC(acc, 0, x + radius);
+            ADD_FAR(bulk, acc, 0, x + radius + 1);
+            SAVE(x, bulk);
+        }
+        /* Subtract previous pixel from "-radius".
+           Add pixels from radius. */
+        for (x = edgeA; x < edgeB; x++) {
+            MOVE_ACC(acc, x - radius - 1, x + radius);
+            ADD_FAR(bulk, acc, x - radius - 1, x + radius + 1);
+            SAVE(x, bulk);
+        }
+        /* Subtract previous pixel from "-radius".
+           Add last pixel. */
+        for (x = edgeB; x <= lastx; x++) {
+            MOVE_ACC(acc, x - radius - 1, lastx);
+            ADD_FAR(bulk, acc, x - radius - 1, lastx);
+            SAVE(x, bulk);
+        }
+    }
+    else
+    {
+        for (x = 0; x < edgeB; x++) {
+            MOVE_ACC(acc, 0, x + radius);
+            ADD_FAR(bulk, acc, 0, x + radius + 1);
+            SAVE(x, bulk);
+        }
+        for (x = edgeB; x < edgeA; x++) {
+            MOVE_ACC(acc, 0, lastx);
+            ADD_FAR(bulk, acc, 0, lastx);
+            SAVE(x, bulk);
+        }
+        for (x = edgeA; x <= lastx; x++) {
+            MOVE_ACC(acc, x - radius - 1, lastx);
+            ADD_FAR(bulk, acc, x - radius - 1, lastx);
+            SAVE(x, bulk);
+        }
+    }
+
+    #undef MOVE_ACC
+    #undef ADD_FAR
+    #undef SAVE
+}
+
+
+void inline
+ImagingLineBoxBlur8(UINT8 *lineOut, UINT8 *lineIn, int lastx, int radius, int edgeA,
+    int edgeB, UINT32 ww, UINT32 fw)
+{
+    int x;
+    UINT32 acc;
+    UINT32 bulk;
+
+    #define MOVE_ACC(acc, subtract, add) \
+        acc += lineIn[add] - lineIn[subtract];
+
+    #define ADD_FAR(bulk, acc, left, right) \
+        bulk = (acc * ww) + (lineIn[left] + lineIn[right]) * fw;
+
+    #define SAVE(x, bulk) \
+        lineOut[x] = (UINT8)((bulk + (1 << 23)) >> 24)
+
+    acc = lineIn[0] * (radius + 1);
+    for (x = 0; x < edgeA - 1; x++) {
+        acc += lineIn[x];
+    }
+    acc += lineIn[lastx] * (radius - edgeA + 1);
+
+    if (edgeA <= edgeB)
+    {
+        for (x = 0; x < edgeA; x++) {
+            MOVE_ACC(acc, 0, x + radius);
+            ADD_FAR(bulk, acc, 0, x + radius + 1);
+            SAVE(x, bulk);
+        }
+        for (x = edgeA; x < edgeB; x++) {
+            MOVE_ACC(acc, x - radius - 1, x + radius);
+            ADD_FAR(bulk, acc, x - radius - 1, x + radius + 1);
+            SAVE(x, bulk);
+        }
+        for (x = edgeB; x <= lastx; x++) {
+            MOVE_ACC(acc, x - radius - 1, lastx);
+            ADD_FAR(bulk, acc, x - radius - 1, lastx);
+            SAVE(x, bulk);
+        }
+    }
+    else
+    {
+        for (x = 0; x < edgeB; x++) {
+            MOVE_ACC(acc, 0, x + radius);
+            ADD_FAR(bulk, acc, 0, x + radius + 1);
+            SAVE(x, bulk);
+        }
+        for (x = edgeB; x < edgeA; x++) {
+            MOVE_ACC(acc, 0, lastx);
+            ADD_FAR(bulk, acc, 0, lastx);
+            SAVE(x, bulk);
+        }
+        for (x = edgeA; x <= lastx; x++) {
+            MOVE_ACC(acc, x - radius - 1, lastx);
+            ADD_FAR(bulk, acc, x - radius - 1, lastx);
+            SAVE(x, bulk);
+        }
+    }
+
+    #undef MOVE_ACC
+    #undef ADD_FAR
+    #undef SAVE
+}
+
+
+
+Imaging
+ImagingHorizontalBoxBlur(Imaging imOut, Imaging imIn, float floatRadius)
+{
+    ImagingSectionCookie cookie;
+
+    int y;
+
+    int radius = (int) floatRadius;
+    UINT32 ww = (UINT32) (1 << 24) / (floatRadius * 2 + 1);
+    UINT32 fw = ((1 << 24) - (radius * 2 + 1) * ww) / 2;
+
+    int edgeA = MIN(radius + 1, imIn->xsize);
+    int edgeB = MAX(imIn->xsize - radius - 1, 0);
+
+    UINT32 *lineOut = calloc(imIn->xsize, sizeof(UINT32));
+    if (lineOut == NULL)
+        return ImagingError_MemoryError();
+
+    // printf(">>> %d %d %d\n", radius, ww, fw);
+
+    ImagingSectionEnter(&cookie);
+
+    if (imIn->image8)
+    {
+        for (y = 0; y < imIn->ysize; y++) {
+            ImagingLineBoxBlur8(
+                (imIn == imOut ? (UINT8 *) lineOut : imOut->image8[y]),
+                imIn->image8[y],
+                imIn->xsize - 1,
+                radius, edgeA, edgeB,
+                ww, fw
+            );
+            if (imIn == imOut) {
+                // Commit.
+                memcpy(imOut->image8[y], lineOut, imIn->xsize);
+            }
+        }
+    }
+    else
+    {
+        for (y = 0; y < imIn->ysize; y++) {
+            ImagingLineBoxBlur32(
+                imIn == imOut ? (pixel *) lineOut : (pixel *) imOut->image32[y],
+                (pixel *) imIn->image32[y],
+                imIn->xsize - 1,
+                radius, edgeA, edgeB,
+                ww, fw
+            );
+            if (imIn == imOut) {
+                // Commit.
+                memcpy(imOut->image32[y], lineOut, imIn->xsize * 4);
+            }
+        }
+    }
+
+    ImagingSectionLeave(&cookie);
+
+    free(lineOut);
+
+    return imOut;
+}
+
+
+Imaging
+ImagingBoxBlur(Imaging imOut, Imaging imIn, float radius, int n)
+{
+    int i;
+
+    if (n < 1) {
+        return ImagingError_ValueError(
+            "number of passes must be greater than zero"
+        );
+    }
+
+    if (strcmp(imIn->mode, imOut->mode) ||
+        imIn->type  != imOut->type  ||
+        imIn->bands != imOut->bands ||
+        imIn->xsize != imOut->xsize ||
+        imIn->ysize != imOut->ysize)
+        return ImagingError_Mismatch();
+
+    if (imIn->type != IMAGING_TYPE_UINT8)
+        return ImagingError_ModeError();
+
+    if (!(strcmp(imIn->mode, "RGB") == 0 ||
+          strcmp(imIn->mode, "RGBA") == 0 ||
+          strcmp(imIn->mode, "RGBX") == 0 ||
+          strcmp(imIn->mode, "CMYK") == 0 ||
+          strcmp(imIn->mode, "L") == 0 ||
+          strcmp(imIn->mode, "LA") == 0))
+        return ImagingError_ModeError();
+
+    Imaging imTransposed = ImagingNew(imIn->mode, imIn->ysize, imIn->xsize);
+    if (!imTransposed)
+        return NULL;
+
+    /* Apply blur in one dimension.
+       Use imOut as a destination at first pass,
+       then use imOut as a source too. */
+    ImagingHorizontalBoxBlur(imOut, imIn, radius);
+    for (i = 1; i < n; i ++) {
+        ImagingHorizontalBoxBlur(imOut, imOut, radius);
+    }
+    /* Transpose result for blur in another direction. */
+    ImagingTranspose(imTransposed, imOut);
+
+    /* Reuse imTransposed as a source and destination there. */
+    for (i = 0; i < n; i ++) {
+        ImagingHorizontalBoxBlur(imTransposed, imTransposed, radius);
+    }
+    /* Restore original orientation. */
+    ImagingTranspose(imOut, imTransposed);
+
+    ImagingDelete(imTransposed);
+
+    return imOut;
+}
+
+
+Imaging ImagingGaussianBlur(Imaging imOut, Imaging imIn, float radius,
+    int passes)
+{
+    float sigma2, L, l, a;
+
+    sigma2 = radius * radius / passes;
+    // from http://www.mia.uni-saarland.de/Publications/gwosdek-ssvm11.pdf
+    // [7] Box length.
+    L = sqrt(12.0 * sigma2 + 1.0);
+    // [11] Integer part of box radius.
+    l = floor((L - 1.0) / 2.0);
+    // [14], [Fig. 2] Fractional part of box radius.
+    a = (2 * l + 1) * (l * (l + 1) - 3 * sigma2);
+    a /= 6 * (sigma2 - (l + 1) * (l + 1));
+
+    return ImagingBoxBlur(imOut, imIn, l + a, passes);
+}

libImaging/Imaging.h

@@ -263,7 +263,8 @@ extern Imaging ImagingFilter(
     FLOAT32 offset, FLOAT32 divisor);
 extern Imaging ImagingFlipLeftRight(Imaging imOut, Imaging imIn);
 extern Imaging ImagingFlipTopBottom(Imaging imOut, Imaging imIn);
-extern Imaging ImagingGaussianBlur(Imaging im, Imaging imOut, float radius);
+extern Imaging ImagingGaussianBlur(Imaging imOut, Imaging imIn, float radius,
+    int passes);
 extern Imaging ImagingGetBand(Imaging im, int band);
 extern int ImagingGetBBox(Imaging im, int bbox[4]);
 typedef struct { int x, y; INT32 count; INT32 pixel; } ImagingColorItem;
@@ -309,7 +310,8 @@ extern Imaging ImagingTransform(
     ImagingTransformFilter filter, void* filter_data,
     int fill);
 extern Imaging ImagingUnsharpMask(
-    Imaging im, Imaging imOut, float radius, int percent, int threshold);
+    Imaging imOut, Imaging im, float radius, int percent, int threshold);
+extern Imaging ImagingBoxBlur(Imaging imOut, Imaging imIn, float radius, int n);
 
 extern Imaging ImagingCopy2(Imaging imOut, Imaging imIn);
 extern Imaging ImagingConvert2(Imaging imOut, Imaging imIn);

libImaging/UnsharpMask.c

@@ -9,385 +9,83 @@
 #include "Python.h"
 #include "Imaging.h"
 
-#define PILUSMVERSION "0.6.1"
 
-/* version history
+typedef UINT8 pixel[4];
 
-0.6.1   converted to C and added to PIL 1.1.7
 
-0.6.0   fixed/improved float radius support (oops!)
+static inline UINT8 clip8(int in)
-        now that radius can be a float (properly), changed radius value to
-            be an actual radius (instead of diameter).  So, you should get
-            similar results from PIL_usm as from other paint programs when
-            using the SAME values (no doubling of radius required any more).
-            Be careful, this may "break" software if you had it set for 2x
-            or 5x the radius as was recommended with earlier versions.
-        made PILusm thread-friendly (release GIL before lengthly operations,
-            and re-acquire it before returning to Python).  This makes a huge
-            difference with multi-threaded applications on dual-processor
-            or "Hyperthreading"-enabled systems (Pentium4, Xeon, etc.)
-
-0.5.0   added support for float radius values!
-
-0.4.0   tweaked gaussian curve calculation to be closer to consistent shape
-            across a wide range of radius values
-
-0.3.0   changed deviation calculation in gausian algorithm to be dynamic
-        _gblur now adds 1 to the user-supplied radius before using it so
-            that a value of "0" returns the original image instead of a
-            black one.
-        fixed handling of alpha channel in RGBX, RGBA images
-        improved speed of gblur by reducing unnecessary checks and assignments
-
-0.2.0   fixed L-mode image support
-
-0.1.0   initial release
-
-*/
-
-static inline UINT8 clip(double in)
 {
-    if (in >= 255.0)
+    if (in >= 255)
-	return (UINT8) 255;
+       return 255;
-    if (in <= 0.0)
+    if (in <= 0)
-	return (UINT8) 0;
+        return 0;
     return (UINT8) in;
 }
 
-static Imaging
-gblur(Imaging im, Imaging imOut, float floatRadius, int channels, int padding)
-{
-    ImagingSectionCookie cookie;
-
-    float *maskData = NULL;
-    int y = 0;
-    int x = 0;
-    float z = 0;
-    float sum = 0.0;
-    float dev = 0.0;
-
-    float *buffer = NULL;
-
-    int *line = NULL;
-    UINT8 *line8 = NULL;
-
-    int pix = 0;
-    float newPixel[4];
-    int channel = 0;
-    int offset = 0;
-    INT32 newPixelFinals;
-
-    int radius = 0;
-    float remainder = 0.0;
-
-    int i;
-
-    /* Do the gaussian blur */
-
-    /* For a symmetrical gaussian blur, instead of doing a radius*radius
-       matrix lookup, you get the EXACT same results by doing a radius*1
-       transform, followed by a 1*radius transform.  This reduces the
-       number of lookups exponentially (10 lookups per pixel for a
-       radius of 5 instead of 25 lookups).  So, we blur the lines first,
-       then we blur the resulting columns. */
-
-    /* first, round radius off to the next higher integer and hold the
-       remainder this is used so we can support float radius values
-       properly. */
-
-    remainder = floatRadius - ((int) floatRadius);
-    floatRadius = ceil(floatRadius);
-
-    /* Next, double the radius and offset by 2.0... that way "0" returns
-       the original image instead of a black one.  We multiply it by 2.0
-       so that it is a true "radius", not a diameter (the results match
-       other paint programs closer that way too). */
-    radius = (int) ((floatRadius * 2.0) + 2.0);
-
-    /* create the maskData for the gaussian curve */
-    maskData = malloc(radius * sizeof(float));
-    /* FIXME: error checking */
-    for (x = 0; x < radius; x++) {
-	z = ((float) (x + 2) / ((float) radius));
-	dev = 0.5 + (((float) (radius * radius)) * 0.001);
-	/* you can adjust this factor to change the shape/center-weighting
-	   of the gaussian */
-	maskData[x] = (float) pow((1.0 / sqrt(2.0 * 3.14159265359 * dev)),
-				  ((-(z - 1.0) * -(x - 1.0)) /
-				   (2.0 * dev)));
-    }
-
-    /* if there's any remainder, multiply the first/last values in
-       MaskData it.  this allows us to support float radius values. */
-    if (remainder > 0.0) {
-	maskData[0] *= remainder;
-	maskData[radius - 1] *= remainder;
-    }
-
-    for (x = 0; x < radius; x++) {
-	/* this is done separately now due to the correction for float
-	   radius values above */
-	sum += maskData[x];
-    }
-
-    for (i = 0; i < radius; i++) {
-	maskData[i] *= (1.0 / sum);
-	/* printf("%f\n", maskData[i]); */
-    }
-
-    /* create a temporary memory buffer for the data for the first pass
-       memset the buffer to 0 so we can use it directly with += */
-
-    /* don't bother about alpha/padding */
-    buffer = calloc((size_t) (im->xsize * im->ysize * channels),
-		    sizeof(float));
-    if (buffer == NULL)
-	return ImagingError_MemoryError();
-
-    /* be nice to other threads while you go off to lala land */
-    ImagingSectionEnter(&cookie);
-
-    /* memset(buffer, 0, sizeof(buffer)); */
-
-    newPixel[0] = newPixel[1] = newPixel[2] = newPixel[3] = 0;
-
-    /* perform a blur on each line, and place in the temporary storage buffer */
-    for (y = 0; y < im->ysize; y++) {
-	if (channels == 1 && im->image8 != NULL) {
-	    line8 = (UINT8 *) im->image8[y];
-	} else {
-	    line = im->image32[y];
-	}
-	for (x = 0; x < im->xsize; x++) {
-	    newPixel[0] = newPixel[1] = newPixel[2] = newPixel[3] = 0;
-	    /* for each neighbor pixel, factor in its value/weighting to the
-	       current pixel */
-	    for (pix = 0; pix < radius; pix++) {
-		/* figure the offset of this neighbor pixel */
-		offset =
-		    (int) ((-((float) radius / 2.0) + (float) pix) + 0.5);
-		if (x + offset < 0)
-		    offset = -x;
-		else if (x + offset >= im->xsize)
-		    offset = im->xsize - x - 1;
-
-		/* add (neighbor pixel value * maskData[pix]) to the current
-		   pixel value */
-		if (channels == 1) {
-		    buffer[(y * im->xsize) + x] +=
-			((float) ((UINT8 *) & line8[x + offset])[0]) *
-			(maskData[pix]);
-		} else {
-		    for (channel = 0; channel < channels; channel++) {
-			buffer[(y * im->xsize * channels) +
-			       (x * channels) + channel] +=
-			    ((float) ((UINT8 *) & line[x + offset])
-			     [channel]) * (maskData[pix]);
-		    }
-		}
-	    }
-	}
-    }
-
-    /* perform a blur on each column in the buffer, and place in the
-       output image */
-    for (x = 0; x < im->xsize; x++) {
-	for (y = 0; y < im->ysize; y++) {
-	    newPixel[0] = newPixel[1] = newPixel[2] = newPixel[3] = 0;
-	    /* for each neighbor pixel, factor in its value/weighting to the
-	       current pixel */
-	    for (pix = 0; pix < radius; pix++) {
-		/* figure the offset of this neighbor pixel */
-		offset =
-		    (int) (-((float) radius / 2.0) + (float) pix + 0.5);
-		if (y + offset < 0)
-		    offset = -y;
-		else if (y + offset >= im->ysize)
-		    offset = im->ysize - y - 1;
-		/* add (neighbor pixel value * maskData[pix]) to the current
-		   pixel value */
-		for (channel = 0; channel < channels; channel++) {
-		    newPixel[channel] +=
-			(buffer
-			 [((y + offset) * im->xsize * channels) +
-			  (x * channels) + channel]) * (maskData[pix]);
-		}
-	    }
-	    /* if the image is RGBX or RGBA, copy the 4th channel data to
-	       newPixel, so it gets put in imOut */
-	    if (strcmp(im->mode, "RGBX") == 0
-		|| strcmp(im->mode, "RGBA") == 0) {
-	      newPixel[3] = (float) ((UINT8 *) & line[x + offset])[3];
-	    }
-
-	    /* pack the channels into an INT32 so we can put them back in
-	       the PIL image */
-	    newPixelFinals = 0;
-	    if (channels == 1) {
-		newPixelFinals = clip(newPixel[0]);
-	    } else {
-		/* for RGB, the fourth channel isn't used anyways, so just
-		   pack a 0 in there, this saves checking the mode for each
-		   pixel. */
-		/* this doesn't work on little-endian machines... fix it! */
-		newPixelFinals =
-		    clip(newPixel[0]) | clip(newPixel[1]) << 8 |
-		    clip(newPixel[2]) << 16 | clip(newPixel[3]) << 24;
-	    }
-	    /* set the resulting pixel in imOut */
-	    if (channels == 1) {
-		imOut->image8[y][x] = (UINT8) newPixelFinals;
-	    } else {
-		imOut->image32[y][x] = newPixelFinals;
-	    }
-	}
-    }
-
-    /* free the buffer */
-    free(buffer);
-
-    /* get the GIL back so Python knows who you are */
-    ImagingSectionLeave(&cookie);
-
-    return imOut;
-}
-
-Imaging ImagingGaussianBlur(Imaging im, Imaging imOut, float radius)
-{
-    int channels = 0;
-    int padding = 0;
-
-    if (strcmp(im->mode, "RGB") == 0) {
-	channels = 3;
-	padding = 1;
-    } else if (strcmp(im->mode, "RGBA") == 0) {
-	channels = 3;
-	padding = 1;
-    } else if (strcmp(im->mode, "RGBX") == 0) {
-	channels = 3;
-	padding = 1;
-    } else if (strcmp(im->mode, "CMYK") == 0) {
-	channels = 4;
-	padding = 0;
-    } else if (strcmp(im->mode, "L") == 0) {
-	channels = 1;
-	padding = 0;
-    } else
-	return ImagingError_ModeError();
-
-    return gblur(im, imOut, radius, channels, padding);
-}
 
 Imaging
-ImagingUnsharpMask(Imaging im, Imaging imOut, float radius, int percent,
+ImagingUnsharpMask(Imaging imOut, Imaging imIn, float radius, int percent,
                    int threshold)
 {
     ImagingSectionCookie cookie;
-
     Imaging result;
-    int channel = 0;
-    int channels = 0;
-    int padding = 0;
 
-    int x = 0;
+    int x, y, diff;
-    int y = 0;
 
-    int *lineIn = NULL;
+    pixel *lineIn = NULL;
-    int *lineOut = NULL;
+    pixel *lineOut = NULL;
     UINT8 *lineIn8 = NULL;
     UINT8 *lineOut8 = NULL;
 
-    int diff = 0;
+    /* First, do a gaussian blur on the image, putting results in imOut
-
+       temporarily. All format checks are in gaussian blur. */
-    INT32 newPixel = 0;
+    result = ImagingGaussianBlur(imOut, imIn, radius, 3);
-
-    if (strcmp(im->mode, "RGB") == 0) {
-	channels = 3;
-	padding = 1;
-    } else if (strcmp(im->mode, "RGBA") == 0) {
-	channels = 3;
-	padding = 1;
-    } else if (strcmp(im->mode, "RGBX") == 0) {
-	channels = 3;
-	padding = 1;
-    } else if (strcmp(im->mode, "CMYK") == 0) {
-	channels = 4;
-	padding = 0;
-    } else if (strcmp(im->mode, "L") == 0) {
-	channels = 1;
-	padding = 0;
-    } else
-	return ImagingError_ModeError();
-
-    /* first, do a gaussian blur on the image, putting results in imOut
-       temporarily */
-    result = gblur(im, imOut, radius, channels, padding);
     if (!result)
         return NULL;
 
-    /* now, go through each pixel, compare "normal" pixel to blurred
+    /* Now, go through each pixel, compare "normal" pixel to blurred
-       pixel.  if the difference is more than threshold values, apply
+       pixel. If the difference is more than threshold values, apply
        the OPPOSITE correction to the amount of blur, multiplied by
        percent. */
 
     ImagingSectionEnter(&cookie);
 
-    for (y = 0; y < im->ysize; y++) {
+    for (y = 0; y < imIn->ysize; y++) {
-	if (channels == 1) {
+        if (imIn->image8)
-	    lineIn8 = im->image8[y];
+        {
+            lineIn8 = imIn->image8[y];
             lineOut8 = imOut->image8[y];
-	} else {
+            for (x = 0; x < imIn->xsize; x++) {
-	    lineIn = im->image32[y];
-	    lineOut = imOut->image32[y];
-	}
-	for (x = 0; x < im->xsize; x++) {
-	    newPixel = 0;
                 /* compare in/out pixels, apply sharpening */
-	    if (channels == 1) {
+                diff = lineIn8[x] - lineOut8[x];
-		diff =
-		    ((UINT8 *) & lineIn8[x])[0] -
-		    ((UINT8 *) & lineOut8[x])[0];
                 if (abs(diff) > threshold) {
                     /* add the diff*percent to the original pixel */
-		    imOut->image8[y][x] =
+                    lineOut8[x] = clip8(lineIn8[x] + diff * percent / 100);
-			clip((((UINT8 *) & lineIn8[x])[0]) +
-			     (diff * ((float) percent) / 100.0));
                 } else {
-		    /* newPixel is the same as imIn */
+                    /* new pixel is the same as imIn */
-		    imOut->image8[y][x] = ((UINT8 *) & lineIn8[x])[0];
+                    lineOut8[x] = lineIn8[x];
                 }
             }
+        } else {
+            lineIn = (pixel *)imIn->image32[y];
+            lineOut = (pixel *)imOut->image32[y];
+            for (x = 0; x < imIn->xsize; x++) {
+                /* compare in/out pixels, apply sharpening */
+                diff = lineIn[x][0] - lineOut[x][0];
+                lineOut[x][0] = abs(diff) > threshold ?
+                    clip8(lineIn[x][0] + diff * percent / 100) : lineIn[x][0];
 
-	    else {
+                diff = lineIn[x][1] - lineOut[x][1];
-		for (channel = 0; channel < channels; channel++) {
+                lineOut[x][1] = abs(diff) > threshold ?
-		    diff = (int) ((((UINT8 *) & lineIn[x])[channel]) -
+                    clip8(lineIn[x][1] + diff * percent / 100) : lineIn[x][1];
-				  (((UINT8 *) & lineOut[x])[channel]));
+
-		    if (abs(diff) > threshold) {
+                diff = lineIn[x][2] - lineOut[x][2];
-			/* add the diff*percent to the original pixel
+                lineOut[x][2] = abs(diff) > threshold ?
-			   this may not work for little-endian systems, fix it! */
+                    clip8(lineIn[x][2] + diff * percent / 100) : lineIn[x][2];
-			newPixel =
+
-			    newPixel |
+                diff = lineIn[x][3] - lineOut[x][3];
-			    clip((float) (((UINT8 *) & lineIn[x])[channel])
+                lineOut[x][3] = abs(diff) > threshold ?
-				 +
+                    clip8(lineIn[x][3] + diff * percent / 100) : lineIn[x][3];
-				 (diff *
-				  (((float) percent /
-				    100.0)))) << (channel * 8);
-		    } else {
-			/* newPixel is the same as imIn
-			   this may not work for little-endian systems, fix it! */
-			newPixel =
-			    newPixel | ((UINT8 *) & lineIn[x])[channel] <<
-			    (channel * 8);
-		    }
-		}
-		if (strcmp(im->mode, "RGBX") == 0
-		    || strcmp(im->mode, "RGBA") == 0) {
-		    /* preserve the alpha channel
-		       this may not work for little-endian systems, fix it! */
-		    newPixel =
-			newPixel | ((UINT8 *) & lineIn[x])[channel] << 24;
-		}
-		imOut->image32[y][x] = newPixel;
             }
         }
     }

setup.py

@@ -37,7 +37,7 @@ _LIB_IMAGING = (
     "RankFilter", "RawDecode", "RawEncode", "Storage", "SunRleDecode",
     "TgaRleDecode", "Unpack", "UnpackYCC", "UnsharpMask", "XbmDecode",
     "XbmEncode", "ZipDecode", "ZipEncode", "TiffDecode", "Incremental",
-    "Jpeg2KDecode", "Jpeg2KEncode")
+    "Jpeg2KDecode", "Jpeg2KEncode", "BoxBlur")
 
 
 def _add_directory(path, dir, where=None):