循环遍历2D数组的有效方法

Question

我有一个2D数组2000x4000，对于该数组中的每个单元，我必须将该单元的值与10个相邻单元（在+/- X和+/- Y中）制作的蒙版的标准偏差进行比较。

例如，这就是我现在正在做的事情：

import numpy as np
from astropy.stats import sigma_clipped_stats

BPmap=[]
N=10
a=np.random.random((2000,4000))
for row in range(N,a.shape[0]-N):
    BPmap_row=[]
    for column in range(N,a.shape[1]-N):
        Bpmap_data=np.array(a[row-N:row+N,column-N:column+N].ravel()) 
        mean, median, std = sigma_clipped_stats(Bpmap_data, sigma=3,iters=5)
        BPmap_Nsigma=float(a[row][column]-median)/std                 
        BPmap_row.append(BPmap_Nsigma)
    BPmap.append(BPmap_row)

这有一个明显的问题，就是我正在执行2000x4000 = 8000000循环，这花费了很长时间。我需要找到一种非常有效的方法来执行这些操作，但是我不知道如何做。

Answer 1

我们通常避免对numpy使用double-for循环；它很慢，而且有了智能索引（array [:, iN] ...），我们可以在一个循环中完成很多工作。
但是对于您的卷积问题，这可能是最简单的~~（并且只有？）~~方式来完成您想要的事情。（编辑：不是。请参见下面的@Masoud的答案）。

Bpmap_data=np.array(a[row-N:row+N,column-N:column+N].ravel())在每个循环中创建一个新数组。
在不创建新数组的情况下（即直接使用numpy视图）计算中值和标准差会更快。

实际上，速度要快40倍（在我的Google colab实例上）

要获得比以前快400倍的算法，请查看对2D数组使用scipy过滤器的@Masoud答案。

import numpy as np
from astropy.stats import sigma_clipped_stats


N=10
a=np.random.random((80,40))


def f():
  """Your original code"""
  BPmap=[]
  for row in range(N,a.shape[0]-N):
      BPmap_row=[]
      for column in range(N,a.shape[1]-N):
          Bpmap_data=np.array(a[row-N:row+N,column-N:column+N].ravel()) 
          mean, median, std = sigma_clipped_stats(Bpmap_data, sigma=3,iters=5)
          BPmap_Nsigma=float(a[row][column]-median)/std                 
          BPmap_row.append(BPmap_Nsigma)
      BPmap.append(BPmap_row)
  return BPmap

def f2():
  """this little guy is improving a lot your work"""
  BPmap=[]
  for row in range(N,a.shape[0]-N):
      BPmap_row=[]
      for column in range(N,a.shape[1]-N):
          # the next 3 lines do not need any more memory
          view = a[row-N:row+N,column-N:column+N]
          std_without_outliers = view[view - view.mean() < 3*view.std()].std()
          median = np.median(view)
          # back to your workflow
          BPmap_Nsigma=float(a[row][column]-median)/std_without_outliers                 
          BPmap_row.append(BPmap_Nsigma)
      BPmap.append(BPmap_row)
  return BPmap

%time _ = f()
%time _ = f2()

f() == f2()

>>>CPU times: user 39.7 s, sys: 14.2 ms, total: 39.7 s
Wall time: 39.7 s
CPU times: user 969 ms, sys: 2.99 ms, total: 972 ms
Wall time: 976 ms
True

编辑
实际上，sigma_clipped_stats(a[row-N:row+N,column-N:column+N])确实会减慢循环速度。我怀疑sigma_clipped_stats会创建其参数的副本。

我从3 sigma切割中消除了异常值后采取了std

我在这里展示了一种使用纯numpy做到这一点的方法；这确实比以前使用的功能要快。

最后，f（）= f2（）那么为什么再使用此astropy函数呢？

Answer 2

代码中存在一些降低性能的问题：

1. 如上所述这里，避免使用for循环。
2. 您实际上是将每个数字平方10 * 10次。

代替for循环，您可以使用Scipy.ndimage和opencv图书馆员执行卷积。尽管这些库用于图像处理，但它们对于处理任何2D数组非常有效。这是一个代码，执行与使用Scipy.ndimage工具所需的任务相同的代码，但是速度提高了1000倍（对于200X400阵列，该速度为23ms vs 27s）。我使用了此处提供的算法来计算标准差：

import numpy as np
from scipy.ndimage.filters import uniform_filter, median_filter

a=np.random.random((200,400))
c1 = uniform_filter(a, size = (10,10))
c2 = uniform_filter(a*a, size = (10,10))
std = ((c2 - c1*c1)**.5)
med = median_filter(a, size=(10, 10))

BPmap = (a - med)/std