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

Answer 1

Ну, это smoothing operation in image processing, что может быть достигнуто с 2D свертка. Вы немного по-другому работаете над приграничными элементами. Итак, если граничные элементы выпущены для точности, вы можете использовать scipy's 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

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')

Все 3 метода возвращают одно и то же значение. Тем не менее, обратите внимание, что 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 кажется быстрее, и это потому, что свертка делится на два одномерных сверточника (аналог gaussian_filter, который также является сепарабельным).

Другие неразделимые фильтры с разными ядрами, скорее всего, будут быстрее использовать signal модуль (тот, что в @Divakar's) решение.

Скорость обоих fftconvolve а также uniform_filter остается постоянным для разных размеров ядра, в то время как convolve2d становится немного медленнее.

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)