Find the largest rectangular bounding box containing only ones in a 2d mask NumPy array

I have a 2d mask where I'd like to find the largest rectangular bounding box containing only ones.

Is there a more efficient way than determining all possibilities for bounding boxes and then comparing their size?

Example code:

import matplotlib.patches as patches
import matplotlib.pyplot as plt
import numpy as np

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

def find_largest_bbox(mask):
    # TODO
    return 3, 1, 4, 4

bbox = find_largest_bbox(mask)

plt.close()
plt.imshow(mask)
x, y = bbox[0] - 0.5, bbox[1] - 0.5
w, h = bbox[2] - bbox[0] + 1, bbox[3] - bbox[1] + 1
rect = patches.Rectangle((x, y), w, h, linewidth=2, edgecolor="red", facecolor="none")
plt.gca().add_patch(rect)
plt.show()

One option would be to use lir (that uses Numba) :

# takes some time for compilation
# required only after installation
import largestinteriorrectangle as lir

def find_largest_bbox(mask):
    return lir.lir(mask.astype("bool"))

x, y, w, h = find_largest_bbox(mask)

im = plt.imshow(mask)

xp, *_, yp = im.get_extent()

plt.gca().add_patch(
    patches.Rectangle(
        (x + xp, y + yp), w, h,
        lw=2, fc="none", ec="r",
    )
)

I think you can do it by dynamic programming (ends up as O(N^3) time complexity).

Let B[i,j,w] be the biggest area that ends (bottom-right corner) at [i,j] and has width w. Then, pseudocode:

if mask[i,j]>0 and ones>=w then
    B[i,j,w] = B[i-1,j,w] + w
else
    B[i,j,w] = 0

where ones is the current tally of consecutive ones for that row.

If you are careful then you only need store B[j,w] and update it at successive rows.

Note that the routine returns -1,-1,-1,-1 for a mask of all zeros, so you'd need to deal outside with that edge case.

import matplotlib.patches as patches
import matplotlib.pyplot as plt
import numpy as np

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

def find_largest_bbox( mask ):
    ni, nj = mask.shape
    bestArea = 0
    x1, y1, x2, y2 = -1, -1, -1, -1
    B = np.zeros( ( nj, nj + 1 ) )
    for i in range( ni ):
        ones = 0
        for j in range( nj ):
            if mask[i,j] > 0: 
                ones += 1
                for w in range( 1, j + 2 ):
                    if ones >= w: 
                        B[j,w] += w
                        if B[j,w] > bestArea:
                            bestArea = B[j,w]
                            x1, y1, x2, y2 = j - w + 1, i - bestArea // w + 1, j, i
                    else:
                        B[j,w] = 0
            else:
                ones = 0
                B[j,:] = 0
    return [ x1, y1, x2, y2 ]


bbox = find_largest_bbox(mask)

plt.close()
plt.imshow(mask)
x, y = bbox[0] - 0.5, bbox[1] - 0.5
w, h = bbox[2] - bbox[0] + 1, bbox[3] - bbox[1] + 1
rect = patches.Rectangle((x, y), w, h, linewidth=2, edgecolor="red", facecolor="none")
plt.gca().add_patch(rect)
plt.show()

Your plot:

Here is my answer that is purely algorithmic:

I iterate over the possible rectangles that could fit in the mask, in decreasing order, area-wise. The first one being populated with only 1, is necessarily one of the largest rectangular bounding box containing you're looking for:

import matplotlib.patches as patches
import matplotlib.pyplot as plt
import numpy as np
from array import array

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

mask = np.ones_like(mask)

def get_divisor(i: int) -> array:
    div_arr = array("I", []) # I is for unsigned int (memory optimization)
    for d in range(1,(i+1)//2 + 1):
        if i%d == 0:
            div_arr.append(d)
    return div_arr


def get_divisors_until(n: int) -> dict[int,array]:
    if n==0 : return None
    divisors = {}
    for i in range(1,n+1):
        divisors[i] = get_divisor(i)
    return divisors

def find_box_of_ones(mask: np.ndarray, height: int, width: int):
    '''Assume that mask contains at least a 1,
       AND mask contains only 0 or 1 (boolean mask)
       AND height <= mask.shape[0] 
       AND width <= mask.shape[1]'''
    H, W = mask.shape
    off_x, off_y = H - height + 1, W - width + 1

    box_area = height * width
    for x in range(off_x):
        for y in range(off_y):
            nb_of_ones = np.sum(mask[x: x + height, y: y + width])
            if nb_of_ones == box_area:
                return True, y, x, y + width - 1, x + height - 1 #order to fit your convention
    
    return False, 0

def fit(h:int, w:int, shape:tuple[int,int]) -> bool:
    return (h <= shape[0]) and (w <= shape[1])

def find_largest_bbox(mask: np.ndarray):
    if np.sum(mask) == 0: raise ValueError("The Array contains no 1")
    area = mask.size
    possible_box_heights = get_divisors_until(area)

    if area == np.sum(mask): return 0,0,*mask.shape #the whole mask is a one.

    for size, pos_heights in reversed(possible_box_heights.items()):
        for height in pos_heights:
            width = size//height #always an int but // is slightly faster
            if not fit(height, width, mask.shape):
                continue
            found, *box_coord = find_box_of_ones(mask, height, width)
            if found:
                return box_coord


bbox = find_largest_bbox(mask)

plt.close()
plt.imshow(mask)
x, y = bbox[0] - 0.5, bbox[1] - 0.5
w, h = bbox[2] - bbox[0] + 1, bbox[3] - bbox[1] + 1
rect = patches.Rectangle((x, y), w, h, linewidth=2, edgecolor="red", facecolor="none")
plt.gca().add_patch(rect)
plt.show()

I used some optimization (using arrays to store integers, for memory usage), but basically, I look for every "possible heights" of rectangles for a certain surface area. Each of theses height (at fixed size) define a single rectangle (the width is deduced by doing width = size//height)

Then the idea is to "scan" your mask with the made up rectangle and check if the window contains only ones.

An optimization I see is to make a clever use of generator, to avoid loading entirely the possible_box_heights dict before iterating on it

If you need further explanations or you found some bug I didn't see, please tell me so I can update this answer