a.smiet - 1 year ago 233
Python Question

# bounding box of numpy array

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

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

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

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

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

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

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.

Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download