18

我正在尝试使用 numpy 在 python 中执行二维卷积

我有一个二维数组,如下所示,行的内核为 H_r,列的为 H_c

data = np.zeros((nr, nc), dtype=np.float32)

#fill array with some data here then convolve

for r in range(nr):
    data[r,:] = np.convolve(data[r,:], H_r, 'same')

for c in range(nc):
    data[:,c] = np.convolve(data[:,c], H_c, 'same')

data = data.astype(np.uint8);

它没有产生我期望的输出,这段代码看起来还可以吗,我认为问题在于从 float32 转换为 8bit。什么是最好的方法来做到这一点

谢谢

4

8 回答 8

21

也许它不是最优化的解决方案,但这是我之前在 Python 的 numpy 库中使用的一个实现:

def convolution2d(image, kernel, bias):
    m, n = kernel.shape
    if (m == n):
        y, x = image.shape
        y = y - m + 1
        x = x - m + 1
        new_image = np.zeros((y,x))
        for i in range(y):
            for j in range(x):
                new_image[i][j] = np.sum(image[i:i+m, j:j+m]*kernel) + bias
    return new_image

我希望这段代码可以帮助其他有同样疑问的人。

问候。

于 2017-03-03T12:47:24.433 回答
6

编辑 [2019 年 1 月]

@Tashus 下面的评论是正确的,因此@dudemeister 的回答可能更准确。他建议的函数也更有效,因为它避免了直接的 2D 卷积和所需的操作数量。

可能的问题

我相信你正在做两个一维卷积,第一个是每列,第二个是每行,然后用第二个的结果替换第一个的结果。

请注意,numpy.convolve使用'same'参数返回一个与提供的最大数组形状相同的数组,因此当您进行第一次卷积时,您已经填充了整个data数组。

在这些步骤中可视化数组的一种好方法是使用Hinton 图,这样您就可以检查哪些元素已经具有值。

可能的解决方案

如果您的卷积矩阵是使用一维矩阵的结果,您可以尝试添加两个卷积的结果(使用data[:,c] += ..而不是data[:,c] =在第二个循环上) :forH_rH_c

卷积核加法

另一种方法是使用scipy.signal.convolve2d2d 卷积数组,这可能是您最初想要做的。

于 2016-10-04T12:48:23.327 回答
5

由于您已经分离了内核,您应该简单地使用 scipy 中的 sepfir2d 函数:

from scipy.signal import sepfir2d
convolved = sepfir2d(data, H_r, H_c)

另一方面,您那里的代码看起来还不错...

于 2010-03-15T15:18:42.397 回答
1

我检查了许多实现,但没有找到适合我的目的,这应该很简单。所以这是一个非常简单的 for 循环实现

def convolution2d(image, kernel, stride, padding):
    image = np.pad(image, [(padding, padding), (padding, padding)], mode='constant', constant_values=0)

    kernel_height, kernel_width = kernel.shape
    padded_height, padded_width = image.shape

    output_height = (padded_height - kernel_height) // stride + 1
    output_width = (padded_width - kernel_width) // stride + 1

    new_image = np.zeros((output_height, output_width)).astype(np.float32)

    for y in range(0, output_height):
        for x in range(0, output_width):
            new_image[y][x] = np.sum(image[y * stride:y * stride + kernel_height, x * stride:x * stride + kernel_width] * kernel).astype(np.float32)
    return new_image
于 2020-12-12T21:20:30.327 回答
1

它也可能不是最优化的解决方案,但它比@omotto 提出的解决方案快大约十倍,并且它只使用基本的 numpy 函数(如 reshape、expand_dims、tile ......)并且没有“for”循环:

def gen_idx_conv1d(in_size, ker_size):
    """
    Generates a list of indices. This indices correspond to the indices
    of a 1D input tensor on which we would like to apply a 1D convolution.

    For instance, with a 1D input array of size 5 and a kernel of size 3, the
    1D convolution product will successively looks at elements of indices [0,1,2],
    [1,2,3] and [2,3,4] in the input array. In this case, the function idx_conv1d(5,3) 
    outputs the following array: array([0,1,2,1,2,3,2,3,4]).

    args:
        in_size: (type: int) size of the input 1d array.
        ker_size: (type: int) kernel size.

    return:
        idx_list: (type: np.array) list of the successive indices of the 1D input array
        access to the 1D convolution algorithm.

    example:
        >>> gen_idx_conv1d(in_size=5, ker_size=3)
        array([0, 1, 2, 1, 2, 3, 2, 3, 4])
    """
    f = lambda dim1, dim2, axis: np.reshape(np.tile(np.expand_dims(np.arange(dim1),axis),dim2),-1)
    out_size = in_size-ker_size+1
    return f(ker_size, out_size, 0)+f(out_size, ker_size, 1)

def repeat_idx_2d(idx_list, nbof_rep, axis):
    """
    Repeats an array of indices (idx_list) a number of time (nbof_rep) "along" an axis
    (axis). This function helps to browse through a 2d array of size
    (len(idx_list),nbof_rep).

    args:
        idx_list: (type: np.array or list) a 1D array of indices.
        nbof_rep: (type: int) number of repetition.
        axis: (type: int) axis "along" which the repetition will be applied.

    return
        idx_list: (type: np.array) a 1D array of indices of size len(idx_list)*nbof_rep.

    example:
        >>> a = np.array([0, 1, 2])
        >>> repeat_idx_2d(a, 3, 0) # repeats array 'a' 3 times along 'axis' 0
        array([0, 0, 0, 1, 1, 1, 2, 2, 2])

        >>> repeat_idx_2d(a, 3, 1) # repeats array 'a' 3 times along 'axis' 1
        array([0, 1, 2, 0, 1, 2, 0, 1, 2])

        >>> b = np.reshape(np.arange(3*4), (3,4))
        >>> b[repeat_idx_2d(np.arange(3), 4, 0), repeat_idx_2d(np.arange(4), 3, 1)]
        array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11])
    """
    assert axis in [0,1], "Axis should be equal to 0 or 1."
    tile_axis = (nbof_rep,1) if axis else (1,nbof_rep)
    return np.reshape(np.tile(np.expand_dims(idx_list, 1),tile_axis),-1)

def conv2d(im, ker):
    """
    Performs a 'valid' 2D convolution on an image. The input image may be
    a 2D or a 3D array.

    The output image first two dimensions will be reduced depending on the 
    convolution size. 

    The kernel may be a 2D or 3D array. If 2D, it will be applied on every
    channel of the input image. If 3D, its last dimension must match the
    image one.

    args:
        im: (type: np.array) image (2D or 3D).
        ker: (type: np.array) convolution kernel (2D or 3D).

    returns:
        im: (type: np.array) convolved image.

    example:
        >>> im = np.reshape(np.arange(10*10*3),(10,10,3))/(10*10*3) # 3D image
        >>> ker = np.array([[0,1,0],[-1,0,1],[0,-1,0]]) # 2D kernel
        >>> conv2d(im, ker) # 3D array of shape (8,8,3)
    """
    if len(im.shape)==2: # it the image is a 2D array, it is reshaped by expanding the last dimension
        im = np.expand_dims(im,-1)

    im_x, im_y, im_w = im.shape

    if len(ker.shape)==2: # if the kernel is a 2D array, it is reshaped so it will be applied to all of the image channels
        ker = np.tile(np.expand_dims(ker,-1),[1,1,im_w]) # the same kernel will be applied to all of the channels 

    assert ker.shape[-1]==im.shape[-1], "Kernel and image last dimension must match."

    ker_x = ker.shape[0]
    ker_y = ker.shape[1]

    # shape of the output image
    out_x = im_x - ker_x + 1 
    out_y = im_y - ker_y + 1

    # reshapes the image to (out_x, ker_x, out_y, ker_y, im_w)
    idx_list_x = gen_idx_conv1d(im_x, ker_x) # computes the indices of a 1D conv (cf. idx_conv1d doc)
    idx_list_y = gen_idx_conv1d(im_y, ker_y)

    idx_reshaped_x = repeat_idx_2d(idx_list_x, len(idx_list_y), 0) # repeats the previous indices to be used in 2D (cf. repeat_idx_2d doc)
    idx_reshaped_y = repeat_idx_2d(idx_list_y, len(idx_list_x), 1)

    im_reshaped = np.reshape(im[idx_reshaped_x, idx_reshaped_y, :], [out_x, ker_x, out_y, ker_y, im_w]) # reshapes

    # reshapes the 2D kernel
    ker = np.reshape(ker,[1, ker_x, 1, ker_y, im_w])

    # applies the kernel to the image and reduces the dimension back to the one of original input image
    return np.squeeze(np.sum(im_reshaped*ker, axis=(1,3)))

我试图添加很多注释来解释该方法,但总体思路是将 3D 输入图像重塑为 5D 形状(output_image_height、kernel_height、output_image_width、kernel_width、output_image_channel),然后使用基本的直接应用内核数组乘法。当然,这种方法会使用更多内存(在执行过程中,图像的大小因此乘以 kernel_height*kernel_width)但速度更快。

为了完成这个重塑步骤,我“过度使用”了 numpy 数组的索引方法,特别是,将 numpy 数组作为索引提供给 numpy 数组的可能性。

这种方法也可用于使用基本数学函数在 Pytorch 或 Tensorflow 中重新编码 2D 卷积乘积,但我毫不怀疑它会比现有的 nn.conv2d 运算符慢...

我真的很喜欢只使用 numpy 基本工具来编写这种方法。

于 2020-11-10T20:13:50.447 回答
0

尝试第一轮然后转换为 uint8:

data = data.round().astype(np.uint8);
于 2013-01-30T14:51:12.240 回答
0

最明显的方法之一是对内核进行硬编码。

img = img.convert('L')
a = np.array(img)
out = np.zeros([a.shape[0]-2, a.shape[1]-2], dtype='float')
out += a[:-2, :-2]
out += a[1:-1, :-2]
out += a[2:, :-2]
out += a[:-2, 1:-1]
out += a[1:-1,1:-1]
out += a[2:, 1:-1]
out += a[:-2, 2:]
out += a[1:-1, 2:]
out += a[2:, 2:]
out /= 9.0
out = out.astype('uint8')
img = Image.fromarray(out)

这个例子做了一个完全展开的盒子模糊 3x3。您可以将具有不同值的值相乘,然后将它们除以不同的量。但是,如果你真的想要最快和最肮脏的方法,就是这样。我认为它比 Guillaume Mougeot 的方法高出 5 倍。他的方法比其他方法高出 10 倍。

如果你正在做类似高斯模糊的事情,它可能会丢失一些步骤。并且需要增加一些东西。

于 2020-11-26T04:50:54.187 回答
-2

此代码不正确:

for r in range(nr):
    data[r,:] = np.convolve(data[r,:], H_r, 'same')

for c in range(nc):
    data[:,c] = np.convolve(data[:,c], H_c, 'same')

请参阅从多维卷积到一维的 Nussbaumer 转换。

于 2015-07-03T08:04:43.040 回答