How to Make a Gaussian Filter?

The 2D Gaussian function is defined as: $$G(x,y) = \frac{1}{2\pi \alpha}e^{-\frac{x^{2}+y^{2}}{2\alpha}}$$

And this is a Gaussian kernel:

If I take the top left corner as the origin and set $$\alpha=1$$, then at $$x$$=4 and $$y$$=0, $$G(4,0)$$ = $$5.3\times10^{-3}$$. But here at $$G(4,0) = 3.66\times10^{-3}$$. Then how they design this Gaussian kernel? I am new on this topic. So, this question may seem too much beginner level.

Cool facts about the Gaussian surface:

It is a rotation:

$$G(x,y) = \frac{1}{2\pi \alpha}e^{-\frac{x^{2}+y^{2}}{2\alpha}} = \frac{1}{2\pi \alpha}e^{-\frac{r^{2}}{2\alpha}} = G(r)$$

where $$r = \sqrt{x^2 + y^2}$$

It is separable:

$$G(x,y) = \frac{1}{2\pi \alpha}e^{-\frac{x^{2}+y^{2}}{2\alpha}} = \frac{1}{2\pi \alpha} e^{-\frac{x^2}{2\alpha} } e^{-\frac{y^2}{2\alpha} } = G_x(x)G_y(y)$$

With proper normalizations

$$G_x(t) = G_y(t) = G_t(t) = \frac{1}{\sqrt{2\pi \alpha}}e^{-\frac{t^{2}}{2\alpha}}$$

This means every slice of a Guassian surface is a Guassian function.

I used to do a lot of smoothing on scatter dot diagrams to make them nice surfaces.

The first is the same as DC. You want the sum of your entries to equal the denominator you are using. Yours is fudged a little bit from the nearest solution found here.

Here is the closest match for a denominatory of 271. You can see where they tweaked to get 273.

[[  1.   4.   6.   4.   1.]
[  4.  16.  26.  16.   4.]
[  6.  26.  43.  26.   6.]
[  4.  16.  26.  16.   4.]
[  1.   4.   6.   4.   1.]]

[[  0.78997163   3.54040722   5.83714469   3.54040722   0.78997163]
[  3.54040722  15.86700434  26.16026756  15.86700434   3.54040722]
[  5.83714469  26.16026756  43.13098958  26.16026756   5.83714469]
[  3.54040722  15.86700434  26.16026756  15.86700434   3.54040722]
[  0.78997163   3.54040722   5.83714469   3.54040722   0.78997163]]



The bottom one is the actual surface rescaled by the denominator.

I am wondering like Royi in the comment if this is supposed to be an exercise in integer approximation. So, I got curious and without further comment on the issues that are raised by this program, here it is:

import numpy as np

#==========================================================
def main():

N = 5

alpha = 1.0

theSurface = np.zeros( ( N, N ) )

L = ( N - 1 ) / 2

C = 1.0 / ( 2.0 * np.pi * alpha )

for x in range( N ):
dx  = x - L
dx2 = dx * dx
for y in range( N ):
dy  = y - L
dy2 = dy * dy
r2  = dx2 + dy2

s = C * np.exp( -r2 / ( 2.0 * alpha ) )

theSurface[x,y] = s

print( theSurface )

theIntegerMatches = []

for theDenominator in range( 100, 10000 ):
theApprox = np.round( theDenominator * theSurface )

theSum = np.sum( theApprox )
if abs( theSum - theDenominator ) < 0.5:
theIntegerMatches.append( ( theDenominator, theApprox ) )

for theMatch in theIntegerMatches:
theDenominator = theMatch[0]
theApprox      = theMatch[1]

theModel = theApprox / theDenominator

theDiff  = theSurface - theModel
theDiff2 = np.multiply( theDiff, theDiff )

theDiffSum  = np.sum( theDiff )
theDiff2Sum = np.sum( theDiff2 )

print( theDenominator, theApprox[0][0], theDiffSum, theDiff2Sum, 1.0 / theDiffSum )

print

theMatch = theIntegerMatches[-1]
theDenominator = theMatch[0]
theApprox      = theMatch[1]

theActual = theDenominator * theSurface

print( theApprox )
print( theActual )

#==========================================================
main()


Any formula based on r^2 is going to be a rotation around an axis.

That being said. In integer implementations it is customary to choose a power of two denominator so that the division can be performed with a simple bit shift.

Keep this in mind when dealing with any filter like this: Your filter is going to work on a surface that can be approximated by the following Taylor series representation. At any point, there is a best fit set of coefficients. In this case, your filter zeros out several of these.

\begin{aligned} f &= C_1 \begin{bmatrix} 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 \\ \end{bmatrix} + C_{x} \begin{bmatrix} -2 & 1 & 0 & 1 & 2 \\ -2 & 1 & 0 & 1 & 2 \\ -2 & 1 & 0 & 1 & 2 \\ -2 & 1 & 0 & 1 & 2 \\ -2 & 1 & 0 & 1 & 2 \\ \end{bmatrix} + C_{y} \begin{bmatrix} -2 & -2 & -2 & -2 & -2 \\ -1 & -1 & -1 & -1 & -1 \\ 0 & 0 & 0 & 0 & 0 \\ 1 & 1 & 1 & 1 & 1 \\ 2 & 2 & 2 & 2 & 2 \\ \end{bmatrix} \\ &+ C_{xx} \begin{bmatrix} 4 & 1 & 0 & 1 & 4 \\ 4 & 1 & 0 & 1 & 4 \\ 4 & 1 & 0 & 1 & 4 \\ 4 & 1 & 0 & 1 & 4 \\ 4 & 1 & 0 & 1 & 4 \\ \end{bmatrix} + C_{yy} \begin{bmatrix} 4 & 4 & 4 & 4 & 4 \\ 1 & 1 & 1 & 1 & 1 \\ 0 & 0 & 0 & 0 & 0 \\ 1 & 1 & 1 & 1 & 1 \\ 4 & 4 & 4 & 4 & 4 \\ \end{bmatrix} + C_{xy} \begin{bmatrix} 4 & 2 & 0 & -2 & -4 \\ 2 & 1 & 0 & -1 & -2 \\ 0 & 0 & 0 & 0 & 0 \\ -2 & -1 & 0 & 1 & 2 \\ -4 & -2 & 0 & 2 & 4 \\ \end{bmatrix} \end{aligned}

These are not the only basis matrices, but they are probably the easiest to understand.

Note that the given Gaussian achieves its maximum at $$x=y=0$$. So that value corresponds to the center of the matrix. The corner values are given by $$G(2,2)$$. Furthermore, the values are quantized. You can try to estimate the chosen value of $$\alpha$$ from the given values.

EDIT: If you assume $$\alpha=1$$ and you evaluate the 2D-Gaussian, multiply it by $$273$$ and then round the values to integers, you arrive at the following matrix:

$$\begin{matrix} 1 & 4 & 6 & 4 & 1\\ 4 & 16 & 26 & 16 & 4\\ 6 & 26 & 43 & 26 & 6\\ 4 & 16 & 26 & 16 & 4\\ 1 & 4 & 6 & 4 & 1 \end{matrix}$$

This is already pretty close to the matrix you got. In practice you won't notice much difference between these two matrices. It could be that the matrix you have is the result of some averaging of the Gaussian across the pixels, i.e., instead of evaluating the Gaussian at one point you average over the pixel surface. But, again, I don't think that this will make much of a difference when using that matrix as a convolution kernel.

• Thanks for your clarification. But still $G(2,2)$ does not match the answer $\frac{1}{273}$. – Swakshar Deb Aug 6 at 6:17
• @SwaksharDeb: Please take a look at my edited answer. – Matt L. Aug 6 at 9:15
• @MattL. I've actually been studying the "Discrete Gaussian" vs a "Sampled Gaussian" as the Eigenfunction of the DFT. I was wondering if the discrepancy might be the quantified version of that. Or just a tweak to get 273 for some reason. 271 gives integers that add up to the denominator. – Cedron Dawg Aug 6 at 12:18
• @CedronDawg: I think they probably averaged the Gaussian over the size of one pixel, as mentioned in my answer. – Matt L. Aug 6 at 13:46
• Olli's link under the question has an excellent work up by LD. In the past, I've always used binomials and shifting when going integer. – Cedron Dawg Aug 6 at 14:03