[SciPy-dev] fast gaussian filtering

Sturla Molden sturla@molden...
Fri Dec 19 09:59:45 CST 2008


Gaussian filtering is a very useful method for smoothing signals and 
images, as well as for kernel density estimation. I have previously used 
it to estimate the firing rates of neurons.

I've noticed that SciPy implements Gaussian filtering (in the ndimage 
package) using convolution. It is implemented by loops in Python, not 
even using a vectorized numpy.convolve or an FFT. For the sake of speed, 
SciPy is using the worst thinkable implementation.

To make it fast:

Gaussian filtering can be implemented using recursion (Young & Vliet, 
1995; Signal Processing, 44: 139-151). The recursive approximation is 
very accurate and does not introduce ringing. It is anti-causal 
(forward-backward) and has zero phase response. The filtering can be 
done in C using scipy's signal.lfilter function.

I found a solution in plain C on the Matlab central. But there is really 
no benefit from doing this in C as lfilter, flipud and fliplr does all 
the hard work. An efficient Matlab version could have been constructed 
similarly, without resorting to C.

http://www.mathworks.com/matlabcentral/fileexchange/5087
http://www.ph.tn.tudelft.nl/~lucas/publications/papersLJvV.html

Note that Young & Vliet (1995) also provides recursive solutions for the 
first, second and third derivatives. They are equally easy to implement 
using lfilter.


Regards,
Sturla Molden




from numpy import array, zeros, ones, flipud, fliplr
from scipy.signal import lfilter
from math import sqrt

def __gausscoeff(s):
     if s < .5: raise ValueError, \
         'Sigma for Gaussian filter must be >0.5 samples'
     q = 0.98711*s - 0.96330 if s > 0.5 else 3.97156 \
         - 4.14554*sqrt(1.0 - 0.26891*s)
     b = zeros(4)
     b[0] = 1.5785 + 2.44413*q + 1.4281*q**2 + 0.422205*q**3
     b[1] = 2.44413*q + 2.85619*q**2 + 1.26661*q**3
     b[2] = -(1.4281*q**2 + 1.26661*q**3)
     b[3] = 0.422205*q**3
     B = 1.0 - ((b[1] + b[2] + b[3])/b[0])
     # convert to a format compatible with lfilter's
     # difference equation
     B = array([B])
     A = ones(4)
     A[1:] = -b[1:]/b[0]
     return B,A

def Gaussian1D(signal, sigma, padding=0):
     n = signal.shape[0]
     tmp = zeros(n + padding)
     if tmp.shape[0] < 4: raise ValueError, \
         'Signal and padding too short'
     tmp[:n] = signal
     B,A = __gausscoeff(sigma)
     tmp = lfilter(B, A, tmp)
     tmp = tmp[::-1]
     tmp = lfilter(B, A, tmp)
     tmp = tmp[::-1]
     return tmp[:n]

def Gaussian2D(image, sigma, padding=0):
     n,m = image.shape[0],image.shape[1]
     tmp = zeros((n + padding, m + padding))
     if tmp.shape[0] < 4: raise ValueError, \
         'Image and padding too small'
     if tmp.shape[1] < 4: raise ValueError, \
         'Image and padding too small'
     B,A = __gausscoeff(sigma)
     tmp[:n,:m] = image
     tmp = lfilter(B, A, tmp, axis=0)
     tmp = flipud(tmp)
     tmp = lfilter(B, A, tmp, axis=0)
     tmp = flipud(tmp)
     tmp = lfilter(B, A, tmp, axis=1)
     tmp = fliplr(tmp)
     tmp = lfilter(B, A, tmp, axis=1)
     tmp = fliplr(tmp)
     return tmp[:n,:m]











More information about the Scipy-dev mailing list