Median filter of masked arrays

Question

I have seen several discussions in this forum about calculating the median of masked arrays, such as images. What I want is slightly more subtle, it is to apply a median filter on my image. I know a way to do it, but is far too slow, and would appreciate ways of speeding the process up.

For example, assume I have a masked array of shape (10,10) and I want to apply a median filter with a box (3,3) not using those elements that are masked. My goal is to substitute the value in each pixel of the image with the value of the masked median of the box.

Assuming a very simple case, we can build the "image" and the mask as:

 im = numpy.random.uniform(size=(10,10))
 mask = numpy.zeros_like(im)
 mask[1:3,:] = 1
 masked_im = numpy.ma.array(im, mask=mask)

Now, to actually make the median filter we can do it on a brute-force way with:

 lx, ly = im.shape
 side = 3
 im_filt = numpy.zeros_like(im)
 for jj in range(ly):
     for ii in range(lx):
         minx, maxx = max([ii-side/2,0]), min([ii+side/2+1,lx])
         miny, maxy = max([jj-side/2,0]), min([jj+side/2+1,ly])
         im_filt[ii,jj] = numpy.ma.median(masked_im[minx:maxx, miny:maxy])

This solves the problem and gives a good result, but as I said, it is painfully slow. One (to me, surprising) way to slightly speed up the process is to use the mask and the image separately, like:

 im_filt2 = numpy.zeros_like(im)
 for jj in range(ly):
     for ii in range(lx):
         minx, maxx = max([ii-side/2,0]), min([ii+side/2+1,lx])
         miny, maxy = max([jj-side/2,0]), min([jj+side/2+1,ly])
         zoom_im = im[minx:maxx, miny:maxy]
         zoom_msk = mask[minx:maxx, miny:maxy]
         im_filt2[ii,jj] = numpy.median(zoom_im[zoom_msk == 0]) 

This brings the execution time from 0.018 to 0.002, which is obviously better (why??) if not by the factor ~50 that I was looking for.

Any input?

Answer 1

I guess that the difference is mainly due to overhead in accessing the MaskedArray object (which is a kind of wrapper around ndarray).

For an efficient median filter in numpy you may also try scikit-image. It also accepts a mask argument.