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
#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;
}