Calculating variance of an image python efficiently

Question

I'm working on a project in which need to get the variance of an image. Currently I'm taking 2 approaches (both work but are very slow):

1. Calculating the variance for each pixel individually:

This is the code using numpy, varianceMatrix is the output

varianceMatrix = np.zeros(im.shape,np.uint8)
w = 1              # the radius of pixels neighbors 
ny = len(im)
nx = len(im[0])


for i in range(w,nx-w):
    for j in range(w,ny-w):

        sampleframe = im[j-w:j+w, i-w:i+w]
        variance    = np.var(sampleframe)
        varianceMatrix[j][i] = int(variance)

return varianceMatrix   

2. Using an existing scipy function:

This is the scipy function:

from scipy import ndimage

varianceMatrix = ndimage.generic_filter(im, np.var, size = 3)

The scipy function is faster, but not so much. I'm looking for a better alternative to calculate the variance.

Any ideas???

Answer 1

Here a fast solution using OpenCV:

import cv2

def winVar(img, wlen):
  wmean, wsqrmean = (cv2.boxFilter(x, -1, (wlen, wlen),
    borderType=cv2.BORDER_REFLECT) for x in (img, img*img))
  return wsqrmean - wmean*wmean

On my machine and for the following example, winVar() is 2915 times faster than ndimage.generic_filter() and 10.8 times faster than sliding_img_var() (see pv.'s answer):

In [66]: img = np.random.randint(0, 256, (500,500)).astype(np.float)

In [67]: %timeit winVar(img, 3)
100 loops, best of 3: 1.76 ms per loop

In [68]: %timeit ndimage.generic_filter(img, np.var, size=3)
1 loops, best of 3: 5.13 s per loop

In [69]: %timeit sliding_img_var(img, 1)
100 loops, best of 3: 19 ms per loop

Result matches that of ndimage.generic_filter():

In [70]: np.allclose(winVar(img, 3), ndimage.generic_filter(img, np.var, size=3))
Out[70]: True

Answer 2

You can use a well-known sliding window stride trick to speed up the computation. It add two "virtual dimensions" to the end of the array without copying the data, and then computes the variance over them.

Note that in your code, im[j-w:j+w, ..] goes over indices j-w,j-w+1,...,j+w-1, the last one is exclusive, which you might not have meant. Also, the variances are larger than the uint8 range, so you end up with integer wraparound.

import numpy as np
import time
np.random.seed(1234)

img = (np.random.rand(200, 200)*256).astype(np.uint8)

def sliding_window(a, window, axis=-1):
    shape = list(a.shape) + [window]
    shape[axis] -= window - 1
    if shape[axis] < 0:
        raise ValueError("Array too small")
    strides = a.strides + (a.strides[axis],)
    return np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides)

def sliding_img_var(img, window):
    if window <= 0:
        raise ValueError("invalid window size")
    buf = sliding_window(img, 2*window, 0)
    buf = sliding_window(buf, 2*window, 1)

    out = np.zeros(img.shape, dtype=np.float32)
    np.var(buf[:-1,:-1], axis=(-1,-2), out=out[window:-window,window:-window])
    return out

def looping_img_var(im, w):
    nx, ny = img.shape
    varianceMatrix = np.zeros(im.shape, np.float32)
    for i in range(w,nx-w):
        for j in range(w,ny-w):
            sampleframe = im[j-w:j+w, i-w:i+w]
            variance    = np.var(sampleframe)
            varianceMatrix[j][i] = variance
    return varianceMatrix

np.set_printoptions(linewidth=1000, edgeitems=5)
start = time.time()
print(sliding_img_var(img, 1))
time_sliding = time.time() - start
start = time.time()
print(looping_img_var(img, 1))
time_looping = time.time() - start
print("duration: sliding: {0} s, looping: {1} s".format(time_sliding, time_looping))

Answer 3

If the method using ndimage.generic_filter is not fast enough, you can write your own optimized implementation for the variance calculation in Cython.