Сглаживание без заполнения пропущенных значений нулями

Я хотел бы сгладить карту, которая не покрывает все небо. Эта карта не является гауссовой и не имеет нулевого значения, поэтому поведение 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 в scipy.gaussian_filter(). healpy.smoothing() неявно использует mode=constant с cval=0, но мне нужно что-то вроде mode=reflect.

Есть ли разумный способ решить эту проблему?


python healpy
person Daniel Lenz    schedule 24.04.2018    source источник

Ответы (2)


arrow_upward
1
arrow_downward

Эта проблема связана со следующим вопросом и ответом (отказ от ответственности: от меня):

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

person David    schedule 01.05.2018

arrow_upward
2
arrow_downward

Самый простой способ справиться с этим — удалить среднее значение карты, выполнить сглаживание с помощью 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 , основная проблема заключается в том, насколько сложна ваша граница. Если это легко, как разрез по широте, то все это легко, потому что общий случай очень сложен и может иметь много угловых случаев, которые вам нужно обработать:

Определите «пограничные слои»

  • получить все недостающие пиксели
  • найдите соседей и найдите, у кого из них есть действительный сосед, и отметьте его как «первую границу»
  • повторите этот алгоритм и найдите пиксели, у которых есть соседний пиксель "первой границы", и пометьте его как "вторую границу"
  • повторяйте, пока у вас не будут все слои, которые вам нужны

Заполните отраженные значения

  • цикл на пограничных слоях
  • цикл на каждом пикселе слоя
  • найдите действительных соседей, вычислите их барицентр, теперь предположим, что линия между центром граничного пикселя и барицентром проходит перпендикулярно границе маски, а граница маски находится на полпути
  • теперь удлините эту линию, удвоив ее в направлении внутри маски, возьмите интерполированное значение карты в этом месте и назначьте его текущему отсутствующему пикселю.
  • повторите это для других слоев, играя с длиной линии.
person Andrea Zonca    schedule 28.04.2018
comment
Спасибо @Andrea Zonca, это действительно хорошие предложения. Я поиграю с этим и, вероятно, приму ответ позже! Маска действительно сложна (думаю, вырезание интенсивности + точечные источники) - person Daniel Lenz; 29.04.2018