26

从二维坐标列表和第三个变量（速度）中，我创建了一个覆盖整个采样区域的二维 numpy 数组。我的目的是创建一个图像，其中每个像素都包含位于其中的点的平均速度。之后用高斯滤波器过滤该图像。

问题是该区域的采样不均匀。Nan因此，我在图像中间有几个没有信息（ ）的像素。当我尝试通过高斯滤波器过滤阵列时，Nan传播会破坏整个图像。

我需要过滤这个图像，但拒绝所有没有信息的像素。换句话说，如果一个像素不包含信息，那么过滤就不应该考虑它。

这是我的平均代码示例：

Mean_V = np.zeros([len(x_bins), len(y_bins)])

for i, x_bin in enumerate(x_bins[:-1]):
    bin_x = (x > x_bins[i]) & (x <= x_bins[i+1])
    for j, y_bin in enumerate(y_bins[:-1]):
        bin_xy = (y[bin_x] > y_bins[j]) & (y[bin_x] <= y_bins[j+1])
        if (sum(x > 0 for x in bin_xy) > 0) :
            Mean_V[i,j]=np.mean(V[bin_x][bin_xy])
        else:
            Mean_V[i,j]=np.nan

编辑：

网上冲浪我在2013年做的这个问题就结束了。这个问题的解决方案可以在astropy库中找到：

http://docs.astropy.org/en/stable/convolution/

Astropy 的卷积将 NaN 像素替换为来自其邻居的核加权插值。

谢谢各位！！

4

3 回答 3

37

用一句话来说：

忽略给定数组U中的 NaN 的高斯滤波器可以通过将标准高斯滤波器应用于两个辅助数组V和W并通过取两者的比率来获得结果Z来轻松获得。

这里，V是原始U的副本，其中 NaN 被零替换，W是一个数组，其中零表示 NaN 在原始U中的位置。

这个想法是，用零替换 NaN 会在过滤后的数组中引入一个误差，但是，可以通过将相同的高斯滤波器应用于另一个辅助数组并将两者结合来进行补偿。

在 Python 中：

import numpy as np
import scipy as sp
import scipy.ndimage

sigma=2.0                  # standard deviation for Gaussian kernel
truncate=4.0               # truncate filter at this many sigmas

U=sp.randn(10,10)          # random array...
U[U>2]=np.nan              # ...with NaNs for testing

V=U.copy()
V[np.isnan(U)]=0
VV=sp.ndimage.gaussian_filter(V,sigma=sigma,truncate=truncate)

W=0*U.copy()+1
W[np.isnan(U)]=0
WW=sp.ndimage.gaussian_filter(W,sigma=sigma,truncate=truncate)

Z=VV/WW

数字：

这里高斯滤波器的系数设置为 [0.25,0.50,0.25] 用于演示目的，它们的总和为 0.25+0.50+0.25=1，不失一般性。

在用零替换 NaN 并应用高斯滤波器（参见下面的 VV）后，很明显零会引入错误，即，由于“缺失”数据，系数 0.25+0.50=0.75 不再加起来为 1因此低估了“真实”价值。

然而，这可以通过使用第二个辅助数组（参见下面的 WW）来补偿，在使用相同的高斯滤波之后，它只包含系数的总和。

因此，将两个滤波后的辅助数组相除会重新调整系数，以使它们总和为 1，而忽略 NaN 位置。

array U         1   2   NaN 1   2    
auxiliary V     1   2   0   1   2    
auxiliary W     1   1   0   1   1
position        a   b   c   d   e

filtered VV_b   = 0.25*V_a  + 0.50*V_b  + 0.25*V_c
                = 0.25*1    + 0.50*2    + 0
                = 1.25

filtered WW_b   = 0.25*W_a  + 0.50*W_b  + 0.25*W_c
                = 0.25*1    + 0.50*1    + 0
                = 0.75

ratio Z         = VV_b / WW_b  
                = (0.25*1 + 0.50*2) / (0.25*1    + 0.50*1)
                = 0.333*1 + 0.666*2
                = 1.666

更新 - 除零：

以下内容包含@AndyL 和@amain 来自下方评论的有用问题和答案，谢谢！

当高斯核的支持范围内只有 NaN 条目时，大面积的 NaN 可能会导致某些位置的分母为零（WW=0）（理论上支持是无限的，但实际上内核通常会被截断，请参阅'上面代码示例中的 truncate' 参数）。在这种情况下，提名者也变为零（VV=0），因此 numpy 会抛出“RuntimeWarning：true_divide 中遇到的无效值”并在相应位置返回 NaN。

这可能是最一致/最有意义的结果，如果您可以忍受 numpy 警告，则无需进一步调整。

于 2016-03-30T11:21:06.540 回答
6

enter image description here

I stepped over this question a while ago and used davids answer (thanks!). As it turned out in the meantime, the task of applying a gaussian filter to a array with nans is not as well defined as I thought.

As descibed in ndimage.gaussian_filter, there are different options to process values at the border of the image (reflection, constant extrapolation, ...). A similar decision has to be made for the nan values in the image.

  • Some idea might be, so linarly interpolate nan values, but the question arrises, what to do with nans at the image borders.
  • filter_nan_gaussian_david: Davids approach is equivalent to assuming some mean-neighborhood-value at each nan-point. This leads to a change in the total intensity (see sum value in colum 3), but does a great job otherwise.
  • filter_nan_gaussian_conserving: This approach is to spead the intesity of each point by a gaussian filter. The intensity, which is mapped to nan-pixels is reshifted back to the origin. If this maskes sense, depends on the application. I have a closed area surronded by nans and want to preseve the total intensity + avoid distortions at the boundaries.
  • filter_nan_gaussian_conserving2: Speads intesity of each point by a gaussian filter. The intensity, which is mapped to nan-pixels is redirected to the other pixels with the same Gaussian weighting. This leads to a relative reduction of the intensity at the origin in the vicinity of many nans / border pixels. This is illustrated in the last row very right.

Code

import numpy as np
from scipy import ndimage
import matplotlib as mpl
import matplotlib.pyplot as plt

def filter_nan_gaussian_conserving(arr, sigma):
    """Apply a gaussian filter to an array with nans.

    Intensity is only shifted between not-nan pixels and is hence conserved.
    The intensity redistribution with respect to each single point
    is done by the weights of available pixels according
    to a gaussian distribution.
    All nans in arr, stay nans in gauss.
    """
    nan_msk = np.isnan(arr)

    loss = np.zeros(arr.shape)
    loss[nan_msk] = 1
    loss = ndimage.gaussian_filter(
            loss, sigma=sigma, mode='constant', cval=1)

    gauss = arr.copy()
    gauss[nan_msk] = 0
    gauss = ndimage.gaussian_filter(
            gauss, sigma=sigma, mode='constant', cval=0)
    gauss[nan_msk] = np.nan

    gauss += loss * arr

    return gauss

def filter_nan_gaussian_conserving2(arr, sigma):
    """Apply a gaussian filter to an array with nans.

    Intensity is only shifted between not-nan pixels and is hence conserved.
    The intensity redistribution with respect to each single point
    is done by the weights of available pixels according
    to a gaussian distribution.
    All nans in arr, stay nans in gauss.
    """
    nan_msk = np.isnan(arr)

    loss = np.zeros(arr.shape)
    loss[nan_msk] = 1
    loss = ndimage.gaussian_filter(
            loss, sigma=sigma, mode='constant', cval=1)

    gauss = arr / (1-loss)
    gauss[nan_msk] = 0
    gauss = ndimage.gaussian_filter(
            gauss, sigma=sigma, mode='constant', cval=0)
    gauss[nan_msk] = np.nan

    return gauss

def filter_nan_gaussian_david(arr, sigma):
    """Allows intensity to leak into the nan area.
    According to Davids answer:
        https://stackoverflow.com/a/36307291/7128154
    """
    gauss = arr.copy()
    gauss[np.isnan(gauss)] = 0
    gauss = ndimage.gaussian_filter(
            gauss, sigma=sigma, mode='constant', cval=0)

    norm = np.ones(shape=arr.shape)
    norm[np.isnan(arr)] = 0
    norm = ndimage.gaussian_filter(
            norm, sigma=sigma, mode='constant', cval=0)

    # avoid RuntimeWarning: invalid value encountered in true_divide
    norm = np.where(norm==0, 1, norm)
    gauss = gauss/norm
    gauss[np.isnan(arr)] = np.nan
    return gauss



fig, axs = plt.subplots(3, 4)
fig.suptitle('black: 0, white: 1, red: nan')
cmap = mpl.cm.get_cmap('gist_yarg_r')
cmap.set_bad('r')
def plot_info(ax, arr, col):
    kws = dict(cmap=cmap, vmin=0, vmax=1)
    if col == 0:
        title = 'input'
    elif col == 1:
        title = 'filter_nan_gaussian_conserving'
    elif col == 2:
        title = 'filter_nan_gaussian_david'
    elif col == 3:
        title = 'filter_nan_gaussian_conserving2'
    ax.set_title(title + '\nsum: {:.4f}'.format(np.nansum(arr)))
    ax.imshow(arr, **kws)

sigma = (1,1)

arr0 = np.zeros(shape=(6, 10))
arr0[2:, :] = np.nan
arr0[2, 1:3] = 1

arr1 = np.zeros(shape=(6, 10))
arr1[2, 1:3] = 1
arr1[3, 2] = np.nan

arr2 = np.ones(shape=(6, 10)) *.5
arr2[3, 2] = np.nan

plot_info(axs[0, 0], arr0, 0)
plot_info(axs[0, 1], filter_nan_gaussian_conserving(arr0, sigma), 1)
plot_info(axs[0, 2], filter_nan_gaussian_david(arr0, sigma), 2)
plot_info(axs[0, 3], filter_nan_gaussian_conserving2(arr0, sigma), 3)

plot_info(axs[1, 0], arr1, 0)
plot_info(axs[1, 1], filter_nan_gaussian_conserving(arr1, sigma), 1)
plot_info(axs[1, 2], filter_nan_gaussian_david(arr1, sigma), 2)
plot_info(axs[1, 3], filter_nan_gaussian_conserving2(arr1, sigma), 3)

plot_info(axs[2, 0], arr2, 0)
plot_info(axs[2, 1], filter_nan_gaussian_conserving(arr2, sigma), 1)
plot_info(axs[2, 2], filter_nan_gaussian_david(arr2, sigma), 2)
plot_info(axs[2, 3], filter_nan_gaussian_conserving2(arr2, sigma), 3)

plt.show()
于 2020-04-28T13:23:54.330 回答
0

如何将 Z=VV/WW 替换为 Z=VV/(WW+epsilon) 与 epsilon=0.000001 以自动处理之前提案中没有任何观察的情况

于 2019-01-29T19:32:25.193 回答