ограничивающая рамка из массива numpy

Question

Предположим, что у вас есть 2D-массив numpy с некоторыми случайными значениями и окружающими нулями.

Пример «наклонена прямоугольник»:

import numpy as np 
from skimage import transform 

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

Теперь я хочу, чтобы найти наименьший прямоугольник для всех ненулевых данных. Например:

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] 

Что бы быстрый способ достичь этого результата? Я уверен, что есть лучший способ, поскольку функция np.where занимает довольно много времени, если я, например. с использованием наборов данных 1000x1000.

Edit: Если работать в 3D ...

Answer 1

Вы можете примерно вдвое сократить время выполнения с помощью np.any, чтобы ограничить количество строк и столбцов, содержащих ненулевые значения векторов 1D, а не поиска индексов все ненулевые значения, используя 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 

Некоторые ориентиры:

%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 

Продление такой подход к 3D случае просто включает в себя проведение сокращения вдоль каждой пару осей:

Легко обобщить это N размеров с помощью itertools.combinations для перебора каждой уникальной комбинации осей, чтобы выполнить редукцию над:

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) 

---

Если вы знаете координаты углов исходного ограничивающего прямоугольника, угла поворота и центра вращения, вы можете получить координаты преобразованных углов прямоугольной рамки непосредственно путем вычисления соответствующего affine transformation matrix и усеивание его с входом координаты:

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 

Это работает очень немного быстрее, чем при использовании np.any для малого например массива:

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

Однако, поскольку скорость этого метода не зависит от размера входного массива, он может быть намного быстрее для больших массивов.

Расширение подхода к трансформации в 3D немного сложнее, поскольку теперь вращение состоит из трех различных компонентов (один относительно оси x, один относительно оси y и один вокруг оси z), но основной метод тот же:

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 

Я по существу только изменил функцию выше, используя матричные выражения вращения от here - у меня не было времени, чтобы написать тест-случай еще, так что используйте с осторожностью.

Answer 2

Ниже приведен алгоритм для вычисления ограничивающего прямоугольника для N одномерных массивов,

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 

, которые могут быть использованы с 2D, 3D и т.д. массивы следующим образом,

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

Несмотря на то, заменив np.diff и np.nonzero звонки одним np.where может быть лучше.

Answer 3

Я смог выжать немного больше производительности, заменив np.where на np.argmax и работая над булевой маской.

 
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

Это было примерно на 10 мкс быстрее для меня, чем решение bbox2 выше, на том же контрольном. Также должен быть способ использовать результат argmax для поиска ненулевых строк и столбцов, избегая дополнительного поиска, сделанного с помощью np.any, но для этого может потребоваться некоторая сложная индексация, с которой я не мог эффективно работать с простой векторный код.