Gaussian filtering with the Image Processing Toolbox

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

4 Responses to “Gaussian filtering with the Image Processing Toolbox”

  1. On November 16th, 2009, at 16:50, A Khan said:

    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

  2. On March 18th, 2011, at 8:52, Madeena Sultana said:

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

  3. On January 15th, 2016, at 23:39, Markus said:

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

  4. On January 18th, 2016, at 17:55, Cris said:

    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 very low 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).

Leave a Reply

You can use these HTML tags: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>

Note: I moderate all comments. Comments without a clear relation to the text above will not be published.