Tengo las siguientes celdas:
cells = np.array([[1, 1, 1], [1, 1, 0], [1, 0, 0], [1, 0, 1], [1, 0, 0], [1, 1, 1]])y quiero calcular adyacencias horizontales y verticales para llegar a este resultado:
# horizontal adjacency array([[3, 2, 1], [2, 1, 0], [1, 0, 0], [1, 0, 1], [1, 0, 0], [3, 2, 1]]) # vertical adjacency array([[6, 2, 1], [5, 1, 0], [4, 0, 0], [3, 0, 1], [2, 0, 0], [1, 1, 1]])La solución real se ve así:
def get_horizontal_adjacency(cells): adjacency_horizontal = np.zeros(cells.shape, dtype=int) for y in range(cells.shape[0]): span = 0 for x in reversed(range(cells.shape[1])): if cells[y, x] > 0: span += 1 else: span = 0 adjacency_horizontal[y, x] = span return adjacency_horizontal def get_vertical_adjacency(cells): adjacency_vertical = np.zeros(cells.shape, dtype=int) for x in range(cells.shape[1]): span = 0 for y in reversed(range(cells.shape[0])): if cells[y, x] > 0: span += 1 else: span = 0 adjacency_vertical[y, x] = span return adjacency_verticalEl algoritmo es básicamente (para la adyacencia horizontal):
Dado que necesito recorrer dos veces todos los elementos de la matriz, esto es lento para matrices más grandes (por ejemplo, imágenes).
¿Hay alguna manera de mejorar el algoritmo usando vectorización o alguna otra magia numpy?
Resumen :
Joni y Mark Setchell han hecho excelentes sugerencias.
Creé un pequeño Repo con una imagen de muestra y un archivo de python con las comparaciones. Los resultados son asombrosos:
Hice un intento muy rápido con Numba, pero no lo he comprobado demasiado, aunque los resultados parecen correctos:
#!/usr/bin/env python3 # https://stackoverflow.com/q/69854335/2836621 # magick -size 1920x1080 xc:black -fill white -draw "circle 960,540 960,1040" -fill black -draw "circle 960,540 960,800" a.png import cv2 import numpy as np import numba as nb def get_horizontal_adjacency(cells): adjacency_horizontal = np.zeros(cells.shape, dtype=int) for y in range(cells.shape[0]): span = 0 for x in reversed(range(cells.shape[1])): if cells[y, x] > 0: span += 1 else: span = 0 adjacency_horizontal[y, x] = span return adjacency_horizontal @nb.jit('void(uint8[:,::1], int32[:,::1])',parallel=True) def nb_get_horizontal_adjacency(cells, result): for y in nb.prange(cells.shape[0]): span = 0 for x in range(cells.shape[1]-1,-1,-1): if cells[y, x] > 0: span += 1 else: span = 0 result[y, x] = span return # Load image im = cv2.imread('a.png', cv2.IMREAD_GRAYSCALE) %timeit get_horizontal_adjacency(im) result = np.zeros((im.shape[0],im.shape[1]),dtype=np.int32) %timeit nb_get_horizontal_adjacency(im, result)Los tiempos son buenos, mostrando una aceleración de 4000x, si funciona correctamente:
In [15]: %timeit nb_get_horizontal_adjacency(im, result) 695 µs ± 9.12 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each) In [17]: %timeit get_horizontal_adjacency(im) 2.78 s ± 44.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)Aporte
La imagen de entrada se creó en dimensiones de 1080p, es decir, 1920x1080, con ImageMagick usando:
magick -size 1920x1080 xc:black -fill white -draw "circle 960,540 960,1040" -fill black -draw "circle 960,540 960,800" a.pngSalida (contraste ajustado)
Como ya se indicó en los comentarios, este es un ejemplo perfecto en el que es más fácil simplemente reescribir la función por medio de Cython o Numba. Dado que Mark ya proporcionó una solución Numba, permítanme proporcionar una solución Cython. Primero, programemos su solución en mi máquina para una comparación justa:
In [5]: %timeit nb_get_horizontal_adjacency(im, result) 836 µs ± 36 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each) Suponiendo que la imagen im es un np.ndarray con dtype=np.uint8 , una solución Cython paralelizada se ve así:
In [6]: %%cython -f -a -c=-O3 -c=-march=native -c=-fopenmp --link-args=-fopenmp from cython import boundscheck, wraparound, initializedcheck from libc.stdint cimport uint8_t, uint32_t from cython.parallel cimport prange import numpy as np @boundscheck(False) @wraparound(False) @initializedcheck(False) def cy_get_horizontal_adjacency(uint8_t[:, ::1] cells): cdef int nrows = cells.shape[0] cdef int ncols = cells.shape[1] cdef uint32_t[:, ::1] adjacency_horizontal = np.zeros((nrows, ncols), dtype=np.uint32) cdef int x, y, span for y in prange(nrows, nogil=True, schedule="static"): span = 0 for x in reversed(range(ncols)): if cells[y, x] > 0: span += 1 else: span = 0 adjacency_horizontal[y, x] = span return np.array(adjacency_horizontal, copy=False)En mi máquina, esto es casi dos veces más rápido:
In [7]: %timeit cy_get_horizontal_adjacency(im) 431 µs ± 4.38 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)