1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
|
/**
IIR Gauss Filter v1.0
By Stephan Soller <stephan.soller@helionweb.de>
Based on the paper "Recursive implementation of the Gaussian filter" by Ian T. Young and Lucas J. van Vliet.
Licensed under the MIT license
QUICK START
#include ...
#include ...
#define IIR_GAUSS_BLUR_IMPLEMENTATION
#include "iir_gauss_blur.h"
...
int width = 0, height = 0, components = 1;
uint8_t* image = stbi_load("foo.png", &width, &height, &components, 0);
float sigma = 10;
iir_gauss_blur(width, height, components, image, sigma);
stbi_write_png("foo.blurred.png", width, height, components, image, 0);
This example uses stb_image.h to load the image, then blurs it and writes the result using stb_image_write.h.
`sigma` controls the strength of the blur. Higher values give you a blurrier image.
DOCUMENTATION
This is a single header file library. You'll have to define IIR_GAUSS_BLUR_IMPLEMENTATION before including this file to
get the implementation. Otherwise just the header will be included.
The library only has a single function: iir_gauss_blur(width, height, components, image, sigma).
- `width` and `height` are the dimensions of the image in pixels.
- `components` is the number of bytes per pixel. 1 for a grayscale image, 3 for RGB and 4 for RGBA.
The function can handle an arbitrary number of channels, so 2 or 7 will work as well.
- `image` is a pointer to the image data with `width * height` pixels, each pixel having `components` bytes (interleaved
8-bit components). There is no padding between the scanlines of the image.
This is the format used by stb_image.h and stb_image_write.h and easy to work with.
- `sigma` is the strength of the blur. It's a number > 0.5 and most people seem to just eyeball it.
Start with e.g. a sigma of 5 and go up or down until you have the blurriness you want.
There are more informed ways to choose this parameter, see CHOOSING SIGMA below.
The function mallocs an internal float buffer with the same dimensions as the image. If that turns out to be a
bottleneck fell free to move that out of the function. The source code is quite short and straight forward (even if the
math isn't).
The function is an implementation of the paper "Recursive implementation of the Gaussian filter" by Ian T. Young and
Lucas J. van Vliet. It has nothing to do with recursive function calls, instead it's a special way to construct a
filter. Other (convolution based) gauss filters apply a kernel for each pixel and the kernel grows as sigma gets larger.
Meaning their performance degrades the more blurry you want your image to be.
Instead The algorithm in the paper gets it done in just a few passes: A horizontal forward and backward pass and a
vertical forward and backward pass. The work done is independent of the blur radius and so you can have ridiculously
large blur radii without any performance impact.
CHOOSING SIGMA
There seem to be several rules of thumb out there to get a sigma for a given "blur radius". Usually this is something
like `radius = 2 * sigma`. So if you want to have a blur radius of 10px you can use `sigma = (1.0 / 2.0) * radius` to
get the sigma for it (5.0). I'm not sure what that "radius" is supposed to mean though.
For my own projects I came up with two different kinds of blur radii and how to get a sigma for them: Given a big white
area on a black background, how far will the white "bleed out" into the surrounding black? How large is the distance
until the white (255) gets blurred down to something barely visible (smaller than 16) or even to nothing (smaller than
1)? There are to estimates to get the sigma for those radii:
sigma = (1.0 / 1.42) * radius16;
sigma = (1.0 / 3.66) * radius1;
Personally I use `radius16` to calculate the sigma when blurring normal images. Think: I want to blur a pixel across a
circle with the radius x so it's impact is barely visible at the edges.
When I need to calculate padding I use `radius1`: When I have a black border of 100px around the image I can use a
`radius1` of 100 and be reasonable sure that I still got black at the edges. So given a `radius1` blur strength I can
use it as a padding width as well.
I created those estimates by applying different sigmas (1 to 100) to a test image and measuring the effects with GIMP.
So take it with a grain of salt (or many). They're reasonable estimates but by no means exact. I tried to solve the
normal distribution to calculate the perfect sigma but gave up after a lot of confusion. If you know an exact solution
let me know. :)
VERSION HISTORY
v1.0 2018-08-30 Initial release
**/
#ifndef IIR_GAUSS_BLUR_HEADER
#define IIR_GAUSS_BLUR_HEADER
template <typename T>
void iir_gauss_blur(unsigned int width, unsigned int height, unsigned char components, T* image, float sigma);
#endif // IIR_GAUSS_BLUR_HEADER
#ifdef IIR_GAUSS_BLUR_IMPLEMENTATION
#include <stdlib.h>
#include <math.h>
template <typename T>
void iir_gauss_blur(unsigned int width, unsigned int height, unsigned char components, T* image, float sigma) {
// Create IDX macro but push any previous definition (and restore it later) so we don't overwrite a macro the user has possibly defined before us
#pragma push_macro("IDX")
#define IDX(x, y, n) ((y)*width*components + (x)*components + n)
// Allocate buffers
float* buffer = (float*)malloc(width * height * components * sizeof(buffer[0]));
// Calculate filter parameters for a specified sigma
// Use Equation 11b to determine q, do nothing if sigma is to small (should have no effect) or negative (doesn't make sense)
float q;
if (sigma >= 2.5)
q = 0.98711 * sigma - 0.96330;
else if (sigma >= 0.5)
q = 3.97156 - 4.14554 * sqrtf(1.0 - 0.26891 * sigma);
else
return;
// Use equation 8c to determine b0, b1, b2 and b3
float b0 = 1.57825 + 2.44413*q + 1.4281*q*q + 0.422205*q*q*q;
float b1 = 2.44413*q + 2.85619*q*q + 1.26661*q*q*q;
float b2 = -( 1.4281*q*q + 1.26661*q*q*q );
float b3 = 0.422205*q*q*q;
// Use equation 10 to determine B
float B = 1.0 - (b1 + b2 + b3) / b0;
// Horizontal forward pass (from paper: Implement the forward filter with equation 9a)
// The data is loaded from the byte image but stored in the float buffer
for(unsigned int y = 0; y < height; y++) {
float prev1[components], prev2[components], prev3[components];
for(unsigned char n = 0; n < components; n++) {
prev1[n] = image[IDX(0, y, n)];
prev2[n] = prev1[n];
prev3[n] = prev2[n];
}
for(unsigned int x = 0; x < width; x++) {
for(unsigned char n = 0; n < components; n++) {
float val = B * image[IDX(x, y, n)] + (b1 * prev1[n] + b2 * prev2[n] + b3 * prev3[n]) / b0;
buffer[IDX(x, y, n)] = val;
prev3[n] = prev2[n];
prev2[n] = prev1[n];
prev1[n] = val;
}
}
}
// Horizontal backward pass (from paper: Implement the backward filter with equation 9b)
for(unsigned int y = height-1; y < height; y--) {
float prev1[components], prev2[components], prev3[components];
for(unsigned char n = 0; n < components; n++) {
prev1[n] = buffer[IDX(width-1, y, n)];
prev2[n] = prev1[n];
prev3[n] = prev2[n];
}
for(unsigned int x = width-1; x < width; x--) {
for(unsigned char n = 0; n < components; n++) {
float val = B * buffer[IDX(x, y, n)] + (b1 * prev1[n] + b2 * prev2[n] + b3 * prev3[n]) / b0;
buffer[IDX(x, y, n)] = val;
prev3[n] = prev2[n];
prev2[n] = prev1[n];
prev1[n] = val;
}
}
}
// Vertical forward pass (from paper: Implement the forward filter with equation 9a)
for(unsigned int x = 0; x < width; x++) {
float prev1[components], prev2[components], prev3[components];
for(unsigned char n = 0; n < components; n++) {
prev1[n] = buffer[IDX(x, 0, n)];
prev2[n] = prev1[n];
prev3[n] = prev2[n];
}
for(unsigned int y = 0; y < height; y++) {
for(unsigned char n = 0; n < components; n++) {
float val = B * buffer[IDX(x, y, n)] + (b1 * prev1[n] + b2 * prev2[n] + b3 * prev3[n]) / b0;
buffer[IDX(x, y, n)] = val;
prev3[n] = prev2[n];
prev2[n] = prev1[n];
prev1[n] = val;
}
}
}
// Vertical backward pass (from paper: Implement the backward filter with equation 9b)
// Also write the result back into the byte image
for(unsigned int x = width-1; x < width; x--) {
float prev1[components], prev2[components], prev3[components];
for(unsigned char n = 0; n < components; n++) {
prev1[n] = buffer[IDX(x, height-1, n)];
prev2[n] = prev1[n];
prev3[n] = prev2[n];
}
for(unsigned int y = height-1; y < height; y--) {
for(unsigned char n = 0; n < components; n++) {
float val = B * buffer[IDX(x, y, n)] + (b1 * prev1[n] + b2 * prev2[n] + b3 * prev3[n]) / b0;
image[IDX(x, y, n)] = val;
prev3[n] = prev2[n];
prev2[n] = prev1[n];
prev1[n] = val;
}
}
}
// Free temporary buffers and restore any potential IDX macro
free(buffer);
#pragma pop_macro("IDX")
}
#endif // IIR_GAUSS_BLUR_IMPLEMENTATION
|
