python - 平滑而不用零填充缺失值

Question

我想平滑一张不能覆盖整个天空的 map 。该映射不是高斯分布，也不是零，因此 healpy 的默认行为(用 0 填充缺失值)会导致该掩码边缘偏向较低值:

import healpy as hp

nside = 128
npix = hp.nside2npix(nside)

arr = np.ones(npix)
mask = np.zeros(npix, dtype=bool)

mask[:mask.size//2] = True

arr[~mask] = hp.UNSEEN
arr_sm = hp.smoothing(arr, fwhm=np.radians(5.))

hp.mollview(arr, title='Input array')
hp.mollview(arr_sm, title='Smoothed array')

我想通过将屏蔽值的权重设置为零来保留锐边，而不是将值设置为零。这似乎很困难，因为 healpy 在调和空间中执行平滑。

更具体地说，我想模仿 scipy.gaussian_filter() 中的 mode 关键字。 。 healpy.smoothing() 隐式使用 mode=constant 和 cval=0，但我需要类似 mode=reflect.

有什么合理的方法可以解决这个问题吗？

Answer 1

处理此问题的最简单方法是删除 map 的平均值，使用 hp.smoothing 执行平滑。 ，然后将偏移量加回来。这可以解决这个问题，因为现在 map 是零均值的，所以零填充不会产生边框效果。

def masked_smoothing(m, fwhm_deg=5.0):
#make sure m is a masked healpy array
m = hp.ma(m)
offset = m.mean()
smoothed=hp.smoothing(m - offset, fwhm=np.radians(fwhm_deg))<br/>
return smoothed + offset


我能想到的另一个选项是一些迭代算法，在平滑之前以“反射”模式填充 map ，可能在 cython 中实现或numba ，主要问题是你的边界有多复杂。如果它像纬度切割一样简单，那么所有这一切都很容易，因为一般情况非常复杂，并且可能有很多需要处理的极端情况:

识别“边界层”

• 获取所有缺失的像素
• 查找邻居并找出哪个有有效邻居并将其标记为“第一边界”
• 重复此算法并找到具有“第一边界”像素邻居的像素并将其标记为“第二边界”
• 重复此操作，直到获得所需的所有图层

填充反射值

• 在边界层上循环
• 在每层像素上循环
• 找到有效的邻居，计算它们的重心，现在假设边界像素中心和重心之间的线垂直穿过掩模边界，并且掩模边界位于中间
• 现在通过在蒙版内的方向加倍来扩展这条线，在该位置获取 map 的插值并将其分配给当前缺失的像素
• 通过调整线条的长度对其他层重复此操作。