Ограничивающий прямоугольник массива
Предположим, у вас есть двумерный массив с некоторыми случайными значениями и окружающими нулями.
Пример "наклоненный прямоугольник":
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.
Редактировать: должен также работать в 3D...
4 ответа
Вы можете примерно вдвое сократить время выполнения, используя np.any уменьшить строки и столбцы, которые содержат ненулевые значения, до одномерных векторов, вместо того, чтобы находить индексы всех ненулевых значений, используя 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
Расширение этого подхода к трехмерному случаю просто включает выполнение сокращения по каждой паре осей:
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
Легко обобщить это на 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)
Если вы знаете координаты углов исходного ограничивающего прямоугольника, угол поворота и центр вращения, вы можете получить координаты преобразованных углов ограничивающего прямоугольника непосредственно, вычислив соответствующую матрицу аффинного преобразования и поставив ее с помощью ввода координаты:
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
По сути, я только что изменил приведенную выше функцию с помощью выражений матрицы вращения отсюда - у меня еще не было времени написать тестовый пример, поэтому используйте его с осторожностью.
Я смог выжать немного больше производительности, заменив np.where с np.argmax и работает над логической маской.
def bbox (img):
img = (img> 0)
row = np.any(img, axis=1)
cols = np.any(img, axis=0)
rmin, rmax = np.argmax(строки), img.shape[0] - 1 - np.argmax(np.flipud(строки))
cmin, cmax = np.argmax (столбцы), img.shape[1] - 1 - np.argmax(np.flipud(столбцы))
возврат rmin, rmax, cmin, cmaxДля меня это было примерно на 10 мкс быстрее, чем решение bbox2, описанное выше, в том же тесте. Также должен быть способ просто использовать результат argmax, чтобы найти ненулевые строки и столбцы, избегая дополнительного поиска, выполняемого с помощью np.any, но это может потребовать некоторой сложной индексации, которую я не смог эффективно использовать с простым векторизованным кодом.
Вот алгоритм для вычисления ограничивающей рамки для 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 может быть лучше.
Я знаю, что этот пост устарел и на него уже ответили, но я считаю, что нашел оптимизированный подход для больших массивов и массивов, загружаемых как np.memmaps.
Я использовал ответ ali_m, который был оптимизирован Алленом Зеленером для небольших ndarray, но этот подход оказался довольно медленным для np.memmaps.
Ниже приведена моя реализация, которая имеет очень схожие скорости производительности с подходом ali_m для массивов, которые помещаются в рабочую память, но намного превосходит их при ограничении больших массивов или np.memmaps.
import numpy as np
from numba import njit, prange
@njit(parallel=True, nogil=True, cache=True)
def bound(volume):
"""
Bounding function to bound large arrays and np.memmaps
volume: A 3D np.array or np.memmap
"""
mins = np.array(volume.shape)
maxes = np.zeros(3)
for z in prange(volume.shape[0]):
for y in range(volume.shape[1]):
for x in range(volume.shape[2]):
if volume[z,y,x]:
if z < mins[0]:
mins[0] = z
elif z > maxes[0]:
maxes[0] = z
if y < mins[1]:
mins[1] = y
elif y > maxes[1]:
maxes[1] = y
if x < mins[2]:
mins[2] = x
elif x > maxes[2]:
maxes[2] = x
return mins, maxes
Мой подход несколько неэффективен в том смысле, что он просто перебирает каждую точку, а не сглаживает массивы по определенным измерениям. Однако я обнаружил, что выравнивание np.memmaps с использованием np.any() с аргументом измерения довольно медленное. Я попытался использовать numba для ускорения выравнивания, но он не поддерживает np.any() с аргументами. Таким образом, я пришел к своему итеративному подходу, который, кажется, работает достаточно хорошо.
На моем компьютере (16-дюймовый MacBook Pro 2019 года, 6-ядерный i7, 16 ГБ, 2667 МГц DDR4) я могу связать np.memmap с формой (1915, 4948, 3227) примерно за 33 секунды , а не к подходу ali_m, который занимает около 250 секунд.
Не уверен, что кто-нибудь когда-нибудь это увидит, но, надеюсь, это поможет в нишевых случаях, когда необходимо связать np.memmaps.