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

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

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

Answer 2

在图像处理中最常用的方法是使用求和区域表,这是1984年本文介绍的一个想法.这个想法是,当你通过添加窗口来计算数量时,移动窗口例如右边一个像素,您不需要在新窗口中添加所有项目,只需从总数中减去最左侧的列,并添加新的最右侧列.因此,如果您在数组的两个维度上创建累积和数组,则可以在窗口上获得具有几个总和和减法的总和.如果为数组及其正方形保留求和区域表,则很容易得到这两者的差异.这是一个实现:

def windowed_sum(a, win):
    table = np.cumsum(np.cumsum(a, axis=0), axis=1)
    win_sum = np.empty(tuple(np.subtract(a.shape, win-1)))
    win_sum[0,0] = table[win-1, win-1]
    win_sum[0, 1:] = table[win-1, win:] - table[win-1, :-win]
    win_sum[1:, 0] = table[win:, win-1] - table[:-win, win-1]
    win_sum[1:, 1:] = (table[win:, win:] + table[:-win, :-win] -
                       table[win:, :-win] - table[:-win, win:])
    return win_sum

def windowed_var(a, win):
    win_a = windowed_sum(a, win)
    win_a2 = windowed_sum(a*a, win)
    return (win_a2 - win_a * win_a / win/ win) / win / win

要看到这个有效:

>>> a = np.arange(25).reshape(5,5)
>>> windowed_var(a, 3)
array([[ 17.33333333,  17.33333333,  17.33333333],
       [ 17.33333333,  17.33333333,  17.33333333],
       [ 17.33333333,  17.33333333,  17.33333333]])
>>> np.var(a[:3, :3])
17.333333333333332
>>> np.var(a[-3:, -3:])
17.333333333333332

这应该比基于卷积的方法更快地运行几个切口.