python/NumPy metodo più veloce per 2g filtraggio kernel rango su array mascherati (e/o posizionamento selettiva)

Question

dato un 2D NumPy gamma

MyArray = np.array([[ 8.02, 9.54, 0.82, 7.56, 2.26, 9.47], 
      [ 2.68, 7.3 , 2.74, 3.03, 2.25, 8.84], 
      [ 2.21, 3.62, 0.55, 2.94, 5.77, 0.21], 
      [ 5.78, 5.72, 8.85, 0.24, 5.37, 9.9 ], 
      [ 9.1 , 7.21, 4.14, 9.95, 6.73, 6.08], 
      [ 1.8 , 5.14, 5.02, 6.52, 0.3 , 6.11]]) 

e una serie maschera

MyMask = np.array([[ 0., 0., 1., 1., 0., 1.], 
      [ 1., 0., 0., 0., 0., 1.], 
      [ 0., 0., 0., 1., 0., 0.], 
      [ 0., 1., 1., 1., 1., 0.], 
      [ 0., 1., 0., 1., 0., 0.], 
      [ 0., 1., 0., 0., 1., 1.]]) 

voglio correre un filtro mediano "bucato" che ignora gli elementi mascherati.

Ad esempio, un filtro rango con un kernel

k = np.array([[ 1, 1, 1], 
       [ 1, 0, 1], 
       [ 1, 1, 1]]); 

porrebbe su : classificare la zona definita dal kernel per ogni elemento di e restituire la mediana dei soli elementi non mascherati (averaging se l'array è un numero pari).

Ora, attualmente lo sto facendo in loop unptonici, usando bottleneck.nanmedian mappando la maschera ai NaN. Questo mi sta dando esattamente quello di cui ho bisogno, ma speravo di fare affidamento su routine di manipolazione di array 2D.

scipy.signal.order_filter e scipy.ndimage.filters.rank_filter bothavailable (rank_filter sembra essere molto più veloce), ma risulta che ordinare NaN e Inf nella parte superiore della matrice prima di restituire il rango e polarizzando il risultato. Sembra che nessuno di questi metodi supporti gli array numpy.ma (masking), né accettano una schiera di ranghi selettivi (quindi potrei riempire tutte le maschere con 0 e compensare il mio rank), né c'è un modo ovvio per variare il kernel per ogni posizione.

Mi chiedo se mi sia mancata una combinazione e/o una funzione Python, o se dovrei cercare di implementare una nuova routine in Cython.

Ignorando movimentazione confine, i punti interni del problema sopra sarebbero

[[ 0.  0.  0.  0.  0.  0. ] 
[ 0.  3.18 3.62 2.26 2.645 0. ] 
[ 0.  2.74 3.325 2.74 2.64 0. ] 
[ 0.  3.88 3.62 4.955 6.08 0. ] 
[ 0.  5.02 5.77 5.77 6.52 0. ] 
[ 0.  0.  0.  0.  0.  0. ]] 

Answer 1

Un modo è quello di sacrificare utilizzo di RAM rinunciare al Python loop. Cioè facciamo esplodere l'array originale in modo da poter applicare il filtro su tutti i sub-array in una sola volta. Che è simile a Numpy broadcasting.

Per un array di 1000x1000 la funzione vettorizzata esegue circa 100 volte più veloce, nei miei test.

Nel mio codice ho usato il mascheramento NaN, ma con alcune linee aggiuntive di codice è possibile utilizzare anche gli array numpy.ma. E non avevo una funzione nanmedian quindi ho usato nanmean, tuttavia le prestazioni dovrebbero essere paragonabili.

import numpy as np 
from numpy.lib.stride_tricks import as_strided 

# test data 
N = 1000 
A = np.random.rand(N, N)*10 
mask = np.random.choice([True, False], size=(N, N)) 

def filter_loop(A, mask): 
    kernel = np.array([[1,1,1],[1,0,1],[1,1,1]], bool) 
    A = A.copy() 
    A[mask] = np.nan 
    N = A.shape[0] - 2 # assuming square matrix 
    out = np.empty((N, N)) 
    for i in xrange(N): 
     for j in xrange(N): 
      out[i,j] = np.nanmean(A[i:i+3, j:j+3][kernel]) 
    return out  

def filter_broadcast(A, mask): 
    A = A.copy() 
    A[mask] = np.nan 
    N = A.shape[0] - 2 
    B = as_strided(A, (N, N, 3, 3), A.strides+A.strides) 
    B = B.copy().reshape((N, N, 3*3)) 
    B[:,:,4] = np.nan 
    return np.nanmean(B, axis=2)