用Python在Nan中高斯滤波图像

Question

从2D坐标列表和第三个变量(速度),我创建了一个覆盖整个采样区域的2D 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并且通过采用两者的比率来获得结果Z,可以容易地获得忽略给定阵列U中的 NaN的高斯滤波器.

这里,V是原始U的副本,其中NaNs被零替换,W是一个零的数组,其中零指示原始U中NaN的位置.

想法是用零替换NaN会在滤波后的阵列中引入误差,然而,可以通过将相同的高斯滤波器应用于另一个辅助阵列并将两者合并来补偿.

在Python中:

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

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=2.0)

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

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

Answer 2

不久前我跨过这个问题并使用了davids 答案（谢谢！）。与此同时，事实证明，将高斯滤波器应用于具有 nans 的数组的任务并不像我想象的那样明确。

正如ndimage.gaussian_filter 中所述，有不同的选项可以处理图像边界处的值（反射、常数外推等）。必须对图像中的 nan 值做出类似的决定。

• 一些想法可能是线性插值 nan 值，但问题出现了，如何处理图像边界处的 nan 。
• filter_nan_gaussian_david：Davids 方法相当于在每个 nan 点假设一些均值邻域值。这会导致总强度发生变化（参见第sum3 列中的值），但在其他方面做得很好。
• filter_nan_gaussian_conserving：这种方法是通过高斯滤波器来扩展每个点的强度。映射到纳米像素的强度被重新移回原点。如果此掩码有意义，则取决于应用程序。我有一个由 nans 包围的封闭区域，并希望保留总强度+避免边界处的扭曲。
• filter_nan_gaussian_conserving2：通过高斯滤波器扩展每个点的强度。映射到纳米像素的强度被重定向到具有相同高斯权重的其他像素。这导致许多纳米/边界像素附近原点处的强度相对降低。这在最后一行非常正确地说明了。

代码

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:
        /sf/answers/2541510401/
    """
    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()