Gaussian Blur - Standard Deviation, Radius and Kernel Size

Question

I've implemented a gaussian blur fragment shader in GLSL. I understand the main concepts behind all of it: convolution, separation of x and y using linearity, multiple passes to increase radius...

I still have a few questions though:

• What's the relationship between sigma and radius?

  I've read that sigma is equivalent to radius, I don't see how sigma is expressed in pixels. Or is "radius" just a name for sigma, not related to pixels?

• How do I choose sigma?

  Considering I use multiple passes to increase sigma, how do I choose a good sigma to obtain the sigma I want at any given pass? If the resulting sigma is equal to the square root of the sum of the squares of the sigmas and sigma is equivalent to radius, what's an easy way to get any desired radius?

• What's the good size for a kernel, and how does it relate to sigma?

  I've seen most implementations use a 5x5 kernel. This is probably a good choice for a fast implementation with decent quality, but is there another reason to choose another kernel size? How does sigma relate to the kernel size? Should I find the best sigma so that coefficients outside my kernel are negligible and just normalize?

Answer 1

What's the relationship between sigma and radius? I've read that sigma is equivalent to radius, I don't see how sigma is expressed in pixels. Or is "radius" just a name for sigma, not related to pixels?

There are three things at play here. The variance, ($\sigma^2$), the radius, and the number of pixels. Since this is a 2-dimensional gaussian function, it makes sense to talk of the covariance matrix $\boldsymbol{\Sigma}$ instead. Be that as it may however, those three concepts are weakly related.

First of all, the 2-D gaussian is given by the equation:

$$ g({\bf z}) = \frac{1}{\sqrt{(2 \pi)^2 |\boldsymbol{\Sigma}|}} e^{-\frac{1}{2} ({\bf z}-\boldsymbol{\mu})^T \boldsymbol{\Sigma}^{-1} \ ({\bf z}-\boldsymbol{\mu})} $$

Where ${\bf z}$ is a column vector containing the $x$ and $y$ coordinate in your image. So, ${\bf z} = \begin{bmatrix} x \\ y\end{bmatrix}$, and $\boldsymbol{\mu}$ is a column vector codifying the mean of your gaussian function, in the $x$ and $y$ directions $\boldsymbol{\mu} = \begin{bmatrix} \mu_x \\ \mu_y\end{bmatrix}$.

Example:

Now, let us say that we set the covariance matrix $\boldsymbol{\Sigma} = \begin{bmatrix} 1 & 0 \\ 0 & 1\end{bmatrix}$, and $\boldsymbol{\mu} = \begin{bmatrix} 0 \\ 0\end{bmatrix}$. I will also set the number of pixels to be $100$ x $100$. Furthermore, my 'grid', where I evaluate this PDF, is going to be going from $-10$ to $10$, in both $x$ and $y$. This means I have a grid resolution of $\frac{10 - (-10)}{100} = 0.2$. But this is completely arbitrary. With those settings, I will get the probability density function image on the left. Now, if I change the 'variance', (really, the covariance), such that $\boldsymbol{\Sigma} = \begin{bmatrix} 9 & 0 \\ 0 & 9\end{bmatrix}$ and keep everything else the same, I get the image on the right.

The number of pixels are still the same for both, $100$ x $100$, but we changed the variance. Suppose instead we do the same experiment, but use $20$ x $20$ pixels instead, but I still ran from $-10$ to $10$. Then, my grid has a resolution of $\frac{10-(-10)}{20} = 1$. If I use the same covariances as before, I get this:

These are how you must understand the interplay between those variables. If you would like the code, I can post that here as well.

How do I choose sigma?

The choice of the variance/covariance-matrix of your gaussian filter is extremely application dependent. There is no 'right' answer. That is like asking what bandwidth should one choose for a filter. Again, it depends on your application. Typically, you want to choose a gaussian filter such that you are nulling out a considerable amount of high frequency components in your image. One thing you can do to get a good measure, is compute the 2D DFT of your image, and overlay its co-efficients with your 2D gaussian image. This will tell you what co-efficients are being heavily penalized.

For example, if your gaussian image has a covariance so wide that it is encompassing many high frequency coefficients of your image, then you need to make its covariance elements smaller.

Answer 2

The parameter sigma is enough to define the Gaussian blur from a continuous point of view. In practice however, images and convolution kernels are discrete. How to choose an optimal discrete approximation of the continuous Gaussian kernel?

The discrete approximation will be closer to the continuous Gaussian kernel when using a larger radius. But this may come at the cost of added computation duration.

Ideally, one would select a value for sigma, then compute a radius that allows to represent faithfully the corresponding continuous Gaussian kernel. For a given approximation error, the bigger sigma is, the larger the radius must be.

Interestingly, this can become very complicated to get it right. When constructing the Gaussian matrix, is the best solution to sample the continuous kernel or are there better approximations? How to normalize the computed discrete kernel to account for truncation? etc.

As a reference, in Mathematica the function GaussianMatrix features several ways to compute a Gaussian discrete matrix, e.g. using discrete Bessel approximation. By default, radius = 2 * sigma, which means that with sigma = 1, the matrix will be 5x5.

Answer 3

It turns out that the rows of Pascal's Triangle approximate a Gaussian quite nicely and have the practical advantage of having integer values whose sum is a power of 2 (we can store these values exactly as integers, fixed point values, or floats). For example, say we wish to construct a 7x7 Gaussian Kernel we can do so using the 7th row of Pascal's triangle as follows:

Note that this filter has the minimum influence at the corners while remaining integer valued. You can use the middle value 20/64 to determine the corresponding standard deviation sigma which is 64/(20 * sqrt(2*pi)) = 1.276 for the approximated Gaussian in this case. You can graph the Gaussian to see this is an excellent fit.

So a good starting point for determining a reasonable standard deviation for a Gaussian Kernel comes from Pascal's Triangle (aka Binomial Coefficients) -- for a (N+1)x(N+1) filter corresponding to the above construction use

Wolfram Alpha's GaussianMatrix[3] just uses r/2 = 1.5. Strangely enough, GaussianMatrix[{3,1.276}] does not yield the same 2D filter as mine and is not the following for x,y between -3 and 3:

I am not sure why not? My 2D filter is an excellent fit.

Answer 4

All of the answers above contained valuable insight, but didn't succinctly answer the practical question: how do I choose sigma? Often the practical question is posed in the form that a trade-off needs to be found between quality of results (Low-pass filtering, noise suppression) and implementation complexity. Suppose a particular HW architecture supports 2D convolution kernel of a particular size - how do I best initialize the kernel for optimal results? @Tarin is correct, this is a complicated subject, the kernel could be optimized once we define "quality of results" in terms of the properties of the resulting filter. But radius = 2 * sigma (@Matthias, above) results to kernels where the maximum (central) coefficient is about 54x the magnitude of the minimum (corner) coefficient, which for 8 bit pixels mean that the long tails of the Gaussians that we quantize to zero by neglecting to include them in our 2D convolution will result to < 4 LSB errors. In other words:

$\sigma = \frac{s-1}{4}$, where s is the kernel size,

is a good practical trade-off between minimizing quantization issues and lowering the cut-off frequency of the resulting low-pass filter.

Also, for practical purposes, define the kernel just with the x,y meshgrid over $e^{-\frac{x^2+y^2}{2\sigma^2}}$ (without the normalization factor), then calculate the 2D sum and normalize the kernel. This prevents image darkening due to the finite discrete kernel size.