Saturday, March 2, 2013

Gaussian blurring using separable kernel in C++

Gaussian blurring is used to reduce the noise and details of the image. It reduces the image's high frequency components and thus it is type of low pass filter. Gaussian blurring is obtained by convolving the image with Gaussian function. For more information about Gaussian function see the Wikipedia page.  Since 2D Gaussian function can be obtained by multiplying two 1D Gaussian functions, the blurring can be obtained by using separable kernel.

Use of Separable Kernel

Convolution is very expensive computationally. For an n x n kernel requires n2 multiplication and the same number of additions per pixel, and there are typically 105 - 106 pixels per image. So a much more efficient algorithm can be used for convolution in the small number of cases where a kernel is separable. A separable n x n kernel can be represented as a vector product of two orthogonal, one-dimensional kernels, each of width n. Convolution with the kernel can be carried out using there one-dimensional components. One of them is applied down the columns of an image, generating intermediate result. The other kernel is then applied along the rows of the intermediate image, producing the final result.

Complexity of standard kernel O(n2)
Complexity of separable kernel O(n)

It is estimated that convolution with a separable 15 x 15 kernel requires just 13 percent of the computation necessary when a non-separable kernel is used.

How to determine if a Kernel is separable?
The easiest way to test if a Kernel is separable is to calculate the rank of the matrix. If the rank is 1, then it is separable (i.e. you can break down the matrix in two vectors 1D vectors:). The rank of the Gaussian Kernel is therefore 1.

In the code below I have used 1D Gaussian function. The discrete value of 1D Gaussian function is calculated using this method and is given by

double coeff[] = {0.0545, 0.2442, 0.4026, 0.2442, 0.0545}


Here sigma = 1 and centre of distribution is assumed to be at x = 0

Source Code : C++

#include<iostream>
#include<opencv2/imgproc/imgproc.hpp>
#include<opencv2/highgui/highgui.hpp>

using namespace std;
using namespace cv;

// reflected indexing for border processing
int reflect(int M, int x)
{
if(x < 0)
{
return -x - 1;
}
if(x >= M)
{
return 2*M - x - 1;
}
    return x;
}

int main()
{

Mat src, dst, temp;
float sum, x1, y1;

/// Load an image

if( !src.data )
{ return -1; }

// coefficients of 1D gaussian kernel with sigma = 1
double coeffs[] = {0.0545, 0.2442, 0.4026, 0.2442, 0.0545};

dst = src.clone();
temp = src.clone();

// along y - direction
for(int y = 0; y < src.rows; y++){
for(int x = 0; x < src.cols; x++){
sum = 0.0;
for(int i = -2; i <= 2; i++){
y1 = reflect(src.rows, y - i);
sum = sum + coeffs[i + 2]*src.at<uchar>(y1, x);
}
temp.at<uchar>(y,x) = sum;
}
}

// along x - direction
for(int y = 0; y < src.rows; y++){
for(int x = 0; x < src.cols; x++){
sum = 0.0;
for(int i = -2; i <= 2; i++){
x1 = reflect(src.cols, x - i);
sum = sum + coeffs[i + 2]*temp.at<uchar>(y, x1);
}
dst.at<uchar>(y,x) = sum;
}
}

namedWindow("final");
imshow("final", dst);

namedWindow("initial");
imshow("initial", src);

waitKey();

return 0;
}


1. Hello!
Nice tutorial, can you please reffer to some literature that explains the principle behind a separable kernel?, I've been trying to find it myself, but I haven't been able to.

For any help you can provide, Thankyou very much!

2. @Roberto, the book Real-Time Rendering is a great book all kind of rendering, and has information on gaussian blur, both "regular" and with separable kernel.

@Bibek, your coefficients are wrong! It adds up to 0.96. Your two 0.224 values should be 0.2442.

1. @Eldon nice catch up. Thank you very much. I think I had written this post in sleepy mind :P

3. in reflect you need add return x at the end

1. Thank you so much. It is corrected.

4. Good morning Mr. Bibek.
I want to ask, why would you use a loop from -2 to 2?
thank you

1. Hello Indra
I have used the Gaussian Kernel of size 5x5. So starting from -2 to 2 will give the expected result. If your query is why not to use from 1 to 5, you must look at the convolution operation first.

5. Hello Bibek,
I tested this codes, but not show source and dst image.
The built did not include any errors.
Thank you !

6. Hello Bibek,
I just tested this codes, but it did not show the src and dst image. Although the program built success, not any error.
Can you talk the causal ?
Thanks !

7. Wrong code. You'll need a change loops in the vertical - along y block. Also it must be combined. Very basic poor code. Shukra jazilan. Inwala