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

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

首先，有不止一种方法可以做到这一点。

这不是最高效的速度，但使用scipy.ndimage.generic_filter 可以让您轻松地在移动窗口上应用任意 python 函数。

举个简单的例子：

result = scipy.ndimage.generic_filter(data, np.std, size=2*radius)

注意边界条件可以通过mode kwarg 来控制。

---

另一种方法是使用一些不同的跨步技巧来查看实际上是一个移动窗口的数组，然后沿最后一个轴应用np.std。 （注意：这取自我之前的答案之一：https://stackoverflow.com/a/4947453/325565）

def strided_sliding_std_dev(data, radius=5):
    windowed = rolling_window(data, (2*radius, 2*radius))
    shape = windowed.shape
    windowed = windowed.reshape(shape[0], shape[1], -1)
    return windowed.std(axis=-1)

def rolling_window(a, window):
    """Takes a numpy array *a* and a sequence of (or single) *window* lengths
    and returns a view of *a* that represents a moving window."""
    if not hasattr(window, '__iter__'):
        return rolling_window_lastaxis(a, window)
    for i, win in enumerate(window):
        if win > 1:
            a = a.swapaxes(i, -1)
            a = rolling_window_lastaxis(a, win)
            a = a.swapaxes(-2, i)
    return a

def rolling_window_lastaxis(a, window):
    """Directly taken from Erik Rigtorp's post to numpy-discussion.
    <http://www.mail-archive.com/numpy-discussion@scipy.org/msg29450.html>"""
    if window < 1:
       raise ValueError, "`window` must be at least 1."
    if window > a.shape[-1]:
       raise ValueError, "`window` is too long."
    shape = a.shape[:-1] + (a.shape[-1] - window + 1, window)
    strides = a.strides + (a.strides[-1],)
    return np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides)

乍一看有点难以理解这里发生了什么。不要插入我自己的答案之一，但我不想重新输入解释，所以请看这里：https://stackoverflow.com/a/4924433/325565，如果您以前没有见过这些“跨步”技巧。

如果我们将时序与 100x100 随机浮点数组与 radius 为 5 进行比较，它比原始版本或 generic_filter 版本快约 10 倍。但是，您在此版本的边界条件上没有灵活性。 （这与您当前所做的相同，而 generic_filter 版本为您提供了很大的灵活性，但牺牲了速度。）

# Your original function with nested loops
In [21]: %timeit sliding_std_dev(data)
1 loops, best of 3: 237 ms per loop

# Using scipy.ndimage.generic_filter
In [22]: %timeit ndimage_std_dev(data)
1 loops, best of 3: 244 ms per loop

# The "stride-tricks" version above
In [23]: %timeit strided_sliding_std_dev(data)
100 loops, best of 3: 15.4 ms per loop

# Ophion's version that uses `np.take`
In [24]: %timeit new_std_dev(data)
100 loops, best of 3: 19.3 ms per loop

“stride-tricks”版本的缺点是，与“普通”跨步滚动窗口技巧不同，此版本确实制作副本，而且它远大于原始数组。如果您在大型阵列上使用它，您会遇到内存问题！ （附带说明一下，就内存使用和速度而言，它基本上等同于 @Ophion 的答案。这只是做同样事情的不同方法。）

Answer 2

您可以先获取索引，然后使用np.take 组成新数组：

def new_std_dev(image_original,radius=5):
    cols,rows=image_original.shape

    #First obtain the indices for the top left position
    diameter=np.arange(radius*2)
    x,y=np.meshgrid(diameter,diameter)
    index=np.ravel_multi_index((y,x),(cols,rows)).ravel()

    #Cast this in two dimesions and take the stdev
    index=index+np.arange(rows-radius*2)[:,None]+np.arange(cols-radius*2)[:,None,None]*(rows)
    data=np.std(np.take(image_original,index),-1)

    #Add the zeros back to the output array
    top=np.zeros((radius,rows-radius*2))
    sides=np.zeros((cols,radius))

    data=np.vstack((top,data,top))
    data=np.hstack((sides,data,sides))
    return data

首先生成一些随机数据并检查时序：

a=np.random.rand(50,20)

print np.allclose(new_std_dev(a),sliding_std_dev(a))
True

%timeit sliding_std_dev(a)
100 loops, best of 3: 18 ms per loop

%timeit new_std_dev(a)
1000 loops, best of 3: 472 us per loop

对于较大的数组，只要您有足够的内存，它总是更快：

a=np.random.rand(200,200)

print np.allclose(new_std_dev(a),sliding_std_dev(a))
True

%timeit sliding_std_dev(a)
1 loops, best of 3: 1.58 s per loop

%timeit new_std_dev(a)
10 loops, best of 3: 52.3 ms per loop

原始函数对于非常小的数组更快，看起来收支平衡点是hgt*wdt >50。需要注意的是，您的函数采用方形框架并将 std dev 放在右下角的索引中，而不是在索引周围采样。这是故意的吗？

Answer 3

很酷的技巧：您可以只给定平方值的总和和窗口中的值的总和来计算标准偏差。

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

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 4

在图像处理中最常用的方法是使用求和面积表，这是 1984 年在this paper 中引入的一个想法。这个想法是，当您通过在窗口上添加来计算数量时，并且移动窗口，例如右移一个像素，不需要把新窗口中的所有项都加进去，只需要从总数中减去最左边的一列，再加上新的最右边一列。因此，如果您在数组的两个维度上创建一个累加和数组，您可以通过一个窗口获得总和，其中包含几个和和一个减法。如果您为数组及其正方形保留总面积表，则很容易从这两个中获得方差。这是一个实现：

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

这应该比基于卷积的方法快几个档次。