2017-08-16 16 views
1

在matlab中有函数用于获取每个条目的n个最近条目的平均值。在python中创建matlab平均过滤器

x=np.array([[0.1,0.8,.2], 
      [0.5,0.2,np.nan], 
      [0.7,0.2,0.9], 
      [0.4,0.7,1], 
      [np.nan,0.14,1]]) 

在MATLAB创建平均滤波器:

x=[[0.1,0.8,.2], 
    [0.5,0.2,nan], 
    [0.7,0.2,0.9], 
    [0.4,0.7,1], 
    [nan,0.14,1]] 

fspecial('average',[3,3]) 
filter2(ave1,x) 

[[ 0.17777778   nan   nan] 
[ 0.27777778   nan   nan] 
[ 0.3    nan   nan] 
[  nan   nan 0.43777778] 
[  nan   nan 0.31555556]] 

我想将其转换为蟒蛇。 我发现这一点: uniform filter 和: skimage.filters.rank.mean

,但结果不一样MATLAB。统一的过滤器:

x=np.array([[0.1,0.8,.2], 
      [0.5,0.2,np.nan], 
      [0.7,0.2,0.9], 
      [0.4,0.7,1], 
      [np.nan,0.14,1]]) 

print(uniform_filter(x, size=3, mode='constant')) 

[[ 0.17777778   nan   nan] 
[ 0.27777778   nan   nan] 
[ 0.3    nan   nan] 
[  nan   nan   nan] 
[  nan   nan   nan]] 

skimage过滤器:

from skimage.filters.rank import mean 
from skimage.morphology import square 
from skimage import img_as_float 
x=np.array([[0.1,0.8,.2], 
     [0.5,0.2,np.nan], 
     [0.7,0.2,0.9], 
     [0.4,0.7,1], 
     [np.nan,0.14,1]]) 

print(mean(x, square(3))) 

[[102 76 76] 
[106 102 97] 
[114 130 127] 
[ 90 142 167] 
[ 79 137 181]] 

print(img_as_float(mean(x, square(3)))) 

[[ 0.4   0.29803922 0.29803922] 
[ 0.41568627 0.4   0.38039216] 
[ 0.44705882 0.50980392 0.49803922] 
[ 0.35294118 0.55686275 0.65490196] 
[ 0.30980392 0.5372549 0.70980392]] 

最后我不得不自己,但在性能上还不成熟:

x=np.array([[0.1,0.8,.2], 
      [0.5,0.2,np.nan], 
      [0.7,0.2,0.9], 
      [0.4,0.7,1], 
      [np.nan,0.14,1]]) 



Winsize=3 
adder=int(Winsize/2) 
result=np.zeros_like(x) 
nan_window_index=np.array([]) 
for i in range(x.shape[0]): 
    for j in range(x.shape[1]): 
     top_left_r= int(i-adder) 
     top_left_c= int(j-adder) 
     bottom_right_r=int(i+adder) 
     bottom_right_c=int(j+adder) 
     sum_list=np.array([]) 

     for r_counter in range(top_left_r, bottom_right_r+1): 
      if r_counter<0 or r_counter > x.shape[0]-1: 
       continue 
      for c_counter in range(top_left_c, bottom_right_c+1): 
       if c_counter<0 or c_counter > x.shape[1]-1: 
        continue 
       if not np.isnan(x[r_counter, c_counter]): 
        sum_list=np.append(sum_list, x[r_counter, c_counter]) 
       else: 
        nan_window_index=np.append(nan_window_index, [[r_counter, c_counter]]) 

     result[i,j]= np.sum(sum_list)/(Winsize*Winsize) 

nan_window_index=np.unique(nan_window_index.reshape(int(len(nan_window_index)/2),2), axis=0) 

for i,j in nan_window_index: 
    top_left_r= int(i-adder) 
    top_left_c= int(j-adder) 
    bottom_right_r=int(i+adder) 
    bottom_right_c=int(j+adder) 

    for r_counter in range(top_left_r, bottom_right_r+1): 
     if r_counter<0 or r_counter > x.shape[0]-1: 
      continue 
     for c_counter in range(top_left_c, bottom_right_c+1): 
      if c_counter<0 or c_counter > x.shape[1]-1: 
       continue 
      result[r_counter, c_counter]=np.nan 

print(result) 

,结果是一样的MATLAB:

[[ 0.17777778   nan   nan] 
[ 0.27777778   nan   nan] 
[ 0.3    nan   nan] 
[  nan   nan 0.43777778] 
[  nan   nan 0.31555556]] 

any suggestio ñ更好的表现?

回答

3

可以使用scipy.signal.convolve(或者scipy.signal.convolve2d,因为它可能会更快)为:

import numpy as np 
# from scipy.signal import convolve 
from scipy.signal import convolve2d 

x=np.array([[0.1,0.8,.2], 
      [0.5,0.2,np.nan], 
      [0.7,0.2,0.9], 
      [0.4,0.7,1], 
      [np.nan,0.14,1]]) 

core = np.full((3,3),1/3**2) 

# convolve(x, core, mode='same') 
convolve2d(x, core, mode='same') 

卷积均匀值相同,统一的过滤器。请注意,这会自动“假设”矩阵外的零,但这符合您所要求的内容,所以它将在您当前的设置中起作用。