1 | // This Source Code Form is subject to the terms of the Mozilla Public |
2 | // License, v. 2.0. If a copy of the MPL was not distributed with this |
3 | // file, You can obtain one at http://mozilla.org/MPL/2.0/. |
4 | |
5 | // An IIR blur. |
6 | // |
7 | // Based on http://www.getreuer.info/home/gaussianiir |
8 | // |
9 | // Licensed under 'Simplified BSD License'. |
10 | // |
11 | // |
12 | // Implements the fast Gaussian convolution algorithm of Alvarez and Mazorra, |
13 | // where the Gaussian is approximated by a cascade of first-order infinite |
14 | // impulsive response (IIR) filters. Boundaries are handled with half-sample |
15 | // symmetric extension. |
16 | // |
17 | // Gaussian convolution is approached as approximating the heat equation and |
18 | // each timestep is performed with an efficient recursive computation. Using |
19 | // more steps yields a more accurate approximation of the Gaussian. A |
20 | // reasonable default value for `numsteps` is 4. |
21 | // |
22 | // Reference: |
23 | // Alvarez, Mazorra, "Signal and Image Restoration using Shock Filters and |
24 | // Anisotropic Diffusion," SIAM J. on Numerical Analysis, vol. 31, no. 2, |
25 | // pp. 590-605, 1994. |
26 | |
27 | // TODO: Blurs right and bottom sides twice for some reason. |
28 | |
29 | use super::ImageRefMut; |
30 | use rgb::ComponentSlice; |
31 | |
32 | struct BlurData { |
33 | width: usize, |
34 | height: usize, |
35 | sigma_x: f64, |
36 | sigma_y: f64, |
37 | steps: usize, |
38 | } |
39 | |
40 | /// Applies an IIR blur. |
41 | /// |
42 | /// Input image pixels should have a **premultiplied alpha**. |
43 | /// |
44 | /// A negative or zero `sigma_x`/`sigma_y` will disable the blur along that axis. |
45 | /// |
46 | /// # Allocations |
47 | /// |
48 | /// This method will allocate a 2x `src` buffer. |
49 | pub fn apply(sigma_x: f64, sigma_y: f64, src: ImageRefMut) { |
50 | let buf_size: usize = (src.width * src.height) as usize; |
51 | let mut buf: Vec = vec![0.0; buf_size]; |
52 | let buf: &mut Vec = &mut buf; |
53 | |
54 | let d: BlurData = BlurData { |
55 | width: src.width as usize, |
56 | height: src.height as usize, |
57 | sigma_x, |
58 | sigma_y, |
59 | steps: 4, |
60 | }; |
61 | |
62 | let data: &mut [u8] = src.data.as_mut_slice(); |
63 | gaussian_channel(data, &d, channel:0, buf); |
64 | gaussian_channel(data, &d, channel:1, buf); |
65 | gaussian_channel(data, &d, channel:2, buf); |
66 | gaussian_channel(data, &d, channel:3, buf); |
67 | } |
68 | |
69 | fn gaussian_channel(data: &mut [u8], d: &BlurData, channel: usize, buf: &mut Vec<f64>) { |
70 | for i: usize in 0..data.len() / 4 { |
71 | buf[i] = data[i * 4 + channel] as f64 / 255.0; |
72 | } |
73 | |
74 | gaussianiir2d(d, buf); |
75 | |
76 | for i: usize in 0..data.len() / 4 { |
77 | data[i * 4 + channel] = (buf[i] * 255.0) as u8; |
78 | } |
79 | } |
80 | |
81 | fn gaussianiir2d(d: &BlurData, buf: &mut Vec<f64>) { |
82 | // Filter horizontally along each row. |
83 | let (lambda_x, dnu_x) = if d.sigma_x > 0.0 { |
84 | let (lambda, dnu) = gen_coefficients(d.sigma_x, d.steps); |
85 | |
86 | for y in 0..d.height { |
87 | for _ in 0..d.steps { |
88 | let idx = d.width * y; |
89 | |
90 | // Filter rightwards. |
91 | for x in 1..d.width { |
92 | buf[idx + x] += dnu * buf[idx + x - 1]; |
93 | } |
94 | |
95 | let mut x = d.width - 1; |
96 | |
97 | // Filter leftwards. |
98 | while x > 0 { |
99 | buf[idx + x - 1] += dnu * buf[idx + x]; |
100 | x -= 1; |
101 | } |
102 | } |
103 | } |
104 | |
105 | (lambda, dnu) |
106 | } else { |
107 | (1.0, 1.0) |
108 | }; |
109 | |
110 | // Filter vertically along each column. |
111 | let (lambda_y, dnu_y) = if d.sigma_y > 0.0 { |
112 | let (lambda, dnu) = gen_coefficients(d.sigma_y, d.steps); |
113 | for x in 0..d.width { |
114 | for _ in 0..d.steps { |
115 | let idx = x; |
116 | |
117 | // Filter downwards. |
118 | let mut y = d.width; |
119 | while y < buf.len() { |
120 | buf[idx + y] += dnu * buf[idx + y - d.width]; |
121 | y += d.width; |
122 | } |
123 | |
124 | y = buf.len() - d.width; |
125 | |
126 | // Filter upwards. |
127 | while y > 0 { |
128 | buf[idx + y - d.width] += dnu * buf[idx + y]; |
129 | y -= d.width; |
130 | } |
131 | } |
132 | } |
133 | |
134 | (lambda, dnu) |
135 | } else { |
136 | (1.0, 1.0) |
137 | }; |
138 | |
139 | let post_scale = |
140 | ((dnu_x * dnu_y).sqrt() / (lambda_x * lambda_y).sqrt()).powi(2 * d.steps as i32); |
141 | buf.iter_mut().for_each(|v| *v *= post_scale); |
142 | } |
143 | |
144 | fn gen_coefficients(sigma: f64, steps: usize) -> (f64, f64) { |
145 | let lambda: f64 = (sigma * sigma) / (2.0 * steps as f64); |
146 | let dnu: f64 = (1.0 + 2.0 * lambda - (1.0 + 4.0 * lambda).sqrt()) / (2.0 * lambda); |
147 | (lambda, dnu) |
148 | } |
149 | |