我有一个代表灰度图像的 2D numpy 数组。我需要提取该数组中的每个 N x N 子数组，子数组之间具有指定的重叠，并计算平均值、标准差或中位数等属性。
下面的代码执行此任务，但速度很慢，因为它使用 Python for 循环。有关如何向量化此计算或以其他方式加快计算的任何想法？
import numpy as np
img = np.random.randn(100, 100)
N = 4
step = 2
h, w = img.shape
out = []
for i in range(0, h - N, step):
outr = []
for j in range(0, w - N, step):
outr.append(np.mean(img[i:i+N, j:j+N]))
out.append(outr)
out = np.array(out)
【问题讨论】：
标签： python python-3.x numpy vectorization
对于均值和标准差，有一个基于cumsum 的快速解决方案。
以下是 500x200 图像、30x20 窗口和步长 5 和 3 的时间安排。为了比较，我使用 skimage.util.view_as_windows 与 numpy mean 和 std。
mn + sd using cumsum 1.1531693299184553 ms
mn using view_as_windows 3.495307120028883 ms
sd using view_as_windows 21.855629019846674 ms
代码：
import numpy as np
from math import gcd
from timeit import timeit
def wsum2d(A, winsz, stepsz, canoverwriteA=False):
M, N = A.shape
m, n = winsz
i, j = stepsz
for X, x, s in ((M, m, i), (N, n, j)):
g = gcd(x, s)
if g > 1:
X //= g
x //= g
s //= g
A = A[:X*g].reshape(X, g, -1).sum(axis=1)
elif not canoverwriteA:
A = A.copy()
canoverwriteA = True
A[x:] -= A[:-x]
A = A.cumsum(axis=0)[x-1::s]
A = A.T
return A
def w2dmnsd(A, winsz, stepsz):
# combine A and A*A into a complex, so overheads apply only once
M21 = wsum2d(A*(A+1j), winsz, stepsz, True)
M2, mean_ = M21.real / np.prod(winsz), M21.imag / np.prod(winsz)
sd = np.sqrt(M2 - mean_*mean_)
return mean_, sd
# test
np.random.seed(0)
A = np.random.random((500, 200))
wsz = (30, 20)
stpsz = (5, 3)
mn, sd = w2dmnsd(A, wsz, stpsz)
from skimage.util import view_as_windows
Av = view_as_windows(A, wsz, stpsz) # this emits a warning on my system
assert np.allclose(mn, np.mean(Av, axis=(2, 3)))
assert np.allclose(sd, np.std(Av, axis=(2, 3)))
from timeit import repeat
print('mn + sd using cumsum ', min(repeat(lambda: w2dmnsd(A, wsz, stpsz), number=100))*10, 'ms')
print('mn using view_as_windows', min(repeat(lambda: np.mean(Av, axis=(2, 3)), number=100))*10, 'ms')
print('sd using view_as_windows', min(repeat(lambda: np.std(Av, axis=(2, 3)), number=100))*10, 'ms')
【讨论】：
如果 Numba 是一个选项，唯一要做的就是避免列表追加（它也适用于列表追加，但 slower. 为了也使用并行化，稍微重写了实现以避免在范围内的步骤，使用 parfor 时不支持。
示例
@nb.njit(error_model='numpy',parallel=True)
def calc_p(img,N,step):
h,w=img.shape
i_w=(h - N)//step
j_w=(w - N)//step
out = np.empty((i_w,j_w))
for i in nb.prange(0, i_w):
for j in range(0, j_w):
out[i,j]=np.std(img[i*step:i*step+N, j*step:j*step+N])
return out
def calc_n(img,N,step):
h, w = img.shape
out = []
for i in range(0, h - N, step):
outr = []
for j in range(0, w - N, step):
outr.append(np.std(img[i:i+N, j:j+N]))
out.append(outr)
return(np.array(out))
时间
所有计时都没有大约0.5s的编译开销（对函数的第一次调用不计入计时）。
#Data
img = np.random.randn(100, 100)
N = 4
step = 2
calc_n :17ms
calc_p :0.033ms
因为这实际上是一个滚动平均值，如果N 变大，还有进一步改进的空间。
【讨论】：
你可以使用 scikit-image block_reduce:
所以你的代码变成了：
import numpy as np
import skimage.measure
N = 4
# Your main array
a = np.arange(9).reshape(3,3)
mean = skimage.measure.block_reduce(a, (N,N), np.mean)
std_dev = skimage.measure.block_reduce(a, (N,N), np.std)
median = skimage.measure.block_reduce(a, (N,N), np.median)
但是，上述代码仅适用于大小为 1 的步幅/步数。
对于均值，您可以使用任何现代 ML 软件包中都可用的均值池。至于中位数和标准差，这似乎是正确的方法。
【讨论】：
N，而不是数组下采样了step。
一般情况可以使用scipy.ndimage.generic_filter解决：
import numpy as np
from scipy.ndimage import generic_filter
img = np.random.randn(100, 100)
N = 4
filtered = generic_filter(img.astype(np.float), np.std, size=N)
step = 2
output = filtered[::step, ::step]
但是，这实际上可能不会比简单的 for 循环快多少。
要应用均值和中值滤波器，您可以分别使用skimage.rank.mean 和skimage.rank.median，这应该会更快。还有scipy.ndimage.median_filter。否则，也可以通过与值为 1./N^2 的 (N, N) 数组进行简单卷积来有效计算均值。对于标准偏差，您可能不得不硬着头皮使用generic_filter，除非您的步长大于或等于 N。
【讨论】：