Estoy buscando un buen enfoque para la eficiente dividir una imagen en regiones más pequeñas, el procesamiento de cada región por separado y, a continuación, volver a montar los resultados de cada proceso en una sola imagen procesada. Matlab había una herramienta para este llamado blkproc (sustituido por blockproc en las nuevas versiones de Matlab).

En un mundo ideal, la función o clase de apoyo de la superposición entre las divisiones de la matriz de entrada demasiado. En la ayuda de Matlab, blkproc se define como:

B = blkproc(A,[m, n],[mborder nborder],diversión,…)

  • A es su matriz de entrada,
  • [m, n] es el tamaño de bloque
  • [mborder, nborder] es el tamaño de la región de la frontera (opcional)
  • la diversión es una función de la aplicación, para cada bloque de

He kluged juntos un enfoque, pero me parece torpe y apuesto a que hay una manera mucho mejor. A riesgo de mi propia vergüenza, aquí está mi código:


import numpy as np

def segmented_process(M, blk_size=(16,16), overlap=(0,0), fun=None):
    rows = []
    for i in range(0, M.shape[0], blk_size[0]):
        cols = []
        for j in range(0, M.shape[1], blk_size[1]):
            cols.append(fun(M[i:i+blk_size[0], j:j+blk_size[1]]))
        rows.append(np.concatenate(cols, axis=1))
    return np.concatenate(rows, axis=0)

R = np.random.rand(128,128)
passthrough = lambda(x):x
Rprime = segmented_process(R, blk_size=(16,16), 
                           overlap=(0,0), 
                           fun=passthrough)

np.all(R==Rprime)
InformationsquelleAutor Carl F. | 2011-02-22

6 Comentarios

  1. 21

    Aquí están algunos ejemplos de las diferentes (loop gratis) manera de trabajar con los bloques:

    import numpy as np
    from numpy.lib.stride_tricks import as_strided as ast
    
    A= np.arange(36).reshape(6, 6)
    print A
    #[[ 0  1  2  3  4  5]
    # [ 6  7  8  9 10 11]
    # ...
    # [30 31 32 33 34 35]]
    
    # 2x2 block view
    B= ast(A, shape= (3, 3, 2, 2), strides= (48, 8, 24, 4))
    print B[1, 1]
    #[[14 15]
    # [20 21]]
    
    # for preserving original shape
    B[:, :]= np.dot(B[:, :], np.array([[0, 1], [1, 0]]))
    print A
    #[[ 1  0  3  2  5  4]
    # [ 7  6  9  8 11 10]
    # ...
    # [31 30 33 32 35 34]]
    print B[1, 1]
    #[[15 14]
    # [21 20]]
    
    # for reducing shape, processing in 3D is enough
    C= B.reshape(3, 3, -1)
    print C.sum(-1)
    #[[ 14  22  30]
    # [ 62  70  78]
    # [110 118 126]]

    Así que intento copiar simplemente la matlab funcionalidad para numpy no es de todas maneras la mejor manera de proceder. A veces un ‘off the hat» pensar es necesario.

    Advertencia:

    En general, las implementaciones basadas en la zancada trucos puede (pero no necesario necesario) sufren penalizaciones de rendimiento. Así que esté preparado para todas las formas de medir su desempeño. En cualquier caso, es aconsejable comprobar primero si la funcionalidad necesaria (o similares suficiente, con el fin de adaptarse fácilmente para) tiene todo listo sido implementado en numpy o scipy.

    Actualización:

    Por favor, tenga en cuenta que no existe un verdadero magic involucrados en el strides, así que voy a dar una simple función para obtener una block_view de cualquier 2D numpy-matriz. Así que aquí vamos:

    from numpy.lib.stride_tricks import as_strided as ast
    
    def block_view(A, block= (3, 3)):
        """Provide a 2D block view to 2D array. No error checking made.
        Therefore meaningful (as implemented) only for blocks strictly
        compatible with the shape of A."""
        # simple shape and strides computations may seem at first strange
        # unless one is able to recognize the 'tuple additions' involved ;-)
        shape= (A.shape[0]/ block[0], A.shape[1]/ block[1])+ block
        strides= (block[0]* A.strides[0], block[1]* A.strides[1])+ A.strides
        return ast(A, shape= shape, strides= strides)
    
    if __name__ == '__main__':
        from numpy import arange
        A= arange(144).reshape(12, 12)
        print block_view(A)[0, 0]
        #[[ 0  1  2]
        # [12 13 14]
        # [24 25 26]]
        print block_view(A, (2, 6))[0, 0]
        #[[ 0  1  2  3  4  5]
        # [12 13 14 15 16 17]]
        print block_view(A, (3, 12))[0, 0]
        #[[ 0  1  2  3  4  5  6  7  8  9 10 11]
        # [12 13 14 15 16 17 18 19 20 21 22 23]
        # [24 25 26 27 28 29 30 31 32 33 34 35]]
    • Este es un truco y he aprendido algo acerca de la numpy estructura interna, pero se ve muy difícil generalizar arbitraria del tamaño de las matrices.
    • F.: ¿Qué entiende usted por arbitrary sized matrices? Si usted está trabajando con matrices 2D, y la necesidad de dividir a cualquier compatible bloques, el strides y shapes cálculos son triviales. Usted puede pedir más relacionado con este aspecto. De todos modos voy a sugerir que usted dé stride tricks una oportunidad, se puede pagar en el largo plazo. Gracias
    • Me gustaría una función que puede tomar cualquier tamaño de bloque como de entrada y no veo cómo calcular los valores de los avances. Puede usted proporciona un algoritmo para el cálculo de ellos basados en el tamaño de bloque? compute_strides(16,16) frente a compute_strides(2,2)? He leído a través de los documentos, pero encontró un duro para seguir. Para ser honesto, yo estoy más cerca de la física de la ciencia de la computación en estos días.
    • F: he actualizado mi respuesta. ¿Te resulta más útil para sus propósitos, ahora? Gracias
    • Creo que lo voy a comparar las dos respuestas en algún momento y votar por el más rápido. Acabo de recibir distraído por el trabajo de otros.
    • Véase también github.com/scikit-image/scikit-image/blob/master/skimage/util/…. Es exactamente el mismo, salvo un adicional de A = np.ascontiguousarray(A). Es esta falta aquí? ¿Qué sucede si la matriz está en Fortran orden?
    • Podría tal vez nos señalan la dirección correcta en cuanto a cómo unir los bloques individuales en una sola matriz?
    • Pavlovic: por Favor, pon más de lo que quieres decir por «cómo unir los bloques individuales en una sola matriz»! Por favor, tenga en cuenta que los «bloques» son sólo una vista a la matriz original, así que no hay necesidad de unirse a cualquiera de los bloques, la única matriz es simplemente la matriz original 😉
    • Ahora entiendo, gracias!
    • Revisando este ejemplo, observe que los pasos son calculados para ‘int32’ tipos. Para ‘int64’ no va a funcionar

  2. 11

    Proceso por cortes o vistas. La concatenación es muy caro.

    for x in xrange(0, 160, 16):
        for y in xrange(0, 160, 16):
            view = A[x:x+16, y:y+16]
            view[:,:] = fun(view)
    • Esto sin duda ayuda con mis problemas de rendimiento. Voy a tener que adaptar un poco para general el tamaño de las matrices. También, sospecho que podría utilizar un poco de itertools para hacer un solo bucle for. Entonces tal vez habría que apoyar la paralelización.
    • Este enfoque es, sin duda el más fácil de entender y proporciona mejoras sustanciales (2-4x), pero he fumado por la forma y la zancada método. El otro enfoque ha demostrado un orden de magnitud de mejora para las pequeñas regiones de la imagen!
  3. 7

    Me tomó ambas entradas, así como mi enfoque original y la comparación de los resultados. Como @comer señala correctamente, los resultados dependen de la naturaleza de los datos de entrada. Sorprendentemente, concatenar beats ver procesando en unos pocos casos. Cada método tiene un sweet spot. Aquí es mi punto de referencia código:

    import numpy as np
    from itertools import product
    def segment_and_concatenate(M, fun=None, blk_size=(16,16), overlap=(0,0)):
    # truncate M to a multiple of blk_size
    M = M[:M.shape[0]-M.shape[0]%blk_size[0], 
    :M.shape[1]-M.shape[1]%blk_size[1]]
    rows = []
    for i in range(0, M.shape[0], blk_size[0]):
    cols = []
    for j in range(0, M.shape[1], blk_size[1]):
    max_ndx = (min(i+blk_size[0], M.shape[0]),
    min(j+blk_size[1], M.shape[1]))
    cols.append(fun(M[i:max_ndx[0], j:max_ndx[1]]))
    rows.append(np.concatenate(cols, axis=1))
    return np.concatenate(rows, axis=0)
    from numpy.lib.stride_tricks import as_strided
    def block_view(A, block= (3, 3)):
    """Provide a 2D block view to 2D array. No error checking made.
    Therefore meaningful (as implemented) only for blocks strictly
    compatible with the shape of A."""
    # simple shape and strides computations may seem at first strange
    # unless one is able to recognize the 'tuple additions' involved ;-)
    shape= (A.shape[0]/ block[0], A.shape[1]/ block[1])+ block
    strides= (block[0]* A.strides[0], block[1]* A.strides[1])+ A.strides
    return as_strided(A, shape= shape, strides= strides)
    def segmented_stride(M, fun, blk_size=(3,3), overlap=(0,0)):
    # This is some complex function of blk_size and M.shape
    stride = blk_size
    output = np.zeros(M.shape)
    B = block_view(M, block=blk_size)
    O = block_view(output, block=blk_size)
    for b,o in zip(B, O):
    o[:,:] = fun(b);
    return output
    def view_process(M, fun=None, blk_size=(16,16), overlap=None):
    # truncate M to a multiple of blk_size
    from itertools import product
    output = np.zeros(M.shape)
    dz = np.asarray(blk_size)
    shape = M.shape - (np.mod(np.asarray(M.shape), 
    blk_size))
    for indices in product(*[range(0, stop, step) 
    for stop,step in zip(shape, blk_size)]):
    # Don't overrun the end of the array.
    #max_ndx = np.min((np.asarray(indices) + dz, M.shape), axis=0)
    #slices = [slice(s, s + f, None) for s,f in zip(indices, dz)]
    output[indices[0]:indices[0]+dz[0], 
    indices[1]:indices[1]+dz[1]][:,:] = fun(M[indices[0]:indices[0]+dz[0], 
    indices[1]:indices[1]+dz[1]])
    return output
    if __name__ == "__main__":
    R = np.random.rand(128,128)
    squareit = lambda(x):x*2
    from timeit import timeit
    t ={}
    kn = np.array(list(product((8,16,64,128), 
    (128, 512, 2048, 4096))  ) )
    methods = ("segment_and_concatenate", 
    "view_process", 
    "segmented_stride")    
    t = np.zeros((kn.shape[0], len(methods)))
    for i, (k, N) in enumerate(kn):
    for j, method in enumerate(methods):
    t[i,j] = timeit("""Rprime = %s(R, blk_size=(%d,%d), 
    overlap = (0,0), 
    fun = squareit)""" % (method, k, k),
    setup="""
    from segmented_processing import %s
    import numpy as np
    R = np.random.rand(%d,%d)
    squareit = lambda(x):x**2""" % (method, N, N),
    number=5
    )
    print "k =", k, "N =", N #, "time:", t[i]
    print ("    Speed up (view vs. concat, stride vs. concat): %0.4f, %0.4f" % (
    t[i][0]/t[i][1], 
    t[i][0]/t[i][2]))

    Y aquí están los resultados:

    ¿Cómo puedo procesar de manera eficiente una colección de matriz en bloques similar a la de Matlab blkproc (blockproc) función
    Tenga en cuenta que la segmentación de los zancada método gana por 3-4x para los pequeños tamaños de bloque. Sólo en grandes tamaños de bloque (128 x 128) y muy grandes matrices (2048 x 2048 y más grande) se hace la vista del enfoque de procesamiento de ganar, y, a continuación, sólo por un pequeño porcentaje. Basado en el bake-off, parece que @de comer se pone la marca de verificación! Gracias a los dos por buenos ejemplos!

    • Buen punto de referencia! Creo que muestra muy bien algunas de las dificultades que uno se enfrentan al tratar de hacer de propósito general código. En realidad, es altamente no trivial problema para producir código que siempre se realiza de manera óptima. Gracias
  4. 3

    Poco tarde para el juego, pero esto haría bloques superpuestos. No he hecho aquí, pero se puede adaptar fácilmente para tamaños de paso para el desplazamiento de la ventana, creo que:

    from numpy.lib.stride_tricks import as_strided
    def rolling_block(A, block=(3, 3)):
    shape = (A.shape[0] - block[0] + 1, A.shape[1] - block[1] + 1) + block
    strides = (A.strides[0], A.strides[1]) + A.strides
    return as_strided(A, shape=shape, strides=strides)
  5. 3

    Incluso más tarde en el juego.
    Hay un Suizo de procesamiento de Imágenes paquete llamado Bob disponible en:
    https://www.idiap.ch/software/bob/docs/releases/last/sphinx/html/index.html
    Tiene un comando python bob.ip.bloque descrito en:
    https://www.idiap.ch/software/bob/docs/releases/last/sphinx/html/ip/generated/bob.ip.block.html#bob.ip.block
    que parece hacer todo lo que el comando de Matlab ‘blockproc no.
    Yo no lo he probado.

    También hay una interesante comando bob.ip.DCTFeatures que incorpora el bloque de la capacidad de extraer o modificar los coeficientes de la DCT de la imagen.

  6. 0

    He encontrado este tutorial – El último código fuente proporciona exactamente la funcionalidad deseada!
    Aún debe trabajar para cualquier dimensionalidad (yo no probarlo)
    http://www.johnvinyard.com/blog/?p=268

    A pesar de la «aplanar» la opción al final del código fuente parece ser un poco buggy. Sin embargo, muy bonita pieza de software!

Dejar respuesta

Please enter your comment!
Please enter your name here