中值过滤器在 FITS 文件上产生意外结果
Median filter produces unexpected result on FITS file
这是基于其他几个尚未完全回答的问题，所以我开始了一个新的 post。我正在努力寻找 50 像素补丁中的掩码数组的中值。图片和蒙版都是901x877telescope张图片。
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
# Use the fits files as input image and mask
hdulist = fits.open('xbulge-w1.fits')
w1data = hdulist[0].data
hdulist3 = fits.open('xbulge-mask.fits')
mask = 1 - hdulist3[0].data
w1masked = np.ma.array(w1data, mask = mask)
# Use general arrays as input image and mask
#w1data = np.arange(790177).reshape(901,877)
#w1masked = np.ma.masked_inside(w1data, 30000, 60000)
side = 50
w, h = w1data.shape
width_index = np.array(range(w//side)) * side
height_index = np.array(range(h//side)) * side
def assign_patch(patch, median, side):
"""Break this loop out to prevent 4 nested 'for' loops"""
for j in range(side):
for i in range(side):
patch[i,j] = median
return patch
for width in width_index:
for height in height_index:
patch = w1masked[width:width+side, height:height+side]
median = np.median(patch)
assign_patch(patch, median, side)
plt.imshow(w1masked)
plt.show()
问题是，当我使用通用数组作为输入图像和蒙版（注释掉的部分）时，它工作正常，但是当我使用 FITS 文件时，它会在上生成 'side' 大小的补丁输出图像。我不知道这是怎么回事。
我不知道你的 FITS 文件是什么样的，但有几件事很突出：
np.median 没有考虑 mask。事实上，在最近的 NumPy 版本中，如果尝试这样做（正确地）会打印警告。您应该改用 np.ma.median。如果您要更新 NumPy，您可能会看到：
UserWarning: Warning: 'partition' will ignore the 'mask' of the MaskedArray.
当你知道可以使用切片赋值时，assign_patch函数是不必要的：
w1masked[width:width+side, height:height+side] = median
# instead of "assign_patch(patch, median, side)"
这也比用双循环替换每个值快得多。
我认为问题实际上是因为您使用 np.median 而不是 np.ma.median。蒙版像素可能有很多值，包括 nan、0、inf、...因此，如果将这些值考虑在内（当它们应该被忽略时）可能会产生任何类型问题，特别是如果 median 开始返回 nans 或类似的。
更一般地说，如果你真的想要一个中值过滤器，你不能只计算一个补丁的中值并将补丁中的所有值替换为该中值。您应该使用考虑掩码的中值滤波器。不幸的是，我从未见过在任何广泛使用的 Python 包中实现这样的过滤器。但是，如果你有 numba，你可以查看我的一个（非常实验性的！）包 numbamisc which contains a median_filter，它考虑到了面具。
这是基于其他几个尚未完全回答的问题，所以我开始了一个新的 post。我正在努力寻找 50 像素补丁中的掩码数组的中值。图片和蒙版都是901x877telescope张图片。
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
# Use the fits files as input image and mask
hdulist = fits.open('xbulge-w1.fits')
w1data = hdulist[0].data
hdulist3 = fits.open('xbulge-mask.fits')
mask = 1 - hdulist3[0].data
w1masked = np.ma.array(w1data, mask = mask)
# Use general arrays as input image and mask
#w1data = np.arange(790177).reshape(901,877)
#w1masked = np.ma.masked_inside(w1data, 30000, 60000)
side = 50
w, h = w1data.shape
width_index = np.array(range(w//side)) * side
height_index = np.array(range(h//side)) * side
def assign_patch(patch, median, side):
"""Break this loop out to prevent 4 nested 'for' loops"""
for j in range(side):
for i in range(side):
patch[i,j] = median
return patch
for width in width_index:
for height in height_index:
patch = w1masked[width:width+side, height:height+side]
median = np.median(patch)
assign_patch(patch, median, side)
plt.imshow(w1masked)
plt.show()
问题是，当我使用通用数组作为输入图像和蒙版（注释掉的部分）时，它工作正常，但是当我使用 FITS 文件时，它会在上生成 'side' 大小的补丁输出图像。我不知道这是怎么回事。
我不知道你的 FITS 文件是什么样的，但有几件事很突出：
np.median没有考虑mask。事实上，在最近的 NumPy 版本中，如果尝试这样做（正确地）会打印警告。您应该改用np.ma.median。如果您要更新 NumPy，您可能会看到：
UserWarning: Warning: 'partition' will ignore the 'mask' of the MaskedArray.
当你知道可以使用切片赋值时，
assign_patch函数是不必要的：w1masked[width:width+side, height:height+side] = median # instead of "assign_patch(patch, median, side)"
这也比用双循环替换每个值快得多。
我认为问题实际上是因为您使用 np.median 而不是 np.ma.median。蒙版像素可能有很多值，包括 nan、0、inf、...因此，如果将这些值考虑在内（当它们应该被忽略时）可能会产生任何类型问题，特别是如果 median 开始返回 nans 或类似的。
更一般地说，如果你真的想要一个中值过滤器，你不能只计算一个补丁的中值并将补丁中的所有值替换为该中值。您应该使用考虑掩码的中值滤波器。不幸的是，我从未见过在任何广泛使用的 Python 包中实现这样的过滤器。但是，如果你有 numba，你可以查看我的一个（非常实验性的！）包 numbamisc which contains a median_filter，它考虑到了面具。