Create a two-dimensional Gaussian kernel

Question

In my code I have a function that generate a 2D gaussian function given sigma like so:

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

#include <iostream>
int main(int argc, char* argv[])
{
  cv::Mat kernel = getGaussianKernel(rows, cols, 50, 50 );
  //Then I use that kernel for some processing

  //return 0;
}
cv::Mat getGaussianKernel(int rows, int cols, double sigmax, double sigmay )
    {
        cv::Mat kernel = cv::Mat::zeros(rows, cols, CV_32FC1); 

        float meanj = (kernel.rows-1)/2, 
              meani = (kernel.cols-1)/2,
              sum = 0,
              temp= 0;

        int sigma=2*sigmay*sigmax;
        for(unsigned j=0;j<kernel.rows;j++)
            for(unsigned i=0;i<kernel.cols;i++)
            {
                temp = exp( -((j-meanj)*(j-meanj) + (i-meani)*(i-meani))  / (sigma));
                if (temp > eps)
                    kernel.at<float>(j,i) = temp;

                sum += kernel.at<float>(j,i);
            }

        if(sum != 0)
            return kernel /= sum;
        else return cv::Mat();
    }

That function takes, after profiling, ~14% of my time and sigma parameter for most of the cases is constant, also rows and cols don't change because I'm using a camera and its resolution won't change.

I want to know if there is a way to use that kernel without wasting 14% of the time or to optimize its function?

I tried to save all the data generated in a file (.yml .txt 13mb) and read them but it takes the same time or little bit more

Answer 1

The two-dimensional Gaussian function can be obtained by composing two one-dimensional Gaussians.

I changed your code slightly so that it would compile (and not optimize away the unused kernel):

#include <iostream>
int main()
{
    int rows = 20000, cols = 20000;
    const auto kernel = getGaussianKernel(rows, cols, 50, 50 );
    std::cout << kernel.total() << std::endl;
}

This runs in

3.90user 0.38system 0:04.29elapsed

Changing to precompute separate x and y Gaussian curves gives a 10x improvement in speed:

#include <opencv2/core/core.hpp>
#include <cmath>

cv::Mat getGaussianKernel(int rows, int cols, double sigmax, double sigmay)
{
    const auto y_mid = (rows-1) / 2.0;
    const auto x_mid = (cols-1) / 2.0;

    const auto x_spread = 1. / (sigmax*sigmax*2);
    const auto y_spread = 1. / (sigmay*sigmay*2);

    const auto denominator = 8 * std::atan(1) * sigmax * sigmay;

    std::vector<double> gauss_x, gauss_y;

    gauss_x.reserve(cols);
    for (auto i = 0;  i < cols;  ++i) {
        auto x = i - x_mid;
        gauss_x.push_back(std::exp(-x*x * x_spread));
    }

    gauss_y.reserve(rows);
    for (auto i = 0;  i < rows;  ++i) {
        auto y = i - y_mid;
        gauss_y.push_back(std::exp(-y*y * y_spread));
    }

    cv::Mat kernel = cv::Mat::zeros(rows, cols, CV_32FC1);
    for (auto j = 0;  j < rows;  ++j)
        for (auto i = 0;  i < cols;  ++i) {
            kernel.at<float>(j,i) = gauss_x[i] * gauss_y[j] / denominator;
        }

    return kernel;
}

0.32user 0.36system 0:00.69elapsed

(both versions compiled with GCC 7.1.0 with -O3 -march=native on an Intel i7-Q6700)

The calculations are reasonably independent, so OpenMP may help. But when I tried, it always worked out slower - probably too much false sharing for it to be useful.

Note that my code doesn't initialize the kernel with zeros - it leaves it uninitialized, and instead writes to every element when it reaches it. This saves two passes over the memory. I precompute the denominator for the same reason, though that's only truly valid for a kernel of infinite size.

Also note that I changed the calculation of mid-point to divide by 2.0 rather than by 2 - I wasn't sure whether the truncation was intentional or not.

---

A further speed increase is gained by using OpenCV to work on a row at a time:

cv::Mat getGaussianKernel(int rows, int cols, double sigmax, double sigmay)
{
    auto gauss_x = cv::Mat_<float>(cols, 1);

    const auto x_mid = (cols-1) / 2.0;
    const auto y_mid = (rows-1) / 2.0;

    const auto x_spread = 1. / (sigmax*sigmax*2);
    const auto y_spread = 1. / (sigmay*sigmay*2);

    for (auto i = 0;  i < cols;  ++i) {
        auto const x = i - x_mid;
        gauss_x[i][0] = std::exp(-x*x * x_spread);
    }

    auto kernel = cv::Mat_<float>(rows, cols);
    for (auto i = 0;  i < rows;  ++i) {
        auto const y = i - y_mid;
        kernel.row(i) = gauss_x * std::exp(-y*y * y_spread);
    }

    const auto denominator = std::accumulate(kernel.begin(), kernel.end(), 0);

    return kernel / denominator;
}

0.17user 0.00system 0:00.17elapsed

This is also shorter, and arguably easier on the eye. In this version, I show an alternative to precomputing the denominator - use std::accumulate to add it together afterwards.

---

The simplest and fastest version, though, is to use the Gaussian distribution provided by OpenCV, and just matrix-multiply the two together:

#include <opencv2/core/core.hpp>
#include <opencv2/imgproc/imgproc.hpp>

cv::Mat getGaussianKernel(int rows, int cols, double sigmax, double sigmay)
{
    auto gauss_x = cv::getGaussianKernel(cols, sigmax, CV_32F);
    auto gauss_y = cv::getGaussianKernel(rows, sigmay, CV_32F);
    return gauss_x * gauss_y.t();
}

---

Do note that if you intend to use the kernel for filtering, it's best to keep it separated, and perform horizontal and vertical passes each with a one-dimensional kernel.

Answer 2

1. Factor out common subexpression. Compiler is great at optimization, but not so good at math, and doesn't know that \$e^{x+y} = e^xe^y\$

   for(unsigned j=0;j<kernel.rows;j++)
       double j_factor = exp(-(j-meanj)*(j-meanj) / (sigma));
       for(unsigned i=0;i<kernel.cols;i++)
       {
           temp = exp(-(i-meani)*(i-meani) / (sigma)) * j_factor );
           ....
   

2. Eliminate the repeated exponentiation. Most of the exponents computed in the outer loop can be reused in the inner. I recommend to precompute a linear array of exponents (in case one of means is integer, and another is half-integer, use an additional factor of \$\sqrt{e}\$).

3. I see no reason to test for temp > eps. Conditionals in the loop kill performance.

Answer 3

First thanks to @toby-speight because I'm using his solution to answer my own question and also to @vnp because I'm using his optimization idea to Eliminate the repeated exponentiation:

So the fastest solution is the first one suggested by @toby-speight. So I used his solution and also eliminated the repeated exponentiation by @vnp, which was a cool idea. Here's my implementation:

PS: I've sigmax=sigmay

#include <iostream>
int main()
{
    int rows = 20000, cols = 20000;
    const auto kernel = getGaussianKernel(rows, cols, 50, 50 );
}

#include <opencv2/highgui.hpp>

cv::Mat getGaussianKernel(int rows, int cols, float sigmax, float sigmay)
{
    const int y_mid = rows / 2; //643 482
    const int x_mid = cols / 2;

    const auto x_spread = 1. / (sigmax*sigmax*2);
    //const auto y_spread = 1. / (sigmay*sigmay*2);

    std::vector<float> gauss_x, gauss_y;

    gauss_x.reserve(cols);
    for (auto i = 0;  i < cols;  ++i) {
        auto x = i - x_mid;
        gauss_x.push_back(std::exp(-x*x*x_spread));

    }

    auto n_rows=x_mid-y_mid;
    float sum=0;
    cv::Mat kernel = cv::Mat::zeros(rows, cols, CV_32FC1);
    for (auto j = 0;  j < rows;  ++j) {
        float temp = gauss_x[n_rows+j];
        for (auto i = 0;  i < cols;  ++i) {
            kernel.at<float>(j,i) = gauss_x[i] * temp ;
            sum += kernel.at<float>(j,i);
        }
    }
    return kernel/sum;
}

This is the version that works faster in my machine i7-2630QM using gcc 5.4.0 -Ofast -march=native