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.

Solution 1

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

Solution 2

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;

Solution 3

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']

Solution 4

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.

Solution 5

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

View more solutions

Comments

[gsgx]

I typed [h1,h2] = meshgrid(-(m-1)/2:(m-1)/2, -(n-1)/2:(n-1)/2) and got a h1 ranging from -2 to 2, not from -1.5 to 1.5. Same problem with h2. But my result is the same. Also, why did you use the values of mesh grid as the values in the formula? What does this represent if you were calculating this for an image?

[petrichor]

You're right! I changed m and n to 4 in order to see whether the code works and then copied the values for this case instead of giving them for value 5. I've fixed it, thanks.

[petrichor]

The values are computed on a grid where the one in the center is the origin, which is h1==0 and h2==0 in our case. All the other pairs represent the other coordinates when you look at the h1,h2 values element by element. During filtering, you can think that this grid will be put on a pixel of the image where the origin of the grid is fit exactly on the pixel. You can read Goz's answer in the link you gave in your question for the details.

[Rookie]

This is flawed: you need to normalize the kernel too, or the image becomes darker depending on W and sigma. simply put: get the sum of kernel values and divide each kernel value by that sum.

[rayryeng]

@Rookie - I've decided to modify this post and add the normalization. This is to allow for those who want a C / C++ solution to use this straight up. Good catch!

[Panfeng Li]

There is a small mistake. We should use [h1, h2] = meshgrid(-(n-1)/2:(n-1)/2, -(m-1)/2:(m-1)/2); So that the answer will the same with fspecial when m not equal n.

[Panfeng Li]

It seems not correct when m, n are even, comparing to the result of fspecial.

[HDJEMAI]

Add some comments to your code, it will be helpful for other people.

[Grant Miller]

It might be helpful to others to add some context to your code in your solution.

[moooeeeep]

I think you can save even more lines if you use scipy or np.convolve as suggested here: stackoverflow.com/q/29920114/1025391

[Tobias Knauss]

I think it's not sufficient to use just one sample value of the kernel function (let it be gaussian or another kernel). The values in the center would be too small, and the values on the borders would be too large. You should sample in a much smaller grid (e.g. 0.001) in the range -3...+3, split the range into the desired number of kernel elements, and integrate all values in each range. This is also suggested by homepages.inf.ed.ac.uk/rbf/HIPR2/gsmooth.htm. I have added a C# example below: stackoverflow.com/a/71968356/2505186

[Cris Luengo]

@TobiasKnauss Your suggested process results in the convolution of a Gaussian with a small square window, and therefore is not a Gaussian any more. Point sampling is the correct procedure.

[Tobias Knauss]

This solution is definitely wrong, because it just calculates the kernel for a larger range, but it does not 'spread' the kernel, i.e. it does not adjust the bandwidth (sigma).