平滑而不用零填充缺失值

Question

我想平滑未覆盖整个天空的地图。该映射不是高斯映射，也不是零，因此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

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

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))
return smoothed + offset

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

识别“边界层”

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

填充反射值

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

Answer 2

此问题与以下问答有关（免责声明：来自我）：

https://stackoverflow.com/a/36307291/5350621

它可以按如下方式转移到您的案例中：

import numpy as np
import healpy as hp

nside = 128
npix = hp.nside2npix(nside)

# using random numbers here to see the actual smoothing
arr = np.random.rand(npix) 
mask = np.zeros(npix, dtype=bool)
mask[:mask.size//2] = True

def masked_smoothing(U, rad=5.0):     
    V=U.copy()
    V[U!=U]=0
    VV=hp.smoothing(V, fwhm=np.radians(rad))    
    W=0*U.copy()+1
    W[U!=U]=0
    WW=hp.smoothing(W, fwhm=np.radians(rad))    
    return VV/WW

# setting array to np.nan to exclude the pixels in the smoothing
arr[~mask] = np.nan
arr_sm = masked_smoothing(arr)
arr_sm[~mask] = hp.UNSEEN

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