# [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

n = signal.shape[0]
if tmp.shape[0] < 4: raise ValueError, \
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]

n,m = image.shape[0],image.shape[1]
if tmp.shape[0] < 4: raise ValueError, \
if tmp.shape[1] < 4: raise ValueError, \
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]

```