How do I gaussian blur an image without using any in-built gaussian functions?

I want to blur my image using the native Gaussian blur formula. I read the Wikipedia article, but I am not sure how to implement this.

How do I use the formula to decide weights?

I do not want to use any built in functions like what MATLAB has

Writing a naive gaussian blur is actually pretty easy. It is done in exactly the same way as any other convolution filter. The only difference between a box and a gaussian filter is the matrix you use.

Imagine you have an image defined as follows:

 0  1  2  3  4  5  6  7  8  9
10 11 12 13 14 15 16 17 18 19
20 21 22 23 24 25 26 27 28 29
30 31 32 33 34 35 36 37 38 39
40 41 42 43 44 45 46 47 48 49
50 51 52 53 54 55 56 57 58 59
60 61 62 63 64 65 66 67 68 69
70 71 72 73 74 75 76 77 78 79
80 81 82 83 84 85 86 87 88 89
90 91 92 93 94 95 96 97 98 99

A 3x3 box filter matrix is defined as follows:

0.111 0.111 0.111
0.111 0.111 0.111
0.111 0.111 0.111

To apply the gaussian blur you would do the following:

For pixel 11 you would need to load pixels 0, 1, 2, 10, 11, 12, 20, 21, 22.

you would then multiply pixel 0 by the upper left portion of the 3x3 blur filter. Pixel 1 by the top middle, pixel 2, pixel 3 by top right, pixel 10 by middle left and so on.

Then add them altogether and write the result to pixel 11. As you can see Pixel 11 is now the average of itself and the surrounding pixels.

Edge cases do get a bit more complex. What values do you use for the values of the edge of the texture? One way can be to wrap round to the other side. This looks good for an image that is later tiled. Another way is to push the pixel into the surrounding places.

So for upper left you might place the samples as follows:

 0  0  1
 0  0  1
10 10 11

I hope you can see how this can easily be extended to large filter kernels (ie 5x5 or 9x9 etc).

The difference between a gaussian filter and a box filter is the numbers that go in the matrix. A gaussian filter uses a gaussian distribution across a row and column.

e.g for a filter defined arbitrarily as (ie this isn't a gaussian, but probably not far off)

0.1 0.8 0.1

the first column would be the same but multiplied into the first item of the row above.

0.01 0.8 0.1
0.08 
0.01 

The second column would be the same but the values would be multiplied by the 0.8 in the row above (and so on).

0.01 0.08 0.01
0.08 0.64 0.08
0.01 0.08 0.01

The result of adding all of the above together should equal 1. The difference between the above filter and the original box filter would be that the end pixel written would have a much heavier weighting towards the central pixel (ie the one that is in that position already). The blur occurs because the surrounding pixels do blur into that pixel, though not as much. Using this sort of filter you get a blur but one that doesn't destroy as much of the high frequency (ie rapid changing of colour from pixel to pixel) information.

These sort of filters can do lots of interesting things. You can do an edge detect using this sort of filter by subtracting the surrounding pixels from the current pixel. This will leave only the really big changes in colour (high frequencies) behind.

Edit: A 5x5 filter kernel is define exactly as above.

e.g if your row is 0.1 0.2 0.4 0.2 0.1 then if you multiply each value in their by the first item to form a column and then multiply each by the second item to form the second column and so on you'll end up with a filter of

0.01 0.02 0.04 0.02 0.01
0.02 0.04 0.08 0.04 0.02
0.04 0.08 0.16 0.08 0.04
0.02 0.04 0.08 0.04 0.02
0.01 0.02 0.04 0.02 0.01

taking some arbitrary positions you can see that position 0, 0 is simple 0.1 * 0.1. Position 0, 2 is 0.1 * 0.4, position 2, 2 is 0.4 * 0.4 and position 1, 2 is 0.2 * 0.4.

I hope that gives you a good enough explanation.

Here's the pseudo-code for the code I used in C# to calculate the kernel. I do not dare say that I treat the end-conditions correctly, though:

double[] kernel = new double[radius * 2 + 1];
double twoRadiusSquaredRecip = 1.0 / (2.0 * radius * radius);
double sqrtTwoPiTimesRadiusRecip = 1.0 / (sqrt(2.0 * Math.PI) * radius);
double radiusModifier = 1.0;

int r = -radius;
for (int i = 0; i < kernel.Length; i++)
{
    double x = r * radiusModifier;
    x *= x;
    kernel[i] = sqrtTwoPiTimesRadiusRecip * Exp(-x * twoRadiusSquaredRecip);
    r++;
}

double div = Sum(kernel);
for (int i = 0; i < kernel.Length; i++)
{
    kernel[i] /= div;
}

Hope this helps.

To use the filter kernel discussed in the Wikipedia article you need to implement (discrete) convolution. The idea is that you have a small matrix of values (the kernel), you move this kernel from pixel to pixel in the image (i.e. so that the center of the matrix is on the pixel), multiply the matrix elements with the overlapped image elements, sum all the values in the result and replace the old pixel value with this sum.

Gaussian blur can be separated into two 1D convolutions (one vertical and one horizontal) instead of a 2D convolution, which also speeds things up a bit.