如何实现NymPy卷积的多维版本? [英] How to implement a multidimensional version of NymPy convolve?
问题描述
此帖子建议我根据 MCVE 提出另一个问题.我的目标是实现 NumPy的卷积形状的输入数组.请考虑我在多变量幂级数的 Cauchy产品的上下文中使用了卷积一词(多元多项式乘法).我已经解释此处.
Following this post I was advised to ask a different question based on MCVE. My objective is to implement the NumPy's convolve for arbitrary shaped input arrays. Please consider that I'm using the term convolution in the context of Cauchy product of multivariate power series (multivariable polynomial multiplication). SciPy functions such as signal.convolve, ndimage.convolve or ndimage.filters.convolve do not work for me as I have explained here.
考虑两个非正方形2D NumPy数组A和B:
Consider two non-square 2D NumPy arrays A and B:
D1=np.array([4,5])
D2=np.array([2,3])
A=np.random.randint(10,size=D1)
B=np.random.randint(10,size=D2)
例如:
[[1 4 4 2 7]
[6 1 7 5 3]
[1 4 3 4 8]
[7 5 8 3 3]]
[[2 2 3]
[5 2 9]]
现在我可以用conv(A,B,K)计算C=convolve(A,B)的元素了:
Now I'm able to calculate the elements of the C=convolve(A,B) with conv(A,B,K):
def crop(A,D1,D2):
return A[tuple(slice(D1[i], D2[i]) for i in range(A.ndim))]
def sumall(A):
sum1=A
for k in range(A.ndim):
sum1 = np.sum(sum1,axis=0)
return sum1
def flipall(A):
return A[[slice(None, None, -1)] * A.ndim]
def conv(A,B,K):
D0=np.zeros(K.shape,dtype=K.dtype)
return sumall(np.multiply(crop(A,np.maximum(D0,np.minimum(A.shape,K-B.shape)) \
,np.minimum(A.shape,K)), \
flipall(crop(B,np.maximum(D0,np.minimum(B.shape,K-A.shape)) \
,np.minimum(B.shape,K)))))
K=np.array([0,0])+1，conve(A,B,K)结果1*2=2，K=np.array([1,0])+1结果5*1+2*6=17，K=np.array([0,1])+1为2*4+1*2=10，K=np.array([1,1])+1给出4*5+6*2+1*1+1*2=36的示例如下:
Fow example for K=np.array([0,0])+1, conve(A,B,K) results 1*2=2, for K=np.array([1,0])+1 results 5*1+2*6=17, for K=np.array([0,1])+1 is 2*4+1*2=10 and for K=np.array([1,1])+1 gives 4*5+6*2+1*1+1*2=36:
[[2 10 ...]
[17 36 ...]
... ]]
现在，如果我知道A和B的尺寸，则可以嵌套一些for循环以填充C，但对我而言并非如此.如何使用conv函数以形状为C.shape=A.shape+B.shape-1的C ndarray填充而不使用for循环?
now if I knew the dimension of A and B I could nest some for loops to populate the C, but that's not the case for me. How can I use the conv function to populate the C ndarray with a shape of C.shape=A.shape+B.shape-1 without using for loops?
推荐答案
在scipy.ndimage和astropy中存在大量n维卷积函数.让我们看看是否可以使用它们中的任何一个.
There exist a multitude of n-dimensional convolution functions in scipy.ndimage and astropy. Let's see if we can use any of them.
首先，我们需要一些数据进行比较.因此，我们来扩展输入空间:
First we need some data to compare against. So let's span up the input space:
d0, d1 = np.array(A.shape) + np.array(B.shape) - 1
input_space = np.array(np.meshgrid(np.arange(d0), np.arange(d1))).T.reshape(-1, 2)
# array([[0, 0],
# [0, 1],
# [0, 2],
# [0, 3],
# [0, 4],
# [0, 5],
# [0, 6],
# [1, 0],
# [1, 1],
# ...
# [4, 5],
# [4, 6]])
并计算在该空间上的卷积:
and calculate your convolution over this space:
out = np.zeros((d0, d1))
for K in input_space:
out[tuple(K)] = conv(A, B, K + 1)
out
# array([[ 2., 10., 19., 24., 30., 20., 21.],
# [ 17., 36., 71., 81., 112., 53., 72.],
# [ 32., 27., 108., 74., 121., 79., 51.],
# [ 19., 46., 79., 99., 111., 67., 81.],
# [ 35., 39., 113., 76., 93., 33., 27.]])
好吧，既然我们知道期望值，那么让我们看看是否可以让scipy和astropy赋予我们相同的值:
Okay, now that we know what values to expect, lets see if we can get scipy and astropy to give us the same values:
import scipy.signal
scipy.signal.convolve2d(A, B) # only 2D!
# array([[ 2., 10., 19., 24., 30., 20., 21.],
# [ 17., 36., 71., 81., 112., 53., 72.],
# [ 32., 27., 108., 74., 121., 79., 51.],
# [ 19., 46., 79., 99., 111., 67., 81.],
# [ 35., 39., 113., 76., 93., 33., 27.]])
import astropy.convolution
astropy.convolution.convolve_fft(
np.pad(A, pad_width=((1, 0), (1, 1)), mode='constant'),
B,
normalize_kernel=False
)
# array([[ 2., 10., 19., 24., 30., 20., 21.],
# [ 17., 36., 71., 81., 112., 53., 72.],
# [ 32., 27., 108., 74., 121., 79., 51.],
# [ 19., 46., 79., 99., 111., 67., 81.],
# [ 35., 39., 113., 76., 93., 33., 27.]])
astropy.convolution.convolve(
np.pad(A, pad_width=((1, 0), (1, 1)), mode='constant'),
np.pad(B, pad_width=((0, 1), (0, 0)), mode='constant'),
normalize_kernel=False
)
# array([[ 2., 10., 19., 24., 30., 20., 21.],
# [ 17., 36., 71., 81., 112., 53., 72.],
# [ 32., 27., 108., 74., 121., 79., 51.],
# [ 19., 46., 79., 99., 111., 67., 81.],
# [ 35., 39., 113., 76., 93., 33., 27.]])
import scipy
scipy.ndimage.filters.convolve(
np.pad(A, pad_width=((0, 1), (0, 2)), mode='constant'),
B,
mode='constant',
cval=0.0,
origin=-1
)
# array([[ 2., 10., 19., 24., 30., 20., 21.],
# [ 17., 36., 71., 81., 112., 53., 72.],
# [ 32., 27., 108., 74., 121., 79., 51.],
# [ 19., 46., 79., 99., 111., 67., 81.],
# [ 35., 39., 113., 76., 93., 33., 27.]])
scipy.ndimage.filters.convolve(
np.pad(A, pad_width=((1, 0), (1, 1)), mode='constant'),
B,
mode='constant',
cval=0.0
)
# array([[ 2., 10., 19., 24., 30., 20., 21.],
# [ 17., 36., 71., 81., 112., 53., 72.],
# [ 32., 27., 108., 74., 121., 79., 51.],
# [ 19., 46., 79., 99., 111., 67., 81.],
# [ 35., 39., 113., 76., 93., 33., 27.]])
如您所见，这只是选择正确的规范化和填充的问题，您可以简单地使用这些库中的任何一个.
As you can see, it is just a matter of chosing the right normalization and padding, and you can simply use any of these libraries.
我建议使用astropy.convolution.convolve_fft，因为它(基于FFT)可能是最快的.
I recommend using astropy.convolution.convolve_fft, as it (being FFT-based) is probably the fastest.
这篇关于如何实现NymPy卷积的多维版本?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
