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

Question

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

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 个内核仍然会花费很多时间。

干杯！

这是一张图片，显示了使用我在问题中发布的代码与使用 @bla 建议的 bwconncomp 时图像大小的函数所花费的时间：

Answer 1

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

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

Answer 2

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

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

从代码中可以看出，相同class 的片段不相碰。但不清楚来自不同 classes 的片段是否不接触（即，如果它们接触，代码会产生相同的输出）。 假设它们不存在，则可以避免这两个循环。

想法是一次性标记所有片段，而不考虑 class。因此，不是每个 class 调用一次 bwlabeln，而是只调用一次。我不知道有多少 class，但是 这可能会大大减少计算时间。

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

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

这可能是代码：

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 找到新 class 开始的索引来改进它是微不足道的.

可以使用@bla bwconncomp 的建议重写上面的代码。此函数 returns 包含一个元胞数组的结构，该元胞数组具有每个标签的所有像素的索引。然后就没有必要使用 regionprops，可以直接找到大小（如@bla 所示）并使用每个标签的第一个索引找到 class（通过索引到 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(:));

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

Huge broadcast variable, optimizing code without parfor?

memory

optimization

matlab

matrix

parfor