Uso de zancadas para un filtro de media móvil eficiente

Recientemente aprendí sobre avances en la respuesta a esta publicación , y me preguntaba cómo podría usarlos para calcular un filtro de promedio móvil de manera más eficiente que lo que propuse en esta publicación (utilizando filtros de convolución).

Esto es lo que tengo hasta ahora. Toma una vista de la matriz original, luego la enrolla por la cantidad necesaria y sum los valores del kernel para calcular el promedio. Soy consciente de que los bordes no se manejan correctamente, pero puedo ocuparme de eso después … ¿Existe una forma mejor y más rápida? El objective es filtrar grandes matrices de punto flotante de hasta 5000×5000 x 16 capas, una tarea en la que scipy.ndimage.filters.convolve es bastante lento.

Tenga en cuenta que estoy buscando conectividad de 8 vecinos, es decir, un filtro de 3×3 toma el promedio de 9 píxeles (8 alrededor del píxel focal) y asigna ese valor al píxel en la nueva imagen.

 import numpy, scipy filtsize = 3 a = numpy.arange(100).reshape((10,10)) b = numpy.lib.stride_tricks.as_strided(a, shape=(a.size,filtsize), strides=(a.itemsize, a.itemsize)) for i in range(0, filtsize-1): if i > 0: b += numpy.roll(b, -(pow(filtsize,2)+1)*i, 0) filtered = (numpy.sum(b, 1) / pow(filtsize,2)).reshape((a.shape[0],a.shape[1])) scipy.misc.imsave("average.jpg", filtered) 

EDITAR Aclaración sobre cómo veo este trabajo:

Código actual:

  1. use stride_tricks para generar una matriz como [[0,1,2], [1,2,3], [2,3,4] …] que corresponde a la fila superior del núcleo de filtro.
  2. Desplácese a lo largo del eje vertical para obtener la fila central del kernel [[10,11,12], [11,12,13], [13,14,15] …] y agréguela a la matriz que obtuve 1)
  3. Repita para obtener la fila inferior del kernel [[20,21,22], [21,22,23], [22,23,24] …]. En este punto, tomo la sum de cada fila y la divido por el número de elementos en el filtro, dándome el promedio de cada píxel, (desplazado por 1 fila y 1 col, y con algunas rarezas alrededor de los bordes, pero puedo cuidar de eso mas tarde)

Lo que esperaba era un mejor uso de stride_tricks para obtener los 9 valores o la sum de los elementos del núcleo directamente, para toda la matriz, o que alguien pueda convencerme de otro método más eficiente …

Para lo que vale, aquí es cómo lo harías usando trucos de “fancy”. Iba a publicar esto ayer, ¡pero me distraje por el trabajo real! 🙂

@Paul y @eat tienen implementaciones agradables que utilizan otras formas de hacerlo. Solo para continuar con las cosas de la pregunta anterior, pensé que publicaría el equivalente en N-dimensional.

Sin embargo, no podrá vencer de forma significativa las funciones de scipy.ndimage para matrices> 1D. ( scipy.ndimage.uniform_filter debería vencer a scipy.ndimage.convolve , sin embargo)

Además, si está intentando obtener una ventana móvil multidimensional, corre el riesgo de explotar el uso de la memoria cada vez que realiza una copia de su matriz sin darse cuenta. Si bien la matriz inicial “rodante” es solo una vista en la memoria de su matriz original, cualquier paso intermedio que copie la matriz hará que una copia sea de órdenes de magnitud mayor que la matriz original (es decir, digamos que está trabajando con una matriz original de 100×100 … La vista en ella (para un tamaño de filtro de (3,3)) será 98x98x3x3 pero usará la misma memoria que la original. Sin embargo, cualquier copia usará la cantidad de memoria que una matriz completa de 98x98x3x3 ¡¡haría!!)

Básicamente, el uso de trucos de zancada loca es ideal cuando se desea vectorizar operaciones de ventanas en movimiento en un solo eje de un ndarray. Hace que sea realmente fácil calcular cosas como una desviación estándar móvil, etc. con muy poca sobrecarga. Cuando quieres comenzar a hacer esto a lo largo de múltiples ejes, es posible, pero normalmente estás mejor con funciones más especializadas. (Como scipy.ndimage , etc.)

En cualquier caso, así es como lo haces:

 import numpy as np def rolling_window_lastaxis(a, window): """Directly taken from Erik Rigtorp's post to numpy-discussion. """ if window < 1: raise ValueError, "`window` must be at least 1." if window > a.shape[-1]: raise ValueError, "`window` is too long." shape = a.shape[:-1] + (a.shape[-1] - window + 1, window) strides = a.strides + (a.strides[-1],) return np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides) def rolling_window(a, window): if not hasattr(window, '__iter__'): return rolling_window_lastaxis(a, window) for i, win in enumerate(window): if win > 1: a = a.swapaxes(i, -1) a = rolling_window_lastaxis(a, win) a = a.swapaxes(-2, i) return a filtsize = (3, 3) a = np.zeros((10,10), dtype=np.float) a[5:7,5] = 1 b = rolling_window(a, filtsize) blurred = b.mean(axis=-1).mean(axis=-1) 

Entonces, lo que obtenemos cuando hacemos b = rolling_window(a, filtsize) es una matriz de 8x8x3x3, que en realidad es una vista en la misma memoria que la matriz original de 10×10. Podríamos haber utilizado tan fácilmente diferentes tamaños de filtro a lo largo de diferentes ejes u operado solo a lo largo de ejes seleccionados de una matriz N-dimensional (es decir, filtsize = (0,3,0,3) en una matriz de 4 dimensiones nos daría una 6 ver).

Luego podemos aplicar una función arbitraria al último eje repetidamente para calcular efectivamente las cosas en una ventana móvil.

Sin embargo, debido a que estamos almacenando arreglos temporales que son mucho más grandes que nuestro arreglo original en cada paso de mean (o std o lo que sea), ¡esto no es en absoluto eficiente en memoria! Tampoco va a ser terriblemente rápido, tampoco.

El equivalente para ndimage es solo:

 blurred = scipy.ndimage.uniform_filter(a, filtsize, output=a) 

Esto manejará una variedad de condiciones de contorno, hará el “desenfoque” en el lugar sin requerir una copia temporal de la matriz y será muy rápido. Los trucos para caminar son una buena manera de aplicar una función a una ventana en movimiento a lo largo de un eje, pero no son una buena manera de hacerlo a lo largo de varios ejes, por lo general …

Solo mis $ 0.02, al menos …

No estoy lo suficientemente familiarizado con Python para escribir código para eso, pero las dos mejores formas de acelerar las circunvoluciones son separar el filtro o usar la transformada de Fourier.

Filtro separado : la convolución es O (M * N), donde M y N son el número de píxeles en la imagen y el filtro, respectivamente. Como el filtrado promedio con un kernel de 3 por 3 equivale a filtrar primero con un kernel de 3 por 1 y luego un kernel de 1 por 3, puede obtener (3+3)/(3*3) = ~ 30% de mejora de la velocidad por convolución consecutiva con dos núcleos 1-d (esto obviamente mejora a medida que el núcleo se hace más grande). Es posible que aún puedas usar trucos de zancada aquí, por supuesto.

Transformada de Fourier : conv(A,B) es equivalente a ifft(fft(A)*fft(B)) , es decir, una convolución en el espacio directo se convierte en una multiplicación en el espacio de Fourier, donde A es su imagen y B es su filtro. Dado que la multiplicación (en cuanto a los elementos) de las transformadas de Fourier requiere que A y B tengan el mismo tamaño, B es una matriz de size(A) con su kernel en el centro de la imagen y ceros en cualquier otra parte. Para colocar un kernel de 3 por 3 en el centro de una matriz, es posible que tenga que rellenar A a tamaño impar. Dependiendo de su implementación de la transformada de Fourier, esto puede ser mucho más rápido que la convolución (y si aplica el mismo filtro varias veces, puede pre-calcular fft(B) , ahorrando otro 30% del tiempo de cálculo).

Una cosa en la que estoy seguro de que debe arreglarse es su matriz de vistas b .

Tiene algunos elementos de memoria no asignada, por lo que obtendrá lockings.

Dada su nueva descripción de su algoritmo, lo primero que debe solucionarse es el hecho de que está avanzando fuera de la asignación de a :

 bshape = (a.size-filtsize+1, filtsize) bstrides = (a.itemsize, a.itemsize) b = numpy.lib.stride_tricks.as_strided(a, shape=bshape, strides=bstrides) 

Actualizar

Debido a que todavía no entiendo el método y parece que hay formas más simples de resolver el problema, solo voy a poner esto aquí:

 A = numpy.arange(100).reshape((10,10)) shifts = [(-1,-1),(-1,0),(-1,1),(0,-1),(0,1),(1,-1),(1,0),(1,1)] B = A[1:-1, 1:-1].copy() for dx,dy in shifts: xstop = -1+dx or None ystop = -1+dy or None B += A[1+dx:xstop, 1+dy:ystop] B /= 9 

… lo que parece ser el enfoque directo. La única operación extraña es que ha asignado y rellenar B solo una vez. Toda la adición, división e indexación debe hacerse independientemente. Si estás haciendo 16 bandas, solo necesitas asignar B una vez si tu intención es guardar una imagen. Incluso si esto no sirve de ayuda, podría aclarar por qué no entiendo el problema, o al menos servir como punto de referencia para cronometrar las aceleraciones de otros métodos. Esto se ejecuta en 2.6 segundos en mi computadora portátil en una matriz de 5k x 5k de float64, de los cuales 0.5 es la creación de B

Veamos:

No es tan claro su pregunta, pero supongo que le gustaría mejorar significativamente este tipo de promedios.

 import numpy as np from numpy.lib import stride_tricks as st def mf(A, k_shape= (3, 3)): m= A.shape[0]- 2 n= A.shape[1]- 2 strides= A.strides+ A.strides new_shape= (m, n, k_shape[0], k_shape[1]) A= st.as_strided(A, shape= new_shape, strides= strides) return np.sum(np.sum(A, -1), -1)/ np.prod(k_shape) if __name__ == '__main__': A= np.arange(100).reshape((10, 10)) print mf(A) 

Ahora, ¿qué tipo de mejoras de rendimiento realmente esperarías?

Actualizar:
En primer lugar, una advertencia: el código en su estado actual no se adapta correctamente a la forma del ‘kernel’. Sin embargo, esa no es mi principal preocupación en este momento (de todos modos, la idea ya está en cómo adaptarnos adecuadamente).

Acabo de elegir la nueva forma de un 4D A intuitivamente, para mí realmente tiene sentido pensar que el centro del “núcleo” 2D se centre en cada posición de la cuadrícula del 2D original.

Pero esa configuración 4D puede no ser realmente la “mejor”. Creo que el verdadero problema aquí es el rendimiento de la sum. Uno debe poder encontrar el “mejor orden” (de 4D A) para utilizar completamente la architecture de caché de sus máquinas. Sin embargo, ese orden puede no ser el mismo para los arreglos ‘pequeños’ que tipo de ‘coopera’ con el caché de sus máquinas y los más grandes, que no lo hacen (al menos no de manera tan directa).

Actualización 2:
Aquí hay una versión ligeramente modificada de mf . Claramente, es mejor cambiar la forma a una matriz 3D primero y luego, en lugar de sumr solo el producto de puntos (esto tiene la ventaja de que el kernel puede ser arbitrario). Sin embargo, sigue siendo 3 veces más lento (en mi máquina) que la función actualizada de Pauls.

 def mf(A): k_shape= (3, 3) k= np.prod(k_shape) m= A.shape[0]- 2 n= A.shape[1]- 2 strides= A.strides* 2 new_shape= (m, n)+ k_shape A= st.as_strided(A, shape= new_shape, strides= strides) w= np.ones(k)/ k return np.dot(A.reshape((m, n, -1)), w)