使用NumPy进行二维卷积

Question

使用NumPy进行二维卷积

pythonnumpyimage-processingmatrixconvolution

50

我正在使用NumPy学习图像处理，并在用卷积进行滤波时遇到问题。

我想对灰度图像进行卷积。 (将一个二维数组与一个较小的二维数组进行卷积)

有没有人有办法改进我的方法？

我知道SciPy支持convolve2d，但我想只使用NumPy来实现卷积。

我已经做了什么

首先，我制作了一个由子矩阵组成的二维数组。

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:]]])

子矩阵看起来很复杂，但我正在做的内容如下图所示。

接下来，我将每个子矩阵与一个滤波器相乘。

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

并将它们相加。

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

因此，这个过程可以被称为我的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)

然而，我发现我的my_convolve2d有3个令人困扰的问题:

生成子矩阵太棘手了，难以阅读，只能在过滤器为3*3时使用。
变量子矩阵的大小似乎太大了，因为它大约比原始矩阵大9倍。
求和似乎有点不直观。简单来说，很丑陋。

感谢您阅读到这里。

更新一下，我为自己编写了一个conv3d。我将把它留作公共领域。

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

- Allosteric

2

谢谢您提供矩阵的图纸 :) 如果我理解正确，您想要一些关于如何使您的解决方案更加优雅的建议？ - Marijn van Vliet

很高兴能帮到你！是的，如果您能提供我克服最后几行中写的三个问题的提示，我将不胜感激。 - Allosteric

我应该补充说明这3个要点是按优先顺序排列的。第一个对我来说非常重要，而最后一个似乎有点琐碎。如果还有其他问题和改进意见，我也会很高兴听取。 - Allosteric

2

第二个图（等号后面）不对吗？每个子矩阵不应该与过滤器逐元素相乘，然后将每个结果子矩阵的元素相加吗？ - Andreas K.

1

@AndyK 它们将产生相同的结果。 - Allosteric

4个回答

19

使用上面提到的as_strided和@Crispin的einsum技巧进行了清理。强制将滤波器大小应用于扩展形状。如果索引兼容，应该甚至允许非方形输入。

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

请进一步简化...请看我的评论...如果您的答案确实是array([[6, 7, 8], [11, 12, 13], [16, 17, 18]])，则np.sum(a_s * filter，axis=(2,3))...其中a_s是分步数组，filter是3x3过滤器。 - NaN

嗨@DanielF，是否有适用于RGB的概括方法？对于einsum符号不是特别熟悉，因此如何推广的想法将是很棒的。 - sirgogo

@sirgogo，我建议您提出一个单独的问题，并链接到这个答案。有些人在numpy上对多维卷积的掌握比我更好。 - Daniel F

@DanielF，谢谢，这是我尝试做的更一般化的版本：[https://stackoverflow.com/questions/50239641/multi-dimensional-batch-image-convolution-using-numpy] - sirgogo

这里内核是否翻转了，还是没有翻转？ - OuttaSpaceTime

显示剩余3条评论

17

你还可以使用fft（一种更快的卷积方法之一）

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
你必须将滤波器填充为与图像相同的大小（放置在一个零矩阵中心）。

干杯，丹

- Dan Erez

如果内核只有9个元素，使用FFT进行卷积肯定不是高效的。这仅适用于大内核。 - Cris Luengo

真的。这取决于内核和图像大小，但fft优于其他算法的门槛非常低。这也取决于您的计算机架构，但对于大多数通用用例来说，fft是高性能库中的首选。 - Dan Erez

1

以下文章比较了各种图像卷积方法的性能，因此我们可以得出结论：Numpy的FFT方法是执行卷积最快的方法。https://laurentperrinet.github.io/sciblog/posts/2017-09-20-the-fastest-2d-convolution-in-the-world.html - Naman Bansal

关于填充：只需在fr2中添加s=x.shape，这样做就可以了吗？就像这样fr2 = fft2(np.flipud(np.fliplr(y)), s=x.shape) # 填充滤波器 - undefined

0

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

在这里查看所有卷积方法及其相应的性能。此外，我发现下面的代码片段更简单。

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

- Naman Bansal

2

这段代码不正确。它会使输出图像发生偏移。核越大，偏移量就越大。 - Cris Luengo

你能详细说明一下它将如何移动输出图像吗？ - Naman Bansal

https://dev59.com/xrLma4cB1Zd3GeqPVx6f#54977551 - Cris Luengo

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Crispin · Accepted Answer

您可以使用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)

为了消除第二个“ugly” sum，修改您的einsum，使输出数组仅包含j和k。这意味着您第二次求和。

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]]