48

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.

submatrices

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)

multiplied_subs

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 procedure 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 variant 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 submatrices
    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 submatrices
    submatrices = strd(img,kernel.shape + sub_shape,img.strides * 2)

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

    return convolved_matrix
greybeard
  • 2,249
  • 8
  • 30
  • 66
Allosteric
  • 871
  • 1
  • 9
  • 14
  • 2
    thank you for providing drawings of the matrices :) If I understand correctly, you want tips on how to make your solution more elegant? – Marijn van Vliet Mar 29 '17 at 07:24
  • Glad it helps! Yes. I would be grateful if you can provide me tips to overcome the 3 problems written in the very last lines. – Allosteric Mar 29 '17 at 07:26
  • I should add that the 3 points are rather arranged in a priority order. The first one is quite important for me and the last one seems kinda trivial. I will also be glad if there are other problems and refinements about it. – Allosteric Mar 29 '17 at 07:28
  • 2
    Isn't the second drawing (after the equality sign) wrong? Shouldn't each submatrix be multiplied (element-wise) with the filter, and then the elements of each of the resulting submatrices summed? – Andreas K. Nov 09 '18 at 13:34
  • 1
    @AndyK They will produce the same result. – Allosteric Nov 10 '18 at 11:38

4 Answers4

36

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]]
Gulzar
  • 23,452
  • 27
  • 113
  • 201
Crispin
  • 2,070
  • 1
  • 13
  • 16
  • if a_s is the strided array and filter is your laplacian like filter, then try... np.sum(a_s*filter, axis=(2,3)) if indeed your answer is array([[ 6, 7, 8], [11, 12, 13], [16, 17, 18]]) – NaN Mar 29 '17 at 08:03
  • Thank you for the tip. I am trying this myself right now. Maybe trivial but I believe that the name filter is not appropriate because it is a built-in function of python. – Allosteric Mar 29 '17 at 08:14
  • 1
    You can do the summation directly in the Einstein summation. See answer – Crispin Mar 29 '17 at 08:16
  • In this question, isn't it better to write, `sub_shape = conv_filter.shape` ? – Allosteric Mar 29 '17 at 08:22
  • It doesn't seem that it would necessarily be the case. `conv_filter.shape` should be `a.shape - sub_shape + 1`. It just so happens in your example `a.shape = [5, 5]` and `sub_shape = (3, 3)` so `conv_filter.shape == sub_shape` is true but not strictly enforced. If you *want* it to be enforced, then you should set `sub_shape = a.shape - conv_filter.shape + 1` – Daniel F Mar 29 '17 at 08:51
  • Got it! So `view_shape = conv_filter.shape`? – Allosteric Mar 29 '17 at 09:00
  • well `view_shape` should be a 4D tuple, so `view.shape = tuple(conv_filter.shape) + tuple(sub_shape)`, so you still need to find `sub_shape` – Daniel F Mar 29 '17 at 09:04
  • Follow-up question here: https://stackoverflow.com/questions/45580013/pure-numpy-2d-mean-convolution-derivative-of-input-image – michaelsnowden Aug 09 '17 at 03:19
  • I have a similar issue for higher dimensions. I was wondering if you could be so kind to take a look at my question [here](https://stackoverflow.com/questions/51794274/convolution-of-numpy-arrays-of-arbitrary-dimension-for-cauchy-product-of-multiva) – Foad S. Farimani Aug 13 '18 at 08:53
  • 1
    Is using `einsum("ij,klij", ...)` sufficient/better? Using the code from this answer with input of different shape, I got `operands could not be broadcast together with remapped shapes`, using `"ij,klij"` seems to work without issues. – Michael Kopp Dec 20 '21 at 21:16
19

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)
Daniel F
  • 13,620
  • 2
  • 29
  • 55
  • simplify further... see my comment below... np.sum(a_s * filter, axis=(2,3)) if indeed your answer is array([[ 6, 7, 8], [11, 12, 13], [16, 17, 18]]) ... where a_s is strided array and filter is the 3x3 filter – NaN Mar 29 '17 at 08:07
  • Not sure why that works @NaN since it is certainly *not* doing what the problem asked - but it does, even for arbitrary `a` matrices – Daniel F Mar 29 '17 at 08:14
  • at least at numpy v12, a.shape and f.shape is a tuple so `s` should be `tuple(np.subtract(a.shape, f.shape)+1)`, I think. – Allosteric Mar 29 '17 at 09:06
  • True! Still used to `shape` giving an `ndarray`. Fixed that. – Daniel F Mar 29 '17 at 09:09
  • Hi @DanielF, is there a generalization of this that works for RGB? Not super familiar with einsum notation, so an idea how to generalize would be awesome. – sirgogo May 06 '18 at 06:14
  • @sirgogo, I suggest asking a seperate question linking this answer. There are people watching `numpy` way better at multidimensional convlcutions than I am. – Daniel F May 07 '18 at 11:03
  • @DanielF, thanks, here is the more general version of what I'm trying to do: [https://stackoverflow.com/questions/50239641/multi-dimensional-batch-image-convolution-using-numpy] – sirgogo May 08 '18 at 18:03
  • is the kernel flipped here or is it without flipping? – OuttaSpaceTime Jan 29 '22 at 13:23
17

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

cheers, Dan

Dan Erez
  • 1,364
  • 15
  • 16
  • It is certainly not efficient to use the FFT for convolution if the kernel has only 9 elements. This is only useful for large kernels. – Cris Luengo Feb 25 '20 at 01:43
  • True. It depends on the kernel and image size but the threshold for fft to outperform is very low. It also depends on your computer architecture but for most generic use-cases fft is the go to in high performance libraries. – Dan Erez Feb 26 '20 at 17:25
  • 1
    The following article compares the performance of various image convolution methods, and thus, we can conclude that the FFT method of Numpy is the fastest method to perform convolution. https://laurentperrinet.github.io/sciblog/posts/2017-09-20-the-fastest-2d-convolution-in-the-world.html – Naman Bansal Jan 30 '21 at 20:24
0

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)))
Naman Bansal
  • 191
  • 2
  • 13