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

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

0 голосов
/ 24 апреля 2018

Я хотел бы сгладить карту, которая не покрывает полное небо.Эта карта не является гауссовой и не имеет нулевого значения, поэтому поведение по умолчанию 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')

enter image description here enter image description here

Я бы хотел сохранить острый край, установив вес маскируемых значений на ноль, вместо того, чтобы устанавливать значения на ноль,Это кажется трудным, потому что healpy выполняет сглаживание в гармоническом пространстве.

Чтобы быть более точным, я хотел бы имитировать ключевое слово mode в scipy.gaussian_filter().healpy.smoothing() неявно использует mode=constant с cval=0, но мне потребуется что-то вроде mode=reflect.

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

Ответы [ 2 ]

0 голосов
/ 01 мая 2018

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

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

0 голосов
/ 28 апреля 2018

Самый простой способ справиться с этим - удалить среднее значение карты, выполнить сглаживание с помощью 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))<br> return smoothed + offset

Другой вариант, который я могу придумать, - это некоторый итерационный алгоритмчтобы заполнить карту в режиме «отражения» перед сглаживанием, возможно, для реализации в cython или numba, основная проблема заключается в том, насколько сложна ваша граница.Если это просто, как сокращение по широте, то все это легко, поскольку общий случай очень сложен и может иметь много угловых случаев, которые вам нужно обработать:

Идентификация «пограничных слоев»

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

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

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