Find the largest rectangular block satisfying some

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

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

Answer 2

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 3

Do you consider convolution as iterating explicitly? If not then it can be used to do what you want. With a simple kernel, say 3x3 1s, you can quickly zero out those non-contiguous 1s.

Edit:

Mathematica has built-in convolution function, you can use that, or brew your own:

Here's the pseudo code (untested of course :)

kernel = [ [1,1,1], [1,1,1], [1,1,1] ]

for row = 1, row <= image_height - 1, row++
  for col = 1, col <= image_width - 1, col++
    compare kernel with the 3x3 matrix at image(row, col):
      if there is 0 on left AND right of the center column, OR
      if there is 0 on top AND bottom of center row, THEN
         zero out whole area from image(row-1, col-1) to image(row+1, col+1)
         # The above may need refinement
  end
end

After that what's left are the contiguous square area of 1s. You can do connected area analysis and determine the biggest area from there.

Answer 4

I cannot compete with Heike's logic, but I can refactor her code a little.

maxBlock[mat_] := Module[{table, maxSeq, pos, i},
  maxSeq = Max[0, Length /@ Split@# ~Cases~ {1 ..}] &;
  table = Join @@
    Table[
       {i++, maxSeq@j},
       {k, Length@mat},
       {j, i = 1; FoldList[BitAnd, mat[[k]], mat~Drop~k]}
    ];
  pos = Ordering[Times @@@ table, -1][[1]];
  {# #2, {##}} & @@ table[[pos]]
]

I believe this is cleaner, and it runs about 20% faster.

Answer 5

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}