Algorithm for finding the fewest rectangles to cover a set of rectangles without overlapping

I have a set of rectangles and I would like to "reduce" the set so I have the fewest number of rectangles to describe the same area as the original set. If possible, I would like it to also be fast, but I am more concerned with getting the number of rectangles as low as possible. I have an approach now which works most of the time.

Currently, I start at the top-left most rectangle and see if I can expand it out right and down while keeping it a rectangle. I do that until it can't expand anymore, remove and split all intersecting rectangles, and add the expanded rectangle back in the list. Then I start the process again with the next top-left most rectangle, and so on. But in some cases, it doesn't work. For example:

With this set of three rectangles, the correct solution would end up with two rectangles, like this:

However, in this case, my algorithm starts by processing the blue rectangle. This expand downwards and splits the yellow rectangle (correctly). But then when the remainder of the yellow rectangle is processed, instead of expanding downwards, it expands right first and takes back the portion that was previously split off. Then the last rectangle is processed and it can't expand right or down, so the original set of rectangles is left. I could tweak the algorithm to expand down first and then right. That would fix this case, but it would cause the same problem in a similar scenario that was flipped.

Edit: Just to clarify, the original set of rectangles do not overlap and do not have to be connected. And if a subset of rectangles are connected, the polygon which completely covers them can have holes in it.

Despite the title to your question, I think you’re actually looking for the minimum dissection into rectangles of a rectilinear polygon. (Jason’s links are about minimum covers by rectangles, which is quite a different problem.)

David Eppstein discusses this problem in section 3 of his 2010 survey article Graph-Theoretic Solutions to Computational Geometry Problems, and he gives a nice summary in this answer on mathoverflow.net:

> The idea is to find the maximum number of disjoint axis-parallel diagonals that have two concave vertices as endpoints, split along those, and then form one more split for each remaining concave vertex. To find the maximum number of disjoint axis-parallel diagonals, form the intersection graph of the diagonals; this graph is bipartite so its maximum independent set can be found in polynomial time by graph matching techniques.

Here’s my gloss on this admirably terse description, using figure 2 from Eppstein’s article. Suppose we have a rectilinear polygon, possibly with holes.

When the polygon is dissected into rectangles, each of the concave vertices must be met by at least one edge of the dissection. So we get the minimum dissection if as many of these edges as possible do double-duty, that is, they connect two of the concave vertices.

So let’s draw the axis-parallel diagonals between two concave vertices that are contained entirely within the polygon. (‘Axis-parallel’ means ‘horizontal or vertical’ here, and a diagonal of a polygon is a line connecting two non-adjacent vertices.) We want to use as many of these lines as possible in the dissection as long as they don’t intersect.

(If there are no axis-parallel diagonals, the dissection is trivial—just make a cut from each concave vertex. Or if there are no intersections between the axis-parallel diagonals then we use them all, plus a cut from each remaining concave vertex. Otherwise, read on.)

The intersection graph of a set of line segments has a node for every line segment, and an edge joins two nodes if the lines cross. Here’s the intersection graph for the axis-parallel diagonals:

It’s bipartite with the vertical diagonals in one part, and the horizontal diagonals in the other part. Now, we want to pick as many of the diagonals as possible as long as they don’t intersect. This corresponds to finding the maximum independent set in the intersection graph.

Finding the maximum independent set in a general graph is an NP-hard problem, but in the special case of a bipartite graph, König’s theorem shows that it’s equivalent to the problem of finding a maximum matching, which can be solved in polynomial time, for example by the Hopcroft–Karp algorithm. A given graph can have several maximum matchings, but any of them will do, as they all have the same size. In the example, all the maximum matchings have three pairs of vertices, for example {(2, 4), (6, 3), (7, 8)}:

(Other maximum matchings in this graph include {(1, 3), (2, 5), (7, 8)}; {(2, 4), (3, 6), (5, 7)}; and {(1, 3), (2, 4), (7, 8)}.)

To get from a maximum matching to the corresponding minimum vertex cover, apply the proof of König’s theorem. In the matching shown above, the left set is L = {1, 2, 6, 7}, the right set is R = {3, 4, 5, 8}, and the set of unmatched vertices in L is U = {1}. There is only one alternating path starting in U, namely 1–3–6, so the set of vertices in alternating paths is Z = {1, 3, 6} and the minimum vertex cover is thus K = (L  Z) ∪ (R ∩ Z) = {2, 3, 7}, shown in red below, with the maximum independent set in green:

Translating this back into the dissection problem, this means that we can use five axis-parallel diagonals in the dissection:

Finally, make a cut from each remaining concave vertex to complete the dissection:

Today I found O(N^5) solution for this problem, and I will share it here.

For the first step, you need to find a way to get the sum of any rectangle in a matrix, with complexity O(1). It's pretty easy to do.

Now for the second step, you need to know dynamic programming. The idea is to store a rectangle and break it into smaller pieces. If the rectangle is empty, you can return 0. And if it's filled, return 1.

There are N^4 states to store the rectangle, plus the O(N) complexity for each state... So you will get an O(N^5) algorithm.

Here's my code. I think it will help.

The input is simple. N, M (size of matrix) After that, the following N lines will have 1s and 0s. Example:

4 9
010000010
111010111
101111101
000101000

#include <bits/stdc++.h>
#define MAX 51
int tab[MAX][MAX];
int N,M;
int sumed[MAX][MAX];
int t(int x,int y) {
    if(x<0||y<0)return 0;
    return sumed[x][y];
}
int subrec(int x1,int y1,int x2,int y2) {
    return t(x2,y2)-t(x2,y1-1)-t(x1-1,y2)+t(x1-1,y1-1);
}
int resp[MAX][MAX][MAX][MAX];
bool exist[MAX][MAX][MAX][MAX];
int dp(int x1,int y1,int x2,int y2) {
    if(exist[x1][y1][x2][y2])return resp[x1][y1][x2][y2];
    exist[x1][y1][x2][y2]=true;
    int soma = subrec(x1,y1,x2,y2);
    int area = (x2-x1+1)*(y2-y1+1);
    if(soma==area){return resp[x1][y1][x2][y2]=1;}
    if(!soma) {return 0;}
    int best = 1000000;
    for(int i = x1;i!=x2;++i) {
        best = std::min(best,dp(x1,y1,i,y2)+dp(i+1,y1,x2,y2));
    }
    for(int i = y1;i!=y2;++i) {
        best = std::min(best,dp(x1,y1,x2,i)+dp(x1,i+1,x2,y2));
    }
    return resp[x1][y1][x2][y2]=best;
}
void backtracking(int x1,int y1,int x2,int y2) {
    int soma = subrec(x1,y1,x2,y2);
    int area = (x2-x1+1)*(y2-y1+1);
    if(soma==area){std::cout<<x1+1<<" "<<y1+1<<" "<<x2+1<<" "<<y2+1<<"\n";return;}
    if(!soma) {return;}
    int best = 1000000;
    int obj = resp[x1][y1][x2][y2];
    for(int i = x1;i!=x2;++i) {
        int ans = dp(x1,y1,i,y2)+dp(i+1,y1,x2,y2);
        if(ans==obj){
            backtracking(x1,y1,i,y2);
            backtracking(i+1,y1,x2,y2);
            return;
        }
    }
    for(int i = y1;i!=y2;++i) {
        int ans = dp(x1,y1,x2,i)+dp(x1,i+1,x2,y2);
        if(ans==obj){
            backtracking(x1,y1,x2,i);
            backtracking(x1,i+1,x2,y2);
            return;
        }
    }
}
int main()
{
    std::cin >> N >> M;
    for(int i = 0; i != N;++i) {
        std::string s;
        std::cin >> s;
        for(int j = 0; j != M;++j) {
            if(s[j]=='1')tab[i][j]++;
        }
    }
    for(int i = 0; i != N;++i) {
        int val = 0;
        for(int j = 0; j != M;++j) {
            val += tab[i][j];
            sumed[i][j]=val;
            if(i)sumed[i][j]+=sumed[i-1][j];
        }
    }
    std::cout << dp(0,0,N-1,M-1) << std::endl;
    backtracking(0,0,N-1,M-1);
}