python - 在 Python 中使用 Nan 对图像进行高斯滤波

Question

根据二维坐标列表和第三个变量(速度)，我创建了一个覆盖整个采样区域的二维 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 像素。

谢谢大家!!

Answer 1

文字:

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

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

想法是用零替换 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 不会求和到一个了，因此低估了“真实”值(value)。

然而，这可以通过使用第二个辅助数组(参见下面的 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)(理论上支持是无限的，但实际上内核通常被截断，请参阅上面代码示例中的“截断”参数)。在这种情况下，提名符也变为零 (VV=0)，因此 numpy 会抛出“RuntimeWarning:在 true_divide 中遇到无效值”并在相应位置返回 NaN。

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