source: libpipi/trunk/pipi/filter/blur.c @ 2902

Last change on this file since 2902 was 2902, checked in by Sam Hocevar, 12 years ago

Support C99 types on Win32 through the same hacks as in libcaca.

File size: 8.7 KB
Line 
1/*
2 *  libpipi       Pathetic image processing interface library
3 *  Copyright (c) 2004-2008 Sam Hocevar <sam@zoy.org>
4 *                All Rights Reserved
5 *
6 *  $Id$
7 *
8 *  This library is free software. It comes without any warranty, to
9 *  the extent permitted by applicable law. You can redistribute it
10 *  and/or modify it under the terms of the Do What The Fuck You Want
11 *  To Public License, Version 2, as published by Sam Hocevar. See
12 *  http://sam.zoy.org/wtfpl/COPYING for more details.
13 */
14
15/*
16 * blur.c: blur functions
17 */
18
19#include "config.h"
20
21#include <stdlib.h>
22#include <stdio.h>
23#include <string.h>
24#include <math.h>
25
26#include "pipi.h"
27#include "pipi_internals.h"
28
29#if !defined TEMPLATE_FILE /* This file uses the template system */
30#define TEMPLATE_FLAGS SET_FLAG_GRAY | SET_FLAG_WRAP
31#define TEMPLATE_FILE "filter/blur.c"
32#include "pipi_template.h"
33
34/* Any standard deviation below this value will be rounded up, in order
35 * to avoid ridiculously low values. exp(-1/(2*0.2*0.2)) is < 10^-5 so
36 * there is little chance that any value below 0.2 will be useful. */
37#define BLUR_EPSILON 0.2
38
39pipi_image_t *pipi_gaussian_blur(pipi_image_t *src, float radius)
40{
41    return pipi_gaussian_blur_ext(src, radius, radius, 0.0, 0.0, 0.0);
42}
43
44pipi_image_t *pipi_gaussian_blur_ext(pipi_image_t *src, float rx, float ry,
45                                     float angle, float dx, float dy)
46{
47    pipi_image_t *ret;
48    double *kernel;
49    double Kx, Ky, t = 0.0, sint, cost, bbx, bby;
50    int i, j, krx, kry, m, n;
51
52    if(rx < BLUR_EPSILON) rx = BLUR_EPSILON;
53    if(ry < BLUR_EPSILON) ry = BLUR_EPSILON;
54
55    sint = sin(angle * M_PI / 180.);
56    cost = cos(angle * M_PI / 180.);
57
58    /* Compute the final ellipse's bounding box */
59    bbx = sqrt(rx * rx * cost * cost + ry * ry * sint * sint);
60    bby = sqrt(ry * ry * cost * cost + rx * rx * sint * sint);
61
62    /* FIXME: the kernel becomes far too big with large values of dx, because
63     * we grow both left and right. Fix the growing direction. */
64    krx = (int)(3. * bbx + .99999 + ceil(abs(dx)));
65    m = 2 * krx + 1;
66    Kx = -1. / (2. * rx * rx);
67
68    kry = (int)(3. * bby + .99999 + ceil(abs(dy)));
69    n = 2 * kry + 1;
70    Ky = -1. / (2. * ry * ry);
71
72    kernel = malloc(m * n * sizeof(double));
73
74    for(j = -kry; j <= kry; j++)
75    {
76        for(i = -krx; i <= krx; i++)
77        {
78            /* FIXME: this level of interpolation sucks. We should
79             * interpolate on the full NxN grid for better quality. */
80            static double const samples[] =
81            {
82                 .0,  .0, 1,
83                -.40, -.40, 0.8,
84                -.30,  .0,  0.9,
85                -.40,  .40, 0.8,
86                 .0,   .30, 0.9,
87                 .40,  .40, 0.8,
88                 .30,  .0,  0.9,
89                 .40, -.40, 0.8,
90                 .0,  -.30, 0.9,
91            };
92            double u, v, ex, ey, d = 0.;
93            unsigned int k;
94
95            for(k = 0; k < sizeof(samples) / sizeof(*samples) / 3; k++)
96            {
97                u = ((double)i + samples[k * 3]) * cost
98                     - ((double)j + samples[k * 3 + 1]) * sint + dx;
99                v = ((double)i + samples[k * 3]) * sint
100                     + ((double)j + samples[k * 3 + 1]) * cost + dy;
101                ex = Kx * u * u;
102                ey = Ky * v * v;
103                d += samples[k * 3 + 2] * exp(ex + ey);
104
105                /* Do not interpolate if this is a standard gaussian. */
106                if(!dx && !dy && !angle)
107                    break;
108            }
109
110            kernel[(j + kry) * m + i + krx] = d;
111            t += d;
112        }
113    }
114
115    for(j = 0; j < n; j++)
116        for(i = 0; i < m; i++)
117            kernel[j * m + i] /= t;
118
119    ret = pipi_convolution(src, m, n, kernel);
120
121    free(kernel);
122
123    return ret;
124}
125
126pipi_image_t *pipi_box_blur(pipi_image_t *src, int size)
127{
128    return pipi_box_blur_ext(src, size, size);
129}
130
131pipi_image_t *pipi_box_blur_ext(pipi_image_t *src, int m, int n)
132{
133    if(src->wrap)
134    {
135        if(src->last_modified == PIPI_PIXELS_Y_F)
136            return boxblur_gray_wrap(src, m, n);
137
138        return boxblur_wrap(src, m, n);
139    }
140    else
141    {
142        if(src->last_modified == PIPI_PIXELS_Y_F)
143            return boxblur_gray(src, m, n);
144
145        return boxblur(src, m, n);
146    }
147}
148
149#else /* XXX: the following functions use the template system */
150
151static pipi_image_t *T(boxblur)(pipi_image_t *src, int m, int n)
152{
153    pipi_image_t *dst;
154    pipi_pixels_t *srcp, *dstp;
155    float *srcdata, *dstdata;
156    double *acc;
157    int x, y, w, h, i, j, i2, j2, size;
158
159    w = src->w;
160    h = src->h;
161    size = (2 * m + 1) * (2 * n + 1);
162
163    srcp = FLAG_GRAY ? pipi_getpixels(src, PIPI_PIXELS_Y_F)
164                     : pipi_getpixels(src, PIPI_PIXELS_RGBA_F);
165    srcdata = (float *)srcp->pixels;
166
167    dst = pipi_new(w, h);
168    dstp = FLAG_GRAY ? pipi_getpixels(dst, PIPI_PIXELS_Y_F)
169                     : pipi_getpixels(dst, PIPI_PIXELS_RGBA_F);
170    dstdata = (float *)dstp->pixels;
171
172    acc = malloc(w * (FLAG_GRAY ? 1 : 4) * sizeof(double));
173
174    /* Step 1: fill the accumulator */
175    for(x = 0; x < w; x++)
176    {
177        double r = 0., g = 0., b = 0., a = 0.;
178        double t = 0.;
179
180        for(j = -n; j <= n; j++)
181        {
182            if(FLAG_WRAP)
183                j2 = (j < 0) ? h - 1 - ((-j - 1) % h) : j % h;
184            else
185                j2 = (j < 0) ? 0 : (j >= h) ? h - 1 : j;
186
187            if(FLAG_GRAY)
188                t += srcdata[j2 * w + x];
189            else
190            {
191                r += srcdata[4 * (j2 * w + x)];
192                g += srcdata[4 * (j2 * w + x) + 1];
193                b += srcdata[4 * (j2 * w + x) + 2];
194                a += srcdata[4 * (j2 * w + x) + 3];
195            }
196        }
197
198        if(FLAG_GRAY)
199            acc[x] = t;
200        else
201        {
202            acc[4 * x] = r;
203            acc[4 * x + 1] = g;
204            acc[4 * x + 2] = b;
205            acc[4 * x + 3] = a;
206        }
207    }
208
209    /* Step 2: blur the image, line by line */
210    for(y = 0; y < h; y++)
211    {
212        double r = 0., g = 0., b = 0., a = 0.;
213        double t = 0.;
214
215        /* 2.1: compute the first pixel */
216        for(i = -m; i <= m; i++)
217        {
218            if(FLAG_WRAP)
219                i2 = (i < 0) ? w - 1 - ((-i - 1) % w) : i % w;
220            else
221                i2 = (i < 0) ? 0 : (i >= w) ? w - 1 : i;
222
223            if(FLAG_GRAY)
224                t += acc[i2];
225            else
226            {
227                r += acc[4 * i2];
228                g += acc[4 * i2 + 1];
229                b += acc[4 * i2 + 2];
230                a += acc[4 * i2 + 3];
231            }
232        }
233
234        /* 2.2: iterate on the whole line */
235        for(x = 0; x < w; x++)
236        {
237            int u, u2, v, v2;
238
239            if(FLAG_GRAY)
240            {
241                dstdata[y * w + x] = t / size;
242            }
243            else
244            {
245                dstdata[4 * (y * w + x)] = r / size;
246                dstdata[4 * (y * w + x) + 1] = g / size;
247                dstdata[4 * (y * w + x) + 2] = b / size;
248                dstdata[4 * (y * w + x) + 3] = a / size;
249            }
250
251            u = x - m;
252            if(FLAG_WRAP)
253                u2 = (u < 0) ? w - 1 - ((-u - 1) % w) : u % w;
254            else
255                u2 = (u < 0) ? 0 : (u >= w) ? w - 1 : u;
256            v = x + m + 1;
257            if(FLAG_WRAP)
258                v2 = (v < 0) ? w - 1 - ((-v - 1) % w) : v % w;
259            else
260                v2 = (v < 0) ? 0 : (v >= w) ? w - 1 : v;
261            if(FLAG_GRAY)
262            {
263                t = t - acc[u2] + acc[v2];
264            }
265            else
266            {
267                r = r - acc[4 * u2] + acc[4 * v2];
268                g = g - acc[4 * u2 + 1] + acc[4 * v2 + 1];
269                b = b - acc[4 * u2 + 2] + acc[4 * v2 + 2];
270                a = a - acc[4 * u2 + 3] + acc[4 * v2 + 3];
271            }
272        }
273
274        /* 2.3: update the accumulator */
275        for(x = 0; x < w; x++)
276        {
277            int u, u2, v, v2;
278
279            u = y - n;
280            if(FLAG_WRAP)
281                u2 = (u < 0) ? w - 1 - ((-u - 1) % w) : u % w;
282            else
283                u2 = (u < 0) ? 0 : (u >= w) ? w - 1 : u;
284            v = y + n + 1;
285            if(FLAG_WRAP)
286                v2 = (v < 0) ? w - 1 - ((-v - 1) % w) : v % w;
287            else
288                v2 = (v < 0) ? 0 : (v >= w) ? w - 1 : v;
289            if(FLAG_GRAY)
290            {
291                acc[x] += srcdata[v2 * w + x] - srcdata[u2 * w + x];
292            }
293            else
294            {
295                int uoff = 4 * (u2 * w + x);
296                int voff = 4 * (v2 * w + x);
297
298                acc[4 * x] += srcdata[voff] - srcdata[uoff];
299                acc[4 * x + 1] += srcdata[voff + 1] - srcdata[uoff + 1];
300                acc[4 * x + 2] += srcdata[voff + 2] - srcdata[uoff + 2];
301                acc[4 * x + 3] += srcdata[voff + 3] - srcdata[uoff + 3];
302            }
303        }
304    }
305
306    free(acc);
307
308    return dst;
309}
310
311#endif
312
Note: See TracBrowser for help on using the repository browser.