# Create a two-dimensional Gaussian kernel

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

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_y.reserve(rows);
for (auto i = 0;  i < rows;  ++i) {
auto y = i - y_mid;
}

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

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.

• Can you re execute your codes with the changes I made, please? and also I didn't succeed to work with the second short version.. it doesn't give the same result as I want! I still working to see why! – Ja_cpp Jul 19 '17 at 20:34
• How much different? It's possible I made a mistake with the arithmetic of denominator (or elsewhere). – Toby Speight Jul 19 '17 at 21:04
• denominator is correct, the first version is correct. I'll check that! – Ja_cpp Jul 20 '17 at 7:14
• I've edited your answer because of some mistakes but my suggestions have been rejected, why? – Ja_cpp Jul 20 '17 at 7:14
• Thanks for your edit suggestion. It seems that the reviewers rejected because the edit deviates from the original intent of the post. I've incorporated the parts that were correct, and fixed the other issues you identified. The denominator value is the normalizing scale factor; I had missed the division of x by sigmax and y by sigmay. – Toby Speight Jul 20 '17 at 8:15
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.

• Thanks it makes the program little bit faster. I'll try 2. The test for temp > eps is to eliminate very small values smaller than eps which I define as #define eps 2.2204e-16 – Ja_cpp Jul 19 '17 at 20:26
• Arguably, it may be worth testing to eliminate subnormals (which can be expensive later). If that's the intent, it may be better to use std::fpclassify() rather than <. See also Avoiding denormal values in C++. – Toby Speight Jul 20 '17 at 7:47

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;

}

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

• Actually, the core of my answer and vnp's boils down to mostly the same suggestion (break down the e^x and e^y), but expressed in different words. I'm glad we were both helpful to you! – Toby Speight Aug 7 '17 at 9:26