Announce

Puki contents have been moved into SONOTS Plugin (20070703)

Article: The 2-D Gaussain Filter

Note: If rendering of math equations (especially, square root and fraction) look wierd, please try reloading. It would fix problems.

Abstract

The Gaussian filter is a smoothing filter used to blur images to suppress noises. The effect of the Gaussian filter is similar to the average filter in this sense, however, the Gaussian filter is more ideal low-pass filter than the average filter.

In this report, I describe properties or practical issues of the Gaussian filter which we have to care when we implement a Gaussian filter.

 Input Output

The Gaussian Distribution

The 1-D (Univariate) Gaussian

The PDF (Probability Density Function) of the Gaussian has a form:

f(x) = \frac{1}{\sqrt{2\pi} \sigma} \exp \left(-\frac{(x-\mu)^2}{2\sigma ^2} \right)

where \mu is the mean and \sigma is the standard deviation of the distribution.

With \mu = 0,

f(x) = \frac{1}{\sqrt{2\pi} \sigma} \exp \left(-\frac{x^2}{2\sigma ^2} \right).

When we create a Gaussian filter function, we usually let \mu be 0 because we will shift the mean point as the form of the convolution.

The Multivariate Gaussain

The PDF of the multivariate gaussian has a form:

f_X(x_1, \dots, x_N) = f_X(\mathbf{x}) = \frac{1}{(2\pi)^{N/2}|\Sigma|^{1/2}} \exp \left( -\frac{1}{2} ( \mathbf{x} - \mathbf{\mu})^\top \Sigma^{-1} (\mathbf{x} - \mathbf{\mu}) \right)

(if $\Sigma$ is nonsingular) where \mathbf{\mu} is the mean vector and \Sigma is the covariance matrix, and $\ \left| \Sigma \right|$ is the determinant of $\ \Sigma$.

PS. This form is a generalization form of the univariate Gaussian. Try to replace $\Sigma$ with $\sigma^2$ and compare with the univariate form.

The 2-D (Bivariate) Gaussian

An image is a 2-D matrix. The Gaussian filter for image is required to be 2-D.

In the 2-dimensional nonsingular case, that is, when

\Sigma =\pmatrix{ \sigma^2_x & \sigma_{xy} \cr \sigma_{xy} & \sigma^2_y \cr }, \mathbf{\mu} = \pmatrix{ \mu_x \cr \mu_y \cr },

by expanding the Multivariate Gaussian form, the PDF of the Bivariate Gaussian becomes

f(x,y) = \frac{1}{2 \pi \sqrt{\sigma^2_x \sigma^2_y - \sigma^2_{xy}}} \exp \left( -\frac{1}{2 (\sigma^2_x \sigma^2_y - \sigma^2_{xy})} \left( (x-\mu_x)^2 \sigma^2_y - 2(x-\mu_x)(y-\mu_y)\sigma_{xy} + (y-\mu_y)^2 \sigma^2_x \right) \right)

We can change mean simply by shifting (x,y) later, therefore, we let \mu be \mathbf{0} for now, then, the PDF becomes as

f(x,y) = \frac{1}{2 \pi \sqrt{\sigma^2_x \sigma^2_y - \sigma^2_{xy}}} \exp \left( -\frac{1}{2 (\sigma^2_x \sigma^2_y - \sigma^2_{xy})} \left( x^2 \sigma^2_y - 2 x y \sigma_{xy} + y^2 \sigma^2_x \right) \right)

Correlation form

We can rewrite \Sigma as

\Sigma =\pmatrix{ \sigma^2_x & \rho \sigma_x \sigma_y \cr \rho \sigma_x \sigma_y & \sigma^2_y \cr }

so that \sigma_{xy} = \rho \sigma_x \sigma_y where $\rho$ is the correlation between $x$ and $y$. In this case,

f(x,y) = \frac{1}{2 \pi \sigma_x \sigma_y \sqrt{1-\rho^2}} \exp \left( -\frac{1}{2 (1-\rho^2)} \left( \frac{x^2}{\sigma_x^2} + \frac{y^2}{\sigma_y^2} - \frac{2 \rho x y}{ (\sigma_x \sigma_y)} \right) \right)

with mean 0.

(However, I rather use the Multivariate Gaussian form when I implement the Gaussian filter because of its better reusability.)

How It Works

(What I want to argue in this report is the practical issues written in the following section.)

Practical Issues

Window size

The Gaussian distribution is a continuous distribution and its support range is infinite. However, it is not practical to create a Gaussian filter of infinite size in the computer world because the larger size of filter requires the higher computational cost. We need a criterion to choose a good size of a window.

The Gaussian distribution has a good property called as 68-95-99.7 rule.

 Picture from wikipedia

From this property, the window size with 4 \lceil \sigma \rceil + 1 (1 is for the center pixel) to support 95% range (equivalent to set rest 5% to be zero probability) would be a sufficient and reasonable choice as a window size. You may use 5 \lceil \sigma \rceil + 1 to support 99% range.

Window shape

The bivariate Gaussian looks as follows:

You may think that the 2-D Gaussian filter window should have an ellipsoidal shape rather than a rectangular shape.

The ellipsoidal shape window can be generated as following matlab code:

>> F = F .* (F > threshold);

where F is the 2-D Gaussian filter (whose values are probabilities). Although the F is still a matrix (rectangular), values smaller than 'threshold' are suppressed to be 0. Therefore, it forms virtually a ellipsoidal shape. The good threshold should be able to be obtained using the 68-95-99.7 rule such as

>> threshold = F(center_x + 2 * ceil(sigma_x), center_y);

However, does this help? Try the following matlab benchmark code

A = rand(1000,1000); B = zeros(1000, 1000); C = rand(1000, 1000);
tic
for i = 1:100, D = A .* B; end
toc
tic
for i = 1:100, D = A .* C; end
toc

My results were

Elapsed time is 2.331364 seconds.
Elapsed time is 2.307469 seconds.

This tells that changing some values in matrix to 0 does not reduce computation time in matlab, but removes details. This is obvious because the multiplication is done in rectangular form.

To reduce computational cost, we have to re-create the multiplication loop so that the multiplication is not performed for pixels whose values are 0. One way to do it is to add "if" statement to check whether values are 0 or not. But, this is a stupid way that makes rather computationally expensive because we have to check all pixels. Another way to create a smart loop is as

for x = 1:size_x
for y = window_y(x,:)
D(x,y) = A(x,y) * B(x,y)
end
end

where window_y shapes as an ellipsoid. However, creating such window_y is still bothersome.

Let us use a rectangular window shape for simplicity and details.

Normalization

Remember an average filter,

\pmatrix{ 1/9 & 1/9 & 1/9 \cr 1/9 & 1/9 & 1/9 \cr 1/9 & 1/9 & 1/9 \cr }.

It is normalized so that the sum of values in the filter becomes 1. The normalization makes sure that the pixel values do not exceed the maximum of the dynamic range (255) and makes sure that the all same value inputs output the same value.

The sum (integral) of Gaussian distribution becomes 1.0 only when we support infinite window size and when we treat the continuity, but the Gaussian filter is discretized and the window size is limited. Therefore, we have to normalize the Gaussian filter so that the sum becomes 1.0.

This can easily be done by the following matlab code:

F % this is the Gaussain filter matrix
F = F ./ sum(sum(F));

Now, notice that the Gaussain distribution itself has a normalization term

\frac{1}{(2\pi)^{N/2}|\Sigma|^{1/2}}

so that the sum (integral) of Gaussian distribution in infinite range becomes 1.0. This term can be omitted because we normalize the Gaussian filter later.

Note on fspecial('gaussian')

Matlab Image Processing toolbox has a function to generate a Gaussian filter, fspecial('gaussian') to be used as

wsize = [3 5]; % window size [rows cols] or [y x]
sigma = 1;
F = fspecial('guassian', wsize, sigma);
O = conv2(I, F, 'same');

Because sigma is a scalar value, the PDF function of the Bivariate Gaussain uses

\Sigma = \pmatrix{ \sigma^2 & 0 \cr 0 & \sigma^2 \cr }.

Thus, the matlab built-in function does not support to use different \sigmas for x and y coordinates and correlation too although it supports different window size for x and y coordinates.

The window shape is a rectangular.

The fspecial normalizes the filter so that the sum of filter values be 1.0.

Implementation

The matlab codes are provided at this section.

The Multivariate Gaussian distribution

The 2-D convolution (An extension of conv2.m)

The 2-D Gaussian filter

If some files are missing, please look cvpr:..

Demonstration

cvGaussFilter2Demo.m

% cvGaussFilter2Demo - Demo of cvGaussFilter2
function cvGaussFilter2Demo
end