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.