提高代码效率：滑动窗口上的标准偏差

Question

我正在尝试改善为图像的每个像素计算位于像素附近的像素的标准偏差的功能。我的函数使用两个嵌入式循环在整个矩阵中运行，这是我程序的瓶颈。我猜有可能通过摆脱numpy摆脱循环来改善它，但是我不知道如何进行。欢迎任何建议！

问候

def sliding_std_dev(image_original,radius=5) :
    height, width = image_original.shape
    result = np.zeros_like(image_original) # initialize the output matrix
    hgt = range(radius,height-radius)
    wdt = range(radius,width-radius)
    for i in hgt:
        for j in wdt:
            result[i,j] = np.std(image_original[i-radius:i+radius,j-radius:j+radius])
    return result

Answer 1

绝妙的技巧：您可以仅在平方值和窗口中的值总和下计算标准偏差。

因此，您可以对数据使用统一滤波器来非常快速地计算标准偏差：

from scipy.ndimage.filters import uniform_filter

def window_stdev(arr, radius):
    c1 = uniform_filter(arr, radius*2, mode='constant', origin=-radius)
    c2 = uniform_filter(arr*arr, radius*2, mode='constant', origin=-radius)
    return ((c2 - c1*c1)**.5)[:-radius*2+1,:-radius*2+1]

这是 可笑的 比原来的功能更快。对于1024x1024阵列且半径为20的情况，旧功能需要34.11秒，新功能需要 0.11秒 ，加速300倍。

这在数学上如何工作？它计算sqrt(mean(x^2) - mean(x)^2)每个窗口的数量。我们可以从标准偏差中得出以下数量sqrt(mean((x - mean(x))^2))：

让我们E成为期望操作符（基本上是mean()），并X成为数据的随机变量。然后：

E[(X - E[X])^2]
= E[X^2 - 2X*E[X] + E[X]^2]
= E[X^2] - E[2X*E[X]] + E[E[X]^2]（通过期望运算符的线性）
= E[X^2] - 2E[X]*E[X] + E[X]^2（再次通过线性和E[X]为常数的事实）
= E[X^2] - E[X]^2

证明使用此技术计算的数量在数学上等同于标准偏差。