Numpy: Sostituisce ogni valore dell'array con la media degli elementi adiacenti

Question

Ho un narray e voglio sostituire ogni valore dell'array con la media degli elementi adiacenti. Il codice qui sotto può fare il lavoro, ma è molto lento quando ho 700 array tutti con forma (7000, 7000), quindi mi chiedo se ci sono modi migliori per farlo. Grazie!

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

Beh, questa è una smoothing operation in image processing, che può essere raggiunto con 2D convoluzione. Stai lavorando in modo un po 'diverso sugli elementi vicini. Quindi, se gli elementi di contorno sono lasciato fuori per la precisione, è possibile utilizzare scipy's convolve2d in questo modo -

from scipy.signal import convolve2d as conv2 

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

Questa operazione specifica è un built-in nel modulo OpenCV come cv2.blur ed è molto efficiente a esso. Il nome descrive fondamentalmente la sua operazione di sfocatura degli array di input che rappresentano le immagini. Credo che l'efficienza derivi dal fatto che internamente è interamente implementato in C per prestazioni con un sottile wrapper Python per gestire gli array NumPy.

Così, l'uscita potrebbe essere alternativamente calcolato con esso, in questo modo -

import cv2 # Import OpenCV module 

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

Ecco uno show-down rapido sul timing su un decentemente grande immagine/array -

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

In seguito alla discussione con @Divakar, trova qui sotto un confronto tra i diversi metodi di convoluzione presenti in 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') 

Tutti e 3 i metodi restituiscono esattamente lo stesso valore. Tuttavia, notare che uniform_filter ha un parametro mode='constant', che indica le condizioni al contorno del filtro e constant == 0 è la stessa condizione al contorno che viene applicata il dominio di Fourier (negli altri 2 metodi). Per diversi casi d'uso è possibile modificare le condizioni al contorno.

Ora, alcune matrici di prova:

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

e alcune temporizzazioni:

%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 

Così, in breve, uniform_filter sembra più veloce, e perché le convolution is separable in due convolutons 1D (simili a gaussian_filter quali è anche separabile).

Altri filtri non separabili con diversi kernel hanno maggiori probabilità di essere più veloci utilizzando la soluzione signal (quella in @Divakar).

La velocità di entrambi fftconvolve e uniform_filter rimane costante per le diverse dimensioni del kernel, mentre convolve2d ottiene leggermente più lento.

Answer 3

Ho avuto un problema simile di recente e ho dovuto trovare una soluzione diversa poiché non posso usare 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)