嵌套循环转换为numpy卷积

Question

我该如何提高这个函数的速度呢？

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，中间的值比如0.767。

Answer 1

你正在做的事情叫做 binary_dilation，但是你的代码里有一个小错误。具体来说，当 x 和 y 的值小于半径时，你会得到负数的索引。这些负数在使用 numpy 的索引规则时会被解释成其他的意思，而这并不是你想要的结果，这会导致你图像的两边出现错误的结果。关于索引的更多信息可以在这里找到.

下面是一些代码，它使用二进制膨胀来实现相同的功能，并修复了上面提到的错误。

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)

Answer 2

这是一个关于卷积的滥用案例。虽然我不会这样做，但如果不这样处理，边界问题会变得很麻烦...

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

其实应该是做同样的事情...