Find all rectangles in a 2D array, trying to cover as much area with as little amount of rectangles as possible. http://stackoverflow.com/questions/5810649/finding-rectangles-in-a-2d-block-grid

find_rects.cpp

#include <stdlib.h>
#include <vector>
#include <utility>
#include <algorithm>

#include "find_rects.hpp"

namespace FindRects {

/* Algorithm was taken from here: http://stackoverflow.com/a/20039017/624766
   http://www.drdobbs.com/database/the-maximal-rectangle-problem/184410529
*/

struct Pair {
	int x;
	int y;
	Pair(int a, int b): x(a), y(b) {}
};

static void update_cache(std::vector<int> *c, int n, int M, int rowWidth, const InputType* input) {
	for (int m = 0; m != M; ++m) {
		if (input[n * rowWidth + m] != 0) {
			(*c)[m]++;
		} else {
			(*c)[m] = 0;
		}
	}
}

static Rect findBiggest(const InputType* input, int M, int N, int rowWidth) {

	Pair best_ll(0, 0); /* Lower-left corner */
	Pair best_ur(-1, -1); /* Upper-right corner */
	int best_area = 0;
	int best_perimeter = 0;

	std::vector<int> c(M + 1, 0); /* Cache */
	std::vector<Pair> s; /* Stack */
	s.reserve(M + 1);

	int m, n;

	/* Main algorithm: */
	for (n = 0; n != N; ++n) {
		int open_width = 0;
		update_cache(&c, n, M, rowWidth, input);
		for (m = 0; m != M + 1; ++m) {
			if (c[m] > open_width) { /* Open new rectangle? */
				s.push_back(Pair(m, open_width));
				open_width = c[m];
			} else if (c[m] < open_width) { /* Close rectangle(s)? */
				int m0, w0, area, perimeter;
				do {
					m0 = s.back().x;
					w0 = s.back().y;
					s.pop_back();
					area = open_width * (m - m0);
					perimeter = open_width + (m - m0);
					/* If the area is the same, prefer squares over long narrow rectangles,
						it finds more rectangles this way when calling findAll() with minLength == 2 or more */
					if (area > best_area || (area == best_area && perimeter < best_perimeter)) {
						best_area = area;
						best_perimeter = perimeter;
						best_ll.x = m0;
						best_ll.y = n;
						best_ur.x = m - 1;
						best_ur.y = n - open_width + 1;
					}
					open_width = w0;
				} while (c[m] < open_width);
				open_width = c[m];
				if (open_width != 0) {
					s.push_back(Pair(m0, w0));
				}
			}
		}
	}
	return Rect(best_ll.x, std::max(0, best_ur.y), 1 + best_ur.x - best_ll.x, 1 + best_ll.y - best_ur.y);
}

Rect findBiggest(const InputType* input, int width, int height) {
	return findBiggest(input, width, height, width);
}

//#define CHECK_BOTH_WAYS 1 /* This will make the algorithm terribly slow, with a factorial complexity */

/** Find biggest rectangle, then recursively search area to the left/right and to the top/bottom of that rectangle
	for smaller rectangles, and choose the one that covers biggest area.
	@return biggest area size, covered by rectangles with side length = maxLength or bigger,
	@param search: limit searching for the following area
	@param output: may be NULL, then the function will only return the area size
*/

static long long findRectsInArea (const InputType* input, int rowWidth, int minLength, OutputType* output, Rect search) {
	if (search.w < minLength || search.h < minLength) {
		return 0; // We reached a size limit
	}
	Rect biggest = findBiggest(input + search.y * rowWidth + search.x, search.w, search.h, rowWidth);

	if (biggest.w < minLength || biggest.h < minLength) {
		return 0; // No rectangles here
	}
	biggest.x += search.x;
	biggest.y += search.y;
	if (biggest.w > OUTPUT_MAX_LENGTH) {
		biggest.w = OUTPUT_MAX_LENGTH;
	}
	if (biggest.h > OUTPUT_MAX_LENGTH) {
		biggest.h = OUTPUT_MAX_LENGTH;
	}
	/* We got two choices to split remaining area into four rectangular regions, where (B) is the biggest rectangle:
		****|***|***				************
		****|***|***				************
		****BBBBB***				----BBBBB---
		****BBBBB***		vs		****BBBBB***
		****BBBBB***				----BBBBB---
		****|***|***				************
		We are not filling the output array in the first recursive call, it's just for determining the resulting area size
	*/

	if (output != NULL) {
		for (int y = biggest.y; y < biggest.y + biggest.h; y++) {
			for (int x = biggest.x; x < biggest.x + biggest.w; x++) {
				output[(y * rowWidth + x) * 2] = 0;
				output[(y * rowWidth + x) * 2 + 1] = 0;
			}
		}
		output[(biggest.y * rowWidth + biggest.x) * 2] = biggest.w;
		output[(biggest.y * rowWidth + biggest.x) * 2 + 1] = biggest.h;
	}

#ifdef CHECK_BOTH_WAYS
	long long splitHorizArea =
		findRectsInArea(input, rowWidth, minLength, NULL,
			Rect(search.x, search.y, biggest.x - search.x, search.h)) +
		findRectsInArea(input, rowWidth, minLength, NULL,
			Rect(biggest.x + biggest.w, search.y, search.x + search.w - biggest.x - biggest.w, search.h)) +
		findRectsInArea(input, rowWidth, minLength, NULL,
			Rect(biggest.x, search.y, biggest.w, biggest.y - search.y)) +
		findRectsInArea(input, rowWidth, minLength, NULL,
			Rect(biggest.x, biggest.y + biggest.h, biggest.w, search.y + search.h - biggest.y - biggest.h));

	long long splitVertArea =
		findRectsInArea(input, rowWidth, minLength, NULL,
			Rect(search.x, search.y, search.w, biggest.y - search.y)) +
		findRectsInArea(input, rowWidth, minLength, NULL,
			Rect(search.x, biggest.y + biggest.h, search.w, search.y + search.h - biggest.y - biggest.h)) +
		findRectsInArea(input, rowWidth, minLength, NULL,
			Rect(search.x, biggest.y, biggest.x - search.x, biggest.h)) +
		findRectsInArea(input, rowWidth, minLength, NULL,
			Rect(biggest.x + biggest.w, biggest.y, search.x + search.w - biggest.x - biggest.w, biggest.h));

	/* Inefficiently perform the recursive call again, this time with non-NULL output array */
	if (splitHorizArea > splitVertArea) {
		if (output != NULL) {
#endif
			return (long long)biggest.w * biggest.h +
			findRectsInArea(input, rowWidth, minLength, output,
				Rect(search.x, search.y, biggest.x - search.x, search.h)) +
			findRectsInArea(input, rowWidth, minLength, output,
				Rect(biggest.x + biggest.w, search.y, search.x + search.w - biggest.x - biggest.w, search.h)) +
			findRectsInArea(input, rowWidth, minLength, output,
				Rect(biggest.x, search.y, biggest.w, biggest.y - search.y)) +
			findRectsInArea(input, rowWidth, minLength, output,
				Rect(biggest.x, biggest.y + biggest.h, biggest.w, search.y + search.h - biggest.y - biggest.h));
#ifdef CHECK_BOTH_WAYS
		}
		return splitHorizArea + (long long)biggest.w * biggest.h;
	} else {
		if (output != NULL) {
			findRectsInArea(input, rowWidth, minLength, output,
				Rect(search.x, search.y, search.w, biggest.y - search.y));
			findRectsInArea(input, rowWidth, minLength, output,
				Rect(search.x, biggest.y + biggest.h, search.w, search.y + search.h - biggest.y - biggest.h));
			findRectsInArea(input, rowWidth, minLength, output,
				Rect(search.x, biggest.y, biggest.x - search.x, biggest.h));
			findRectsInArea(input, rowWidth, minLength, output,
				Rect(biggest.x + biggest.w, biggest.y, search.x + search.w - biggest.x - biggest.w, biggest.h));
		}
		return splitVertArea + (long long)biggest.w * biggest.h;
	}
#endif
}

long long findAll (const InputType* input, int width, int height, int minLength, OutputType* output) {
	return findRectsInArea(input, width, minLength, output, Rect(0, 0, width, height));
}

} /* namespace */

find_rects.hpp

#ifndef __FIND_RECTS_HPP__
#define __FIND_RECTS_HPP__

#include <climits>

namespace FindRects {

	/** Input aray type, typically the smallest integral type */
	typedef unsigned char InputType;
	/** Output array type, which consists of the width and height of any given rectangle at this point */
	typedef unsigned short OutputType;
	enum { OUTPUT_MAX_LENGTH = USHRT_MAX };

	/** A rectangle, (x,y) is the upper-left corner, coordinates start at (0,0) */
	struct Rect {
		int x;
		int y;
		int w;
		int h;
		Rect(int x, int y, int w, int h): x(x), y(y), w(w), h(h) {}
		/** Returns true if a point (x,y) is inside this rectangle */
		bool isPointInside(int x, int y) { return x >= this->x && x < this->x + this->w && y >= this->y && y < this->y + this->h; }
	};

	/** Find a biggest rectangle in a given two-dimensional array
		@param data: two-dimensional array, value 0 is empty area, any non-zero value is considered for searching
		@param width: a width of data array
		@param height: a height of data array
	*/
	Rect findBiggest(const InputType* data, int width, int height);

	/** Split an area into rectangles, trying to cover as much area as possible
		@param minLength: minimal length of the rectangle side, can be 1 to count every element in the array
		@param output: output array of width x height * 2, where each two elements are width/height of a rectangle's top-left corner, or 0 if the area is inside a rectangle
		@return total area size of all rectangles, covered by rectangles with side length = maxLength or bigger,
	*/
	long long findAll (const InputType* input, int width, int height, int minLength, OutputType* output);

}

#endif

find_rects_test.cpp

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "find_rects.hpp"

enum { MM = 16, NN = 12 };

const unsigned char testdata[NN * MM] = {
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0,
0, 0, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0,
0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0,
0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0,
0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0,
0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0,
0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0,
};

int main() {
	FindRects::Rect ret = FindRects::findBiggest(testdata, MM, NN);
	printf("The maximal rectangle has area %d.\n", ret.w * ret.h);
	printf("Location: x=%d, y=%d, w=%d, h=%d\n",
									ret.x, ret.y, ret.w, ret.h);
	for (int y = 0; y < NN; y++) {
		for (int x = 0; x < MM; x++) {
			printf(" %d%s", testdata[y * MM + x], ret.isPointInside(x, y) ? "*" : " ");
		}
		printf("\n");
	}
	unsigned short *output = new unsigned short[MM * NN * 2];
	long long area;
	int numrects;

	/* The algorithm cannot find two 3x3 rectangles in the test data,
		because it finds 6x2 rectangle first, then discards it, because the rectangle height is too small. */
	memset(output, 0, sizeof(unsigned short) * MM * NN * 2);
	area = FindRects::findAll(testdata, MM, NN, 3, output);
	numrects = 0;
	printf("Searching for rects with side length 3 or more:\n", area);
	for (int y = 0; y < NN; y++) {
		for (int x = 0; x < MM; x++) {
			printf(" %d,%d", output[(y * MM + x) * 2], output[(y * MM + x) * 2 + 1]);
			if (output[(y * MM + x) * 2] > 0) {
				numrects++;
			}
		}
		printf("\n");
	}
	printf("Found %d rects with total area %lld\n", numrects, area);
	fflush(stdout);

	memset(output, 0, sizeof(unsigned short) * MM * NN * 2);
	area = FindRects::findAll(testdata, MM, NN, 2, output);
	numrects = 0;
	printf("Searching for rects with side length 2 or more:\n", area);
	for (int y = 0; y < NN; y++) {
		for (int x = 0; x < MM; x++) {
			printf(" %d,%d", output[(y * MM + x) * 2], output[(y * MM + x) * 2 + 1]);
			if (output[(y * MM + x) * 2] > 0) {
				numrects++;
			}
		}
		printf("\n");
	}
	printf("Found %d rects with total area %lld\n", numrects, area);
	fflush(stdout);

	memset(output, 0, sizeof(unsigned short) * MM * NN * 2);
	area = FindRects::findAll(testdata, MM, NN, 1, output);
	printf("Searching for rects with side length 1 or more:\n", area);
	for (int y = 0; y < NN; y++) {
		for (int x = 0; x < MM; x++) {
			printf(" %d,%d", output[(y * MM + x) * 2], output[(y * MM + x) * 2 + 1]);
			if (output[(y * MM + x) * 2] > 0) {
				numrects++;
			}
		}
		printf("\n");
	}
	printf("Found %d rects with total area %lld\n", numrects, area);
	fflush(stdout);

	unsigned short *verifydata = new unsigned short[MM * NN];
	memset(verifydata, 0, sizeof(unsigned short) * MM * NN);

	for (int y = 0; y < NN; y++) {
		for (int x = 0; x < MM; x++) {
			int w = output[(y * MM + x) * 2];
			int h = output[(y * MM + x) * 2 + 1];
			for (int yy = y; yy < y + h; yy++) {
				for (int xx = x; xx < x + w; xx++) {
					verifydata[yy * MM + xx] = 1;
				}
			}
		}
	}

	printf("\n");

	for (int y = 0; y < NN; y++) {
		for (int x = 0; x < MM; x++) {
			if ((verifydata[y * MM + x] != 0) != (testdata[y * MM + x] != 0)) {
				printf("\nERROR: test data does not match calculated rectangle data at x=%d, y=%d\n", x, y);
				return 1;
			}
		}
	}

	delete[] verifydata;
	delete[] output;

	return 0;
}