【问题标题】:Vectorized extraction of submatrices in numpy-arraynumpy-array中子矩阵的矢量化提取
【发布时间】:2020-05-02 09:48:50
【问题描述】:

我的目标是实现一个中值滤波器,它是一个函数,它将(大部分)二维数组中的每个像素替换为其周围像素的中值。可用于图像去噪。

我的实现从原始矩阵中提取子矩阵,其中包含像素本身及其邻居。这种提取目前是通过 for 循环完成的,正如您可以想象的那样,for 循环占用了大约 95% 的执行时间。

这是我当前的实现:

def median_blur(img, fsize=3):
    img_intermediate = np.zeros((img.shape[0] + 2, img.shape[1] + 2), np.uint8)  # First intermediate img, used for padding original image
    img_intermediate[1:img.shape[0]+1, 1:img.shape[1]+1] = img
    img_result = np.empty((*img.shape, fsize, fsize), np.uint8)  # Will contain result, first receives kernel-submatrices

    # Extract submatrices from image
    for i in range(img.shape[0]):
        for j in range(img.shape[1]):
            img_result[i, j] = img_intermediate[i:i+fsize, j:j+fsize]

    img_result = img_result.reshape(*img.shape, fsize**2)  # Reshape Result-Matrix
    img_result = np.median(img_result, axis=2)  # Calculate median in one go
    return img_result.astype('uint8')

如何使用矢量化操作提取这些子矩阵?

作为奖励,如果有计算机视觉经验的人正在阅读本文:除了将中值滤波器应用于中间零填充矩阵之外,还有更好的方法来实现中值滤波器吗?

非常感谢。

【问题讨论】:

  • median filter in scipy 涵盖了您的用例。
  • @edoput 它被多个库(scipy、scikit-image、opencv...)所覆盖,但我的目标是以最高效的方式自己实现它。
  • 使用view-as_windows : stackoverflow.com/questions/42831022 然后使用np.median?
  • @Divakar 感谢您的提及,我不知道该功能!它不会是唯一的 numpy 解决方案,但肯定会做我想做的事。
  • view_as_windows 是一个独立的函数代码库。因此,您可以从源代码中获取它。

标签: python numpy image-processing computer-vision vectorization


【解决方案1】:

这是一个矢量化解决方案。但是,您可以通过注意图像数组的内存顺序来提出更快的解决方案:

from numpy.lib.stride_tricks import as_strided

img_padded = np.pad(img, 1, mode='constant')
sub_shape = (fsize, fsize)
view_shape = tuple(np.subtract(img_padded.shape, sub_shape) + 1) + sub_shape
sub_img = as_strided(img_padded, view_shape, img_padded.strides * 2)
sub_img = sub_img.reshape(img.shape + (fsize**2,))
result = np.median(sub_img, axis=2).astype(int)

先使用pad 函数填充,然后使用as_strided 获取子矩阵,最后将median 应用于您的步幅。

更新:在 cmets 中使用 @Divakar 建议的 view_as_windows

from skimage.util.shape import view_as_windows

img_padded = np.pad(img, 1, mode='constant')
sub_shape = (fsize, fsize)
view_shape = view_as_windows(img_padded, sub_shape, 1)
sub_img = view_shape.reshape(img.shape + (fsize**2,))
result = np.median(sub_img, axis=2).astype(int)

view_as_windows 还提供了类似于 strides 的子矩阵。

示例图像和输出:

img: 
[[ 1  2  3  4  5  6]
 [ 7  8  9 10 11 12]
 [13 14 15 16 17 18]
 [19 20 21 22 23 24]]

median_filtered: 
[[ 0  2  3  4  5  0]
 [ 2  8  9 10 11  6]
 [ 8 14 15 16 17 12]
 [ 0 14 15 16 17  0]]

【讨论】:

  • 漂亮!这使我现有的解决方案大约。快 7 倍!您能否详细说明可能的内存顺序增益?
  • 很高兴它有帮助。如果您知道您的数据类型,可能有一种更快的方法可以使用 numba 或跨步循环来实现这一点。如果您的数组以C-contiguous 方式保存(即行在前),您可能会从中受益。如果您的数组不是C-contiguousview_as_windows 将创建一个新数组。我希望as_strides 也会这样做,但我不确定。一般来说,任何比这更快的东西可能需要进入C-contiguous/F-contiguous 内存分配并了解您的数据类型。
猜你喜欢
  • 2013-10-10
  • 2018-05-24
  • 1970-01-01
  • 2018-02-11
  • 1970-01-01
  • 1970-01-01
  • 2019-10-18
  • 2015-09-04
  • 1970-01-01
相关资源
最近更新 更多