* blur.c: proper box blur implementation that runs linearly with

pixel count.

git-svn-id: file:///srv/caca.zoy.org/var/lib/svn/libpipi/trunk@2743 92316355-f0b4-4df1-b90c-862c8a59935f

pipi/filter/blur.c

@@ -115,27 +115,161 @@ pipi_image_t *pipi_gaussian_blur_ext(pipi_image_t *src, float rx, float ry,
     return ret;
 }
 
-/* FIXME: box blur would be incredibly faster using an accumulator instead
- * of a convolution filter... */
 pipi_image_t *pipi_box_blur(pipi_image_t *src, int size)
 {
     return pipi_box_blur_ext(src, size, size);
 }
 
+/* FIXME: split this into templates for wrap-around and proper gray support */
 pipi_image_t *pipi_box_blur_ext(pipi_image_t *src, int m, int n)
 {
-    pipi_image_t *ret;
-    double *kernel;
-    int i;
+    pipi_image_t *dst;
+    pipi_pixels_t *srcp, *dstp;
+    float *srcdata, *dstdata;
+    double *acc;
+    int x, y, w, h, i, j, size, gray;
 
-    kernel = malloc(m * n * sizeof(double));
-    for(i = 0; i < m * n; i++)
-        kernel[i] = 1. / (m * n);
+    w = src->w;
+    h = src->h;
+    size = (2 * m + 1) * (2 * n + 1);
 
-    ret = pipi_convolution(src, m, n, kernel);
+    gray = (src->last_modified == PIPI_PIXELS_Y_F);
 
-    free(kernel);
+    srcp = gray ? pipi_getpixels(src, PIPI_PIXELS_Y_F)
+                : pipi_getpixels(src, PIPI_PIXELS_RGBA_F);
+    srcdata = (float *)srcp->pixels;
 
-    return ret;
+    dst = pipi_new(w, h);
+    dstp = gray ? pipi_getpixels(dst, PIPI_PIXELS_Y_F)
+                : pipi_getpixels(dst, PIPI_PIXELS_RGBA_F);
+    dstdata = (float *)dstp->pixels;
+
+    acc = malloc(w * (gray ? 1 : 4) * sizeof(double));
+
+    /* Step 1: fill the accumulator */
+    for(x = 0; x < w; x++)
+    {
+        if(gray)
+        {
+            double t = 0.;
+
+            for(j = -n; j <= n; j++)
+            {
+                int j2 = (j < 0) ? h - 1 - ((-j - 1) % h) : j % h;
+                t += srcdata[j2 * w + x];
+            }
+
+            acc[x] = t;
+        }
+        else
+        {
+            double r = 0., g = 0., b = 0., a = 0.;
+
+            for(j = -n; j <= n; j++)
+            {
+                int j2 = (j < 0) ? h - 1 - ((-j - 1) % h) : j % h;
+                r += srcdata[4 * (j2 * w + x)];
+                g += srcdata[4 * (j2 * w + x) + 1];
+                b += srcdata[4 * (j2 * w + x) + 2];
+                a += srcdata[4 * (j2 * w + x) + 3];
+            }
+
+            acc[4 * x] = r;
+            acc[4 * x + 1] = g;
+            acc[4 * x + 2] = b;
+            acc[4 * x + 3] = a;
+        }
+    }
+
+    /* Step 2: blur the image, line by line */
+    for(y = 0; y < h; y++)
+    {
+        double r = 0., g = 0., b = 0., a = 0.;
+        double t = 0.;
+
+        /* 2.1: compute the first pixel */
+        if(gray)
+        {
+            for(i = -m; i <= m; i++)
+            {
+                int i2 = (i < 0) ? w - 1 - ((-i - 1) % w) : i % w;
+                t += acc[i2];
+            }
+        }
+        else
+        {
+            for(i = -m; i <= m; i++)
+            {
+                int i2 = (i < 0) ? w - 1 - ((-i - 1) % w) : i % w;
+                r += acc[4 * i2];
+                g += acc[4 * i2 + 1];
+                b += acc[4 * i2 + 2];
+                a += acc[4 * i2 + 3];
+            }
+        }
+
+        /* 2.2: iterate on the whole line */
+        for(x = 0; x < w; x++)
+        {
+            int u, u2, v, v2;
+
+            if(gray)
+            {
+                dstdata[y * w + x] = t / size;
+            }
+            else
+            {
+                dstdata[4 * (y * w + x)] = r / size;
+                dstdata[4 * (y * w + x) + 1] = g / size;
+                dstdata[4 * (y * w + x) + 2] = b / size;
+                dstdata[4 * (y * w + x) + 3] = a / size;
+            }
+
+            u = x - m;
+            u2 = (u < 0) ? w - 1 - ((-u - 1) % w) : u % w;
+            v = x + m + 1;
+            v2 = (v < 0) ? w - 1 - ((-v - 1) % w) : v % w;
+            if(gray)
+            {
+                t = t - acc[u2] + acc[v2];
+            }
+            else
+            {
+                r = r - acc[4 * u2] + acc[4 * v2];
+                g = g - acc[4 * u2 + 1] + acc[4 * v2 + 1];
+                b = b - acc[4 * u2 + 2] + acc[4 * v2 + 2];
+                a = a - acc[4 * u2 + 3] + acc[4 * v2 + 3];
+            }
+        }
+
+        /* 2.3: update the accumulator */
+        for(x = 0; x < w; x++)
+        {
+            int u, u2, v, v2;
+
+            u = y - n;
+            u2 = (u < 0) ? w - 1 - ((-u - 1) % w) : u % w;
+            v = y + n + 1;
+            v2 = (v < 0) ? w - 1 - ((-v - 1) % w) : v % w;
+            if(gray)
+            {
+                acc[x] += srcdata[v2 * w + x] - srcdata[u2 * w + x];
+            }
+            else
+            {
+                int uoff = 4 * (u2 * w + x);
+                int voff = 4 * (v2 * w + x);
+
+                acc[4 * x] += srcdata[voff] - srcdata[uoff];
+                acc[4 * x + 1] += srcdata[voff + 1] - srcdata[uoff + 1];
+                acc[4 * x + 2] += srcdata[voff + 2] - srcdata[uoff + 2];
+                acc[4 * x + 3] += srcdata[voff + 3] - srcdata[uoff + 3];
+            }
+        }
+    }
+
+    free(acc);
+
+    return dst;
 }