巨大的广播变量，优化代码没有parfor？

Question

我有一个40000乘80000矩阵，我从中获得“簇”的数量（具有相同值的元素组彼此相邻），然后计算每个簇的大小。这是代码的一大块。

FRAGMENTSIZESCLASS = struct([]);  %We store the data in a structure
for class=1:NumberOfClasses
  %-First we create a binary image for each class-%
  BWclass = foto==class;
  %-Second we calculate the number of connected components (fragments)-%
  L = bwlabeln(BWclass);          %returns a label matrix, L, containing labels for the connected components in BWclass
  clear BWclass
  NumberFragments=max(max(L));
  %-Third we calculate the size of each fragment-%
  FragmentSize=zeros(NumberFragments,1);
  for f=1:NumberFragments      % potential improvement: using parfor while saring the memory between workers
    FragmentSize(f,1) = sum(L(:) == f);
  end
  FRAGMENTSIZESCLASS{class}=FragmentSize;
  clear L
end

问题是矩阵L是如此之大，如果我使用parfor循环它变成一个广播变量，然后内存成倍增加，我的内存耗尽。

有关如何解决这个问题的任何想法？我已经看过这个文件：https://ch.mathworks.com/matlabcentral/fileexchange/28572-sharedmatrix但不是一个简单的解决方案，即使我有24个内核仍然会花费很多时间。

干杯!

---

这是一张图片，显示了使用我在问题中发布的代码与使用bwconncomp的代码所需的时间，如@bla所示：

Answer 1

而不是bwlabeln使用内置函数bwconncomp，例如：

...
s=bwconncomp(BWClass);
fragmentsize=sum(cellfun(@numel,s.PixelIdxList));
....

Answer 2

请注意，内存不足的原因可能是因为您使用parfor替换代码中的两个循环之一。在这两种情况下，每个工作线程都会在处理过程中创建一个与foto大小相同的数组。请注意，在内部循环中，sum(L(:) == f)创建一个大小为L的逻辑数组，然后对其值求和（我不认为JIT足够聪明以优化该中间数组）。

简而言之，以你所做的方式在这么大的图像上并行操作是不可行的。并行化的正确方法是将图像切割成图块，并在不同的线程上处理每个图块。如果片段很小（我敢于给出名称的假设），应该可以仅使用小的重叠来处理切片（切片需要重叠，以便每个碎片完全在内部至少在切片上）。在这种情况下删除重复项有点复杂，因此并行化并非易事。但是，我希望下面的建议使得根本无需并行化代码。

从代码中可以清楚地看到，同一类的片段不会触及。但是不清楚来自不同类的片段是否不接触（即如果它们的话，代码将产生相同的输出）。假设它们没有，可以避免两个循环。

我们的想法是一次性标记所有片段，无论课程如何。因此，不是每个类调用一次bwlabeln，而是只调用一次。我不知道有多少类，但这可能会大大减少计算时间。

接下来，使用regionprops为每个片段确定其大小和类别。原则上，该操作也可以通过仅在图像上迭代一次来执行。请注意，您的代码FragmentSize(f,1) = sum(L(:) == f)会在每个片段上迭代一次图像。鉴于图像的大小，可能有数百万个碎片。这可以使时间减少6个数量级。

从这一点开始，我们只处理regionprops的输出，它可能包含（数量级）一百万个元素，一个微不足道的数量（比像素数少三个数量级）。

这可能是代码：

L = bwlabeln(foto>0);
cls  = regionprops(L,foto,'MinIntensity','Area');
clear L
sz = [cls.Area];
cls = [cls.MinIntensity];
NumberOfClasses = max(cls);
FRAGMENTSIZESCLASS = cell(NumberOfClasses,1);
for ii=1:NumberOfClasses
   FRAGMENTSIZESCLASS{ii} = sz(cls==ii);
end

最后一个循环可能没有必要，我没有找到一个快速的选择。我无法想象它的价格昂贵，但是如果是这样的话，通过对cls进行排序并使用diff来找到新类开始的索引来并行化它或者改进它是微不足道的。

可以使用@ bla的bwconncomp建议重写上面的代码。此函数返回一个包含单元数组的结构，该数组具有每个标签的所有像素的索引。然后没有必要使用regionprops，可以直接找到大小（如@bla所示）并使用每个标签的第一个索引来查找类（通过索引到foto）：

cc = bwconncomp(foto>0);
sz = cellfun(@numel,cc.PixelIdxList);
cls = cellfun(@(indx)foto(indx(1)),cc.PixelIdxList);
NumberOfClasses = max(cls);
FRAGMENTSIZESCLASS2 = cell(NumberOfClasses,1);
for ii=1:NumberOfClasses
   FRAGMENTSIZESCLASS2{ii} = sz(cls==ii);
end

对于具有63个片段的256x256的小测试图像，这快3到4倍。但是，考虑到你正在处理的图像的大小，我担心这可能实际上是非常低效的。要知道的唯一方法是尝试两种方法并计算时间！

---

关于您的代码的一些注意事项：

FRAGMENTSIZESCLASS = struct([]);

您将其初始化为空结构数组，但随后使用{}对其进行索引，将其转换为单元数组。如上所述，预分配数组总是好的：

FRAGMENTSIZESCLASS = cell(NumberOfClasses,1);

NumberFragments=max(max(L));

这将创建L在水平轴上的最大投影（80k元素），然后在其中找到最大值。像在其他地方那样重塑矩阵更有效：

NumberFragments = max(L(:));