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!
regards
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 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]
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.
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:
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] (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.