[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