在 Python 2.7 中复制 IDL 'smooth'

Question

我一直在尝试研究如何在 Python 中复制 IDL 的平滑功能，但我无法得到类似的结果。 （免责声明：自从我接触这种数学问题以来可能已经有 10 年了，所以它已经被丢弃，以便为诸如在哪里可以找到最便宜的当地燃料之类的信息让路）。我正在尝试对此进行编码：

smooth(b,w,/nan)

其中 b 是包含 NAN 的 2D 浮点数组（零 - 缺失数据 - 也已转换为 NAN）。

从 IDL 文档来看，使用 boxcar 似乎很顺利，所以从 scipy.ndimage.filters 我尝试过：

bsmooth = uniform_filter(b, w)

我知道这里有一些根本区别：

1. IDL 的默认边缘行为是“复制端点 从原始数组到没有平滑的结果”而我 似乎没有使用统一过滤器执行此操作的选项。
2. NaN 元素的处理。在 IDL 中，/nan 关键字似乎 意味着在可能的情况下，NaN 值将由结果填充 窗口中的其他点。如果没有有效的积分 通过 MISSING 关键字生成结果。我以为我可以 在使用平滑之后近似此行为 scipy.interpolate 的 NearestNDInterpolator （感谢辉煌 亚历克斯在这里的解释： filling gaps on an image using numpy and scipy)

这是我的测试数组：

 >>>b                                                                                                                array([[ 0.97599638,  0.93114936,  0.87070072,  0.5379253 ],                                                              
       [ 0.34873217,         nan,  0.40985891,  0.22407863],                                                              
       [        nan,         nan,         nan,  0.67532134],                                                              
       [        nan,         nan,  0.85441768,         nan]])  

无论我是否使用 /nan 关键字，我的回答都与 IDL 没有一点相似之处。

IDL> smooth(b,2,/nan)
      0.97599638      0.93114936      0.87070072      0.53792530
      0.34873217      0.70728749      0.60817236      0.22407863
             NaN      0.53766960      0.54091913      0.67532134
             NaN             NaN      0.85441768             NaN

IDL> smooth(b,2)
      0.97599638      0.93114936      0.87070072      0.53792530
      0.34873217            -NaN            -NaN      0.22407863
            -NaN            -NaN            -NaN      0.67532134
            -NaN            -NaN      0.85441768             NaN

我承认我发现 scipy 文档在细节上相当稀疏，所以我不知道我是否真的在做我认为我在做的事情。我认为这两种 python 方法都能平滑图像给出不同答案的事实表明，事情并不是我理解的那样。

>>>uniform_filter(b, 2)
array([[ 0.97599638,  0.95357287,  0.90092504,  0.70431301],
       [ 0.66236428,         nan,         nan,         nan],
       [        nan,         nan,         nan,         nan],
       [        nan,         nan,         nan,         nan]])    

我觉得它太空了有点奇怪，所以我尝试了一个包含 100 个元素的数组（仍然使用 2 个窗口）并输出图像。结果（第一张图片是“b”第二张是“bsmooth”）并不是我所希望的：

回到较小的数组并按照以下示例进行操作：http://scipy.github.io/old-wiki/pages/Cookbook/SignalSmooth，我认为它会提供与 uniform_filter 相同的输出，我尝试了：

>>> box = np.array([1,1,1,1])
>>> box = box.reshape(2,2)
>>> box
array([[1, 1],
       [1, 1]])
>>> bsmooth = scipy.signal.convolve2d(b,box,mode='same')
>>> print bsmooth
[[ 0.97599638  1.90714574  1.80185008  1.40862602]
 [ 1.32472855         nan         nan  2.04256356]
 [        nan         nan         nan         nan]
 [        nan         nan         nan         nan]]

显然我完全误解了 scipy 函数，甚至可能是 IDL 函数。如果有人能帮助我尽可能地复制 IDL 平滑功能，我将不胜感激。我面临着相当大的时间压力来获得一个不依赖于 IDL 的解决方案，我正在掷硬币来决定是从头开始编写函数还是开发一种传染性很强的疾病。

如何在 python 中进行同样的平滑处理？

Answer 1

首先：请使用matplotlib.pyplot.imshow 和interpolation="none"，这样看起来更好看，可能还带有灰度。

因此，对于您的示例：在 scipy 和 numpy 中实际上没有将 NaN 视为缺失值的卷积（过滤器）（它们在卷积中传播它们）。至少到目前为止我还没有发现，而且（据我所知）您的边界处理也没有实施。但是边界可以在之后被替换。

如果您想使用NaN 进行卷积，例如可以使用astropy.convolution.convolve。使用过滤器的内核对 NaN 进行插值。但是他们的卷积也有一些缺点：那里也没有实现你想要的边界处理，你的内核必须是奇数形状，你的内核之和不能为零（或非常接近它）

例如：

from astropy.convolution import convolve
import numpy as np
array = np.random.uniform(10,100, (4,4))
array[1,1] = np.nan
kernel = np.ones((3,3))
convolve(array, kernel)

例如一个初始数组

array([[ 97.19514587,  62.36979751,  93.54811286,  30.23567842],
       [ 51.02184613,          nan,  46.14769821,  60.08088041],
       [ 20.86482452,  42.39661484,  36.96961278,  96.89180175],
       [ 45.54453509,  76.61274347,  46.44485141,  25.40985372]])

会变成：

array([[ 266.9009961 ,  406.59680717,  348.69637399,  230.01236989],
       [ 330.16243546,  506.82785931,  524.95440336,  363.87378443],
       [ 292.75477064,  422.31693304,  487.26826319,  311.94469828],
       [ 185.41871792,  268.83318211,  324.72547798,  205.71611967]])

如果你想“标准化”它，astropy 提供了normalize_kernel 参数：

convolved = convolve(array, kernel, normalize_kernel=True)
array([[ 29.58753936,  42.09982189,  49.31793529,  33.00203873],
       [ 49.87040638,  65.67695002,  66.10447436,  40.44026448],
       [ 52.51126383,  63.03914444,  60.85474739,  35.88011742],
       [ 39.40188443,  46.82350749,  40.1380926 ,  22.46090152]])

如果你想用原始数组中的值替换“边缘”值，只需替换它们：

convolved[0,:] = array[0,:]
convolved[-1,:] = array[-1,:]
convolved[:,0] = array[:,0]
convolved[:,-1] = array[:,-1]

这就是现有软件包所提供的（据我所知）。如果您想学习一点Cython 或numba，您可以轻松编写自己的卷积，该卷积不会比 numpy/scipy 慢多少（仅 2-10 倍），但完全符合您的要求而不会弄乱周围。

Answer 2

由于这不是 python 包中可用的东西，而且因为我在研究过程中多次看到这个问题没有令人满意的答案，所以这就是我解决问题的方法。

提供的是我要整理的功能的测试版本。我相信会有更好的方法来完成我所做的事情，因为我对 Python 还是很陌生 - 请推荐任何适当的更改。

绘图使用秋季颜色图只是因为它让我可以清楚地看到 NaN。

我的结果：

IDL propagate
     0.033369284     0.067915268      0.96602046      0.85623550
      0.30435592             NaN             NaN       100.00000
      0.94065958             NaN             NaN      0.90966976
     0.018516513     0.044460904     0.051047217             NaN

python propagate
[[  3.33692829e-02   6.79152655e-02   9.66020487e-01   8.56235492e-01]
 [  3.04355923e-01              nan              nan   1.00000000e+02]
 [  9.40659566e-01              nan              nan   9.09669768e-01]
 [  1.85165123e-02   4.44609040e-02   5.10472165e-02              nan]]


IDL replace
     0.033369284     0.067915268      0.96602046      0.85623550
      0.30435592      0.47452110       14.829881       100.00000
      0.94065958      0.33833817       17.002417      0.90966976
     0.018516513     0.044460904     0.051047217             NaN

python replace
[[  3.33692829e-02   6.79152655e-02   9.66020487e-01   8.56235492e-01]
 [  3.04355923e-01   4.74521092e-01   1.48298812e+01   1.00000000e+02]
 [  9.40659566e-01   3.38338177e-01   1.70024175e+01   9.09669768e-01]
 [  1.85165123e-02   4.44609040e-02   5.10472165e-02              nan]]

我的功能：

#!/usr/bin/env python

# smooth.py  
__version__ = 0.1

# Version 0.1    29 Feb 2016    ELH    Test release

import numpy as np
import matplotlib.pyplot as mp

def Smooth(v1, w, nanopt): 

    # v1 is the input 2D numpy array.
    # w is the width of the square window along one dimension
    # nanopt can be replace or propagate 

    '''
    v1 = np.array(
    [[3.33692829e-02, 6.79152655e-02, 9.66020487e-01, 8.56235492e-01], 
    [3.04355923e-01, np.nan        , 4.86013025e-01, 1.00000000e+02],
    [9.40659566e-01, 5.23314093e-01, np.nan        , 9.09669768e-01],
    [1.85165123e-02, 4.44609040e-02, 5.10472165e-02, np.nan ]])
    w = 2
    '''

    mp.imshow(v1, interpolation='None', cmap='autumn')
    mp.show()

    # make a copy of the array for the output:
    vout=np.copy(v1)

    # If w is even, add one
    if w % 2 == 0:
        w = w + 1

    # get the size of each dim of the input:
    r,c = v1.shape

    # Assume that w, the width of the window is always square.
    startrc = (w - 1)/2
    stopr = r - ((w + 1)/2) + 1
    stopc = c - ((w + 1)/2) + 1

    # For all pixels within the border defined by the box size, calculate the average in the window.
    # There are two options:
    # Ignore NaNs and replace the value where possible.
    # Propagate the NaNs

    for col in range(startrc,stopc):
        # Calculate the window start and stop columns
        startwc = col - (w/2) 
        stopwc = col + (w/2) + 1
        for row in range (startrc,stopr):
            # Calculate the window start and stop rows
            startwr = row - (w/2)
            stopwr = row + (w/2) + 1
            # Extract the window
            window = v1[startwr:stopwr, startwc:stopwc]
            if nanopt == 'replace':
                # If we're replacing Nans, then select only the finite elements
                window = window[np.isfinite(window)]
            # Calculate the mean of the window
            vout[row,col] = np.mean(window)

    mp.imshow(vout, interpolation='None', cmap='autumn')
    mp.show()

    return vout