D = bwdistgeodesic(logical(YourImage), ~YourImage);
Now iterate. if there are no finite values in D, break the loop. Find the largest value in D. It is at the corner of a square of white pixels of that dimension, but you will need to do a little work to figure out which corner it is. When you have figured out which corner it is, record that square location and then set all the D locations inside that square to be inf or nan. Now return back to the beginning of the loop (looking for the largest remaining value)
This is a "greedy" algorithm. It makes no attempt at all to be consistent in block sizes. If there were a situation in which two squares of equal size plus one small square would cover a section, then the algorithm would instead take the larger square leaving the rest to be broken up into a number of smaller squares.