Gaussian blurring using separable kernel in C++
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#include #include using namespace std; using namespace cv; // reflected indexing for border processing int reflect(int M, int x) { if(x = M) { return 2*M - x - 1; }
return x; } int main() { Mat src, dst, temp; float sum, x1, y1; /// Load an image src = imread("salt.jpg", CV_LOAD_IMAGE_GRAYSCALE); 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 (y1, x); } temp.at(y,x) = sum; } } // along x - direction for(int y = 0; y (y, x1); } dst.at (y,x) = sum; } } namedWindow("final"); imshow("final", dst); namedWindow("initial"); imshow("initial", src); waitKey(); return 0; }