提高代码效率:滑动窗口上的标准偏差 [英] 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]

这比原始功能要快 .对于半径为20的1024x1024数组，旧功能需要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屋！