hộp giới hạn của mảng numpy

Question

Giả sử bạn có một mảng numpy 2D với một số giá trị ngẫu nhiên và số 0 xung quanh.

Ví dụ "hình chữ nhật nghiêng":

import numpy as np 
from skimage import transform 

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

Bây giờ tôi muốn tìm hình chữ nhật bounding nhỏ nhất cho tất cả các dữ liệu khác không. Ví dụ:

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] 

Điều gì sẽ là cách nhanh nhất để đạt được kết quả này? Tôi chắc chắn có một cách tốt hơn vì hàm np.where mất khá nhiều thời gian nếu tôi hiểu sử dụng bộ dữ liệu 1000x1000.

Edit: cũng nên làm việc trong không gian 3D ...

Answer 1

Bạn gần như có thể giảm một nửa thời gian thực hiện bằng cách sử dụng np.any để giảm các hàng và cột có chứa giá trị khác không để vector 1D, chứ không phải là việc tìm kiếm các chỉ số của tất cả các giá trị khác không sử dụng 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 

Một số tiêu chuẩn:

%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 

Mở rộng phương pháp này đối với trường hợp 3D chỉ liên quan đến việc thực hiện giảm cùng mỗi cặp trục:

Thật dễ dàng để khái quát này để N kích thước bằng cách sử dụng itertools.combinations để lặp qua mỗi sự kết hợp độc đáo của các trục để thực hiện việc giảm hơn:

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) 

---

Nếu bạn biết tọa độ của các góc của hộp giới hạn ban đầu, góc xoay và tâm xoay, bạn có thể nhận được tọa độ của các góc hộp bao phủ được chuyển đổi trực tiếp bằng cách tính affine transformation matrix tương ứng và nằm rải rác nó với đầu vào tọa độ:

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 

này hoạt động ra là rất nhẹ nhanh hơn so với sử dụng np.any ví dụ nhỏ của bạn mảng:

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

Tuy nhiên, do tốc độ của phương thức này độc lập với kích thước của mảng đầu vào, nó có thể nhanh hơn rất nhiều đối với các mảng lớn hơn.

Mở rộng cách tiếp cận chuyển đổi sang 3D phức tạp hơn một chút, trong đó xoay bây giờ có ba thành phần khác nhau (một về trục x, một về trục y và một về trục z), nhưng cơ bản phương pháp là như nhau:

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 

tôi đã về cơ bản chỉ cần sửa đổi các chức năng trên bằng cách sử dụng biểu thức ma trận xoay từ here - tôi đã không có thời gian để viết một bài kiểm tra hợp cụ thể, vì vậy sử dụng một cách thận trọng.

Answer 2

Dưới đây là một thuật toán để tính toán khung giới hạn cho mảng N chiều,

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 

mà có thể được sử dụng với 2D, 3D, vv mảng như sau,

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

Mặc dù thay thế np.diff và np.nonzero cuộc gọi theo số np.where có thể tốt hơn.

Answer 3

Tôi đã có thể ép thêm một chút hiệu suất bằng cách thay thế np.where bằng np.argmax và làm việc trên mặt nạ boolean.

 
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

Điều này nhanh hơn khoảng 10µs so với giải pháp bbox2 trên cùng một điểm chuẩn. Cũng nên có cách để chỉ sử dụng kết quả của argmax để tìm các hàng và cột khác 0, tránh tìm kiếm bổ sung được thực hiện bằng cách sử dụng np.any, nhưng điều này có thể yêu cầu một số chỉ mục phức tạp mà tôi không thể làm việc hiệu quả mã vectơ đơn giản.