Search code examples
c++convolutiongaussianvariance

How to create a Gaussian kernel of arbitrary width?


How to create a Gaussian kernel by only specifying its width w (3,5,7,9...), and without specifying its variance sigma?

In other word, how to adapt sigma so that the Gaussian distribution 'fits well' w?

I would be interested in a C++ implementation:

void create_gaussian_kernel(int w, std::vector<std::vector<float>>& kernel)
{
    kernel = std::vector<std::vector<float>>(w, std::vector<float>(w, 0.f)); // 2D array of size w x w 
    const Scalar sigma = 1.0; // how to adapt sigma to w ???
    const int hw = (w-1)/2; // half width

    for(int di = -hw; di <= +hw; ++di)
    {
        const int i = hw + di;
        for(int dj = -hw; dj <= +hw; ++dj)
        {
            const int j = hw + dj;
            kernel[i][j] = gauss2D(di, dj, sigma);
        }
    } 
}

Everything I see on the Internet use a fixed size w and a fixed variance sigma :


Solution

  • I found a simple (arbitrary) relation between sigma and w.

    I want the next value outside the kernel (along one axis) below a very small value epsilon:

    exp( - (half_width + 1)^2 / (2 * sigma^2) ) < epsilon
    

    with half_width the kernel 'half width'.

    The result is

    sigma^2 = - (half_width + 1)^2 / (2 * log(epsilon))
    

    I use the following c++ code:

    #include <vector>
    #include <cmath>
    #include <cassert>
    
    using Matrix = std::vector<std::vector<float>>;
    
    // compute sigma^2 that 'fit' the kernel half width 
    float compute_squared_variance(int half_width, float epsilon = 0.001)
    {
        assert(0 < epsilon && epsilon < 1); // small value required
        return - (half_width + 1.0) * (half_width + 1.0) / 2.0 / std::log(epsilon);
    }
    
    float gaussian_exp(float y, float x, float sigma2)
    {
        assert(0 < sigma2);
        return std::exp( - (x*x + y*y) / (2 * sigma2) );
    }
    
    // create a Gaussian kernel of size 2*half_width+1 x 2*half_width+1
    Matrix make_gaussian_kernel(int half_width)
    {
        if(half_width <= 0)
        {
            // kernel of size 1 x 1
            Matrix kernel(1, std::vector<float>(1, 1.0));
            return kernel;
        }
    
        Matrix kernel(2*half_width+1, std::vector<float>(2*half_width+1, 0.0));
    
        const float sigma2 = compute_squared_variance(half_width, 0.1);
    
        float sum = 0;
        for(int di = -half_width; di <= +half_width; ++di)
        {
            const int i = half_width + di;
            for(int dj = -half_width; dj <= +half_width; ++dj)
            {
                const int j = half_width + dj;
                kernel[i][j] = gaussian_exp(di, dj, sigma2);
                sum += kernel[i][j];
            }
        }
    
        assert(0 < sum);
    
        // normalize 
        for(int i=0; i<2*half_width+1; ++i)
        {
            for(int j=0; j<2*half_width+1; ++j)
            {
                kernel[i][j] /= sum;
            }
        }
    
        return kernel;
    }