Gaussian filtering with the Image Processing Toolbox
Edit May 2018
Since publishing this post, the MATLAB Image Processing Toolbox has added the
function imgaussfilt that correctly applies a
Gaussian smoothing filter. For Gaussian derivatives, the recommendations here still apply.
Edit July 2024
I just learned that imgaussfilt uses a truncation of 2 by default (i.e. it computes cutoff = ceil(2*sigma)).
This is not sufficiently large for most applications, see this other post.
There’s a 'FilterSize' name-value argument that should be used with this function to set the correct filter
size, computed as shown in this post.
If you don’t use DIPimage, you probably use
MATLAB’s Image Processing Toolbox. This toolbox makes it really easy to do
convolutions with a Gaussian in the wrong way. On three accounts. The
function fspecial is used to create a convolution kernel for
a Gaussian filter. This kernel is 2D. That’s the first problem. The other two problems are given by the default values
of its parameters. The default value for the kernel size is [3 3]. The default value for the σ (sigma) is 0.5. σ=0.5
is too small for a Gaussian kernel. When a kernel this small is sampled, too much aliasing occurs. Depending on who you
listen to, the smallest σ you can use is either 0.8 or 1.0. The value of σ dictates the shape of the kernel, but it also
determines its size: truncating the Gaussian kernel at [3 3] pixels creates a sharp border which will cause unwanted
artifacts in the output image. All the good properties of the Gaussian filter will be lost. You have to cut off the
Gaussian at a point where it is close to zero. This happens around 3σ. This means that the size of the Gaussian kernel
should be 6σ (cut off at 3σ on ether size of the origin). Or if you need really high numerical accuracy, you can make
this larger, say 10σ. These last two issues are easy to take care of:
img = single(imread('cameraman.tif'));
sigma = 3;
cutoff = ceil(3*sigma);
h = fspecial('gaussian',2*cutoff+1,sigma);
out = conv2(img,h,'same');
However, this does not solve the first issue: the filter h is two-dimensional, but the Gaussian is separable! That
it’s separable means that we can save a lot of computation by convolving with a one-dimensional Gaussian along the
x-axis and then with another along the y-axis. You can accomplish this with:
h = fspecial('gaussian',[1,2*cutoff+1],sigma); % 1D filter
out = conv2(h,h,img,'same');
The separable Gaussian filter is always faster to compute when implemented as two one-dimensional filters than when implemented as one two-dimensional filter. Taking σ=1, this means we do 7+7=14 multiplications per pixel instead of 7×7=49 multiplications per pixel. The savings are increasingly important for larger kernels. With σ=5 this translates to 62 vs 961 multiplications.
Finally, a Gaussian derivative should not be implemented like so:
h = fspecial('gaussian',[1,2*cutoff+1],sigma);
out = conv2(h,h,img,'same');
out = diff(out,1);
Nor like this:
h = fspecial('gaussian',[1,2*cutoff+1],sigma);
dh = diff(h);
out = conv2(dh,h,img,'same');
Computing the finite difference approximation to the derivative of a smoothed image is the same as convolving with the finite difference approximation to the derivative of the Gaussian. As discussed earlier in this blog, the advantage of computing Gaussian derivatives is lost when implemented this way. Instead, compute the exact derivative of the Gaussian, and convolve with that:
h = fspecial('gaussian',[1,2*cutoff+1],sigma);
dh = h .* (-cutoff:cutoff) / (-sigma^2);
out = conv2(dh,h,img,'same');
Also note that, because the derivative of the Gaussian goes to zero slightly slower than the Gaussian itself, it is
worthwhile to increase the value of cutoff to at least ceil(3.5*sigma).