Fill Irregular Shape with Squares without Overlap

I have irregular shapes (think geographic boundaries) that I need to fill with a set of differently sized squares, without overlap. Priority should be placed on using the largest square possible. My approach thus far is to loop over the length/width of the image at an increment and check if a square at that (x,y) is valid. If so, save the square there and mark the region as invalid. I then repeat the process for each square size I have. It does what I want but is far too slow for practically small increments. If I try to do it in one pass using else/if on each square dimension, the smallest square will dominate. I also have no guarantee of optimality regarding square size.

nim[mask] = .5
nim[~mask] = 0
cim = nim.copy()

inc = 500
dims = [4000, 2000, 1000, 500]

regions = []
for d in dims:
    for x in range(0, nim.shape[0], inc):
        for y in range(0, nim.shape[1], inc):
            if (np.count_nonzero(cim[x:x+d, y:y+d])/(d*d) < .0001):
                nim = rect(nim, y, x, d, d, 1)
                cim = blocc(cim, y, x, d, d, 1)
                regions.append([x,y,d])
                
plt.imshow(nim)
plt.show()

Where rect() and blocc() draw the sides of the region and the filled region, respectively.

One approach to the problem is to:
A) create a grid of the largest squares
B) find a good position for the grid, and
C) subdivide the squares that aren't fully inside the shape.

The first step is to overlay the shape with a grid of the largest squares. An example of this is shown below. Green squares are fully within the shape, yellow squares are partially outside the shape, and red squares are fully outside the shape.

The second step is to find the optimal placement for the grid. Defining the smallest square as 1x1 units, and the largest square as 8x8 units, there are 64 possible placements for the grid. The best placement for the grid maximizes the number of green squares. Given two placements that have an equal number of green squares, the criteria for choosing the better placement would be

• minimize the number of yellow squares, or
• put the green squares roughly at the center of the shape

Once the grid placement has been settled, green squares are kept as is, red squares are discarded, and yellow squares are subdivided into four smaller squares. Here's an example of a yellow square being subdivided:

The same rules apply: keep the green, discard red, subdivide yellow. The process of subdividing squares continues recursively until the smallest size square is reached.

The final result looks like this: