Find the largest rectangular block satisfying some condition without iterating explicitly

Question

I have a few large 2D arrays like:

   1   2  3  4  5
   --------------
1 | 0  1  1  1  0
2 | 0  1  1  1  0
3 | 0  1  0  1  1
4 | 0  1  0  1  1

So, the largest rectangular block (by area) satisfying ==1 starts at (1,2) and its dimensions are (2,3).

How to find it with Mathematica without iterating explicitly?

---

NB:

Just to ease your testing here is one of my samples:

matrix = ImageData@Binarize@Import@"http://i.stack.imgur.com/ux7tA.png"

Answer 1

This is my attempt using BitAnd

maxBlock[mat_] := Block[{table, maxSeq, pos},

  maxSeq[list_] := 
   Max[Length[#] & /@ Append[Cases[Split[list], {1 ..}], {}]];

  table = 
   Flatten[Table[
     MapIndexed[{#2[[1]], maxSeq[#1]} &, 
      FoldList[BitAnd[#1, #2] &, mat[[k]], Drop[mat, k]]], {k, 1, 
      Length[mat]}], 1];

  pos = Ordering[(Times @@@ table), -1][[1]];

  {Times[##], {##}} & @@ table[[pos]]]

Result for belisarius' picture:

Timing[maxBlock[Unitize[matrix, 1.]]]

(* {1.13253, {23433, {219, 107}}} *)

On the plus side this code seems faster than David's and Sjoerd's code, but for some reason it returns a rectangle which is one smaller in both dimensions than their result. Since the difference is exactly one I suspect a counting error somewhere but I can't find it at the moment.

Answer 2

Well, just to prove it's possible using functional programming here's my terribly, terribly inefficient brute force approach:

First, I generate a list of all possible squares, sorted in order of descending area:

rectangles = Flatten[
               Table[{i j, i, j}, 
                     {i, Length[matrix]}, 
                     {j, Length[matrix[[1]]]}
               ],1 
             ] // Sort // Reverse;

For a given rectangle I do a ListCorrelate. If a free rectangle of this size can be found in the matrix there should be at least one number in the result that corresponds to the area of that rectangle (assuming the matrix contains only 1's and 0's). We check that using Max. As long as we don't find a match we look for smaller rectangles (LengthWhile takes care of that). We end up with the largest rectangle number that fits in the matrix:

LengthWhile[
   rectangles, 
   Max[ListCorrelate[ConstantArray[1, {#[[2]], #[[3]]}], matrix]] != #[[1]] &
]

On my laptop, using belisarius' image, it took 156 seconds to find that the 11774+1th rectangle (+1 because the LengthWhile returns the number of the last rectangle that doesn't fit) is the largest one that will fit

In[70]:= rectangles[[11774 + 1]]

Out[70]= {23760, 220, 108}

Answer 3

A viable option is to ignore the dictum to avoid iteration.

First a routine to find the largest length given a fixed width. Use it on the transposed matrix to reverse those dimensions. It works by divide and conquer, so is reasonably fast.

maxLength[mat_, width_, min_, max_] := Module[
  {len = Floor[(min + max)/2], top = max, bottom = min, conv},
  While[bottom <= len <= top,
   conv = ListConvolve[ConstantArray[1, {len, width}], mat];
   If[Length[Position[conv, len*width]] >= 1,
    bottom = len;
    len = Ceiling[(len + top)/2],
    top = len;
    len = Floor[(len + bottom)/2]];
   If[len == bottom || len == top, Return[bottom]]
   ];
  bottom
  ]

Here is the slower sweep code. We find the maximal dimensions and for one of them we sweep downward, maximizing the other dimension, until we know we cannot improve on the maximal area. The only efficiency I came up with was to increase the lower bounds based on prior lower bounds, so as to make the maxLength calls slightly faster.

maxRectangle[mat_] := Module[
  {min, dims = Dimensions[mat], tmat = Transpose[mat], maxl, maxw, 
   len, wid, best},
  maxl = Max[Map[Length, Cases[Map[Split, mat], {1 ..}, 2]]];
  maxw = Max[Map[Length, Cases[Map[Split, tmat], {1 ..}, 2]]];
  len = maxLength[tmat, maxw, 1, maxl];
  best = {len, maxw};
  min = maxw*len;
  wid = maxw - 1;
  While[wid*maxl >= min,
   len = maxLength[tmat, wid, len, maxl];
   If[len*wid > min, best = {len, wid}; min = len*wid];
   wid--;
   ];
  {min, best}
  ]

This is better than Sjoerd's by an order of magnitude, being only terrible and not terrible^2.

In[364]:= Timing[maxRectangle[matrix]]

Out[364]= {11.8, {23760, {108, 220}}}

Daniel Lichtblau