caja delimitadora de matriz numpy

Supongamos que tiene una matriz numpy 2D con algunos valores aleatorios y ceros circundantes.

Ejemplo “rectángulo inclinado”:

import numpy as np from skimage import transform img1 = np.zeros((100,100)) img1[25:75,25:75] = 1. img2 = transform.rotate(img1, 45) 

Ahora quiero encontrar el rectángulo delimitador más pequeño para todos los datos distintos de cero. Por ejemplo:

 a = np.where(img2 != 0) bbox = img2[np.min(a[0]):np.max(a[0])+1, np.min(a[1]):np.max(a[1])+1] 

¿Cuál sería la forma más rápida de lograr este resultado? Estoy seguro de que hay una mejor manera ya que la función np.where lleva bastante tiempo si, por ejemplo, estoy usando conjuntos de datos de 1000×1000.

Edición: También debería funcionar en 3D …

Puede reducir aproximadamente a la mitad el tiempo de ejecución utilizando np.any para reducir las filas y columnas que contienen valores distintos de cero a vectores 1D, en lugar de encontrar los índices de todos los valores distintos de cero utilizando np.where :

 def bbox1(img): a = np.where(img != 0) bbox = np.min(a[0]), np.max(a[0]), np.min(a[1]), np.max(a[1]) return bbox def bbox2(img): rows = np.any(img, axis=1) cols = np.any(img, axis=0) rmin, rmax = np.where(rows)[0][[0, -1]] cmin, cmax = np.where(cols)[0][[0, -1]] return rmin, rmax, cmin, cmax 

Algunos puntos de referencia:

 %timeit bbox1(img2) 10000 loops, best of 3: 63.5 µs per loop %timeit bbox2(img2) 10000 loops, best of 3: 37.1 µs per loop 

Extender este enfoque al caso 3D solo implica realizar la reducción a lo largo de cada par de ejes:

 def bbox2_3D(img): r = np.any(img, axis=(1, 2)) c = np.any(img, axis=(0, 2)) z = np.any(img, axis=(0, 1)) rmin, rmax = np.where(r)[0][[0, -1]] cmin, cmax = np.where(c)[0][[0, -1]] zmin, zmax = np.where(z)[0][[0, -1]] return rmin, rmax, cmin, cmax, zmin, zmax 

Es fácil generalizar esto a N dimensiones usando itertools.combinations para iterar sobre cada combinación única de ejes para realizar la reducción sobre:

 import itertools def bbox2_ND(img): N = img.ndim out = [] for ax in itertools.combinations(range(N), N - 1): nonzero = np.any(img, axis=ax) out.extend(np.where(nonzero)[0][[0, -1]]) return tuple(out) 

Si conoce las coordenadas de las esquinas del cuadro delimitador original, el ángulo de rotación y el centro de rotación, puede obtener las coordenadas de las esquinas del cuadro delimitador transformadas directamente calculando la matriz de transformación afín correspondiente y punteando con la entrada coordenadas:

 def bbox_rotate(bbox_in, angle, centre): rmin, rmax, cmin, cmax = bbox_in # bounding box corners in homogeneous coordinates xyz_in = np.array(([[cmin, cmin, cmax, cmax], [rmin, rmax, rmin, rmax], [ 1, 1, 1, 1]])) # translate centre to origin cr, cc = centre cent2ori = np.eye(3) cent2ori[:2, 2] = -cr, -cc # rotate about the origin theta = np.deg2rad(angle) rmat = np.eye(3) rmat[:2, :2] = np.array([[ np.cos(theta),-np.sin(theta)], [ np.sin(theta), np.cos(theta)]]) # translate from origin back to centre ori2cent = np.eye(3) ori2cent[:2, 2] = cr, cc # combine transformations (rightmost matrix is applied first) xyz_out = ori2cent.dot(rmat).dot(cent2ori).dot(xyz_in) r, c = xyz_out[:2] rmin = int(r.min()) rmax = int(r.max()) cmin = int(c.min()) cmax = int(c.max()) return rmin, rmax, cmin, cmax 

Esto resulta ser un poco más rápido que usar np.any para su pequeña matriz de ejemplo:

 %timeit bbox_rotate([25, 75, 25, 75], 45, (50, 50)) 10000 loops, best of 3: 33 µs per loop 

Sin embargo, dado que la velocidad de este método es independiente del tamaño de la matriz de entrada, puede ser bastante más rápida para matrices más grandes.

Extender el enfoque de transformación a 3D es un poco más complicado, ya que la rotación ahora tiene tres componentes diferentes (uno sobre el eje x, uno sobre el eje y y otro sobre el eje z), pero el método básico es el mismo :

 def bbox_rotate_3d(bbox_in, angle_x, angle_y, angle_z, centre): rmin, rmax, cmin, cmax, zmin, zmax = bbox_in # bounding box corners in homogeneous coordinates xyzu_in = np.array(([[cmin, cmin, cmin, cmin, cmax, cmax, cmax, cmax], [rmin, rmin, rmax, rmax, rmin, rmin, rmax, rmax], [zmin, zmax, zmin, zmax, zmin, zmax, zmin, zmax], [ 1, 1, 1, 1, 1, 1, 1, 1]])) # translate centre to origin cr, cc, cz = centre cent2ori = np.eye(4) cent2ori[:3, 3] = -cr, -cc -cz # rotation about the x-axis theta = np.deg2rad(angle_x) rmat_x = np.eye(4) rmat_x[1:3, 1:3] = np.array([[ np.cos(theta),-np.sin(theta)], [ np.sin(theta), np.cos(theta)]]) # rotation about the y-axis theta = np.deg2rad(angle_y) rmat_y = np.eye(4) rmat_y[[0, 0, 2, 2], [0, 2, 0, 2]] = ( np.cos(theta), np.sin(theta), -np.sin(theta), np.cos(theta)) # rotation about the z-axis theta = np.deg2rad(angle_z) rmat_z = np.eye(4) rmat_z[:2, :2] = np.array([[ np.cos(theta),-np.sin(theta)], [ np.sin(theta), np.cos(theta)]]) # translate from origin back to centre ori2cent = np.eye(4) ori2cent[:3, 3] = cr, cc, cz # combine transformations (rightmost matrix is applied first) tform = ori2cent.dot(rmat_z).dot(rmat_y).dot(rmat_x).dot(cent2ori) xyzu_out = tform.dot(xyzu_in) r, c, z = xyzu_out[:3] rmin = int(r.min()) rmax = int(r.max()) cmin = int(c.min()) cmax = int(c.max()) zmin = int(z.min()) zmax = int(z.max()) return rmin, rmax, cmin, cmax, zmin, zmax 

Básicamente, acabo de modificar la función anterior usando las expresiones de matriz de rotación de aquí . No he tenido tiempo de escribir un caso de prueba todavía, así que use con precaución.

Aquí hay un algoritmo para calcular el cuadro delimitador para matrices tridimensionales,

 def get_bounding_box(x): """ Calculates the bounding box of a ndarray""" mask = x == 0 bbox = [] all_axis = np.arange(x.ndim) for kdim in all_axis: nk_dim = np.delete(all_axis, kdim) mask_i = mask.all(axis=tuple(nk_dim)) dmask_i = np.diff(mask_i) idx_i = np.nonzero(dmask_i)[0] if len(idx_i) != 2: raise ValueError('Algorithm failed, {} does not have 2 elements!'.format(idx_i)) bbox.append(slice(idx_i[0]+1, idx_i[1]+1)) return bbox 

que se puede utilizar con matrices 2D, 3D, etc. de la siguiente manera,

 In [1]: print((img2!=0).astype(int)) ...: bbox = get_bounding_box(img2) ...: print((img2[bbox]!=0).astype(int)) ...: [[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0] [0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0] [0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0] [0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0] [0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0] [0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0] [0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0] [0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0] [0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0] [0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]] [[0 0 0 0 0 0 1 1 0 0 0 0 0 0] [0 0 0 0 0 1 1 1 1 0 0 0 0 0] [0 0 0 0 1 1 1 1 1 1 0 0 0 0] [0 0 0 1 1 1 1 1 1 1 1 0 0 0] [0 0 1 1 1 1 1 1 1 1 1 1 0 0] [0 1 1 1 1 1 1 1 1 1 1 1 1 0] [1 1 1 1 1 1 1 1 1 1 1 1 1 1] [1 1 1 1 1 1 1 1 1 1 1 1 1 1] [0 1 1 1 1 1 1 1 1 1 1 1 1 0] [0 0 1 1 1 1 1 1 1 1 1 1 0 0] [0 0 0 1 1 1 1 1 1 1 1 0 0 0] [0 0 0 0 1 1 1 1 1 1 0 0 0 0] [0 0 0 0 0 1 1 1 1 0 0 0 0 0] [0 0 0 0 0 0 1 1 0 0 0 0 0 0]] 

Aunque reemplazar las llamadas np.diff y np.nonzero por un np.where podría ser mejor.

Pude exprimir un poco más el rendimiento al reemplazar np.where con np.argmax y trabajar en una máscara booleana.

 def bbox (img):
     img = (img> 0)
     filas = np.any (img, axis = 1)
     cols = np.any (img, axis = 0)
     rmin, rmax = np.argmax (filas), img.shape [0] - 1 - np.argmax (np.flipud (filas))
     cmin, cmax = np.argmax (cols), img.shape [1] - 1 - np.argmax (np.flipud (cols))
     devuelve rmin, rmax, cmin, cmax 

Esto fue unos 10 µs más rápido para mí que la solución bbox2 anterior en el mismo punto de referencia. También debería haber una manera de usar el resultado de argmax para encontrar las filas y columnas que no sean cero, evitando la búsqueda adicional que se realiza mediante el uso de np.any , pero esto puede requerir un índice complicado que no pude hacer funcionar. Eficientemente con código simple vectorizado.