嵌套for循环到numpy卷积答案

【问题标题】：Nested for loop to numpy convolve嵌套for循环到numpy卷积
【发布时间】：2013-01-23 09:30:26
【问题描述】：

如何提高此功能的速度？

def foo(mri_data, radius):

    mask = mri_data.copy()

    ny = len(mri_data[0,:])
    nx = len(mri_data[:])

    for y in xrange(0, ny):
        for x in xrange(0, nx):
            if (mri_data[x-radius:x+radius,y-radius:y+radius] != 1.0).all():
                mask[x,y] = 0.0                    
    return mask.copy()

它以 numpy 数组的形式接收图像切片。遍历每个像素并测试围绕该像素的边界框。如果框中没有值等于 1，则我们通过将其设置为 0 来丢弃该像素。

有人告诉我可以使用numpy.convolve，但我不确定这有什么关系。

编辑：图像值在二进制范围内，因此最小值为 0.0，最大值为 1.0。介于 ex 之间的值：0.767。

【问题讨论】：

小心边缘。如果你有例如radius=3 和 mri_data = np.arange(8) 那么你的第一个窗口是 mri_data[-3:3] 它返回一个空数组...

标签： python numpy convolution

【解决方案1】：

您可以滥用卷积的一种情况。我不会使用它，但是边界很繁琐......

from scipy.ndimage import convolve

not_one = (mri_data != 1.0) # are you sure you want to compare with float like that?!

conv = convolve(not_one, np.ones((2*radius, 2*radius)))
all_not_one = (conv == (2*radius)**2)

mask[all_not_one] = 0

真的应该做同样的事情......

【讨论】：

【解决方案2】：

您正在做的事情称为binary_dilation，但您的代码中有一个小错误。具体来说，当 x, y 小于半径时，您会得到负索引。这些负数是使用 numpy 索引规则解释的，这不是你想要的 more on indexing here，在图像的两个边缘上给你错误的结果。

这里有一些代码使用二进制膨胀来完成同样的事情，并修复了上面提到的错误。

import numpy as np
from scipy.ndimage import binary_dilation

def foo(mri_data, radius):
    structure = np.ones((2*radius, 2*radius))
    # I set the origin here to match your code
    mask = binary_dilation(mri_data == 1, structure, origin=-1)
    return np.where(mask, mri_data, 0)

【讨论】：

感谢 Bi Rico。我不知道我设置的这种阈值已经是一个过滤器。我也将不得不尝试您的解决方案。