resvg/filter/
iir_blur.rs

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
29use super::ImageRefMut;
30use rgb::ComponentSlice;
31
32struct 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.
49pub fn apply(sigma_x: f64, sigma_y: f64, src: ImageRefMut) {
50    let buf_size = (src.width * src.height) as usize;
51    let mut buf = vec![0.0; buf_size];
52    let buf = &mut buf;
53
54    let d = 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 = src.data.as_mut_slice();
63    gaussian_channel(data, &d, 0, buf);
64    gaussian_channel(data, &d, 1, buf);
65    gaussian_channel(data, &d, 2, buf);
66    gaussian_channel(data, &d, 3, buf);
67}
68
69fn gaussian_channel(data: &mut [u8], d: &BlurData, channel: usize, buf: &mut Vec<f64>) {
70    for i 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 in 0..data.len() / 4 {
77        data[i * 4 + channel] = (buf[i] * 255.0) as u8;
78    }
79}
80
81fn 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
144fn gen_coefficients(sigma: f64, steps: usize) -> (f64, f64) {
145    let lambda = (sigma * sigma) / (2.0 * steps as f64);
146    let dnu = (1.0 + 2.0 * lambda - (1.0 + 4.0 * lambda).sqrt()) / (2.0 * lambda);
147    (lambda, dnu)
148}