在Python 2.7中複製IDL'smooth'

Question

我一直在努力研究如何在Python中複製IDL的平滑函數，並且我無法得到類似於相同結果的任何內容。 （免責聲明：從我觸及這種數學問題開始，這可能已經過去了10年，所以它已經被傾銷，讓信息從哪裏找到最便宜的本地燃料）。我試圖代碼此：

smooth(b,w,/nan) 

其中b是含有NANS二維浮陣列（零 - 缺失的數據 - 也被轉換爲NAN）。

從IDL文件，它看起來更平滑使用棚車，所以從scipy.ndimage.filters我曾嘗試：

bsmooth = uniform_filter(b, w) 

我知道這裏有一些根本性的差異：

1. IDL的默認邊緣行爲是「端點從原始數組複製到沒有平滑的結果 」，而I 似乎沒有選擇使用統一過濾器執行此操作。
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]]) 

我的答案與IDL沒有最小的相似之處，不管我是否使用/ nan關鍵字。

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

由於這是不是在Python包中可用，因爲我看到在我的研究過程中多次提出的問題沒有得到滿意的答案，這裏是我如何解決這個問題。

提供了我的功能的測試版本，我要收拾整理。我確信有更好的方法來完成我已經完成的工作，因爲我對Python還是比較陌生 - 請建議進行適當的更改。

情節使用秋季色彩地圖，只是因爲它允許我清楚地看到NaNs。

我的結果：

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 

Answer 2

第一：請使用matplotlib.pyplot.imshow與interpolation="none"這是更好看，也許與灰度。

所以對於你的例子：在scipy和numpy中實際上沒有卷積（filter）把NaN當作缺失值（它們在卷積中傳播它們）。至少我迄今沒有發現任何東西，並且你的邊界處理也是（據我所知）沒有實現。但是之後這個邊界可能會被取代。

如果你想與NaN進行卷積，你可以使用例如astropy.convolution.convolve。有NaNs使用你的過濾器的內核進行插值。但是它們的卷積也有一些缺點：像你想要的邊框處理也沒有在那裏實現，並且你的內核必須是奇形狀，並且你的內核的總和不能爲零（或者非常接近它）

For例如：

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.], 
     [ 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的倍數），但是它確實是你想要的，而且不會搞亂。