Implementing Gaussian Blur - How to calculate convolution matrix (kernel)

Question

My question is very close to this question: How do I gaussian blur an image without using any in-built gaussian functions?

The answer to this question is very good, but it doesn't give an example of actually calculating a real Gaussian filter kernel. The answer gives an arbitrary kernel and shows how to apply the filter using that kernel but not how to calculate a real kernel itself. I am trying to implement a Gaussian blur in C++ or Matlab from scratch, so I need to know how to calculate the kernel from scratch.

I'd appreciate it if someone could calculate a real Gaussian filter kernel using any small example image matrix.

Answer 1

Here's a calculation in C#, which does not take single samples from the gaussian (or another kernel) function, but it calculates a large number of samples in a small grid and integrates the samples in the desired number of sections.
The calculation is for 1D, but it may easily be extended to 2D.

This calculation uses some other functions, which I did not add here, but I have added the function signatures so that you will know what they do.

This calculation produces the following discrete values for the limits +/- 3 (sum areaSum of integral is 0.997300):

kernel size: normalized kernel values, rounded to 6 decimals  
3:                     0.157731, 0.684538, 0.157731  
5:           0.034674, 0.238968, 0.452716, 0.238968, 0.034674  
7: 0.014752, 0.083434, 0.235482, 0.332663, 0.235482, 0.083434, 0.014752

This calculation produces the following discrete values for the limits +/- 2 (sum areaSum of integral is 0.954500):

kernel size: normalized kernel values, rounded to 6 decimals  
3:                     0.240694, 0.518612, 0.240694
5:           0.096720, 0.240449, 0.325661, 0.240449, 0.096720
7: 0.056379, 0.124798, 0.201012, 0.235624, 0.201012, 0.124798, 0.056379

Code:

using System.Linq;

private static void Main ()
{
  int    positionCount    = 1024;  // Number of samples in the range 0..1.
  double positionStepSize = 1.0 / positionCount;
  double limit            = 3;  // The calculation range of the kernel. +/- 3 is sufficient for gaussian.
  int    sectionCount     = 3;  // The number of elements in the kernel. May be 1, 3, 5, 7, ... (n*2+1)

  // calculate the x positions for each kernel value to calculate.
  double[] positions = CreateSeries (-limit, positionStepSize, (int)(limit * 2 * positionCount + 1));
  // calculate the gaussian function value for each position
  double[] values    = positions.Select (pos => Gaussian (pos)).ToArray ();
  // split the values into equal-sized sections and calculate the integral of each section.
  double[] areas     = IntegrateInSections (values, positionStepSize, sectionCount);
  double   areaSum   = areas.Sum ();
  // normalize to 1
  double[] areas1    = areas.Select (a => a / areaSum).ToArray ();
  double   area1Sum  = areas1.Sum ();  // just to check it's 1 now
}

///-------------------------------------------------------------------
/// <summary>
/// Create a series of <paramref name="i_count"/> numbers, starting at <paramref name="i_start"/> and increasing by <paramref name="i_stepSize"/>.
/// </summary>
/// <param name="i_start">The start value of the series.</param>
/// <param name="i_stepSize">The step size between values in the series.</param>
/// <param name="i_count">The number of elements in the series.</param>
///-------------------------------------------------------------------
public static double[] CreateSeries (double i_start,
                                     double i_stepSize,
                                     int    i_count)
{ ... }

private static readonly double s_gaussian_Divisor = Math.Sqrt (Math.PI * 2.0);

/// ------------------------------------------------------------------
/// <summary>
///   Calculate the value for the given position in a Gaussian kernel.
/// </summary>
/// <param name="i_position"> The position in the kernel for which the value will be calculated. </param>
/// <param name="i_bandwidth"> The width factor of the kernel. </param>
/// <returns> The value for the given position in a Gaussian kernel. </returns>
/// ------------------------------------------------------------------
public static double Gaussian (double i_position,
                               double i_bandwidth = 1)
{
  double position = i_position / i_bandwidth;
  return Math.Pow (Math.E, -0.5 * position * position) / s_gaussian_Divisor / i_bandwidth;
}

/// ------------------------------------------------------------------
/// <summary>
///   Calculate the integrals in the given number of sections of all given values with the given distance between the values.
/// </summary>
/// <param name="i_values"> The values for which the integral will be calculated. </param>
/// <param name="i_distance"> The distance between the values. </param>
/// <param name="i_sectionCount"> The number of sections in the integration. </param>
/// ------------------------------------------------------------------
public static double[] IntegrateInSections (IReadOnlyCollection<double> i_values,
                                            double                      i_distance,
                                            int                         i_sectionCount)
{ ... }

Answer 2

OK, a late answer but in case of...

Using the @moooeeeep answer, but with numpy;

import numpy as np
radius = 3
sigma = radius/2.

k = np.arange(2*radius +1)
row = np.exp( -(((k - radius)/(sigma))**2)/2.)
col = row.transpose()
out = np.outer(row, col)
out = out/np.sum(out)
for line in out:
    print(["%.3f" % x for x in line])

Just a bit less of lines.

Answer 3

// my_test.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"

#include <cmath>
#include <vector>
#include <iostream>
#include <iomanip>
#include <string>

//https://stackoverflow.com/questions/8204645/implementing-gaussian-blur-how-to-calculate-convolution-matrix-kernel
//https://docs.opencv.org/2.4/modules/imgproc/doc/filtering.html#getgaussiankernel
//http://dev.theomader.com/gaussian-kernel-calculator/

double gaussian(double x, double mu, double sigma) {
    const double a = (x - mu) / sigma;
    return std::exp(-0.5 * a * a);
}

typedef std::vector<double> kernel_row;
typedef std::vector<kernel_row> kernel_type;

kernel_type produce2dGaussianKernel(int kernelRadius, double sigma) {
    kernel_type kernel2d(2 * kernelRadius + 1, kernel_row(2 * kernelRadius + 1));
    double sum = 0;
    // compute values
    for (int row = 0; row < kernel2d.size(); row++)
        for (int col = 0; col < kernel2d[row].size(); col++) {
            double x = gaussian(row, kernelRadius, sigma)
                * gaussian(col, kernelRadius, sigma);
            kernel2d[row][col] = x;
            sum += x;
        }
    // normalize
    for (int row = 0; row < kernel2d.size(); row++)
        for (int col = 0; col < kernel2d[row].size(); col++)
            kernel2d[row][col] /= sum;
    return kernel2d;
}

char* gMatChar[10] = {
    "          ",
    "         ",
    "        ",
    "       ",
    "      ",
    "     ",
    "    ",
    "   ",
    "  ",
    " "
};

static int countSpace(float aValue)
{
    int count = 0;
    int value = (int)aValue;
    while (value > 9)
    {
        count++;
        value /= 10;
    }
    return count;
}

int main() {
    while (1)
    {
        char str1[80]; // window size
        char str2[80]; // sigma
        char str3[80]; // coefficient
        int space;

        int i, ch;
        printf("\n-----------------------------------------------------------------------------\n");
        printf("Start generate Gaussian matrix\n");
        printf("-----------------------------------------------------------------------------\n");
        // input window size
        printf("\nPlease enter window size (from 3 to 10) It should be odd (ksize/mod 2 = 1 ) and positive: Exit enter q \n");
        for (i = 0; (i < 80) && ((ch = getchar()) != EOF)
            && (ch != '\n'); i++)
        {
            str1[i] = (char)ch;
        }

        // Terminate string with a null character
        str1[i] = '\0';
        if (str1[0] == 'q')
        {
            break;
        }
        int input1 = atoi(str1);
        int window_size = input1 / 2;
        printf("Input window_size was: %d\n", input1);

        // input sigma
        printf("Please enter sigma. Use default press Enter . Exit enter q \n");
        str2[0] = '0';
        for (i = 0; (i < 80) && ((ch = getchar()) != EOF)
            && (ch != '\n'); i++)
        {
            str2[i] = (char)ch;
        }

        // Terminate string with a null character
        str2[i] = '\0';
        if (str2[0] == 'q')
        {
            break;
        }
        float input2 = atof(str2);
        float sigma;
        if (input2 == 0)
        {
            // Open-CV sigma � Gaussian standard deviation. If it is non-positive, it is computed from ksize as sigma = 0.3*((ksize-1)*0.5 - 1) + 0.8 .
            sigma = 0.3*((input1 - 1)*0.5 - 1) + 0.8;
        }
        else
        {
            sigma = input2;
        }
        printf("Input sigma was: %f\n", sigma);

        // input Coefficient K
        printf("Please enter Coefficient K. Use default press Enter . Exit enter q \n");
        str3[0] = '0';
        for (i = 0; (i < 80) && ((ch = getchar()) != EOF)
            && (ch != '\n'); i++)
        {
            str3[i] = (char)ch;
        }

        // Terminate string with a null character
        str3[i] = '\0';
        if (str3[0] == 'q')
        {
            break;
        }
        int input3 = atoi(str3);
        int cK;
        if (input3 == 0)
        {
            cK = 1;
        }
        else
        {
            cK = input3;
        }
        float sum_f = 0;
        float temp_f;
        int sum = 0;
        int temp;
        printf("Input Coefficient K was: %d\n", cK);

        printf("\nwindow size=%d | Sigma = %f Coefficient K = %d\n\n\n", input1, sigma, cK);

        kernel_type kernel2d = produce2dGaussianKernel(window_size, sigma);
        std::cout << std::setprecision(input1) << std::fixed;
        for (int row = 0; row < kernel2d.size(); row++) {
            for (int col = 0; col < kernel2d[row].size(); col++)
            {
                temp_f = cK* kernel2d[row][col];
                sum_f += temp_f;
                space = countSpace(temp_f);
                std::cout << gMatChar[space] << temp_f << ' ';
            }
            std::cout << '\n';
        }
        printf("\n Sum array = %f | delta = %f", sum_f, sum_f - cK);

        // rounding
        printf("\nRecommend use round(): window size=%d | Sigma = %f Coefficient K = %d\n\n\n", input1, sigma, cK);
        sum = 0;
        std::cout << std::setprecision(0) << std::fixed;
        for (int row = 0; row < kernel2d.size(); row++) {
            for (int col = 0; col < kernel2d[row].size(); col++)
            {
                temp = round(cK* kernel2d[row][col]);
                sum += temp;
                space = countSpace((float)temp);
                std::cout << gMatChar[space] << temp << ' ';
            }
            std::cout << '\n';
        }
        printf("\n Sum array = %d | delta = %d", sum, sum - cK);

        // recommented
        sum_f = 0;
        int cK_d = 1 / kernel2d[0][0];
        cK_d = cK_d / 2 * 2;
        printf("\nRecommend: window size=%d | Sigma = %f Coefficient K = %d\n\n\n", input1, sigma, cK_d);
        std::cout << std::setprecision(input1) << std::fixed;
        for (int row = 0; row < kernel2d.size(); row++) {
            for (int col = 0; col < kernel2d[row].size(); col++)
            {
                temp_f = cK_d* kernel2d[row][col];
                sum_f += temp_f;
                space = countSpace(temp_f);
                std::cout << gMatChar[space] << temp_f << ' ';
            }
            std::cout << '\n';
        }
        printf("\n Sum array = %f | delta = %f", sum_f, sum_f - cK_d);

        // rounding
        printf("\nRecommend use round(): window size=%d | Sigma = %f Coefficient K = %d\n\n\n", input1, sigma, cK_d);
        sum = 0;
        std::cout << std::setprecision(0) << std::fixed;
        for (int row = 0; row < kernel2d.size(); row++) {
            for (int col = 0; col < kernel2d[row].size(); col++)
            {
                temp = round(cK_d* kernel2d[row][col]);
                sum += temp;
                space = countSpace((float)temp);
                std::cout << gMatChar[space] << temp << ' ';
            }
            std::cout << '\n';
        }
        printf("\n Sum array = %d | delta = %d", sum, sum - cK_d);

    }
}

Answer 4

To implement the gaussian blur you simply take the gaussian function and compute one value for each of the elements in your kernel.

Usually you want to assign the maximum weight to the central element in your kernel and values close to zero for the elements at the kernel borders. This implies that the kernel should have an odd height (resp. width) to ensure that there actually is a central element.

To compute the actual kernel elements you may scale the gaussian bell to the kernel grid (choose an arbitrary e.g. sigma = 1 and an arbitrary range e.g. -2*sigma ... 2*sigma) and normalize it, s.t. the elements sum to one. To achieve this, if you want to support arbitrary kernel sizes, you might want to adapt the sigma to the required kernel size.

Here's a C++ example:

#include <cmath>
#include <vector>
#include <iostream>
#include <iomanip>

double gaussian( double x, double mu, double sigma ) {
    const double a = ( x - mu ) / sigma;
    return std::exp( -0.5 * a * a );
}

typedef std::vector<double> kernel_row;
typedef std::vector<kernel_row> kernel_type;

kernel_type produce2dGaussianKernel (int kernelRadius) {
  double sigma = kernelRadius/2.;
  kernel_type kernel2d(2*kernelRadius+1, kernel_row(2*kernelRadius+1));
  double sum = 0;
  // compute values
  for (int row = 0; row < kernel2d.size(); row++)
    for (int col = 0; col < kernel2d[row].size(); col++) {
      double x = gaussian(row, kernelRadius, sigma)
               * gaussian(col, kernelRadius, sigma);
      kernel2d[row][col] = x;
      sum += x;
    }
  // normalize
  for (int row = 0; row < kernel2d.size(); row++)
    for (int col = 0; col < kernel2d[row].size(); col++)
      kernel2d[row][col] /= sum;
  return kernel2d;
}

int main() {
  kernel_type kernel2d = produce2dGaussianKernel(3);
  std::cout << std::setprecision(5) << std::fixed;
  for (int row = 0; row < kernel2d.size(); row++) {
    for (int col = 0; col < kernel2d[row].size(); col++)
      std::cout << kernel2d[row][col] << ' ';
    std::cout << '\n';
  }
}

The output is:

$ g++ test.cc && ./a.out
0.00134 0.00408 0.00794 0.00992 0.00794 0.00408 0.00134 
0.00408 0.01238 0.02412 0.03012 0.02412 0.01238 0.00408 
0.00794 0.02412 0.04698 0.05867 0.04698 0.02412 0.00794 
0.00992 0.03012 0.05867 0.07327 0.05867 0.03012 0.00992 
0.00794 0.02412 0.04698 0.05867 0.04698 0.02412 0.00794 
0.00408 0.01238 0.02412 0.03012 0.02412 0.01238 0.00408 
0.00134 0.00408 0.00794 0.00992 0.00794 0.00408 0.00134 

As a simplification you don't need to use a 2d-kernel. Easier to implement and also more efficient to compute is to use two orthogonal 1d-kernels. This is possible due to the associativity of this type of a linear convolution (linear separability). You may also want to see this section of the corresponding wikipedia article.

---

Here's the same in Python (in the hope someone might find it useful):

from math import exp

def gaussian(x, mu, sigma):
  return exp( -(((x-mu)/(sigma))**2)/2.0 )

#kernel_height, kernel_width = 7, 7
kernel_radius = 3 # for an 7x7 filter
sigma = kernel_radius/2. # for [-2*sigma, 2*sigma]

# compute the actual kernel elements
hkernel = [gaussian(x, kernel_radius, sigma) for x in range(2*kernel_radius+1)]
vkernel = [x for x in hkernel]
kernel2d = [[xh*xv for xh in hkernel] for xv in vkernel]

# normalize the kernel elements
kernelsum = sum([sum(row) for row in kernel2d])
kernel2d = [[x/kernelsum for x in row] for row in kernel2d]

for line in kernel2d:
  print ["%.3f" % x for x in line]

produces the kernel:

['0.001', '0.004', '0.008', '0.010', '0.008', '0.004', '0.001']
['0.004', '0.012', '0.024', '0.030', '0.024', '0.012', '0.004']
['0.008', '0.024', '0.047', '0.059', '0.047', '0.024', '0.008']
['0.010', '0.030', '0.059', '0.073', '0.059', '0.030', '0.010']
['0.008', '0.024', '0.047', '0.059', '0.047', '0.024', '0.008']
['0.004', '0.012', '0.024', '0.030', '0.024', '0.012', '0.004']
['0.001', '0.004', '0.008', '0.010', '0.008', '0.004', '0.001']

Answer 5

 function kernel = gauss_kernel(m, n, sigma)
 % Generating Gauss Kernel

 x = -(m-1)/2 : (m-1)/2;
 y = -(n-1)/2 : (n-1)/2;

 for i = 1:m
     for j = 1:n
         xx(i,j) = x(i);
         yy(i,j) = y(j);
     end
 end

 kernel = exp(-(xx.*xx + yy.*yy)/(2*sigma*sigma));

 % Normalize the kernel
 kernel  = kernel/sum(kernel(:));

 % Corresponding function in MATLAB
 % fspecial('gaussian', [m n], sigma)

Answer 6

Gaussian blur in python using PIL image library. For more info read this: http://blog.ivank.net/fastest-gaussian-blur.html

from PIL import Image
import math

# img = Image.open('input.jpg').convert('L')
# r = radiuss
def gauss_blur(img, r):
    imgData = list(img.getdata())

    bluredImg = Image.new(img.mode, img.size)
    bluredImgData = list(bluredImg.getdata())

    rs = int(math.ceil(r * 2.57))

    for i in range(0, img.height):
        for j in range(0, img.width):
            val = 0
            wsum = 0
            for iy in range(i - rs, i + rs + 1):
                for ix in range(j - rs, j + rs + 1):
                    x = min(img.width - 1, max(0, ix))
                    y = min(img.height - 1, max(0, iy))
                    dsq = (ix - j) * (ix - j) + (iy - i) * (iy - i)
                    weight = math.exp(-dsq / (2 * r * r)) / (math.pi * 2 * r * r)
                    val += imgData[y * img.width + x] * weight
                    wsum += weight 
            bluredImgData[i * img.width + j] = round(val / wsum)

    bluredImg.putdata(bluredImgData)
    return bluredImg

Answer 7

It's as simple as it sounds:

double sigma = 1;
int W = 5;
double kernel[W][W];
double mean = W/2;
double sum = 0.0; // For accumulating the kernel values
for (int x = 0; x < W; ++x) 
    for (int y = 0; y < W; ++y) {
        kernel[x][y] = exp( -0.5 * (pow((x-mean)/sigma, 2.0) + pow((y-mean)/sigma,2.0)) )
                         / (2 * M_PI * sigma * sigma);

        // Accumulate the kernel values
        sum += kernel[x][y];
    }

// Normalize the kernel
for (int x = 0; x < W; ++x) 
    for (int y = 0; y < W; ++y)
        kernel[x][y] /= sum;

Answer 8

You can create a Gaussian kernel from scratch as noted in MATLAB documentation of fspecial. Please read the Gaussian kernel creation formula in the algorithms part in that page and follow the code below. The code is to create an m-by-n matrix with sigma = 1.

m = 5; n = 5;
sigma = 1;
[h1, h2] = meshgrid(-(m-1)/2:(m-1)/2, -(n-1)/2:(n-1)/2);
hg = exp(- (h1.^2+h2.^2) / (2*sigma^2));
h = hg ./ sum(hg(:));

h =

    0.0030    0.0133    0.0219    0.0133    0.0030
    0.0133    0.0596    0.0983    0.0596    0.0133
    0.0219    0.0983    0.1621    0.0983    0.0219
    0.0133    0.0596    0.0983    0.0596    0.0133
    0.0030    0.0133    0.0219    0.0133    0.0030

Observe that this can be done by the built-in fspecial as follows:

fspecial('gaussian', [m n], sigma)
ans =

    0.0030    0.0133    0.0219    0.0133    0.0030
    0.0133    0.0596    0.0983    0.0596    0.0133
    0.0219    0.0983    0.1621    0.0983    0.0219
    0.0133    0.0596    0.0983    0.0596    0.0133
    0.0030    0.0133    0.0219    0.0133    0.0030

I think it is straightforward to implement this in any language you like.

EDIT: Let me also add the values of h1 and h2 for the given case, since you may be unfamiliar with meshgrid if you code in C++.

h1 =

    -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

h2 =

    -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