Calculating Variance Of An Image Python Efficiently

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): Calculating the variance for each pixel in

Solution 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

Solution 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])
    returnout

def looping_img_var(im, w):
    nx, ny = img.shape
    varianceMatrix = np.zeros(im.shape, np.float32)
    for i inrange(w,nx-w):
        for j inrange(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() -startstart= 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))

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