Masked median filters

[Jim421616]

As part of studying towards a MSc in Astronomy, I'm trying to filter an image so that the value of each pixel is equal to the value of the median of the pixels within a 50x50 square around it, excluding any masked pixels.
Here's my latest attempt:

import numpy as np
from astropy.io import fits
from skimage.util.shape import view_as_windows

#   Use the fits files as input image and mask
hdulist = fits.open('xbulge-w1.fits')
image = hdulist[0].data
hdulist3 = fits.open('xbulge-mask.fits')
mask = 1 - hdulist3[0].data
imagemasked = np.ma.masked_array(image, mask = mask)

side = 20
window_shape = (side, side)

Afiltered = view_as_windows(imagemasked, window_shape)

# collapse the last two dimensions in one
flatten_view = Afiltered.reshape(Afiltered.shape[0], Afiltered.shape[1], -1)
print(flatten_view.shape)

# resampling the image by taking median
median_view = np.ma.median(flatten_view, axis=2)

In the code here, I've used some methods from the documentation of skimage.util.view_as_windows to produce the filtered image.
It looks like it's ignoring the mask, but apart from that, I think it works fine. Any advice on processing the masked image?