Update: this project has been revised and published as an article on IPOL and surveys several methods in addition to Alvarez–Mazorra: “A Survey of Gaussian Convolution Algorithms”.
## OverviewThis C/C++/MATLAB source code implements approximate Gaussian convolution in 1D, 2D, and 3D using the efficient recursive filtering algorithm of Alvarez and Mazorra, Alvarez, Mazorra, “Signal and Image Restoration using Shock Filters and Anisotropic Diffusion,” The code is implemented in ANSI C for use in C and C++ programs, and MEX wrappers are included for use in MATLAB. In case MEX is unavailable, (slower) MATLAB M-code implementations are also included. Gaussian convolution is attractive because it is rotationally invariant, satisfies a semi-group property, and it does not create new extrema or enhance existing extrema. Furthermore, the Gaussian decays rapidly in both space and frequency. Thus Gaussian convolution is very popular in signal and image processing. The algorithm of Alvarez and Mazorra provides an efficient approximation of Gaussian convolution. It uses a cascade of first order recursive filters, also called infinite impulse response or IIR filters, to simulate time steps of the heat equation. The speed vs. accuracy of the approximation is controlled by selecting the number of time steps. Useful aspects of this algorithm are - computational complexity is independent of the Gaussian standard deviation
- the computation is in-place, and memory overhead is small and fixed
- converges to exact solution in the limit number of time steps → ∞
Caveats: - the method is approximate: it is not exactly invariant to rotations or axis reflections (but accuracy can be improved by increasing the number of time steps)
- this method is less useful for small sigma, where convolution with a truncated Gaussian (FIR approximation) is more accurate/efficient
- for high accuracy, consider instead Fourier-based Gaussian convolution
The algorithm of Alvarez and Mazorra provides an efficient approximation of Gaussian convolution. A useful property of the algorithm is that the computational complexity is independent of the Gaussian standard deviation. It uses a cascade of first order recursive filters, also called infinite impulse response or IIR filters, to simulate time steps of the heat equation. The speed vs. accuracy of the approximation is controlled by selecting the number of time steps. ## LicenseFiles gaussianiir1d.c, gaussianiir2d.c, and gaussianiir3d.c, and the files of the same name but with extension .h and .m may be used with either the GPLv3 or simplified BSD license. You can redistribute it and/or modify it under, at your option, the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version, or the terms of the simplified BSD license. ## DemoIncluded is a simple program blurdemo to apply 2D Gaussian convolution to color images. Usage: BMP images are supported, and optionally blurdemo can be compiled to use the libjpeg, libpng, and libtiff libraries to support JPEG, PNG, and TIFF images. Options:
Example:
The blurdemo program can be compiled with GCC using makefile.gcc as
or with MSVC using makefile.vc by opening a Visual Studio Command Prompt (under Start Menu → Microsoft Visual Studio → Visual Studio Tools) and entering
## Usage in C/C++The convolution itself is in gaussianiir1d.c, gaussianiir2d.c, and gaussianiir3d.c, implementing respectively 1D, 2D, and 3D. There is no dependency among these files, so you may select the one that matches the dimension of your application and ignore the other two. The calling syntax of the convolution functions are void gaussianiir1d(float *data, long length, float sigma, int numsteps); void gaussianiir2d(float *image, long width, long height, float sigma, int numsteps); void gaussianiir3d(float *volume, long width, long height, long depth, float sigma, int numsteps); In all three cases, the first argument is the input data, a contiguous array
which is modified in-place and is the convolution output upon return.
Boundaries are handled with half-sample symmetric extension. The standard
deviation is specified by In gaussianiir1d, In gaussianiir2d, image is expected to be in row-major order,
In gaussianiir3d, volume is expected to be ordered such that
## Usage in MATLABMATLAB MEX wrappers are included in gaussianiir1d.c, gaussianiir2d.c, and gaussianiir3d.c so that they can be compiled as fast MEX functions, usable from MATLAB. To compile the MEX functions, enter >> mex gaussianiir1d.c >> mex gaussianiir2d.c >> mex gaussianiir3d.c on the MATLAB console. This should compile the MEX functions gaussianiir1d, gaussianiir2d, and gaussianiir3d. Alternatively, native M-code implementations of gaussianiir1d, gaussianiir2d,
and gaussianiir3d are included, but beware that they are significantly slower
than the MEX versions. To determine whether you are using the MEX version or
the M-code version, use MATLAB's >> which gaussianiir1d /home/pascal/projects/gaussianiir/gaussianiir1d.mexa64 If the output ends in .m, you are using the M-code version. Otherwise, you are using the compiled MEX function (the exact extension depends on platform, e.g., .mexw32 for 32-bit Windows). The usage of gaussianiir1d is
where
x is a matrix or multidimensional array, the convolution is performed
independently on each column as if you had done
for j = 1:size(x,2) y(:,j) = gaussianiir1d(x(:,j), sigma, numsteps) end The functions gaussianiir2d and gaussianiir3d are used in the same way, except
that gaussianiir2d performs 2D convolution along the first two dimensions of |

Research >