numpy中2D的卷积卷积
Flo*_*oor 5 python numpy convolution
我尝试使用for循环实现2D数组的跨步卷积
arr = np.array([[2,3,7,4,6,2,9],
[6,6,9,8,7,4,3],
[3,4,8,3,8,9,7],
[7,8,3,6,6,3,4],
[4,2,1,8,3,4,6],
[3,2,4,1,9,8,3],
[0,1,3,9,2,1,4]])
arr2 = np.array([[3,4,4],
[1,0,2],
[-1,0,3]])
def stride_conv(arr1,arr2,s,p):
beg = 0
end = arr2.shape[0]
final = []
for i in range(0,arr1.shape[0]-1,s):
k = []
for j in range(0,arr1.shape[0]-1,s):
k.append(np.sum(arr1[beg+i : end+i, beg+j:end+j] * (arr2)))
final.append(k)
return np.array(final)
stride_conv(arr,arr2,2,0)
这导致3*3阵列:
array([[ 91, 100, 88],
[ 69, 91, 117],
[ 44, 72, 74]])
是否有一个numpy函数或scipy函数来做同样的事情？我的做法并不那么好.我该如何对此进行矢量化？
我认为我们可以做一个“有效”的 fft 卷积并只在跨步位置挑选那些结果，如下所示：
def strideConv(arr,arr2,s):
cc=scipy.signal.fftconvolve(arr,arr2[::-1,::-1],mode='valid')
idx=(np.arange(0,cc.shape[1],s), np.arange(0,cc.shape[0],s))
xidx,yidx=np.meshgrid(*idx)
return cc[yidx,xidx]
这给出了与其他人的答案相同的结果。但我想这仅在内核大小为奇数时才有效。
此外，arr2[::-1,::-1]为了与其他人保持一致，我还翻转了内核，您可能希望根据上下文省略它。
更新：
我们目前有几种不同的方法可以单独使用 numpy 和 scipy 进行 2D 或 3D 卷积，我考虑进行一些比较，以了解哪种方法在不同大小的数据上更快。我希望这不会被视为题外话。
方法一：FFT卷积（使用scipy.signal.fftconvolve）：
def padArray(var,pad,method=1):
if method==1:
var_pad=numpy.zeros(tuple(2*pad+numpy.array(var.shape[:2]))+var.shape[2:])
var_pad[pad:-pad,pad:-pad]=var
else:
var_pad=numpy.pad(var,([pad,pad],[pad,pad])+([0,0],)*(numpy.ndim(var)-2),
mode='constant',constant_values=0)
return var_pad
def conv3D(var,kernel,stride=1,pad=0,pad_method=1):
'''3D convolution using scipy.signal.convolve.
'''
var_ndim=numpy.ndim(var)
kernel_ndim=numpy.ndim(kernel)
stride=int(stride)
if var_ndim<2 or var_ndim>3 or kernel_ndim<2 or kernel_ndim>3:
raise Exception("<var> and <kernel> dimension should be in 2 or 3.")
if var_ndim==2 and kernel_ndim==3:
raise Exception("<kernel> dimension > <var>.")
if var_ndim==3 and kernel_ndim==2:
kernel=numpy.repeat(kernel[:,:,None],var.shape[2],axis=2)
if pad>0:
var_pad=padArray(var,pad,pad_method)
else:
var_pad=var
conv=fftconvolve(var_pad,kernel,mode='valid')
if stride>1:
conv=conv[::stride,::stride,...]
return conv
方法 2：特殊转换（请参阅此 anwser）：
def conv3D2(var,kernel,stride=1,pad=0):
'''3D convolution by sub-matrix summing.
'''
var_ndim=numpy.ndim(var)
ny,nx=var.shape[:2]
ky,kx=kernel.shape[:2]
result=0
if pad>0:
var_pad=padArray(var,pad,1)
else:
var_pad=var
for ii in range(ky*kx):
yi,xi=divmod(ii,kx)
slabii=var_pad[yi:2*pad+ny-ky+yi+1:1, xi:2*pad+nx-kx+xi+1:1,...]*kernel[yi,xi]
if var_ndim==3:
slabii=slabii.sum(axis=-1)
result+=slabii
if stride>1:
result=result[::stride,::stride,...]
return result
方法 3：跨步视图转换，如 Divakar 所建议的：
def asStride(arr,sub_shape,stride):
'''Get a strided sub-matrices view of an ndarray.
<arr>: ndarray of rank 2.
<sub_shape>: tuple of length 2, window size: (ny, nx).
<stride>: int, stride of windows.
Return <subs>: strided window view.
See also skimage.util.shape.view_as_windows()
'''
s0,s1=arr.strides[:2]
m1,n1=arr.shape[:2]
m2,n2=sub_shape[:2]
view_shape=(1+(m1-m2)//stride,1+(n1-n2)//stride,m2,n2)+arr.shape[2:]
strides=(stride*s0,stride*s1,s0,s1)+arr.strides[2:]
subs=numpy.lib.stride_tricks.as_strided(arr,view_shape,strides=strides)
return subs
def conv3D3(var,kernel,stride=1,pad=0):
'''3D convolution by strided view.
'''
var_ndim=numpy.ndim(var)
kernel_ndim=numpy.ndim(kernel)
if var_ndim<2 or var_ndim>3 or kernel_ndim<2 or kernel_ndim>3:
raise Exception("<var> and <kernel> dimension should be in 2 or 3.")
if var_ndim==2 and kernel_ndim==3:
raise Exception("<kernel> dimension > <var>.")
if var_ndim==3 and kernel_ndim==2:
kernel=numpy.repeat(kernel[:,:,None],var.shape[2],axis=2)
if pad>0:
var_pad=padArray(var,pad,1)
else:
var_pad=var
view=asStride(var_pad,kernel.shape,stride)
#return numpy.tensordot(aa,kernel,axes=((2,3),(0,1)))
if numpy.ndim(kernel)==2:
conv=numpy.sum(view*kernel,axis=(2,3))
else:
conv=numpy.sum(view*kernel,axis=(2,3,4))
return conv
我做了3组比较：
- 对二维数据进行卷积，输入大小不同，核大小不同，stride=1，pad=0。结果如下（颜色为用于卷积的时间重复 10 次）：
所以“FFT conv”通常是最快的。“Special conv”和“Stride-view conv”随着内核大小的增加而变慢，但随着接近输入数据的大小而再次减小。最后一个子图显示了最快的方法，所以紫色的大三角形表示 FFT 是赢家，但请注意左侧有一个细的绿色列（可能太小看不到，但它在那里），表明“特殊转换”对于非常小的内核（小于大约 5x5）具有优势。当内核大小接近输入时，“stride-view conv”是最快的（见对角线）。
比较 2：对 3D 数据进行卷积。
设置：pad=0，stride=2，输入维度= nxnx5，核形状= fxfx5。
当内核大小在输入的中间时，我跳过了“Special Conv”和“Stride-view conv”的计算。基本上“Special Conv”现在没有任何优势，“Stride-view”对于小内核和大内核都比 FFT 快。
另一个注意事项：当大小超过 350 时，我注意到“跨步视图转换”的内存使用量达到了可观的峰值。
比较 3：对 3D 数据进行更大步幅的卷积。
设置：pad=0，stride=5，输入维度= nxnx10，核形状= fxfx10。
这次我省略了“特别会议”。对于更大的区域，“Stride-view conv”超过了 FFT，最后的子图显示差异接近 100%。可能是因为随着步幅的增加，FFT 方法会浪费更多的数字，因此“步幅视图”对于小内核和大内核都具有更多优势。
使用signal.convolve2dfromscipy怎么样？
我的方法类似于 Jason 的方法，但使用索引。
def strideConv(arr, arr2, s):
return signal.convolve2d(arr, arr2[::-1, ::-1], mode='valid')[::s, ::s]
请注意，内核必须反转。有关详细信息，请参阅此处和此处的讨论。否则使用signal.correlate2d.
例子：
>>> strideConv(arr, arr2, 1)
array([[ 91, 80, 100, 84, 88],
[ 99, 106, 126, 92, 77],
[ 69, 98, 91, 93, 117],
[ 80, 79, 87, 93, 61],
[ 44, 72, 72, 63, 74]])
>>> strideConv(arr, arr2, 2)
array([[ 91, 100, 88],
[ 69, 91, 117],
[ 44, 72, 74]])
- 有没有办法在不浪费计算的情况下做到这一点？就像这个例子一样，在所需的步幅之间完成了一堆无用的计算 (5认同)
忽略padding参数和尾部窗口没有足够的长度来卷积第二个数组,这里有一种方法np.lib.stride_tricks.as_strided-
def strided4D(arr,arr2,s):
strided = np.lib.stride_tricks.as_strided
s0,s1 = arr.strides
m1,n1 = arr.shape
m2,n2 = arr2.shape
out_shp = (1+(m1-m2)//s, m2, 1+(n1-n2)//s, n2)
return strided(arr, shape=out_shp, strides=(s*s0,s*s1,s0,s1))
def stride_conv_strided(arr,arr2,s):
arr4D = strided4D(arr,arr2,s=s)
return np.tensordot(arr4D, arr2, axes=((2,3),(0,1)))
或者,我们可以使用内置的scikit-image view_as_windows来优雅地获得这些窗口,就像这样 -
from skimage.util.shape import view_as_windows
def strided4D_v2(arr,arr2,s):
return view_as_windows(arr, arr2.shape, step=s)
- @Divakar 但是您的代码不适用于除 `arr=7x7` 以外的形状（这恰好给出了 `1+(m1-m2)//s = 3 = 内核大小`），有时 `arr4D` 甚至会有 ` nan在里面。更改为 `(1+(m1-m2)//s, 1+(n1-n2)//s, m2, n2)` 会正确匹配张量点轴 (`(2,3),(0,1 )`)。我认为@NaN 也想询问这个问题，因为他的 `np.sum()` 方法现在也有效。 (4认同)