无需显式迭代即可找到满足某些条件的最大矩形块

Question

我有一些大型二维数组，例如：

   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

因此，满足 ==1 的最大矩形块（按面积）从 (1,2) 开始，其尺寸为 (2,3)。

如何在不显式迭代的情况下使用 Mathematica 找到它？

---

注意：

为了简化您的测试，这里是我的示例之一：

matrix = ImageData@Binarize@Import@"https://i.sstatic.net/ux7tA.png"

Answer 1

这是我

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]]]

对 belisarius 图片使用 BitAnd 结果的尝试：

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

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

从好的方面来看，这段代码似乎比 David 和 Sjoerd 的代码更快，但由于某种原因，它返回一个矩形，该矩形在两个维度上都比它们的小结果。由于差异正好是一，我怀疑某处存在计数错误，但目前我找不到它。

Answer 2

好吧，只是为了证明使用函数式编程是可能的，这是我非常非常低效的强力方法：

首先，我生成所有可能的正方形的列表，按面积降序排序：

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

对于给定的矩形，我执行 ListCorrelate< /代码>。如果可以在矩阵中找到这一大小的空闲矩形，则结果中应该至少有一个数字与该矩形的面积相对应（假设矩阵仅包含 1 和 0）。我们使用 Max 进行检查。只要我们找不到匹配项，我们就会寻找较小的矩形（LengthWhile 会处理这个问题）。我们最终得到了适合矩阵的最大矩形数：

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

在我的笔记本电脑上，使用贝利撒留的图像，花了 156 秒才找到第 11774+1 个矩形（+1，因为 LengthWhile 返回最后一个不适合的矩形的编号）是适合的最大矩形

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

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

Answer 3

一个可行的选择是忽略该格言以避免迭代。

第一个例程是在给定固定宽度的情况下查找最大长度。在转置矩阵上使用它来反转这些维度。它通过分而治之的方式工作，因此速度相当快。

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
  ]

这是较慢的扫描代码。我们找到最大尺寸，对于其中一个尺寸，我们向下扫描，最大化另一个尺寸，直到我们知道我们无法改进最大面积。我想到的唯一效率是根据先前的下限增加下限，以便使 maxLength 调用稍微快一些。

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}
  ]

这比 Sjoerd 的好一个数量级，只是可怕而不是可怕^2。

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

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

丹尼尔·利希布劳

Answer 4

我无法与 Heike 的逻辑竞争，但我可以稍微重构她的代码。

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]]
]

我相信这更干净，而且运行速度快了大约 20%。

Answer 5

您认为卷积是显式迭代吗？如果没有，那么它可以用来做你想做的事。使用简单的内核（例如 3x3 1），您可以快速将那些不连续的 1 归零。

编辑：

Mathematica 有内置的卷积函数，你可以使用它，或者自己编写：

这是伪代码（当然未经测试:)

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

之后剩下的是 1 的连续正方形区域。您可以进行连通区域分析并从中确定最大区域。