1

我想平滑一张不覆盖整个天空的地图。该映射不是高斯映射，也不是零，因此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在谐波空间中执行平滑。

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

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

4

2 回答 2

2

处理此问题的最简单方法是移除地图的平均值，使用 执行平滑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

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

识别“边界层”

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

填充反射值

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

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

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

于 2018-05-01T11:17:31.057 回答