not - open image python numpy

Usando zancadas para un eficiente filtro de promedio móvil

Question

Hace poco aprendí sobre los strides 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 (usando filtros de convolución).

Esto es lo que tengo hasta ahora. Toma una vista de la matriz original, luego la rueda por la cantidad necesaria y suma los valores del núcleo para calcular el promedio. Soy consciente de que los bordes no se manejan correctamente, pero puedo ocuparme de eso luego ... ¿Hay una manera mejor y más rápida? El objetivo es filtrar grandes matrices de coma flotante de hasta 5000x5000 x 16 capas de tamaño, una tarea que scipy.ndimage.filters.convolve es bastante lenta en.

Tenga en cuenta que estoy buscando conectividad de 8 vecinos, es decir, un filtro de 3x3 toma un 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 esto funcionando:

Código actual:

1. utilice 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 del filtro.
2. Avance por el eje vertical para obtener la fila central del núcleo [[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 suma de cada fila y la divido por la cantidad de elementos en el filtro, obteniendo el promedio de cada píxel, (desplazado por 1 fila y 1 columna, y con algunas rarezas alrededor de los bordes, pero puedo ocúpate de eso más tarde).

Lo que esperaba era un mejor uso de stride_tricks para obtener los 9 valores o la suma de los elementos del kernel directamente, para toda la matriz, o que alguien me puede convencer de otro método más eficiente ...

Answer 1

No estoy lo suficientemente familiarizado con Python para escribir el código para eso, pero las dos mejores formas de acelerar las convoluciones es 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 es equivalente a filtrar primero con un kernel de 3 por 1 y luego un kernel de 1 por 3, puedes obtener (3+3)/(3*3) = ~ Mejora del 30% de velocidad por convolución consecutiva con dos núcleos de 1-d (esto obviamente mejora a medida que el núcleo se hace más grande). Puede que aún pueda usar trucos de zancadas 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 tu imagen y B es tu filtro. Como la multiplicación (a nivel de elemento) de las transformadas de Fourier requiere que A y B tengan el mismo tamaño, B es una matriz de size(A) con su núcleo en el centro de la imagen y ceros en cualquier otro lugar. Para colocar un núcleo de 3 por 3 en el centro de una matriz, es posible que tenga que aplicar A al 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 calcular previamente fft(B) , ahorrando otro 30% del tiempo de cálculo).

Answer 2

Por lo que vale, así es como lo harías usando trucos de "fantasía". ¡Iba a publicar esto ayer, pero me distrajo con el trabajo real! :)

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

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

Además, si intenta obtener una ventana móvil multidimensional, corre el riesgo de explotar el uso de la memoria cada vez que haga una copia inadvertidamente de su matriz. Mientras que la matriz "rodada" inicial es solo una vista en la memoria de la matriz original, cualquier paso intermedio que copie la matriz hará una copia que sea de órdenes de magnitud más grandes que la matriz original (es decir, digamos que está trabajando con una matriz original de 100x100 ... La vista dentro de ella (para un tamaño de filtro de (3,3)) será de 98x98x3x3 pero usará la misma memoria que el original. Sin embargo, cualquier copia usará la cantidad de memoria que una matriz completa de 98x98x3x3 ¡¡haría!!)

Básicamente, usar trucos locos es genial para cuando quieres vectorizar las operaciones de la ventana móvil en un solo eje de un ndarray. Hace que sea muy fácil calcular cosas como una desviación estándar móvil, etc. con muy poca sobrecarga. Cuando desee comenzar a hacer esto a lo largo de múltiples ejes, es posible, pero generalmente estará 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.
    <http://www.mail-archive.com/[email protected]/msg29450.html>"""
    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 lo 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 10x10. Podríamos haber utilizado fácilmente diferentes tamaños de filtro a lo largo de diferentes ejes u operar solo a lo largo de los ejes seleccionados de una matriz N-dimensional (es decir, filtsize = (0,3,0,3) en una matriz de 4 dimensiones nos daría una dimensión 6 ver).

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

Sin embargo, debido a que estamos almacenando matrices temporales que son mucho más grandes que nuestra matriz original en cada paso de la mean (o std o lo que sea), ¡esto no es para nada eficiente desde el punto de vista de la 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, haga el "desenfoque" en el lugar sin requerir una copia temporal de la matriz y sea muy rápido. Los trucos en zancada son una buena forma de aplicar una función a una ventana en movimiento a lo largo de un eje, pero no son una buena forma de hacerlo a lo largo de múltiples ejes, por lo general ....

Solo mi $ 0.02, en cualquier caso ...

Answer 3

Una cosa de la que estoy seguro necesita ser corregida es tu vista array b .

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

Dada su nueva descripción de su algoritmo, lo primero que debe corregir es el hecho de que está caminando 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 aún no entiendo el método y parece que hay formas más simples de resolver el problema, 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

... que parece el enfoque directo. La única operación extraña es que ha asignado y poblado B solo una vez. Toda la adición, división e indexación tiene que hacerse independientemente. Si está haciendo 16 bandas, solo necesita asignar B una vez si su intención es guardar una imagen. Incluso si esto no es de ayuda, podría aclarar por qué no entiendo el problema, o al menos servir como un 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

Answer 4

Veamos:

No está tan claro de su pregunta, pero estoy asumiendo ahora que le gustaría mejorar significativamente este tipo de promediado.

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 ''kernel''. Sin embargo, esa no es mi principal preocupación en este momento (de todos modos, la idea ya está lista para adaptarse adecuadamente).

Acabo de elegir la nueva forma de un 4D A de manera intuitiva, para mí realmente tiene sentido pensar en un centro ''kernel'' 2D que se centre en cada posición de cuadrícula del A 2D original.

Pero esa forma 4D puede no ser realmente la "mejor". Creo que el verdadero problema aquí es el rendimiento de la suma. Uno debe ser capaz de encontrar el ''mejor orden'' (del 4D A) para utilizar completamente la arquitectura de caché de su máquina. Sin embargo, ese orden puede no ser el mismo para arreglos ''pequeños'' que tipo de ''coopera'' con la memoria caché de la máquina y aquellos más grandes, que no (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 de una matriz en 3D primero y luego, en lugar de sumar, solo hacer un producto por puntos (esto tiene la ventaja de que todo eso, ese núcleo puede ser arbitrario). Sin embargo, sigue siendo un 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)