在 Python 2.7 中复制 IDL 'smooth'
Replicate IDL 'smooth' in Python 2.7
我一直在努力研究如何在 Python 中复制 IDL 的平滑函数，但我无法得到类似的结果。 （免责声明：我接触这种数学问题可能已经有 10 年了，所以它已被丢弃，以便为诸如在哪里可以找到最便宜的当地燃料之类的信息让路）。我正在尝试编写代码：
smooth(b,w,/nan)
其中 b 是包含 NAN 的二维浮点数组（零 - 缺失数据 - 也已转换为 NAN）。
从 IDL 文档来看，使用 boxcar 似乎很顺利，所以从 scipy.ndimage.filters 我试过：
bsmooth = uniform_filter(b, w)
我知道这里有一些根本的区别：
- IDL 的默认边行为是“端点被复制
从原始数组到没有平滑的结果”，而我
似乎没有使用统一过滤器执行此操作的选项。
- NaN 元素的处理。在 IDL 中，/nan 关键字似乎
意味着在可能的情况下 NaN 值将由结果填充
window 中的其他点。如果没有有效的点
通过 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 个元素的数组（仍然使用 window 为 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 smooth 函数，我将不胜感激。我在相当大的时间压力下要为此找到一个不依赖 IDL 的解决方案，我正在掷硬币来决定是从头开始编写函数代码还是患上传染性很强的疾病。
如何在 python 中执行相同的平滑处理？
首先：请将 matplotlib.pyplot.imshow 与 interpolation="none" 一起使用，这样看起来更好看，也许还可以使用灰度。
因此对于您的示例：scipy 和 numpy 中实际上没有卷积（过滤器）将 NaN 视为缺失值（它们在卷积中传播它们）。至少到目前为止我发现 none 并且您的边界处理（据我所知）也没有实施。但边界可以在之后被替换。
如果你想用 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]])
如果你想"normalize"它，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]])
如果您想用原始数组中的值替换 "edge" 值，只需替换它们：
convolved[0,:] = array[0,:]
convolved[-1,:] = array[-1,:]
convolved[:,0] = array[:,0]
convolved[:,-1] = array[:,-1]
这就是现有软件包所提供的（据我所知）。如果你想学习一点 Cython 或 numba 你可以很容易地编写你自己的卷积，它不会比 numpy/scipy 慢很多（只有 2-10 倍）但是完全是你想要的而不是乱七八糟的。
因为这不是 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
我一直在努力研究如何在 Python 中复制 IDL 的平滑函数，但我无法得到类似的结果。 （免责声明：我接触这种数学问题可能已经有 10 年了，所以它已被丢弃，以便为诸如在哪里可以找到最便宜的当地燃料之类的信息让路）。我正在尝试编写代码：
smooth(b,w,/nan)
其中 b 是包含 NAN 的二维浮点数组（零 - 缺失数据 - 也已转换为 NAN）。
从 IDL 文档来看，使用 boxcar 似乎很顺利，所以从 scipy.ndimage.filters 我试过：
bsmooth = uniform_filter(b, w)
我知道这里有一些根本的区别：
- IDL 的默认边行为是“端点被复制 从原始数组到没有平滑的结果”，而我 似乎没有使用统一过滤器执行此操作的选项。
- NaN 元素的处理。在 IDL 中，/nan 关键字似乎 意味着在可能的情况下 NaN 值将由结果填充 window 中的其他点。如果没有有效的点 通过 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 个元素的数组（仍然使用 window 为 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 smooth 函数，我将不胜感激。我在相当大的时间压力下要为此找到一个不依赖 IDL 的解决方案，我正在掷硬币来决定是从头开始编写函数代码还是患上传染性很强的疾病。
如何在 python 中执行相同的平滑处理？
首先：请将 matplotlib.pyplot.imshow 与 interpolation="none" 一起使用，这样看起来更好看，也许还可以使用灰度。
因此对于您的示例：scipy 和 numpy 中实际上没有卷积（过滤器）将 NaN 视为缺失值（它们在卷积中传播它们）。至少到目前为止我发现 none 并且您的边界处理（据我所知）也没有实施。但边界可以在之后被替换。
如果你想用 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]])
如果你想"normalize"它，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]])
如果您想用原始数组中的值替换 "edge" 值，只需替换它们：
convolved[0,:] = array[0,:]
convolved[-1,:] = array[-1,:]
convolved[:,0] = array[:,0]
convolved[:,-1] = array[:,-1]
这就是现有软件包所提供的（据我所知）。如果你想学习一点 Cython 或 numba 你可以很容易地编写你自己的卷积，它不会比 numpy/scipy 慢很多（只有 2-10 倍）但是完全是你想要的而不是乱七八糟的。
因为这不是 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