2D Gaussian Derivation

I'm currently working on some image manipulation that requires a Gaussian Point Spread function that isn't uniform in the x and y directions so thought its worth revisiting the derivation from an older blog post along with some thoughts on optimization:

In one dimension the Gaussian function looks like:

\$f(x) = Ae^{- \frac{(x-b)^2} {2\sigma^2} } \$

where \$\sigma\approx2.718281828 \$  which is Euler's Number, and b is the point over which the bell curve will be centred. You should recognize the \$(x-b)^2\$  as the first step in calculating the distance between two points. A is the amplitude of the function. The bigger A is, the higher the peak produced.

As can been seen from this equation \$\sigma\$ controls the spread of the bell shaped curve produced. If its not immediately obvious then keep in mind as you divide by a bigger number then the fraction gets smaller and \$anything^0 = 1 \$

Now that we understand the function in one dimension lets extend it to 2 dimensions, after all that's what I am interested in for my image manipulation.

It is very simple to extend the function to 2 dimensions as we are really looking at the distance of a point from a centre location. For now lets assume our \$\sigma\$ (curve spread) is the same in each direction and that our centre point is \$bx,by\$

\$f(x,y) = Ae^{- \frac{(x-b_x)^2+(y-b_y)^2} {2\sigma^2} } \$

Again you should see the \$(x-b_x)^2+(y-b_y)^2\$ as the distance of the two points from the centre point - we are just missing the square root.

This is the equation we used to calculate our Gaussian PSF kernel, for this example we are going to use the following parameters:

\$A=15\$ and \$\sigma=1.4\$

We then take the points from -2 to 2  = 5  in each direction and plug them into our equation as x,y values. The resultant value is rounded and stored in our matrix.

For Example: x=2, y=2

\$f(2,2) = Ae^{-(\frac{2^2+2^2}{2(1.4^2)})}\$

\$f(2,2) = Ae^{-(\frac{8}{3.92})}\$

\$f(2,2) = Ae^{-2.0408}\$

\$f(2,2) = 15*0.1299\$

\$f(2,2) = 1.9488\$    which can be rounded up to 2 and this is the value we store at 2,2 in our kernel.

Note that as we are squaring our differences from the position to the centre in each direction the values at: -2,-2  ; -2,2  ; 2,-2   are all the same as the 2,2 one calculated above. We can use this as an optimization in calculating our kernel coefficients.

As mentioned above the \$\sigma\$ value can vary for x and y directions. This causes our 2D curve to be stretched/compressed in the x or y direction.

Our equation then becomes:

\$f(x,y) = Ae^{-( \frac{(x-b_x)^2}{2\sigma_x^2}+\frac{(y-b_y)^2} {2\sigma_y^2} )} \$

knowing about our one dimensional Gaussian function we can clearly see how the above function works: the two component directions are calculated first, added together, then the result used as the power to which we are raising e.

Once we have calculated our kernel coefficients we can apply them to our image. Remember that this is a separable convolution so don't implement it in the trivial manner by reading in a 5x5 block around your pixel of interest, multiplying, adding and finally dividing. Rather, apply the x and y convolutions separately which will involve 4 times fewer reads and these reads will be in a more cache friendly manner.