Numpy: Get rectangle area just the size of mask

If you prefer pure Numpy, you can achieve this using np.where and np.meshgrid:

i, j = np.where(mask)
indices = np.meshgrid(np.arange(min(i), max(i) + 1),
                      np.arange(min(j), max(j) + 1),
                      indexing='ij')
sub_image = image[indices]

np.where returns a tuple of arrays specifying, pairwise, the indices in each axis for each non-zero element of mask. We then create arrays of all the row and column indices we will want using np.arange, and use np.meshgrid to generate two grid-shaped arrays that index the part of the image we're interested in. Note that we specify matrix-style indexing using index='ij' to avoid having to transpose the result (the default is Cartesian-style indexing).

Essentially, meshgrid constructs indices so that:

image[indices][a, b] == image[indices[0][a, b], indices[1][a, b]]

Example

Start with the following:

>>> image = np.arange(12).reshape((4, 3))
>>> image
array([[ 0,  1,  2],
       [ 3,  4,  5],
       [ 6,  7,  8],
       [ 9, 10, 11]])

Let's say we want to extract the [[3,4],[6,7]] sub-matrix, which is the bounding rectangle for the the following mask:

>>> mask = np.array([[0,0,0],[0,1,0],[1,0,0],[0,0,0]])
>>> mask
array([[0, 0, 0],
       [0, 1, 0],
       [1, 0, 0],
       [0, 0, 0]])

Then, applying the above method:

>>> i, j = np.where(mask)
>>> indices = np.meshgrid(np.arange(min(i), max(i) + 1), np.arange(min(j), max(j) + 1), indexing='ij')
>>> image[indices]
array([[3, 4],
       [6, 7]])

Here, indices[0] is a matrix of row indices, while indices[1] is the corresponding matrix of column indices:

>>> indices[0]
array([[1, 1],
       [2, 2]])
>>> indices[1]
array([[0, 1],
       [0, 1]])

I think using np.amax and np.amin and cropping the image is much faster.

i, j = np.where(mask)
indices = np.meshgrid(np.arange(min(i), max(i) + 1),
              np.arange(min(j), max(j) + 1),
              indexing='ij')
sub_image = image[indices]

Time taken: 50 msec

where = np.array(np.where(mask))

x1, y1 = np.amin(where, axis=1)
x2, y2 = np.amax(where, axis=1)
sub_image = image[x1:x2, y1:y2]

Time taken: 5.6 msec

I don't get Hans's results when running the two methods (using NumPy 1.18.5). In any case, there is a much more efficient method, where you take the arg-max along each dimension

i, j = np.where(mask)
y, x = np.meshgrid(
    np.arange(min(i), max(i) + 1),
    np.arange(min(j), max(j) + 1),
    indexing="ij",
)

Took 38 ms

where = np.array(np.where(mask))
y1, x1 = np.amin(where, axis=1)
y2, x2 = np.amax(where, axis=1) + 1
sub_image = image[y1:y2, x1:x2]

Took 35 ms

maskx = np.any(mask, axis=0)
masky = np.any(mask, axis=1)
x1 = np.argmax(maskx)
y1 = np.argmax(masky)
x2 = len(maskx) - np.argmax(maskx[::-1])
y2 = len(masky) - np.argmax(masky[::-1])
sub_image = image[y1:y2, x1:x2]

Took 2 ms

Timings script