bounding box of numpy array
Asked Answered
C

4

49

Suppose you have a 2D numpy array with some random values and surrounding zeros.

Example "tilted rectangle":

import numpy as np
from skimage import transform

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

Now I want to find the smallest bounding rectangle for all the nonzero data. For example:

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]

What would be the fastest way to achieve this result? I am sure there is a better way since the np.where function takes quite a time if I am e.g. using 1000x1000 data sets.

Edit: Should also work in 3D...

Convalescence answered 14/7, 2015 at 7:49 Comment(1)
lol github copilot just pasted a link to this question and refused to elaborate furtherCocky
A
103

You can roughly halve the execution time by using np.any to reduce the rows and columns that contain non-zero values to 1D vectors, rather than finding the indices of all non-zero values using 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

Some benchmarks:

%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

Extending this approach to the 3D case just involves performing the reduction along each pair of axes:

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

It's easy to generalize this to N dimensions by using itertools.combinations to iterate over each unique combination of axes to perform the reduction over:

import itertools

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

If you know the coordinates of the corners of the original bounding box, the angle of rotation, and the centre of rotation, you could get the coordinates of the transformed bounding box corners directly by computing the corresponding affine transformation matrix and dotting it with the input coordinates:

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

This works out to be very slightly faster than using np.any for your small example array:

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

However, since the speed of this method is independent of the size of the input array, it can be quite a lot faster for larger arrays.

Extending the transformation approach to 3D is slightly more complicated, in that the rotation now has three different components (one about the x-axis, one about the y-axis and one about the z-axis), but the basic method is the same:

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

I've essentially just modified the function above using the rotation matrix expressions from here - I haven't had time to write a test-case yet, so use with caution.

Asserted answered 14/7, 2015 at 9:8 Comment(7)
Nice! How can I extend this to the 3D case? Can I still use np.any somehow?Convalescence
@ali_m: bbox2 is a very good solution, especially if there are large numbers of empty rows/columns, about an order of magnitude faster than: https://mcmap.net/q/356310/-is-there-a-quot-bounding-box-quot-function-slice-with-non-zero-values-for-a-ndarray-in-numpy, but I'm guessing that performance would be similar or worse in the extreme case where there are no non-zero rows/columns.Raymonraymond
@Raymonraymond I'd be surprised if that solution could beat bbox2, even for very large fully dense arrays. In that solution the input and output arrays for np.argwhere increase quadratically with the size of the array, whereas they only increase linearly for np.where in bbox2. One hack that could make it even faster would be to use np.argmax(rows) and rows.size - 1 - np.argmax(rows[::-1]) rather than np.where to get the first and last non-zero values in rows and cols.Asserted
I found a possible bug in this code. xmin, ymin and zmin should be added -1, and xmax, ymax and zmax shoud be added +1.Disenthral
I think the ND solution requires some reversing as the itertools.combinations does yield the reverse of the required order of axes.Screwdriver
@Screwdriver You're right - I didn't notice because the example image is rotationally symmetric about the origin. Generating the combinations of axes in reverse order yields the expected result.Asserted
This is nice solution! Any idea on how to extend it with periodic boundary conditions (or "wrapping around the borders")?Sirenasirenic
S
5

Here is an algorithm to calculate the bounding box for N dimensional arrays,

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

which can be used with 2D, 3D, etc arrays as follows,

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

Although replacing the np.diff and np.nonzero calls by one np.where might be better.

Sophist answered 14/7, 2015 at 9:13 Comment(1)
It's slower than ali_m's approach but very general, I like it!Convalescence
T
5

I was able to squeeze out a little more performance by replacing np.where with np.argmax and working on a boolean mask.

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

This was about 10µs faster for me than the bbox2 solution above on the same benchmark. There should also be a way to just use the result of argmax to find the non-zero rows and columns, avoiding the extra search done by using np.any, but this may require some tricky indexing that I wasn't able to get working efficiently with simple vectorized code.

Trogon answered 19/5, 2016 at 0:41 Comment(1)
Slightly less efficient for me, with many all-zero rows/cols.Raymonraymond
F
2

I know this post is old and has already been answered, but I believe I've identified an optimized approach for large arrays and arrays loaded as np.memmaps.

I was using ali_m's response that was optimized by Allen Zelener for smaller ndarrays, but this approach turns out to be quite slow for np.memmaps.

Below is my implementation that has extremely similar performance speeds to ali_m's approach approach for arrays that fit in the working memory, but that far outperforms when bounding large arrays or 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

My approach is somewhat inefficient in the sense that it just iterates over every point rather than flattening the arrays over specific dimensions. However, I found flattening np.memmaps using np.any() with a dimension argument to be quite slow. I tried using numba to speed up the flattening, but it doesn't support np.any() with arguments. As such, I came to my iterative approach that seems to perform quite well.

On my computer (2019 16" MacBook Pro, 6-core i7, 16 GB 2667 MHz DDR4), I'm able to bound a np.memmap with a shape of (1915, 4948, 3227) in ~33 seconds, as opposed to the ali_m approach that takes around ~250 seconds.

Not sure if anyone will ever see this, but hopefully it helps in the niche cases of needing to bound np.memmaps.

Formula answered 1/1, 2022 at 22:16 Comment(0)

© 2022 - 2024 — McMap. All rights reserved.