python - 提高代码效率: standard deviation on sliding windows

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

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