How to approximate gaussian kernel for image blur

Question

From wiki, a $3 \times3$ gaussian kernel is approximated as: $$\frac{1}{16}\begin{bmatrix}1&2&1\\2&4&2\\1&2&1 \end{bmatrix}.$$ Applying this kernel in an image equals to applying an one-dimensional kernel in x-direction then again in y-direction, so one-dimensional kernel is $$\frac{1}{4}\begin{bmatrix}1&2&1\end{bmatrix}.\qquad \text{for 3*3 kernel}$$ $$\frac{1}{16}\begin{bmatrix}1&4&6&4&1\end{bmatrix}.\qquad \text{for 5*5 kernel}$$ My question is how to derive this approximation?

Answer 1

The continuous Gaussian, whatever its dimension (1D, 2D), is a very important function in signal and image processing. As most data is discrete, and filtering can be costly, it has been and still is, subject of quantities of optimization and quantification/quantization schemes. In one 1D, the three most direct for a finite-length filter are illustrated below:

In solid blue, the continuous Gaussian. For an $L=7$-length, I choose a Gaussian with parameter $s=\sqrt{L/2}$, in consistence with the first approximation $(P)$ given by the binomial formula or the Pascal triangle, as detailed by @Olli Niemitalo.

Each filter coefficient (red crosses) is of the form $$h_P[k] = \frac{\binom{L}{k}}{2^L}\,.$$ Their sum is naturally equal to unity, and the coefficients are dyadic-rationals, in the form of $a/2^b$. They are practical for integer operations, on integer-valued pixels. As you can see, they are slightly above the Gaussian around the center, and below on the tails. When $L\to \infty, the shape tends to a certain bell function shape, see From Pascal's Triangle to the Bell-shaped Curve.

A second approximation (green circles), often called the truncated Gaussian, consists in sampling ($S$) $x\mapsto G(x)=\frac{1}{s\sqrt{\pi}}e^{-\frac{x^2}{s^2}}$ at integer values of $x$: $[-3,-2,\ldots,0,\ldots,3]$. They need to be normalized to 1-unit sum:

$$h_S[k] = \frac{G(k-(L-1)/2)}{\sum_{l=1}^L G(l-(L-1)/2)}\,.$$

The third approximation ($A$), blue circles and bars) is area dependent: once again after normalization, the coefficients are proportional to the area under the Gaussian, on the interval $[k-1/2,k+1/2]$:

$$h_A[k] \propto \int_{k-1/2}^{k+1/2} G(x)dx\,.$$

The exact value depends on the quadrature formulae.

In our settings (with the chosen $s$), coefficients are, respectively from left to right (pascal, Sampling, Area):

$$ \begin{Bmatrix} 0.0156 & 0.0232 & 0.0255\\ 0.0938 & 0.0968 & 0.0998\\ 0.2344 & 0.2282 & 0.2262\\ 0.3125 & 0.3036 & 0.2970\\ 0.2344 & 0.2282 & 0.2262\\ 0.0938 & 0.0968 & 0.0998\\ 0.0156 & 0.0232 & 0.0255\\ \end{Bmatrix} $$

Only the first column admits dyadic/integer operations. To work with integer arithmetic, one can round to the nearest integer. As $2^7=128$, this quantization choice (another approximation) leads to:

$$\frac{1}{128} \begin{Bmatrix} 2 & 3 & 3\\ 12 & 12 & 13\\ 30 & 29 & 29\\ 40 & 39 & 38\\ 30 & 29 & 29\\ 12 & 12 & 13\\ 2 & 3 & 3\\ \end{Bmatrix} $$

With $2^6=64$, you get:
$$\frac{1}{64} \begin{Bmatrix} 1 & 1 & 2\\ 6 & 6 & 6\\ 15 & 15 & 14\\ 20 & 19 & 19\\ 15 & 15 & 14\\ 6 & 6 & 6\\ 1 & 1 & 2\\ \end{Bmatrix} $$

In the case you mention in comments, the formula for $\sigma$ is a bit different: $\sigma=(3L+1)/20$, equal to $1.4$ for $L=7$. I don't know its origin yet. As said above, those are the main principles. In the Gaussian Kernel Calculator demo, with those values, you get:

$$\begin{Bmatrix} 0.031251\\ 0.106235 \\0.221252\\ 0.282524 \\0.221252\\ 0.106235\\ 0.031251 \end{Bmatrix} $$ which is quite close (up to the integral quadrature precision) to your: $$\begin{Bmatrix} 0.03125\\0.109375\\0.21875\\0.28125\\0.21875\\0.109375\\0.03125 \end{Bmatrix} $$

as the former, multiplied by 64, gives:

$$\begin{Bmatrix} 2.0001 \\ 6.7990 \\ 14.1601 \\ 18.0815 \\ 14.1601 \\ 6.7990 \\ 2.0001 \end{Bmatrix} $$

The ($A$) approximation is possibly more accurate with respect to the integral version of the convolution, with a continuous Gaussian and a discrete or piecewise constant signal.

In 2D, one can tensorize two 1D filters, in which case the 2D filter has rank one. Or one extends the above method to 2D. The 2D coefficients can be obtained by two-dimensional integral estimations:

$$h_A[m,n] \propto \iint_{[m-1/2,m+1/2]\times [n-1/2,n+1/2]} G(x,y)dxdy\,.$$

See for instance Gaussian smoothing, esp. Figure 3, Discrete approximation to Gaussian function with $\sigma=1.0$:

Answer 2

Some possible explanations for the coefficients:

Binomial coefficients

The 1-d kernels are probability mass functions of binomial distributions with probability parameter $p=1/2$ to make them symmetrical. Binomial distributions can be approximated by Gaussian distributions, so it should be true that Gaussian distributions can also be approximated by binomial distributions. You can obtain binomial distributions with $p = 1/2$ by convolving the length 2 kernel:

$$\frac{1}{2}\begin{bmatrix}1&1\end{bmatrix}$$

by the length 2 kernel multiple times. In Octave, without the normalization factors for clarity:

>> f = [1 1];
>> g = f #1/2
g =

   1   1

>> g = conv(g, f) #1/4
g =

   1   2   1

>> g = conv(g, f) #1/8
g =

   1   3   3   1

>> g = conv(g, f) #1/16
g =

   1   4   6   4   1

>> g = conv(g, f) #1/32
g =

    1    5   10   10    5    1

>> g = conv(g, f) #1/64
g =

    1    6   15   20   15    6    1

For large kernels the first and last binomial coefficients become very small and have hardly any effect on the result.

OpenCV

In comments to this answer, you also mention the length 7 kernel: $$\begin{bmatrix}0.03125&0.109375&0.21875&0.28125&0.21875&0.109375&0.03125\end{bmatrix}\\ = \frac{1}{64}\begin{bmatrix}2&7&14&18&14&7&2\end{bmatrix}.$$

This can be found in the source code of OpenCV function cv::getGaussianKernel where the small kernels up to that kernel are hard-coded:

static const float small_gaussian_tab[][SMALL_GAUSSIAN_SIZE] =
    {
        {1.f},
        {0.25f, 0.5f, 0.25f},
        {0.0625f, 0.25f, 0.375f, 0.25f, 0.0625f},
        {0.03125f, 0.109375f, 0.21875f, 0.28125f, 0.21875f, 0.109375f, 0.03125f}
    };

These also match the kernels cited in the question. The function is also documented explaining the general calculation, but not the rounding in the small kernels. For sizes $3,$ $5,$ $7, $ this appears to be to the nearest multiple of $\frac{1}{4},$ $\frac{1}{16},$ $\frac{1}{64},$ respectively, using the general calculation and then rounding. Based on the documentation, the general calculation goes as:

$$\begin{array}{ll}N&\quad\text{filter length}\\ \sigma = 0.3\times((N - 1)\times0.5 - 1) + 0.8&\quad\text{standard deviation}\\ \exp(- (i - (N - 1)/2)^2/(2\sigma^2))&\quad\text{unnormalized coefficients, with }0\le i<N\end{array}$$

They do not document how that was derived. Seems ad hoc / empirical. After the above calculation, the coefficients are normalized so that their sum is 1.