## 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.

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 worth while to increase the value of `cutoff`

to at least `ceil(3.5*sigma)`

.

Hi,

You’ve made some useful points in this article. It might also be useful to know that two convolutions with Gaussian filters applied sequentially (or “in cascade” as they say) is equivalent to a single Gaussian convolution, with variances related as follows:

sigma_eq^2 = sigma_1^2 + sigma_2^2

It helped me a lot, thank you very much Cris.

Can you elaborate on sigma 0.5 and “When a kernel this small is sampled, too much aliasing occurs”?

Markus,

A Gaussian kernel with sigma=0.5 pixels has too much power in the frequencies above the Nyquist frequency. These frequencies are aliased when sampling. That is, their power is assigned to frequencies under the Nyquist frequency, causing low-frequency signals to appear in the sampled signal that were never present in the original signal.

The Nyquist frequency is simply the highest frequency that can be represented by the sampled signal/image, and depends with a simple fraction on the sampling frequency. Thus, the closer together the pixels are, the more detail can be represented in the image.

A Gaussian kernel has an infinite number of frequencies, but frequencies above a certain value have

verylow power, and can be ignored. Depending on where exactly you place this threshold (how large an error you are willing to accept), you need a sampling frequency of about 1*sigma to 1/0.8*sigma. This translates to, when picking a fixed sampling grid, a minimum sigma of 0.8 to 1 the sampling spacing.When you sample a Gaussian with a sigma of 0.5 pixels, the resulting sampled sigma does not have a Guassian shape, and does not present some of the beneficial properties of the Gaussian.

If you want to convolve with such a small Guassian, you should do this through the Fourier domain, or following the recipe of Rein van den Boomgaard (presented at Scale-Space 2001).

Should not the last line of code be:

instead of:

?

Johannes,

No, it is

`(dh,h)`

for a derivative along the x-axis, and`(h,dh)`

for a derivative along the y-axis.Thanks Cris!

Ok, so the goal was to find the derivative in only one direction.

So I then interpret the expression conv2(dh,h,img) as follows:

the derivative (of the Gaussian) along x and the pure Gaussian along y.

Consequently conv2(dh,h,img) yields:

the pure Gaussian along x and the derivative (of the Gaussian) along y.

Would it then be correct to say that conv2(dh,dh,img) yields:

the derivative along x plus the derivative along y, i.e. the magnitude of the gradient?

Johannes,

The combination

`(dh,dh)`

would be a 2nd order derivative, corresponding to ∂²/∂x∂y.To get the magnitude of the gradient, you need to add the square of the two partial derivatives, which you cannot do with a single convolution (the convolution is always a linear operation, and squaring is a non-linear operation).

Thank you Cris, that makes sens, now I finally understood the operation!

Just a note: When playing with Matlab with 3D images (CT-scans), I noticed that when defining a filter as above by

h = fspecial(‘gaussian’,[1,2*cutoff+1],sigma)

it is generally faster and also more flexible to use

J = imfilter(I,h,’replicate’,’conv’);

J = imfilter(J,permute(h,[2,3,1]),’replicate’,’conv’);

J = imfilter(J,permute(h,[3,1,2]),’replicate’,’conv’);

instead of applying the convn() function in all three dimensions.

The problem with convn(), and respectively conv2(), is that the replication pattern at the boarders cannot be defined (it is automatically padded with 0). Matlab’s internal function imgaussfilt3() seems to peroform almost equally fast as the code above, I guess it separates the kernel too.

Cheers!