resvg/filter/
iir_blur.rs

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