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

Last change on this file since 2643 was 2643, checked in by Sam Hocevar, 14 years ago
  • blur.c: adapt the kernel size to large values of dx and/or dy.
  • blur.c: fix the dx/dy sign meaning.
File size: 5.1 KB
Line 
1/*
2 *  libpipi       Proper image processing implementation 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#include "common.h"
21
22#include <stdlib.h>
23#include <stdio.h>
24#include <string.h>
25#include <math.h>
26
27#include "pipi.h"
28#include "pipi_internals.h"
29
30/* Any standard deviation below this value will be rounded up, in order
31 * to avoid ridiculously low values. exp(-1/(2*0.2*0.2)) is < 10^-5 so
32 * there is little chance that any value below 0.2 will be useful. */
33#define BLUR_EPSILON 0.2
34
35pipi_image_t *pipi_gaussian_blur(pipi_image_t *src, float radius)
36{
37    return pipi_gaussian_blur_ext(src, radius, radius, 0.0, 0.0);
38}
39
40pipi_image_t *pipi_gaussian_blur_ext(pipi_image_t *src, float rx, float ry,
41                                     float dx, float dy)
42{
43    pipi_image_t *dst;
44    pipi_pixels_t *srcp, *dstp;
45    float *srcdata, *dstdata;
46    double *kernel, *buffer;
47    double K;
48    int x, y, i, w, h, kr, kw, gray;
49
50    if(rx < BLUR_EPSILON) rx = BLUR_EPSILON;
51    if(ry < BLUR_EPSILON) ry = BLUR_EPSILON;
52
53    w = src->w;
54    h = src->h;
55
56    gray = (src->last_modified == PIPI_PIXELS_Y_F);
57
58    srcp = gray ? pipi_getpixels(src, PIPI_PIXELS_Y_F)
59                : pipi_getpixels(src, PIPI_PIXELS_RGBA_F);
60    srcdata = (float *)srcp->pixels;
61
62    dst = pipi_new(w, h);
63    dstp = gray ? pipi_getpixels(dst, PIPI_PIXELS_Y_F)
64                : pipi_getpixels(dst, PIPI_PIXELS_RGBA_F);
65    dstdata = (float *)dstp->pixels;
66
67    buffer = malloc(w * h * (gray ? 1 : 4) * sizeof(double));
68
69    /* FIXME: the kernel becomes far too big with large values of dx, because
70     * we grow both left and right. Fix the growing direction. */
71    kr = (int)(3. * rx + .99999 + ceil(abs(dx)));
72    kw = 2 * kr + 1;
73    K = -1. / (2. * rx * rx);
74
75    kernel = malloc(kw * sizeof(double));
76    for(i = -kr; i <= kr; i++)
77        kernel[i + kr] = exp(K * ((double)i + dx) * ((double)i + dx));
78
79    for(y = 0; y < h; y++)
80    {
81        for(x = 0; x < w; x++)
82        {
83            if(gray)
84            {
85                double Y = 0., t = 0.;
86                int x2;
87
88                for(i = -kr; i <= kr; i++)
89                {
90                    double f = kernel[i + kr];
91
92                    x2 = x + i;
93                    if(x2 < 0) x2 = 0;
94                    else if(x2 >= w) x2 = w - 1;
95
96                    Y += f * srcdata[y * w + x2];
97                    t += f;
98                }
99
100                buffer[y * w + x] = Y / t;
101            }
102            else
103            {
104                double R = 0., G = 0., B = 0., t = 0.;
105                int x2, off = 4 * (y * w + x);
106
107                for(i = -kr; i <= kr; i++)
108                {
109                    double f = kernel[i + kr];
110
111                    x2 = x + i;
112                    if(x2 < 0) x2 = 0;
113                    else if(x2 >= w) x2 = w - 1;
114
115                    R += f * srcdata[(y * w + x2) * 4];
116                    G += f * srcdata[(y * w + x2) * 4 + 1];
117                    B += f * srcdata[(y * w + x2) * 4 + 2];
118                    t += f;
119                }
120
121                buffer[off] = R / t;
122                buffer[off + 1] = G / t;
123                buffer[off + 2] = B / t;
124            }
125        }
126    }
127
128    free(kernel);
129
130    kr = (int)(3. * ry + .99999 + ceil(abs(dy)));
131    kw = 2 * kr + 1;
132    K = -1. / (2. * ry * ry);
133
134    kernel = malloc(kw * sizeof(double));
135    for(i = -kr; i <= kr; i++)
136        kernel[i + kr] = exp(K * ((double)i + dy) * ((double)i + dy));
137
138    for(y = 0; y < h; y++)
139    {
140        for(x = 0; x < w; x++)
141        {
142            if(gray)
143            {
144                double Y = 0., t = 0.;
145                int y2;
146
147                for(i = -kr; i <= kr; i++)
148                {
149                    double f = kernel[i + kr];
150
151                    y2 = y + i;
152                    if(y2 < 0) y2 = 0;
153                    else if(y2 >= h) y2 = h - 1;
154
155                    Y += f * buffer[y2 * w + x];
156                    t += f;
157                }
158
159                dstdata[y * w + x] = Y / t;
160            }
161            else
162            {
163                double R = 0., G = 0., B = 0., t = 0.;
164                int y2, off = 4 * (y * w + x);
165
166                for(i = -kr; i <= kr; i++)
167                {
168                    double f = kernel[i + kr];
169
170                    y2 = y + i;
171                    if(y2 < 0) y2 = 0;
172                    else if(y2 >= h) y2 = h - 1;
173
174                    R += f * buffer[(y2 * w + x) * 4];
175                    G += f * buffer[(y2 * w + x) * 4 + 1];
176                    B += f * buffer[(y2 * w + x) * 4 + 2];
177                    t += f;
178                }
179
180                dstdata[off] = R / t;
181                dstdata[off + 1] = G / t;
182                dstdata[off + 2] = B / t;
183            }
184        }
185    }
186
187    free(buffer);
188    free(kernel);
189
190    return dst;
191}
192
Note: See TracBrowser for help on using the repository browser.