平滑而不用零填充缺失值

Smoothing without filling missing values with zeros

我想平滑一张没有覆盖整个天空的地图。该图不是高斯分布的,也不是零值,因此 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 以谐波 space.

执行平滑

更具体地说,我想模仿 scipy.gaussian_filter() 中的 mode 关键字。 healpy.smoothing() 隐含地使用 mode=constantcval=0,但我需要像 mode=reflect.

这样的东西

有什么合理的方法可以克服这个问题吗?

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

def masked_smoothing(m, fwhm_deg=5.0): #make 确保 m 是一个屏蔽的 healpy 数组 m = hp.ma(m) 偏移量 = m.mean() 平滑=hp.smoothing(m - 偏移量,fwhm=np.radians(fwhm_deg))<br> return 平滑 + 偏移

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

识别"border layers"

  • 获取所有丢失的像素
  • 查找邻居,找出哪一个有有效的邻居并将其标记为"first border"
  • 重复此算法并找到具有 "first border" 像素邻居的像素并将其标记为 "second border"
  • 重复直到获得所需的所有图层

填充反射值

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

此问题与以下问答相关(免责声明:来自我):

可以通过以下方式转入您的案例:

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')