How to generate 2d gaussian kernel using 2d convolution in python?

Question

From my workout instruction:

A 2D Gaussian can be formed by convolution of a 1D Gaussian with its transpose.

Here is my 1d gaussian function:

def gauss1d(sigma, filter_length=11):
    # INPUTS
    # @ sigma         : sigma of gaussian distribution
    # @ filter_length : integer denoting the filter length
    # OUTPUTS
    # @ gauss_filter  : 1D gaussian filter without normalization

    rng = range(-int(filter_length/2),int(filter_length/2)+1)
    gauss_filter = [np.exp((-x**2) / (2*sigma**2)) for x in rng]

    # The formula used above has been given in the instruction.
    return np.array(gauss_filter)

And 2d convolution function which performs a 2D convolution between image and filt, image being a 2D image.

def myconv2(image, filt):
    # INPUTS
    # @ image         : 2D image, as numpy array of size mxn
    # @ filt          : 1D or 2D filter of size kxl
    # OUTPUTS
    # img_filtered    : 2D filtered image, of size (m+k-1)x(n+l-1)

    m, n = image.shape
    k, l = filt.shape
    offsety = k // 2
    offsetx = l // 2
    img_filtered = np.zeros((m+k-1, n+l-1), "double")
    image = np.pad(image, ((offsety,offsety),(offsetx, offsetx)), mode='constant')
    
    for i in range(offsety, m+offsety):
        for j in range(offsetx, n+offsetx):
            box_vals = image[ i - offsety : i + offsety+1, j-offsetx: j+offsetx+1]
            new_val = np.sum( filt *  box_vals)
            img_filtered[i][j] = np.sum(new_val)
    return img_filtered

A simple presentation of how function works for 5x5 input image and 3x3 filter kernel:

With having following 1d gaussian and its transpose, I call myconv2 function :

sigma = 3
filter_length = 5

gauss = gauss1d(sigma, filter_length).reshape(1,filter_length)
guass
array([[0.18073067, 0.20897821, 0.22058223, 0.20897821, 0.18073067]]) 

gauss_t = np.transpose(gauss)
gauss_t 
array([[0.18073067],
       [0.20897821],
       [0.22058223],
       [0.20897821],
       [0.18073067]])

myconv2(gauss, guass_t)
array([[0.        , 0.        , 0.        , 0.        , 0.        ],
       [0.        , 0.        , 0.        , 0.        , 0.        ],
       [0.03986597, 0.04609688, 0.04865652, 0.04609688, 0.03986597],
       [0.        , 0.        , 0.        , 0.        , 0.        ],
       [0.        , 0.        , 0.        , 0.        , 0.        ]])

As you can see its not actually a 2d gaussian kernel and some values are missing. I don't know what I am missing and what should I consider in my code to reach the goal. Thanks.

Answer 1

You could just do a matrix multiplication. The convolution should also work, just beware of the padding.

gaus2d = gauss.T @ gauss

Your conv2d implementation does not seem to be right. I suggest you to implement a 'valid' convolution (or cross correlation):

simple_valid_cross_correlation(img, filt):
    ih, iw = img.shape
    fh, fw = filt.shape
    result = np.zeros((ih - fh + 1, iw - fw + 1))
    for i in range(result.shape[0]):
        for j in range(result.shape[1]):
            result[i, j] = np.sum(filt * img[i:i+fh, j:j+fw])
    return result

gauss_pad = np.pad(gauss.T, ((0, 0), (gauss.shape[1]-1, gauss.shape[1]-1)))
gauss2d = simple_valid_cross_correlation(gauss_pad, gauss)

There is also scipy.signal.convolve2d if you don't want to implement your own conv. I think it may be faster