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: usize = (src.width * src.height) as usize;
50 let mut buf: Vec = vec![0.0; buf_size];
51 let buf: &mut Vec = &mut buf;
52
53 let d: BlurData = 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: &mut [u8] = src.data.as_mut_slice();
62 gaussian_channel(data, &d, channel:0, buf);
63 gaussian_channel(data, &d, channel:1, buf);
64 gaussian_channel(data, &d, channel:2, buf);
65 gaussian_channel(data, &d, channel:3, buf);
66}
67
68fn gaussian_channel(data: &mut [u8], d: &BlurData, channel: usize, buf: &mut [f64]) {
69 for i: usize 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: usize 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: f64 = (sigma * sigma) / (2.0 * steps as f64);
145 let dnu: f64 = (1.0 + 2.0 * lambda - (1.0 + 4.0 * lambda).sqrt()) / (2.0 * lambda);
146 (lambda, dnu)
147}
148