Numpy：用相邻元素的平均值替换数组中的每个值

Question

我有一个ndarray，我想用相邻元素的平均值替换数组中的每个值。下面的代码可以完成这项工作，但当我有700个形状为（7000，7000）的数组时，它非常缓慢，因此我想知道是否有更好的方法来完成它。谢谢！

a = np.array(([1,2,3,4,5,6,7,8,9],[4,5,6,7,8,9,10,11,12],[3,4,5,6,7,8,9,10,11]))
row,col = a.shape
new_arr = np.ndarray(a.shape)
for x in xrange(row):
    for y in xrange(col):
        min_x = max(0, x-1)
        min_y = max(0, y-1)
        new_arr[x][y] = a[min_x:(x+2),min_y:(y+2)].mean()
print new_arr

- Chiefscreation

---

对我来说，它根本没有运行缓慢，也不像应该运行缓慢。 - Eli Sadoff

@EliSadoff 如果我有700个数组，每个数组的形状都是(7000,7000)，那么就会这样。 - Chiefscreation

我不确定如何在Python中实现，但您是否考虑过多线程或并行处理？我知道在C语言中可以这样做来加速大量数据处理。 - Michael

你没有正确设置 max_x 和 max_y... - Mad Physicist

Answer 1

这是一个关于图像处理中平滑操作的讨论，可以通过2D卷积实现。在处理边缘元素时需要采用不同的方法。因此，如果为了精度而忽略边界元素，可以使用scipy的convolve2d函数，如下所示 -

from scipy.signal import convolve2d as conv2

out = (conv2(a,np.ones((3,3)),'same')/9.0

这个特定操作是OpenCV模块中的一个内置函数，称为cv2.blur，非常高效。其名称基本上描述了对表示图像的输入数组进行模糊处理的操作。我认为其高效性来自于其完全由C实现以提高性能，并使用薄的Python包装器处理NumPy数组。

因此，可以使用它来替代地计算输出，如下所示-

import cv2 # Import OpenCV module

out = cv2.blur(a.astype(float),(3,3))

这里是关于一个相当大的图像/数组计时的快速展示：

In [93]: a = np.random.randint(0,255,(5000,5000)) # Input array

In [94]: %timeit conv2(a,np.ones((3,3)),'same')/9.0
1 loops, best of 3: 2.74 s per loop

In [95]: %timeit cv2.blur(a.astype(float),(3,3))
1 loops, best of 3: 627 ms per loop

- Divakar

---

干得好。我正要写类似的东西。+1。 - rayryeng

@rayryeng 很高兴看到你使用NumPy标签！ ;) - Divakar

@Divakar 它在我的动态中 :) 不过我最近没有回答问题... 我这边有很多事情要处理。 - rayryeng

uniform_filter是这个的快捷方式...但由于signal模块在傅里叶域中操作，所以我猜这会更快？ - Imanol Luengo

@Chiefscreation Scipy 中也有最大值过滤器，链接在这里：http://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.ndimage.filters.maximum_filter.html。 - Divakar

Answer 2

在与@Divakar的讨论后，以下是scipy中不同卷积方法的比较：

import numpy as np
from scipy import signal, ndimage

def conv2(A, size):
    return signal.convolve2d(A, np.ones((size, size)), mode='same') / float(size**2)

def fftconv(A, size):
    return signal.fftconvolve(A, np.ones((size, size)), mode='same') / float(size**2)

def uniform(A, size):
    return ndimage.uniform_filter(A, size, mode='constant')

所有三种方法都返回完全相同的值。但是，请注意，uniform_filter 有一个参数 mode='constant'，它指示过滤器的边界条件，而 constant == 0 是强制傅里叶域（在其他两种方法中）的相同边界条件。对于不同的用例，您可以更改边界条件。

现在是一些测试矩阵：

A = np.random.randn(1000, 1000)

以下是一些时间：

%timeit conv2(A, 3)     # 33.8 ms per loop
%timeit fftconv(A, 3)   # 84.1 ms per loop
%timeit uniform(A, 3)   # 17.1 ms per loop

%timeit conv2(A, 5)     # 68.7 ms per loop
%timeit fftconv(A, 5)   # 92.8 ms per loop
%timeit uniform(A, 5)   # 17.1 ms per loop

%timeit conv2(A, 10)     # 210 ms per loop
%timeit fftconv(A, 10)   # 86 ms per loop
%timeit uniform(A, 10)   # 16.4 ms per loop

%timeit conv2(A, 30)     # 1.75 s per loop
%timeit fftconv(A, 30)   # 102 ms per loop
%timeit uniform(A, 30)   # 16.5 ms per loop

简而言之，uniform_filter 更快，因为在两个1D卷积中（类似于gaussian_filter），卷积是可分离的。

其他具有不同内核的不可分离滤波器更可能使用signal模块（@Divakar的解决方案中的一个）更快。

fftconvolve和uniform_filter的速度对于不同的内核大小保持恒定，而convolve2d会稍微变慢。

- Imanol Luengo

---

好的发现！现在得消化所有这些东西。 - Divakar

@Divakar 刚刚发现了一些非常有趣的事情，对于不同的小内核大小，最后两种方法保持恒定的执行时间。 - Imanol Luengo

我觉得我理解了uniform_filter的双通道实现方式。虽然fftconvolve的内部实现看起来很混乱。再次感谢您提供这些有用的发现！ - Divakar

您IP地址为143.198.54.68，由于运营成本限制，当前对于免费用户的使用频率限制为每个IP每72小时10次对话，如需解除限制，请点击左下角设置图标按钮（手机用户先点击左上角菜单按钮）。 - Imanol Luengo

Answer 3

最近我也遇到了类似的问题，但由于无法使用scipy，我不得不寻找另一种解决方案。

import numpy as np

a = np.random.randint(100, size=(7000,7000)) #Array of 7000 x 7000
row,col = a.shape

column_totals =  a.sum(axis=0) #Dump the sum of all columns into a single array

new_array = np.zeros([row,col]) #Create an receiving array

for i in range(row):
    #Resulting row = the value of all rows minus the orignal row, divided by the row number minus one. 
    new_array[i] = (column_totals - a[i]) / (row - 1)

print(new_array)

- GeorgeLPerkins