ograniczające pole numpy array

Question

Załóżmy, że masz tablicę 2D numpy z pewnymi losowymi wartościami i otaczającymi zerami.

Przykład „przechylił prostokąt”:

import numpy as np 
from skimage import transform 

img1 = np.zeros((100,100)) 
img1[25:75,25:75] = 1. 
img2 = transform.rotate(img1, 45) 

Teraz chcę znaleźć najmniejszy prostokąt, dla wszystkich danych niezerowych. Na przykład:

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] 

Jaki byłby najszybszym sposobem do osiągnięcia tego celu? Jestem pewien, że jest lepszy sposób, ponieważ funkcja np.where zajmuje trochę czasu, jeśli jestem np. przy użyciu zestawów danych 1000x1000.

Edycja: Należy również pracować w trybie 3D ...

Answer 1

Można z grubsza o połowę czas wykonania za pomocą np.any zmniejszyć wiersze i kolumny, które zawierają wartości niezerowe wektory 1D, niż znalezienie wskaźników Wszystkie wartości niezerowe wykorzystujące 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 

Niektóre Benchmarki:

%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 

rozszerzenie tego podejścia do sprawy 3D po prostu polega na przeprowadzeniu redukcji wzdłuż każdej pary osiach:

Łatwo uogólnić to N wymiarów za pomocą itertools.combinations iteracyjne nad każdą unikalną kombinację osi przeprowadzić redukcję nad:

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) 

---

Jeśli znasz współrzędne narożników oryginalnej ramki ograniczającej, kąt obrotu i środek obrotu, możesz uzyskać współrzędne przekształconych narożników obwiedni bezpośrednio poprzez obliczenie odpowiedniej affine transformation matrix i rozsianych go z wejściem współrzędne:

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 

Działa to się bardzo nieznacznie szybciej niż przy użyciu np.any za mały przykład tablicy:

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

Jednak od prędkości tej metody jest niezależna od rozmiaru tablicy wejściowej, może być znacznie szybsza w przypadku większych tablic.

Rozszerzenie podejścia do transformacji na 3D jest nieco bardziej skomplikowane, ponieważ obrót ma teraz trzy różne komponenty (jeden o osi X, jeden o osi Y i jeden o osi Z), ale podstawowy metoda jest taka sama:

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 

ja w zasadzie tylko zmodyfikowane funkcję powyżej używając wyrażeń macierz obrotu z here - nie miałem czasu, aby napisać test-sprawę jeszcze, więc należy używać z rozwagą.

Answer 2

Tutaj jest algorytm do obliczania obwiedni dla N wymiarów macierzy,

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 

, które mogą być stosowane z 2d, 3d itp tablice w następujący

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]] 

Choć zastąpienie np.diff i np.nonzero połączenia przez jeden np.where może być lepiej.

Answer 3

Udało mi się wycisnąć nieco więcej wydajności, zastępując np.where z np.argmax i pracując na masce boolowskiej.

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

To było około 10μs szybsze dla mnie niż powyższe rozwiązanie bbox2 na tym samym benchmarku. Powinien również istnieć sposób użycia wyniku argmax do znalezienia niezerowych wierszy i kolumn, unikając dodatkowego wyszukiwania wykonanego przy użyciu np.any, ale może to wymagać trochę trudnego indeksowania, którego nie mogłem efektywnie pracować z prosty wektorowy kod.