Perform local standard deviation in Python

Question

I am trying to obtain the local standard deviation of each pixel of an image.This means that for each pixel I want to calculate the standard deviation of its value and its neighbors'. I used this library I developed the following code:

def stdd(image, N):
    width = image.shape[0]
    heigth = image.shape[1]
    desv = np.zeros((width,heigth))
    for i in range (width):
        for j in range (heigth):
            if i < N :
                mini = 0
            else :
                mini = i - N
            if (i+N) > width :
                maxi = width
            else : 
                maxi = N + i
            if j < N :
                minj = 0
            else :
                minj = j - N
            if (j+N) > heigth :
                maxj = heigth
            else : 
                maxj = N + j
            window = image[mini:maxi,minj:maxj]
            desv[i,j] = window.std()
    return desv

Where N is the size of the local matrix for each pixel and image is a numpy.array() image The problem of my code is that it takes too much to process it, and I would like to know if there is an already defined function which optimizes it

Answer 1

You can try the following

from scipy.ndimage import generic_filter
import numpy as np
generic_filter(img, np.std, size=3)

Answer 2

You could try to calculate the standard deviations all at once, using the following identity:

To get the sum of all elements in a local neighborhood, you can use a convolution.

def std_convoluted(image, N):
    im = np.array(image, dtype=float)
    im2 = im**2
    ones = np.ones(im.shape)
    
    kernel = np.ones((2*N+1, 2*N+1))
    s = scipy.signal.convolve2d(im, kernel, mode="same")
    s2 = scipy.signal.convolve2d(im2, kernel, mode="same")
    ns = scipy.signal.convolve2d(ones, kernel, mode="same")
    
    return np.sqrt((s2 - s**2 / ns) / ns)

Warning: Please not that while the results look good on a few test images, this function not return the same results as your code, but I can't spot the error right now. (If someone sees it: would you be so kind to point it out or edit it?)

Anyway, the idea is still valid and runs a lot faster (about a factor of 10 on my computer).