代码之家 › 专栏 › 技术社区 › Karpov

巨大的广播变量,没有parfor优化代码?

parfor optimization memory matrix matlab

Karpov · 技术社区 · 7 年前

我有一个40000 x 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核,仍然需要很多时间。

干杯

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

2 回复 | 直到 6 年前

bla 7 年前

而不是 bwlabeln 使用内置功能 bwconncomp ,例如:

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

Cris Luengo 7 年前

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

简而言之,在如此大的映像上以您的方式并行化操作是不可行的。将其并行化的正确方法是将图像切割成块,并在不同的线程上处理每个块。如果碎片很小(我敢假设这个名字),应该可以只使用一个小的重叠来处理瓷砖(瓷砖需要重叠,以便每个碎片都完全在瓷砖内部,至少在瓷砖上)。在这种情况下,删除重复项要复杂一些,因此并行化并不简单。 然而,我希望下面的建议根本没有必要对代码进行并行化。

从代码中可以清楚地看到,同一类的片段不会相互影响。但不清楚不同类的片段是否不接触(即,如果它们接触,代码将产生相同的输出)。 假设他们没有 ,可以避免两个循环。

其想法是一次性标记所有片段,不考虑类别。所以与其打电话 bwlabeln 每堂课一次,你只叫它一次。我不知道有多少课,但是 这可能会大大减少计算时间。

接下来,使用 regionprops 确定每个片段的大小和类别。原则上,也可以通过仅对图像进行一次迭代来执行此操作。请注意,您的代码, 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 查找新类开始的索引。

可以使用@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

对于256x256、63个片段的小测试图像来说,这要快3到4倍。然而,考虑到您正在处理的图像的大小,我担心这实际上可能效率很低。 唯一知道的方法是尝试这两种方法并计时!

关于代码的一些注意事项:

FRAGMENTSIZESCLASS = struct([]);

将其初始化为空结构数组,但使用 {} 将其转换为单元格数组。正如我在上面所做的那样,预先分配数组总是很好的:

FRAGMENTSIZESCLASS = cell(NumberOfClasses,1);

NumberFragments=max(max(L));

这将创建L在水平轴上的最大投影(80k元素),然后在此范围内找到最大投影。与其他地方一样,重塑矩阵更有效:

NumberFragments = max(L(:));