Non-Disjunct Rectangle Edge Covering for 2D Squares on a Grid

2017-02-07

An answer to this question on Stack Overflow.

Question

Even though the title sounds complicated, my actual problem should not be too hard to model. However, I have not been able to find a good algorithm to do the following:

I want to cover a set of squares on a grid with a fixed number n of rectangles. These rectangles may overlap and they only have to cover the outer edges of my shape.

Why Not Brute-Force?

The number of different rectangles on a sqare m x m grid is

$\frac{m^2 (1+m)^2}{4}$ .

Thereby the number of combinations a brute-force-approach would have to try is in

$O((m^4)^n))$

That would be 27,680,640,625 combinations for a 10 x 10 grid and only 3 rectangles.

Example

The initial grid with some squares on it could look like this:

n = 1: The optimal way to cover this shape with a single rectangle would then be:

n = 2: The amount of covered empty squares can be reduced using two rectangles like this:

(Note that the center is now covered by two rectangles)

Valid Cover

I am looking for a solution that does cover at least all squares that are part of the outer edge, i.e. all filled squares that share an edge on the grid width an empty square.

All squares that are not part of the outer edge of the shape may or may not be covered, the covering rectangles may or may not intersect.

Target Function

Given a fixed number of covering rectangles n, I want to cover all filled squares but minimise the the number of covered empty squares outside of the shape. This means that the empty square in the center should not count towards the target function that has to be minimised (I could also fill all the holes before applying the algorithm without it making a difference).

The value of the target function for my example is thereby:

n  |  target function
---|-----------------
1  |  11
2  |   3

Additional Constraints

Note that the original set of squares may not be connected and that the number of non-connected subshapes may even exceed the number of covering rectangles.

Alternate Description

To simplify the problem, you could also just work on a transformed version of the input data:

Then the aim is to cover all blue squares and minimise the number of covered white squares using n rectangles which may intersect.

Answer

Well, I haven't yet thought of a P-class solution, but it did occur to me that this problem may be a good candidate for stochastic solutions.

Notably, there's an easily-defined feasible starting point: just set all of the cover rectangles to the extents of the bounding box of the target squares.

From this initial state, new valid states can be generated by reducing one of the bounds of the cover rectangles and checking to see that all of the target squares are still covered.

Further, the path between any two states is likely to be short (each rectangle can be reduced to its appropriate dimension in O(√n) time, where n is the number of squares in the bounding box), meaning that it's easy to move around the search space. Though this comes with the caveat that some possible solutions are separated by a narrow path back through the initial state, which means that rerunning the algorithm we're about to develop a few times is probably good.

Given the foregoing, simulated annealing is a possible means to address the problem. The following Python script implements it:

#!/usr/bin/env python3
import random
import numpy as np
import copy
import math
import scipy
import scipy.optimize
#Generate a grid
class Grid:
  def __init__(self,grid_array):
    self.grid    = np.array(grid_array)
    self.width   = len(self.grid[0]) #Use inclusive coordinates
    self.height  = len(self.grid)    #Use inclusive coordinates
    #Convert into a list of cells
    self.cells = {}
    for y in range(len(self.grid)):
      for x in range(len(self.grid[y])):
        self.cells[(x,y)] = self.grid[y][x]
    #Find all cells which are border cells (the ones we need covered)
    self.borders = []
    for c in self.cells:
      for dx in [-1,0,1]:                #Loop through neighbors
        for dy in [-1,0,1]:
          n = (c[0]+dx,c[1]+dy)          #This is the neighbor
          if self.cells[c]==1 and self.cells.get(n, 1)==0: #See if this cell has a neighbor with value 0. Use default return to simplify code
            self.borders.append(c)
    #Ensure grid contains only valid target cells
    self.grid = np.zeros((self.height,self.width))
    for b in self.borders:
      self.grid[b[1],b[0]] = 1
    self.ntarget = np.sum(self.grid)
  def copy(self):
    return self.grid.copy()
#A state is valid if the bounds of each rectangle are inside the bounding box of
#the target squares and all the target squares are covered.
def ValidState(rects):
  #Check bounds
  if not (np.all(0<=rects[0::4]) and np.all(rects[0::4]<g.width)): #x
    return False
  if not (np.all(0<=rects[1::4]) and np.all(rects[1::4]<g.height)): #y
    return False
  if not (np.all(0<=rects[2::4]) and np.all(rects[2::4]<=g.width)): #w
    return False
  if not (np.all(0<=rects[3::4]) and np.all(rects[3::4]<=g.height)): #h
    return False
  fullmask = np.zeros((g.height,g.width))
  for r in range(0,len(rects),4):
    fullmask[rects[r+1]:rects[r+3],rects[r+0]:rects[r+2]] = 1
  return np.sum(fullmask * g.grid)==g.ntarget
#Mutate a randomly chosen bound of a rectangle. Keep trying this until we find a
#mutation that leads to a valid state.
def MutateRects(rects):
  current_state = rects.copy()
  while True:
    rects = current_state.copy()
    c = random.randint(0,len(rects)-1)
    rects[c] += random.randint(-1,1)
    if ValidState(rects):
      return rects
#Determine the score of a state. The score is the sum of the number of times
#each empty space is covered by a rectangle. The best solutions will minimize
#this count.
def EvaluateState(rects):
  score   = 0
  invgrid = -(g.grid-1) #Turn zeros into ones, and ones into zeros
  for r in range(0,len(rects),4):
    mask = np.zeros((g.height,g.width))
    mask[rects[r+1]:rects[r+3],rects[r+0]:rects[r+2]] = 1
    score += np.sum(mask * invgrid)
  return score
#Print the list of rectangles (useful for showing output)
def PrintRects(rects):
  for r in range(0,len(rects),4):
    mask = np.zeros((g.height,g.width))
    mask[rects[r+1]:rects[r+3],rects[r+0]:rects[r+2]] = 1
    print(mask)
#Input grid is here
gridi = [[0,0,1,0,0],
         [0,1,1,1,0],
         [1,1,0,1,1],
         [0,1,1,1,0],
         [0,1,0,1,0]]
g = Grid(gridi)
#Number of rectangles we wish to solve with
rect_count = 2
#A rectangle is defined as going from (x,y)-(w,h) where (w,h) is an upper bound
#on the array coordinates. This allows efficient manipulation of rectangles as
#numpy arrays
rects = []
for r in range(rect_count):
  rects += [0,0,g.width,g.height]
rects = np.array(rects)
#Might want to run a few times since the initial state is something of a
#bottleneck on moving around the search space
sols = []
for i in range(10):
  #Use simulated annealing to solve the problem
  sols.append(scipy.optimize.basinhopping(
    func      = EvaluateState,
    take_step = MutateRects,
    x0        = rects,
    disp      = True,
    niter     = 3000
  ))
#Get a minimum solution and display it
PrintRects(min(sols, key=lambda x: x['lowest_optimization_result']['fun'])['x'])

Here's a display of the algorithm's progress for the the ten runs I specify in my example code above as a function of the number of iterations (I've added some jitter so you can see all of the lines):

You'll note that most (8/10) of the runs find the minima at 8 early on. Likewise, of the 6/10 runs that find the minima at 5, most of them do so early on. This suggests that it may be better to run many shorter searches rather than a few long searches. Choosing appropriate lengths and numbers of runs will be a matter of experimentation.

Note that EvaluateState adds points for each time an empty square is covered by a rectangle. This disincentivizes redundant coverage which may be necessary to find a solution or may result in getting to a solution faster. It's pretty common for cost functions to include this sort of thing. Experimenting with a cost function that directly asks for what you want is easy - just replace EvaluateState as follows:

#Determine the score of a state. The score is the sum of the number of times
#each empty space is covered by a rectangle. The best solutions will minimize
#this count.
def EvaluateState(rects):
  score   = 0
  invgrid = -(g.grid-1) #Turn zeros into ones, and ones into zeros
  mask = np.zeros((g.height,g.width))
  for r in range(0,len(rects),4):
    mask[rects[r+1]:rects[r+3],rects[r+0]:rects[r+2]] = 1
  score += np.sum(mask * invgrid)
  return score

Using this cost function does seem to produce better results in this case:

That may be because it provides more transition paths for rectangles between feasible states. But I'd keep the other function in mind if you hit difficulties.