How to Make a Gaussian Filter?

Question

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.

Answer 1

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.

Answer 2

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.