Create a two-dimensional Gaussian kernel
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 &lt;iostream&gt;
int main()
{
    int rows = 20000, cols = 20000;
    const auto kernel = getGaussianKernel(rows, cols, 50, 50 );
    std::cout &lt;&lt; kernel.total() &lt;&lt; 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 &lt;opencv2/core/core.hpp&gt;
#include &lt;cmath&gt;

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&lt;double&gt; gauss_x, gauss_y;

    gauss_x.reserve(cols);
    for (auto i = 0;  i &lt; 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 &lt; 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 &lt; rows;  ++j)
        for (auto i = 0;  i &lt; cols;  ++i) {
            kernel.at&lt;float&gt;(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_&lt;float&gt;(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 &lt; cols;  ++i) {
        auto const x = i - x_mid;
        gauss_x[i][0] = std::exp(-x*x * x_spread);
    }

    auto kernel = cv::Mat_&lt;float&gt;(rows, cols);
    for (auto i = 0;  i &lt; 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 &lt;opencv2/core/core.hpp&gt;
#include &lt;opencv2/imgproc/imgproc.hpp&gt;

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.