2D Convolution in Python similar to Matlab's conv2

Question

I have been trying to do Convolution of a 2D Matrix using SciPy, and Numpy but have failed. For SciPy I tried, sepfir2d and scipy.signal.convolve and Convolve2D for Numpy. Is there a simple function like conv2 in Matlab for Python?

Here is an example:

 A= [ 5     4     5     4;
      3     2     3     2;
      5     4     5     4;
      3     2     3     2 ]

I want to convolve it with [0.707 0.707]

And the result as by conv2 from Matlab is

3.5350    6.3630    6.3630    6.3630    2.8280
2.1210    3.5350    3.5350    3.5350    1.4140
3.5350    6.3630    6.3630    6.3630    2.8280
2.1210    3.5350    3.5350    3.5350    1.4140

Some function to compute this output in Python? I will be grateful for a response.

Solution 1

There are a number of different ways to do it with scipy, but 2D convolution isn't directly included in numpy. (It's also easy to implement with an fft using only numpy, if you need to avoid a scipy dependency.)

scipy.signal.convolve2d, scipy.signal.convolve, scipy.signal.fftconvolve, and scipy.ndimage.convolve will all handle a 2D convolution (the last three are N-d) in different ways.

scipy.signal.fftconvolve does the convolution in the fft domain (where it's a simple multiplication). This is much faster in many cases, but can lead to very small differences in edge effects than the discrete case, and your data will be coerced into floating point with this particular implementation. Additionally, there's unnecessary memory usage when convolving a small array with a much larger array. All in all, fft-based methods can be dramatically faster, but there are some common use cases where scipy.signal.fftconvolve is not an ideal solution.

scipy.signal.convolve2d, scipy.signal.convolve, and scipy.ndimage.convolve all use a discrete convolution implemented in C, however, they implement it in different ways.

scipy.ndimage.convolve keeps the same data type, and gives you control over the location of the output to minimize memory usage. If you're convolving uint8's (e.g. image data), it's often the best option. The output will always be the same shape as the first input array, which makes sense for images, but perhaps not for more general convolution. ndimage.convolve gives you a lot of control over how edge effects are handled through the mode kwarg (which functions completely differently than scipy.signal's mode kwarg).

Avoid scipy.signal.convolve if you're working with 2d arrays. It works for the N-d case, but it's suboptimal for 2d arrays, and scipy.signal.convolve2d exists to do the exact same thing a bit more efficiently. The convolution functions in scipy.signal give you control over the output shape using the mode kwarg. (By default, they'll behave just like matlab's conv2.) This is useful for general mathematical convolution, but less useful for image processing. However, scipy.signal.convolve2d is generally slower than scipy.ndimage.convolve.

There are a lot of different options partly due to duplication in the different submodules of scipy and partly because there are different ways to implement a convolution that have different performance tradeoffs.

If you can give a bit more detail about your use case, we can recommend a better solution. If you're convolving two arrays of roughly the same size, and they're already floats, fftconvolve is an excellent choice. Otherwise, scipy.ndimage.convolve may beat it.

Solution 2

scipy's convolved1d() does what you want, just treats the edges a bit differently:

sp.ndimage.filters.convolve1d(A,[0.707,0.707],axis=1,mode='constant')

will give you:

array([[ 6.363,  6.363,  6.363,  2.828],
       [ 3.535,  3.535,  3.535,  1.414],
       [ 6.363,  6.363,  6.363,  2.828],
       [ 3.535,  3.535,  3.535,  1.414]])

If you want the exact same result, just add a column of zeros to A like this:

sp.ndimage.filters.convolve1d(np.c_[np.zeros((4,1)),A],[0.707,0.707],axis=1,mode='constant')

and you will get:

array([[ 3.535,  6.363,  6.363,  6.363,  2.828],
       [ 2.121,  3.535,  3.535,  3.535,  1.414],
       [ 3.535,  6.363,  6.363,  6.363,  2.828],
       [ 2.121,  3.535,  3.535,  3.535,  1.414]])

From my experience you can do in scipy/numpy most of what you do in Matlab very easily (and more).

Comments

[user1343318]

Thanks a ton Bitwise, I am going to try it.

[user1343318]

Thank you Joe for a detailed response. I am trying to do Stationary Wavelet Transform of an Image to be used in denoising. That is where I want to use it.

[user1343318]

I can of course, and if I hadn't found a way, I must have. But I guess my implementation will not be optimized, so I am looking for a package.

[Joe Kington]

In that case, scipy.ndimage.convolve1d, as Bitwise mentioned is probably your best choice. It's optimized for the particular use case of convolving a 1d array with a 2d array. Hope that helps!

[Joe Kington]

Also, if you don't want to have to pad A with zeros, scipy.signal.convolve2d may be a better option. The shape of the array returned will be just like matlab's conv2, by default. (I forgot to include that information in my answer originally) If you're working with an image, though, I'd assume that you want the returned array to be the same shape as the original. In that case, ndimage is your best bet.

[Daniel Conde Marin]

@downvoter: Please, could you explain why the downvotes

[rayryeng]

@Daniel - Because you are suggesting a rather inefficient solution when there are tools already available to compute the convolution for you. Also, judging from the context of the OP's question, they want to use a package. I didn't downvote you, but if I was feeling in a bad way, I would.

[dtk]

Ftr: I was bitten when trying to use scipy.signal.convolve2d or scipy.ndimage.convolve because my weights weren't 2D... enough -.- Cf. np.array([0.707, 0.707]).shape vs np.array([[0.707, 0.707]]).shape. Stupid mistake when literally translating code.

[Robyc]

Note that while scipy.convolve2d only performs a "direct" convolution (i.e. in the original domain), scipy.convolve (N-dimensional) works either in the original domain or in the FFT-domain, depending on which is computationally better. With relatively big images the time saved thanks to the FFT can be consistent! By using scipy.convolve instead of scipy.convolve2d I got an average of 40x speed boost when convolving 3000x3000 with 41x41 arrays.