[SciPy-User] Rotated, Anisotropic Gaussian Filtering (Kernel Density Estimation)

Patrick Marsh patrickmarshwx@gmail....
Tue Oct 16 13:02:45 CDT 2012


Greetings,


I know that people on this list are way smarter than I, so hopefully
someone can help me out here.

I have a gridded dataset of 1s and 0s with which I'm needing to apply a
rotated, anisotropic Gaussian filter to achieve a kernel density estimate.
Currently I have some Cython code that I wrote to do this. The code for
this can be found here: https://gist.github.com/3900591. In essence, I
search through the grid for a value of 1. When a 1 is found, a weighting
function is applied to all the surrounding grid points. Their are rotation
matrices at work throughout the code to handle the fact the axes of the
anisotropic Gaussian kernel can be off-cartesian axes. The code works, and
is reasonably efficient for small domains or large grid spacing. However,
as grid spacing decreases, the performance takes a substantial hit.

Recently I started playing around with Scipy's ndimage. The Gaussian filter
routines are exactly what I've been looking for -- as long as my sigma
values lie along the cartesian axes. I attempted to rectify this situation
by first rotating the underyling data grid so that it lined up with the
cartesian grid. (In other words, if the anisotropic Gaussian rotated angle
was 45 degrees, I rotated the underlying data by 45 degrees so they
aligned.) I then applied the Gaussian filters and finally rotated the data
grid back to the original position. This works perfectly…sometimes.

The problem with this approach arrises from the rotating of the data grid.
Since I can't really afford to lose the sharpness of my underlying data
(they are actual observations and need to remain 1s and 0s), I chose to use
a spline of order 0. This works some of the time, but, unfortunatley, this
does not always conserve the sum of the total number of points. For example:


import numpy as np
from scipy import ndimage

dist = 21
midpoint = np.floor(dist/2)
hist = np.zeros((dist, dist), dtype=float)
hist[midpoint, midpoint] = 1
hist2 = ndimage.rotate(hist.copy(), 15, order=0, reshape=True,
prefilter=False)
print(hist.sum(), hist2.sum())

>> 1.0, 0.0


results in hist2 being all 0s, whereas hist has a sum of 1. So, here are my
questions:

1. Is there a way to do an image rotation such that grid points aren't lost
when using a spline of order 0?
2. Is there an alternative way to achieve what I'm trying to do?
3. Is there a way to speed up the Cython code linked above?



Thanks for letting me pick your brain…

Patrick
---
Patrick Marsh
Ph.D. Candidate / Liaison to the HWT
School of Meteorology / University of Oklahoma
Cooperative Institute for Mesoscale Meteorological Studies
National Severe Storms Laboratory
http://www.patricktmarsh.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20121016/1c5908a8/attachment.html 


More information about the SciPy-User mailing list