Return to Answer

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

Note that my code doesn't initialize the kernel with zeroes - 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.

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

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.

This is also shorter, and arguably easier on the eye.

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.

#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 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));
 }
 gauss_y.reserve(rows);
 for (auto i = 0; i < rows; ++i) {
 auto y = i - y_mid;
 gauss_y.push_back(std::exp(-y*y));
 }
 cv::Mat kernel = cv::Mat::zeros(rows, cols, CV_32FC1);
 for (auto j = 0; j < rows; ++j)
 for (auto i = 0; i < rows; ++i) {
 kernel.at<float>(j,i) = gauss_x[i] * gauss_y[j] / denominator;
 }
 return kernel;
}

cv::Mat getGaussianKernel(int rows, int cols, double sigmax, double sigmay)
{
 auto gauss_x = cv::Mat_<float>(rows, 1);
 const auto x_mid = (cols-1) / 2.0;
 const auto y_mid = (rows-1) / 2.0;
 const auto denominator = 8 * atan(1) * sigmax * sigmay;
 for (auto i = 0; i < cols; ++i) {
 auto const x = i - x_mid;
 gauss_x[i][0] = std::exp(-x*x);
 }
 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)/denominator);
 }
 return kernel;
}

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

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