[SciPy-user] Ignoring pixels with gaussian_filter

[Thomas Robitaille]

Hello,

I am currently using gaussian_filter to smooth an image stored in a
numpy array. My problem is that some pixels have no defined value, and
are set to NaN. If gaussian_filter stumbles upon a NaN value, it will
set all pixels within a certain radius of that value to NaN. Consider
the following example:

---

import numpy as np
from scipy.ndimage import gaussian_filter

a = np.ones((9,9))
a[4,4] = np.nan

b = gaussian_filter(a,sigma=2)

print "Original:"
print a

print "Smoothed:"
print b

---

Of course, I can simply use

a[np.where(np.isnan(a))] = 0.

to reset all NaN pixels to zero, but then the result depends on the
value I replace NaNs by.

Is there a way to get gaussian_filter to simply ignore such pixels
when smoothing? (apart from writing a gaussian filter algorithm from
scratch!)

Thanks for any help,

Thomas

[Zachary Pincus]

Hi Thomas,

Post by Thomas Robitaille
Hello,
I am currently using gaussian_filter to smooth an image stored in a
numpy array. My problem is that some pixels have no defined value, and
are set to NaN. If gaussian_filter stumbles upon a NaN value, it will
set all pixels within a certain radius of that value to NaN.

...

Post by Thomas Robitaille
Of course, I can simply use
a[np.where(np.isnan(a))] = 0.
to reset all NaN pixels to zero, but then the result depends on the
value I replace NaNs by.
Is there a way to get gaussian_filter to simply ignore such pixels
when smoothing? (apart from writing a gaussian filter algorithm from
scratch!)

No simple way to get gaussian_filter to ignore nan pixels when doing
the convolution.
You could write your own convolution function in cython (not as bad as
it sounds, or you could fill in the missing data with something
reasonable (like the median of the neighbors) before gaussian
filtering. These should give pretty similar results.

You could check out an earlier thread on numpy-discussion titled "Help
with interpolating missing values from a 3D scanner" for some
discussion about filling in missing values.

Zach

[Thomas Robitaille]

Hi Zach,

Post by Zachary Pincus

Post by Thomas Robitaille
Is there a way to get gaussian_filter to simply ignore such pixels
when smoothing? (apart from writing a gaussian filter algorithm from
scratch!)

No simple way to get gaussian_filter to ignore nan pixels when doing
the convolution.
You could write your own convolution function in cython [snip]

I looked into this, and figured out how to write my own filter. It
does the job, but is very slow. The code is below, and takes 11s on a
300x300 array with a filter size of 31x31 on my computer. How could I
use cython to speed things up? (I know that I should really pass the
filter size and the filter sigma as arguments using e.g.
extra_arguments, but I want to try and get it to run as fast as
possible by computing the gaussian filter only once.)

Thanks!

Thomas

---

import numpy as np
from scipy.ndimage import generic_filter

# size of filter (in both directions)
s = 31

# center of filter (in both directions)
c = 16

# sigma of filter (in both directions)
sigma = 5.

# define gaussian function
def gaussian(cx, cy, w):
return lambda x,y: np.exp(-(((cx-x)/w)**2+((cy-y)/w)**2)/2)

# define gaussian filter
x,y = np.mgrid[0:s,0:s]
filt = gaussian(c,c,sigma)(x,y).ravel()

# define custom filter
def custom_filter(values):
mask = np.where(np.isnan(values) == False)
return np.sum(values[mask]*filt[mask])/np.sum(filt[mask])

# the function to test the custom filter
def do():
a = np.ones((300,300))
b =
generic_filter(a,custom_filter,size=s,mode='constant',cval=np.nan)

---

[alex]

Post by Thomas Robitaille
...
# define gaussian function
return lambda x,y: np.exp(-(((cx-x)/w)**2+((cy-y)/w)**2)/2)
...

A neat mathematical property of gaussian blur that is not true of 2d
kernels in general is that it can be applied to the x and y axes
separately. That is, it can be implemented as two 1d passes instead of
one 2d pass. This could speed up your code a lot if you aren't already
doing it.

Alex

[Thomas Robitaille]

Thanks a lot for pointing this out! In the end I've used this property
and implemented a fortran routine to do the smoothing, wrapped using
f2py. The resulting routine is almost as fast as the default
gaussian_filter from scipy.

Cheers,

Thomas

Post by alex

Post by Thomas Robitaille
...
# define gaussian function
return lambda x,y: np.exp(-(((cx-x)/w)**2+((cy-y)/w)**2)/2)
...

A neat mathematical property of gaussian blur that is not true of 2d
kernels in general is that it can be applied to the x and y axes
separately. That is, it can be implemented as two 1d passes instead of
one 2d pass. This could speed up your code a lot if you aren't already
doing it.
Alex
_______________________________________________
SciPy-user mailing list
http://mail.scipy.org/mailman/listinfo/scipy-user