[SciPy-dev] recursive Gaussian filter in C

[Sturla Molden]

I have written a C version of the recursive Gaussian filter. It is more 
accurate than gaussian_filter in ndimage (less truncation error near 
edges, as shown before) and also faster.  Here are some timings on 
filtering the "lenna" test image on my computer:
sigma = 5
elapsed time (iir with openmp): 81.366640 ms
elapsed time (iir): 107.282360 ms
elapsed time (ndimage): 137.548760 ms

sigma = 9
elapsed time (iir with openmp): 44.403760 ms
elapsed time (iir): 75.285720 ms
elapsed time (ndimage): 163.750920 ms

sigma = 21
elapsed time (iir with openmp): 45.063040 ms
elapsed time (iir): 75.941400 ms
elapsed time (ndimage): 313.359080 ms

sigma = 101
elapsed time (iir with openmp): 56.134120 ms
elapsed time (iir): 87.622240 ms
elapsed time (ndimage): 1210.016680 ms


It is still only written for np.float64, but it would be easy to make 
optimized versions for various dtypes, including rgb images.

You don't really see the scaling effect of OpenMP here, as the only 
thing that changes with sigma is the amount of zero-padding. Anyhow, 
this beats ndimage on speed and accuracy, and scales much better with 
sigma.

It is not restricted to 2 dimensions. It can filter along any axis of an 
ndarray. Thus another use case is fast kernel density estimation in nd 
space.


Regards,
Sturla Molden