Rolling Statistics on Image Python

Question

I need to calculate the local statistics of a image depending on a 2D Window block defined by the user. Stats include : Mean, Variance, Skew, Kurtosis. I need to traverse through each pixel of the image and find the neighboring pixels depending on the window size.

The code that I used was:

scipy.ndimage.generic_filter(array,numpy.var,size=3)

But the performance through this is very low. I even tried strides-numpy but that too isn't showing much difference (wasn't able to compute skewness, kurtosis). I'm not familiar with Cython so have not ventured into that option. So is there any other way to accomplish this without Cython?

Answer 1

The reason uniform_filter() is so much faster than generic_filter() is due to Python -- for generic_filter(), Python gets called for each pixel, while for uniform_filter(), the whole image is processed in native code. (I found OpenCV's boxFilter() even faster than uniform_filter(), see my answer to a "window variance" question.)

In the remainder of this answer, I show how to do a skew calculation using uniform_filter(), which dramatically speeds up a generic_filter()-based version such as:

import scipy.ndimage as ndimage, scipy.stats as st
ndimage.generic_filter(img, st.skew, size=(1,5))

SciPy's st.skew() (see, e.g., v0.17.0) appears to calculate the skew as

m3 / m2**1.5

where m3 = E[(X-m)**3] (the third central moment), m2 = E[(X-m)**2] (the variance), and m = E[X] (the mean).

To use uniform_filter(), one has to write this in terms of raw moments such as m3p = E[X**3] and m2p = E[X**2] (a prime symbol is usually used to distinguish the raw moment from the central one):

m3 = E[(X-m)**3] = ... = m3p - 3*m*m2p + 2*m**3
m2 = E[(X-m)**2] = ... = m2p - m*m

(In case my "..." skips too much, this answer has the full derivation for m2.) Then one can implement skew() using uniform_filter() (or boxFilter() for some additional speedup):

def winSkew(img, wsize):
  imgS = img*img
  m, m2p, m3p = (ndimage.uniform_filter(x, wsize) for x in (img, imgS, imgS*img))
  mS = m*m
  return (m3p-3*m*m2p+2*mS*m)/(m2p-mS)**1.5

Compared to generic_filter(), winSkew() gives a 654-fold speedup on the following example on my machine:

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

In [186]: %timeit ndimage.generic_filter(img, st.skew, size=(1,5))
1 loops, best of 3: 14.2 s per loop

In [188]: %timeit winSkew(img, (1,5))
10 loops, best of 3: 21.7 ms per loop

And the two calculations give essentially identical results:

In [190]: np.allclose(winSkew(img, (1,5)), ndimage.generic_filter(img, st.skew, size=(1,5)))
Out[190]: True

The code for a Kurtosis calculation can be derived the same way.

Answer 2

The problem is that generic_filter() cannot assume that your filter is separable along the x or y axes. Thus it must operate as a true 2D filter rather than a series of two 1D filters, so run-time will be much slower.

The mean filter and is equivalent (I think) to the uniform_filter(), which if you read the documentation, is implemented as a series of two 1d uniform filters.

I compared timing via this code block:

import numpy as np
from scipy import ndimage as ndi
from scipy import misc
baboonfile = '/Users/curt/Downloads/BaboonRGB.jpg'  #local download of http://read.pudn.com/downloads169/sourcecode/graph/texture_mapping/776733/Picture/BaboonRGB__.jpg
im = misc.imread(baboonfile)

meanfilt2D = ndi.generic_filter(im, np.mean, size=[3, 3, 1])
%timeit meanfilt2D = ndi.generic_filter(im, np.mean, size=[3, 3, 1])
print meanfilt2D.shape

meanfiltU = ndi.uniform_filter(im, size=[3, 3, 1])
%timeit meanfiltU = ndi.uniform_filter(im, size=[3, 3, 1])
print meanfiltU.shape

The output of that block was:

1 loops, best of 3: 5.22 s per loop
(512, 512, 3)
100 loops, best of 3: 11.8 ms per loop
(512, 512, 3)

so true two-dimensional generic_filter() takes 5 seconds for a small image but the two-pass 1D uniform_filter() takes only milliseconds. (N.B.: The difference image meanfilt2D-meanfiltU was not identically zero, but the maximum element was 2; I think the differences are caused by rounding and the imprecise datatype (uint8) used for im.)

For variance and other filters, you should see this old Stack Overflow post which answers a highly related question.