提高代码效率:滑动窗口的标准偏差 [英] improving code efficiency: standard deviation on sliding windows

问题描述

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

I am trying to improve function which calculate for each pixel of an image the standard deviation of the pixels located in the neighborhood of the pixel. My function uses two embedded loops to run accross the matrix, and it is the bottleneck of my programme. I guess there is likely a way to improve it by getting rid of the loops thanks to numpy, but I don't know how to proceed. Any advice are welcome!

问候

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

推荐答案

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

Cool trick: you can compute the standard deviation given just the sum of squared values and the sum of values in the window.

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

Therefore, you can compute the standard deviation very fast using a uniform filter on the data:

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 倍.

This is ridiculously faster than the original function. For a 1024x1024 array and a radius of 20, the old function takes 34.11 seconds, and the new function takes 0.11 seconds, a speed-up of 300-fold.

这在数学上是如何工作的?它计算每个窗口的数量 sqrt(mean(x^2) - mean(x)^2).我们可以从标准偏差 sqrt(mean((x - mean(x))^2)) 推导出这个量，如下所示:

How does this work mathematically? It computes the quantity sqrt(mean(x^2) - mean(x)^2) for each window. We can derive this quantity from the standard deviation sqrt(mean((x - mean(x))^2)) as follows:

设 E 是期望操作符(基本上是 mean())，X 是数据的随机变量.然后:

Let E be the expectation operator (basically mean()), and X be the random variable of data. Then:

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

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] (by the linearity of the expectation operator)
= E[X^2] - 2E[X]*E[X] + E[X]^2 (again by linearity, and the fact that E[X] is a constant)
= E[X^2] - E[X]^2

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

which proves that the quantity computed using this technique is mathematically equivalent to the standard deviation.

这篇关于提高代码效率:滑动窗口的标准偏差的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！