Convolve2d just by using Numpy

I am studying image-processing using Numpy and facing a problem with filtering with convolution.

I would like to convolve a gray-scale image. (convolve a 2d Array with a smaller 2d Array)

Does anyone have an idea to refine my method ?

I know that scipy supports convolve2d but I want to make a convolve2d only by using Numpy.

What I have done

First, I made a 2d array the submatrices.

a = np.arange(25).reshape(5,5) # original matrix

submatrices = np.array([
     [a[:-2,:-2], a[:-2,1:-1], a[:-2,2:]],
     [a[1:-1,:-2], a[1:-1,1:-1], a[1:-1,2:]],
     [a[2:,:-2], a[2:,1:-1], a[2:,2:]]])

the submatrices seems complicated but what I am doing is shown in the following drawing.

Next, I multiplied each submatrices with a filter.

conv_filter = np.array([[0,-1,0],[-1,4,-1],[0,-1,0]])
multiplied_subs = np.einsum('ij,ijkl->ijkl',conv_filter,submatrices)

and summed them.

np.sum(np.sum(multiplied_subs, axis = -3), axis = -3)
#array([[ 6,  7,  8],
#       [11, 12, 13],
#       [16, 17, 18]])

Thus this procudure can be called my convolve2d.

def my_convolve2d(a, conv_filter):
    submatrices = np.array([
         [a[:-2,:-2], a[:-2,1:-1], a[:-2,2:]],
         [a[1:-1,:-2], a[1:-1,1:-1], a[1:-1,2:]],
         [a[2:,:-2], a[2:,1:-1], a[2:,2:]]])
    multiplied_subs = np.einsum('ij,ijkl->ijkl',conv_filter,submatrices)
    return np.sum(np.sum(multiplied_subs, axis = -3), axis = -3)

However, I find this my_convolve2d troublesome for 3 reasons.

1. Generation of the submatrices is too awkward that is difficult to read and can only be used when the filter is 3*3
2. The size of the varient submatrices seems to be too big, since it is approximately 9 folds bigger than the original matrix.
3. The summing seems a little non intuitive. Simply said, ugly.

Thank you for reading this far.

Kind of update. I wrote a conv3d for myself. I will leave this as a public domain.

def convolve3d(img, kernel):
    # calc the size of the array of submatracies
    sub_shape = tuple(np.subtract(img.shape, kernel.shape) + 1)

    # alias for the function
    strd = np.lib.stride_tricks.as_strided

    # make an array of submatracies
    submatrices = strd(img,kernel.shape + sub_shape,img.strides * 2)

    # sum the submatraces and kernel
    convolved_matrix = np.einsum('hij,hijklm->klm', kernel, submatrices)

    return convolved_matrix

Answers:

Thank you for visiting the Q&A section on Magenaut. Please note that all the answers may not help you solve the issue immediately. So please treat them as advisements. If you found the post helpful (or not), leave a comment & I’ll get back to you as soon as possible.

Method 1

You could generate the subarrays using as_strided:

import numpy as np

a = np.array([[ 0,  1,  2,  3,  4],
       [ 5,  6,  7,  8,  9],
       [10, 11, 12, 13, 14],
       [15, 16, 17, 18, 19],
       [20, 21, 22, 23, 24]])

sub_shape = (3,3)
view_shape = tuple(np.subtract(a.shape, sub_shape) + 1) + sub_shape
strides = a.strides + a.strides

sub_matrices = np.lib.stride_tricks.as_strided(a,view_shape,strides)

To get rid of your second “ugly” sum, alter your einsum so that the output array only has j and k. This implies your second summation.

conv_filter = np.array([[0,-1,0],[-1,5,-1],[0,-1,0]])
m = np.einsum('ij,ijkl->kl',conv_filter,sub_matrices)

# [[ 6  7  8]
#  [11 12 13]
#  [16 17 18]]

Method 2

Cleaned up using as_strided and @Crispin ‘s einsum trick from above. Enforces the filter size into the expanded shape. Should even allow non-square inputs if the indices are compatible.

def conv2d(a, f):
    s = f.shape + tuple(np.subtract(a.shape, f.shape) + 1)
    strd = numpy.lib.stride_tricks.as_strided
    subM = strd(a, shape = s, strides = a.strides * 2)
    return np.einsum('ij,ijkl->kl', f, subM)

Method 3

You can also use fft (one of the faster methods to perform convolutions)

from numpy.fft import fft2, ifft2
import numpy as np

def fft_convolve2d(x,y):
    """ 2D convolution, using FFT"""
    fr = fft2(x)
    fr2 = fft2(np.flipud(np.fliplr(y)))
    m,n = fr.shape
    cc = np.real(ifft2(fr*fr2))
    cc = np.roll(cc, -m/2+1,axis=0)
    cc = np.roll(cc, -n/2+1,axis=1)
    return cc

• https://gist.github.com/thearn/5424195
• you must pad the filter to be the same size as image ( place it in the middle of a zeros_like mat.)

cheers,
Dan

Method 4

https://laurentperrinet.github.io/sciblog/posts/2017-09-20-the-fastest-2d-convolution-in-the-world.html

Check out all the convolution methods and their respective performances here.
Also, I found the below code snippet to be simpler.

from numpy.fft  import fft2, ifft2
def np_fftconvolve(A, B):
    return np.real(ifft2(fft2(A)*fft2(B, s=A.shape)))

All methods was sourced from stackoverflow.com or stackexchange.com, is licensed under cc by-sa 2.5, cc by-sa 3.0 and cc by-sa 4.0