Cover all 1 pixels in a 0-1 image with minimum number of fixed size square

Problem Definition

For a 0-1 image in pixel shape of MxN and a given integer A, a feasible solution of the problem is defined as a certain number of AxA shape of pixel squares whose union may cover all of the 1-pixels in MxN image, in which the edges of the squares are parallel to edge of MxN image and aligned to pixels. For all of the feasible solutions, I need to find any one of it who contains minimum number of squares.

What I want to ask is, is there any algorithm that may solve such problem or similar problems? Or, any suggestions about the region(Computer Vision? or Discrete Optimization?) I can seek for the algorithm I want?

Example of feasible solution

Problem Set

M = N = 25, A = 15
The to-be-covered area is colored with grey in the image below

Feasible Solution

Three squares covering part of the area, colored with higher saturation.
And the union of the squares, covering all of the grey pixels
So the three squares make up a feasible solution of the problem, and actually its also one of the optimal solution.

Other Information

The to-be-cover area is almost definite to be not convex, but connected.
Actually I don't need a policy to find the mathematical optimal solution. I'd prefer a policy that may find a good enough, almost feasible solution in limited time/complexity.

Solution

Sub-optimal solution would be to use sliding window (square), find the biggest intersection area with the original image and remove that part from the image, then repeat util nothing is left. Here how you can do it with python.

Reading the image

from skimage import io, filters
from skimage.transform import resize
import matplotlib.pyplot as plt

image = 1-io.imread("img.png", as_gray=True).astype(np.uint8)
image = resize(image, (25, 25))
t = filters.threshold_otsu(image)
image = (image>=t).astype(np.uint8)
plt.imshow(image, cmap=plt.cm.gray)

Approach

def get_square_masks(image, square_size=15):

    # copy the image to keep it unchanged
    image_c = image.copy()

    # this will create a view of (square_size, square_size) sliding window with stride=1
    image_parts = util.view_as_windows(image_c, (square_size, square_size), 1)
    masks = []

    # repeat until you will get blank image
    while image_c.sum() !=0 :
        
        overlapping_areas = image_parts.sum(axis=(2, 3))
        r, c = np.where(overlapping_areas==overlapping_areas.max())

        # taking the first coordinates of maximum overlapping area,
        # nice to have a check for others
        # since it is view, setting 0 will set overlapped squares with 0 as well
        image_parts[r[0], c[0]] = 0 

        # creating the square mask
        r_start = r[0]
        c_start = c[0]
        r_end = r_start + square_size
        c_end = c_start + square_size
        mask = np.zeros_like(image_c)
        mask[r_start:r_end,c_start:c_end] = 1
        
        masks.append(mask)
        
    return masks

masks = get_square_masks(image, 15)

I got 4 masks with this naive approach. But I think it is good starting point. Further you may try to optimize the masks by considering them as potential candidates. Since the number of masks depends on the first choice, we can iterate through the retrieved masks and pick other than the first one as our first choice and see whether it is more optimal, something like this:

import numpy as np

def try_reduce_masks(image, masks, square_size=15):
    masks_min = masks
    for m in masks[1:]:
        # remove the mask from image
        new_image = image*np.logical_xor(image, m).astype(np.uint8)
        # get new masks with different starting point
        new_masks = get_square_masks(new_image, square_size)
        new_masks.append(m)
        # by the end if get less masks then save it
        if len(new_masks) < len(masks_min):
            masks_min = new_masks
    return masks_min

masks = try_reduce_masks(image, masks, 15)

In general I think we may start with boundary masks, i.e. with masks which has the biggest boundary intersection. The first iteration would be to find those boundary squares, and then continue as shown.

However, if you will have a look to overlapping_areas from get_square_masks function you will get overlapping areas per each step of sliding window (square). I believe from this map you can find smarter way of finding more optimal solution.

array([[177, 181, 181, 179, 175, 166, 153, 139, 124, 109,  94],
       [178, 183, 184, 183, 180, 172, 159, 146, 132, 117, 102],
       [175, 181, 183, 183, 181, 175, 163, 151, 138, 124, 109],
       [170, 177, 180, 182, 182, 178, 168, 157, 145, 132, 117],
       [161, 170, 175, 179, 181, 179, 171, 162, 151, 139, 125],
       [150, 161, 168, 174, 178, 178, 172, 165, 156, 145, 132],
       [137, 150, 159, 167, 173, 175, 171, 166, 159, 149, 137],
       [123, 137, 148, 158, 166, 170, 168, 165, 160, 152, 141],
       [108, 122, 135, 147, 157, 163, 163, 162, 159, 153, 144],
       [ 94, 107, 121, 134, 146, 154, 156, 157, 156, 152, 145],
       [ 82,  94, 107, 120, 133, 143, 147, 150, 151, 149, 144]])

Visualization

plt.imshow(image)
for m in masks:
    plt.imshow(m, alpha=0.3, interpolation='nearest', cmap=plt.cm.gray)
plt.show()

I believe this can be solved by reinforcement learning by picking appropriate reward function. For example we can punish for small overlapping area and reward for less steps.

EDIT

May be slightly better approach would be to prevent cutting the image into 2 ore more parts. For that purpose at each step we can check the number of connected components, and if it is more than 1, than recover the image and pick the next biggest, something like this:

from skimage import measure, util
import numpy as np

def get_square_masks(image, square_size=15):

    # copy the image to keep it unchanged
    image_c = image.copy()

    # this will create a view of (square_size, square_size) sliding window with stride=1
    image_parts = util.view_as_windows(image_c, (square_size, square_size), step=1)
    masks = []

    # repeat until you will get blank image
    while image_c.sum() !=0 :
        
        overlapping_areas = image_parts.sum(axis=(2, 3))
        overlapping_areas_flat = overlapping_areas.ravel()
        sorted_indices = np.argsort(overlapping_areas_flat)[::-1]
        rr, cc = np.unravel_index(sorted_indices, overlapping_areas.shape)

        # itereate by the biggest overlapping areas, but prevent to have 2 or more connected components
        for r, c in zip(rr, cc):
            # copy the part to recover in case it will cut the image into 2 or more parts
            image_part_c = image_parts[r, c].copy()
            image_parts[r, c] = 0
            # check if you still have one connected component
            if measure.label(image_c).max() <= 1:
                break
            image_parts[r, c] = image_part_c
                    
        # creating the square mask
        r_start = r
        c_start = c
        r_end = r_start + square_size
        c_end = c_start + square_size
        mask = np.zeros_like(image_c)
        mask[r_start:r_end,c_start:c_end] = 1
        
        masks.append(mask)
        
    return masks